- Email: [email protected]
Micromagnetic study of statistical switching in magnetic tunnel junctions stabilized by perpendicular shape anisotropy
Journal Pre-proof Micromagnetic study of statistical switching in magnetic tunnel junctions stabilized by perpendicular shape anisotropy M. d'Aquino, S. Perna, C. Serpico PII:
S0921-4526(19)30646-5
DOI:
https://doi.org/10.1016/j.physb.2019.411744
Reference:
PHYSB 411744
To appear in:
Physica B: Physics of Condensed Matter
Received Date: 30 June 2019 Revised Date:
30 September 2019
Accepted Date: 2 October 2019
Please cite this article as: M. d'Aquino, S. Perna, C. Serpico, Micromagnetic study of statistical switching in magnetic tunnel junctions stabilized by perpendicular shape anisotropy, Physica B: Physics of Condensed Matter (2019), doi: https://doi.org/10.1016/j.physb.2019.411744. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier B.V. All rights reserved.
Micromagnetic study of statistical switching in magnetic tunnel junctions stabilized by perpendicular shape anisotropy M. d’Aquino1 , S. Perna2 , C. Serpico2 1 Engineering
Department, University of Naples ”Parthenope”, 80143 Naples, Italy. University of Naples Federico II, 80125 Naples, Italy.
2 DIETI,
Abstract Magnetic tunnel junctions with free layer thermally stabilized by the combined action of perpendicular magnetocrystalline and shape anisotropy are considered. The spin-transfer torque driven switching dynamics is studied. Analytical formulas for switching time statistical distributions and write-error rate are derived under the assumption of ballistic dynamics. The analytical results are compared with full micromagnetic and macrospin simulations based on stochastic Landau-Lifshitz-Gilbert-Slonczewski equation. This comparison reveals quantitative agreement with the proposed analytical theory. Keywords: Magnetic tunnel junction, perpendicular shape anisotropy, switching times statistics, write-error rates, micromagnetic simulations
1. Introduction The continual demand in magnetic recording of memory devices with increasingly higher information density and faster operation is pushing through the minimization of the area associated with the single bit cell, while the read/write process is realized via spintronic technology. One of the most promising single bit cell device in this route 5
is the magnetic tunnel junction (MTJ), where the information is written due to the spin-transfer torque mechanism and is read due to the tunnel magneto-resistance effect (TMR). The non-volatile nature of the bit stored, the high TMR which allows large enough signal to noise ratio (SNR) in the reading process for low current values, small dissipation and high speed during the writing process, are strong indicators of the possibility to achieve a universal memory[1]. However, the scaling process has limits well summarized by the trilemma of the magnetic recording. Indeed, reducing
10
the volume of the bit cell, the thermal stability can be achieved only if the anisotropy is increased. This in turn implies the increase of the threshold value of dc field or spin-polarized current for the switching process (critical switching current), which is limited by the technology. To overcome this issue, intense research activity has been performed in the last decades and several strategies have been proposed. For instance, precessional switching[2, 3] relies on magnetization reversal driven by a strong torque produced by tilting the external field with respect to the equilibrium 1 Corresponding
author tel. +390815476767, fax +390815476777, email. [email protected] (M. d’Aquino)
Preprint submitted to Physica B
October 4, 2019
15
magnetization, but suffers the need of extremely precise design of the field pulse duration[4, 5]. Another possibility is energy assisted switching[6], where the magnetization switching is achieved by the spin-polarized dc currents or dc magnetic fields combined with a supplementary external action. The latter action can be a heat pulse (heat assisted magnetic recording, HAMR[7, 8]) or a microwave magnetic field (microwave assisted magnetic recording MAMR[9, 10]). In this way, despite the volume reduction and the consequent anisotropy increment, the switching
20
process can be achieved with sustainable current (field) values. A more recent approach for enhancing the information density is the realization of perpendicular MTJs (PMTJs) where the free layer has thickness size of the order or even bigger than the diameter[11],[12]. In this situation, the shape anisotropy helps the crystalline anisotropy to enforce the thermal stabilization of the magnetic state. With this technology, it is possible to satisfy the industrial requirement of a thermal stability factor ∆ at the temperature T = 300 K such that ∆300 > 60 for sub-20 nm diameter cells (the thermal
25
stability factor is defined as ∆T = E/(kB T ), with E being the energy barrier which separates the two metastable states and kB the Boltzmann constant). Moreover, the increase of the thickness reduces the damping α due to the weakness of the spin-pumping effect. This is reflected in small dissipation and lower switching current. In this work, the free layer of the PMTJ[11] is modeled as a magnetic cylinder with diameter d = 10 nm and thickness t = 15 nm, as sketched in fig.1, with saturation magnetization such that µ0 M s = 1.5 T (µ0 is the vacuum
30
permeability), magneto-crystalline uniaxial anisotropy constant Ku = 1.10 × 105 J/m3 which favours magnetization orientations along the cylinder symmetry axis z, and damping α = 0.00425 (the values of parameters are taken from Ref.[11]). The device is assumed to be in contact with a thermal bath at fixed temperature T = 300 K. For such a device, the magnetic metastable states are spatially-uniform and aligned along the symmetry axis z. In this case, the combined action of the magneto-crystalline anisotropy and the shape anisotropy produces the thermal stability. It has
35
been shown[11] that, for T = 300 K, the thermal stability factor ∆300 ranges between 60 and 80, which yields typical values required for magnetic memories[13]. For such values of ∆T , in the absence of any external driving force, thermal activation acts on a time scale which is much larger than that associated with the driven switching dynamics of interest. This also occurs during the switching process, if the driving force (field or current) is well above the threshold value to achieve the switching.
40
Thus, in this condition, the magnetization dynamics can be assumed to be purely ballistic. On the other hand, the effect of thermal fluctuations at room temperature, in the absence of external actions, determines the fact that the initial magnetization orientation follows suitable statistical distribution. As a consequence, the resulting spin-torque switching time will obey appropriate statistical distribution arising from the combination of the above arguments. In this paper, the statistical description of switching times will be given in terms of their probability distribution function
45
(pdf), cumulative distribution function (cdf) and write-error-rate (WER). The latter quantities will be analytically derived as function of material and geometrical parameters and compared with the results of macrospin and full micromagnetic simulations.
2
𝑑 𝑡
𝑧 𝑦 𝑥
Figure 1: Sketch of the MTJ device: the metallic nano-contacts are depicted in grey, the spin current polarizer in green, the non conductive layer in purple and the free layer in yellow.
2. Analytical theory of switching times statistical distributions Magnetization dynamics is described by the stochastic Landau-Lifshitz-Gilbert-Slonczewski (LLGS) equation[14], which in normalized form reads as: dm = m × ∇Σ gL − α∇Σ φ − νm × hN (t) , dt
(1)
where m = M/M s (T ) is the magnetization unit vector, gL (m, ha ) the normalized Gibbs-Landau free energy, ∇Σ the nabla operator on the unit sphere Σ : m2 = 1, α the damping coefficient and φ = gL − [β/(αc p )] ln(1 + c p m · p) the potential describing damping effects and Slonczewski spin-transfer torque effects, where β = 2η J/J p is the dimensionless spin-polarized current with η being the polarization factor, J p = (eγM s2 t)/(g µB ) a normalization current density (γ is the absolute value of the gyromagnetic ratio, e is the absolute value of the electron charge, t is the free layer thickness, g is the Land´e factor, µB is the Bohr magneton) and c p quantifies the spin-torque asymmetry. The term hN (t) is a three-dimensional Gaussian white noise where ν is a parameter which measures its intensity. Equation (1) is interpreted in the Stratonovic sense[15], such that the stochastic magnetization dynamics occurs on the unit sphere Σ. This equation is associated with a partial differential equation, governing the evolution of the transition probability density w(m, t|m0 , t0 ) [14], which is given by: ∂t w = −∇Σ · J w ,
3
(2)
where w satisfies the following normalization condition: Z w(m, t|m0 , t0 ) dm = 1 , Σ
(3)
with the initial condition w(m, t0 |m0 , t0 ) = δ(m − m0 ), while J w is the probability density current, expressed by: J w = (m × ∇Σ gL − α∇Σ φ) w −
ν2 ∇Σ w . 2
(4)
The stationary probability distribution weq is responsible for a divergenceless probability current. It can be seen [14], that weq is given by the following expression: weq (m) =
1 exp(−µ φ(m)) , Z(µ)
(5)
where Z is the partition function and µ = 2α/ν2 . In order to guarantee that, before the injection of current, the stationary distribution is consistent with the Boltzmann statistics and the normalization used, the following Einstein relation must hold[16]:
1 ν2 kB T . = = µ 2α µ0 M s2 V
(6)
Let us assume that the applied field and the polarizer of the spin current are oriented along the symmetry axis z, namely: ha = haz ez and e p = ez . This produces the following expression for the free energy: gL =
Dk keff 1 m · D · m − ha · m = (1 − m2z ) − haz mz + , 2 2 2
(7)
where keff = D⊥ − Dk (Dk , D⊥ are effective demagnetizing factors along symmetry and transverse axes, respectively), and for the potential associated with non-conservative effects: φ = gL −
β ln(1 + c p mz ) . α cp
(8)
In the absence of injected current (β = 0) and large energy barrier ∆T = µkeff /2, the magnetization orientation is statistically distributed in one energy well (e.g. surrounding the equilibrium mz = 1) according to the following distribution expressed in terms of the tilting angle θ (such that mz = cos θ) [17]: ! keff 2 θ . weq (θ) = µkeff θ exp −µ 2
(9)
The dynamics is forced by the spin polarized current and is assumed to occur in a purely ballistic regime, namely in a phase flow fashion[18]. The more the switching current exceeds the threshold value, the more this approximation is accurate. This means that the dynamics can be studied with eq.(1) by setting ν = 0. Moreover, due to the cylindrical symmetry, only the dynamics of magnetization component mz is relevant for the switching process, which is governed by the following equation:
dmz ∂φ = −α(1 − m2z ) . dt ∂mz 4
(10)
It is shown in ref.[17] that exact integration of eq. (10) can be performed in both cases of symmetric (c p = 0) and asymmetric (c p 6= 0) spin transfer torques. In general, the obtained expression can be written as: ts = 0 where keff = keff + αβ c p , h˜ = haz −
β α
1 (1±c p )
1 0 q(H, mz0 , mz f ) , α keff
(11)
± c p (the ± sign when c p 6= 0 depends on whether parallel to anti-parallel
˜ 0 and such that |H|> 1, t s is the switching time, mz0 (mz f ) is the initial switching is considered or viceversa), H = h/k eff (final) condition and the function q(H, mz0 , mz f ) is: ! ! ! H + mz f 1 − mz f 1 + mz f 1 1 1 q(H, mz0 , mz f ) = − ln ln + ln − 2 . 2(H + 1) 1 − mz0 2(H − 1) 1 + mz0 H + mz0 H −1
(12)
With this notation, assuming the initial conditions distributed around mz = 1 according to eq.(9), the switching time t s is defined as the time needed to reach mz = mz f . For a given value mz f (e.g. mz f = −0.9), eq.(11) can be written as function of the sole mz0 : ts =
1 q(mz0 , mz f ) = g(mz0 ) . α keff
(13)
It is important, for what will be discussed in the sequel, to perform the inversion of the function g(mz0 ). In general, such an inversion has to be performed numerically. Nevertheless, it turns out that an approximate expression of the function g−1 (t s ) can be obtained with algebraic manipulation of the function q defined in eq.(12). In fact, one can be easily convinced that the first term is, in general, dominant over the remaining two. Thus, the idea is to derive an approximation q˜ of q in the following form: ! 1 − mz f 1 q˜ = (H, mz0 , mz f ) = − ln + Q2 + Q3 , 2(H + 1) 1 − mz0
(14)
where Q2 , Q3 are appropriate constants the values of which are determined in order to minimize the difference between q and q. ˜ We propose two possible choices for Q2 , Q3 . The former is based on computing the following averages of the second and third term in eq.(12):
Q2 =
1 1 + mz f
Z
1 Q3 = − 1 + mz f
1
−mz f
Z
1 + mz f 1 ln dmz0 = 2(H − 1) 1 + mz0 !
1
−mz f
2 ln
mz f +1 2
m +1 + ln − mzzff −1 + 1 (mz f − 1) + 2
, 2(H − 1)(mz f + 1) H+mz f H+mz f ! ln − + 1 (H − mz f ) ln + 1 H−mz f H + mz f 1 H+1 ln dm = − + . z0 H + mz0 (H − 1)(1 + mz f ) H2 − 1 (H 2 − 1)(1 + mz f )
(15)
(16)
A simpler possibility is given by taking Q2 , Q3 equal to the values of second and third term in eq.(12) at mz0 = −mz f : ! 1 + mz f 1 ln , 2(H − 1) 1 − mz f ! H + mz f 1 ln Q3 = − 2 . H − mz f (H − 1) Q2 =
5
(17) (18)
mzf=-0.9
0.0 01
0.0 0.0 1 05
0.001
0.001
0.001 -3 0.001
mz0
0.001
-4.5
0.001 0.95
-5 0.8
1
0.01 0.005
1 0.00
01 0.0
0.9
00.1 .05
0.001
05 0.0
-4
1 0.00
05 0.0 0.85
0. 0.01 00 5
01 0.0
-3.5
-4
-5 0.8
1 0.0 5 0 0.0
0.01 5 -2.5 0.00
-3.5 0.0 05
-4.5
.1 00.05
0.1 0.05
-2
H
H
-1.5
01 0.0
0.01 0.005
-2.5 0.005 0.001 -3 0.001
00.1 .05
01 0.0
00..015 1 0.0 05 0 0.
0.1 0.05
-1.5 0.1 0.05 -2
mzf=-0.9
0.85
0.9
mz0
0.95
1
Figure 2: Error analysis for formula (14). Diagram of Relative error of |q(H, mz0 , mz f )− q(H, ˜ mz0 , mz f )| for: (left) Q2 , Q3 values computed according to (15)-(16); (right) Q2 , Q3 values computed according to (17)-(18).
By using the function q˜ in place of q, it is possible to perform exact inversion of the function g(mz0 ), which ends up with the expression: 0 mz0 = g−1 (t s ) = 1 − (1 − mz f ) exp[2(α keff t s − Q2 − Q3 )(H + 1)] .
(19)
It is worth to remark that the functions q, q˜ depend only on the dimensionless quantities H, mz0 , mz f , therefore they 50
describe a dynamical behaviour that is somehow universal in this kind of switching process regardless of the particular choice of material and geometrical parameters of the considered system. The accuracy of such approximations has been tested and some results of error analysis as function of mz0 , H for a given value mz f = −0.9 are reported in figure 2. It is apparent that the determination of Q2 , Q3 according to eqs.(15)-(16) is generally more accurate than that coming from eqs.(17)-(18). Indeed, in the large energy barrier limit,
55
the initial conditions mz0 are very close to mz = 1, with statistical distribution given by eq.(9) (see ref. [17]). For these values of mz0 , the approximation of Q2 , Q3 according to eqs.(15)-(16) works better, although both of them for |H|> 1.5 give a relative error which is in the order of percent or even less. We proceed now to the determination of statistical distributions for switching times. The assumption of ballistic switching is crucial in deriving the switching times probability distribution function (pdf). Indeed, one has: Z w(mz = mz f , t s ) = w(mz = mz f , t s |mz0 ) weq (mz0 ) dmz0 ,
(20)
X0
where X0 is the ensemble of initial conditions (it is understood that we assumed t0 = 0). Given t s , there exists one and only one mz0 = g−1 (t s ) such that after the time t s from the current pulse, one has mz (t s ) = mz f . For this reason, one has w(mz = mz f , t s |mz0 ) = δ(mz0 − g−1 (ts)) and then eq.(20) can be written as: w(mz = mz f , t s ) = weq (g−1 (t s )) . 6
(21)
Thus, denoting the switching times pdf with fT s (t s ), one has: fT s (t s ) = weq (g−1 (ts))
dg−1 (t s ) . dt s
By using eqs.(9), (10), one ends up with[17]: " # keff 2 −1 0 ˜ − (g−1 (t s ))2 ]| , fT s (t s ) = µkeff exp −µ arccos (g (t s )) |α(keff g−1 (t s ) + h)[1 2
(22)
(23)
where g−1 (t s ) can be obtained either from the numerical inversion of eq.(12) or can be directly approximated by eq.(19). Analogously, one can obtain the cumulative distribution function (cdf) of switching times, namely: # " keff arccos2 (g−1 (t s )) . F(t s ) = exp −µ 2
(24)
Once the analytical expression of F(t s ) is known, it is immediate to get the formula for the WER. Indeed, it is given by the following equation, which expresses the probability of not switching after time t s : " # keff 2 −1 WER = 1 − FT s (t s ) = 1 − exp −µ arccos (g (t s )) . 2
(25)
It can be also useful deriving the asymptotic behavior of the WER for large t s . From eq.(25), one can see that, in log scale and for t s → +∞, WER(t s ) approaches a straight line with negative slope S WER and starting from WER0 for t s = 0. By using the fact that g−1 (t s )|ts →+∞ → ±1 (+ sign applies for parallel to anti-parallel switching P → AP and − sign for AP → P), it is readily seen that the asymptotic slope and the offset of the asymptotic straight line (WER0 ), are given by:
0 2αkeff d log WER = ± (±1 + H) , t s →+∞ dt s γ Ms 0 WER0 = lim (log WER − S WER t s ) = log ±µ keff (±1 − mz f ) exp(−2(1 ± H) q0 ) ,
S WER = lim
t s →+∞
(26) (27)
where H ≤ −1 during the P → AP switching process (H ≥ 1 for AP → P switching), and the function q0 is given by: ! ! 1 ± mz f H + mz f 1 1 q0 (H, mz f ) = ± − 2 ln . (28) ln 2(±H − 1) 2 H± 1 H −1 Equations (26),(27),(28) provide a very simple analytical tool for quickly drawing asymptotic write-error rate diagrams, namely: 0 log(WER) ≈ WER0 + S WER t s
0 t s < − WER S WER 0 t s ≥ − WER S WER
(29)
We remark that formulas (19),(23)-(28) are closed-form expressions which depend on material and geometrical parameters of the considered magnetic particle and contain no adjustable parameters. 60
In the next section, the numerical study of the considered device by means of macrospin[19] and full micromagnetic simulations[20] is addressed. Then, the analytical formulas are compared with the numerical results. Finally, the validity and the predictive power of the proposed analytical theory are discussed. 7
1
z
s
[M ]
0.5
0
-0.5 magnetic macrospin -1 0
2
4
6
8 10 time [ns]
12
14
16
Figure 3: Comparison between magnetization switching dynamics obtained from macrospin simulation (red curve) and full micromagnetic simulation (blue curve).
3. Numerical results Validity of macrospin approximation. Here, the switching dynamics at room temperature T = 300K for the considered 65
system is simulated. We adopt geometrical and material parameters recently used for experimental demonstration of a PMTJ stabilized by shape anisotropy[11]. To this end, the stochastic LLGS equation (1) is integrated starting with initial magnetization mz = +1 and, after 4 ns, the spin-polarized current pulse is applied. The time interval before the current pulse application is chosen to let the magnetization relax in the energy well around mz = 1. Then, magnetization switching driven by symmetric spin-transfer torque (c p = 0) occurs. The spin-polarized current
70
amplitude is chosen equal to twice the critical switching current. In all the simulations, the switching process is considered complete when mz has reached the final value mz f = −0.9. Both macrospin and micromagnetic simulations are performed by using geometrical integration of LLGS based on the midpoint rule[21]. The adopted mesh cell size for micromagnetic simulations is 2 × 2 × 3 nm3 . In figure 3, a comparison of the switching dynamics of the magnetization z-component obtained integrating the macrospin LLGS and the full micromagnetic LLGS is reported
75
for a single realization. Apart from a small shift in the switching time of about 2 ns, the dynamics appears quite similar. Moreover, stronger indicators of the validity of the macrospin approximation can be obtained by looking at the comparison of the respective free energy quantities (anisotropy, demagnetizing and exchange) and at the average magnetization amplitude normalized by the saturation magnetization. These information are shown in figures 4 and 5, respectively.
80
The range of values for the anisotropy and magnetostatic free energies (G AN and G D in figure 4) is almost identical 8
0.3
GD (magnetic) GEX (magnetic) GAN (magnetic)
0.2
GD (macrospin)
0
s
Free Energy [ M2V]
0.25
GEX (macrospin)
0.15
GAN (macrospin)
0.1 0.05 0 0
2
4
6
8 10 time [ns]
12
14
16
Figure 4: Comparison among different free energy terms (see legend for details) during magnetization switching obtained from macrospin and full micromagnetic simulations.
1
s
|| [M ]
0.99 0.98 0.97 0.96 0.95 0
2
4
6
8 10 time [ns]
12
14
16
Figure 5: Time evolution of spatially-averaged magnetization amplitude, normalized by the saturation magnetization at T = 300 K such that µ0 M s = 1.5 T .
for macrospin and micromagnetic simulations. A substantial energy difference is present on the exchange free energy value which, in the macrospin case is zero, while it is obviously nonzero in the micromagnetic framework. This is due to the effect of the thermal field, which tries to create disuniformities against the exchange field. In this context, the cell size plays an important role. By choosing a smaller cell size, for example a 1 × 1 × 1 nm3 cell, we checked 9
85
that, while the anisotropy and magnetostatic energy are almost identical to those computed for larger mesh size, the exchange energy increases by a large amount (10 times larger). This fact can be reasonably explained from the expression of ν used for micromagnetic integration of the LLGS. The only difference with respect to eq.(6), is the substitution ∆V → V, where ∆V is the computational cell volume[20]. The smaller the cell size, the larger the intensity of thermal fluctuations, and therefore the larger the exchange energy produced by spatial inhomogeneities.
90
Moreover, the increase of ν produces non-uniformities which disappear for bigger cell size. This problem is wellknown in computational micromagnetics[22] and is still relatively open. In figure 5, it is also shown that the average magnetization amplitude computed from the micromagnetic simulation is such that |< M > |≈ 1 throughout the whole switching process. This is a strong numerical evidence for the validity of the macrospin approximation. Similar results have been obtained for higher current values. In addition, it has been
95
checked that, up to J sw ∼ 10 MA/cm2 , the Oersted magnetic field produced by the electric current is µ0 HOe ∼ 10−3 T and, therefore, is negligible with respect to anisotropy and demagnetizing fields which are µ0 HAN , µ0 H M ∼ 10−1 T. Switching times statistical distributions. In order to determine the statistical description of the switching process, the stochastic LLGS equation (1) is repeatedly integrated for a large number N = 1000 of PMTJ replicas, as prescribed by the approach referred to as Langevin dynamics[19]. Each realization of the simulated switching process starts with
100
initial magnetization mz = +1 and, after 4 ns, the spin-polarized current pulse is applied. As mentioned before, the time interval before the current pulse application is chosen to let the magnetization relax in the energy well around mz = 1. Then, magnetization switching driven by symmetric spin-transfer torque (c p = 0) occurs. The spin-polarized current values are chosen J sw = 2 J sw,crit , J sw = 3 J sw,crit , J sw = 4 J sw,crit , where J sw,crit is the critical switching current provided by the analytical prediction βcrit = αkeff obtained by imposing the condition
105
H = −1 = −βcrit /(αkeff ) for P → AP switching (see section 2). In dimensional units, assuming a spin-polarization factor η ∼ 0.1 it follows that J sw,crit ∼ 25 MA/cm2 , which is of the same order of magnitude as the switching current measured in ref.[11]. Figure 6 shows the comparison of the switching time probability distributions obtained from micromagnetic and macrospin simulations. As expected from the aforementioned discussion, the agreement is quantitative. Next, we
110
compare the statistical distributions provided by the analytical theory with the results of macrospin simulations. The function g−1 (t s ) in eq.(23) has been computed both from numerical inversion of eq.(12) and from the approximate expression (19). The comparison is reported in figure 7. One can see that for J sw = 2 J sw,crit the magnetization switching is slightly faster than that predicted by the analytical model. This discrepancy is due to the thermal fluctuations neglected during the switching dynamics, which become relatively more important when the switching current value
115
approaches its threshold value. As expected, the agreement dramatically improves for higher current values where magnetization dynamics is fully compliant with the assumption of ballistic switching. We also observe that the approximate closed-form formula (19) for g−1 (t s ) leads to statistical distributions that are almost indistinguishable from those computed by performing numerical inversion of eq.(12). 10
0.7
Jsw = Jsw = Jsw = Jsw = Jsw = Jsw =
Switching Times pdf
0.6 0.5 0.4
2*Jsw,crit 3*Jsw,crit 4*Jsw,crit 2*Jsw,crit 3*Jsw,crit 4*Jsw,crit
(magnetic) (magnetic) (magnetic) (macrospin) (macrospin) (macrospin)
0.3 0.2 0.1 0 0
5
10 15 switching time [ns]
20
Figure 6: Switching times statistical distributions obtained from micromagnetic simulations (colored dots), and macrospin simulations (dashed lines).
0.7
Jsw = Jsw = Jsw = Jsw = Jsw = Jsw = Jsw = Jsw = Jsw =
Switching Times pdf
0.6 0.5 0.4 0.3
2*Jsw,crit 3*Jsw,crit 4*Jsw,crit 2*Jsw,crit 3*Jsw,crit 4*Jsw,crit 2*Jsw,crit 3*Jsw,crit 4*Jsw,crit
(macrospin) (macrospin) (macrospin) (analytic) (analytic) (analytic) (approx.) (approx.) (approx.)
0.2 0.1 0 0
5
10 switching time [ns]
15
20
Figure 7: Comparison between analytical and numerical (macrospin) calculations of the switching times statistical distributions. Colored dots represent macrospin simulations, colored dashed lines represent the result of the analytic theory where the function g−1 (t s ) is obtained from the numerical inversion of eq. (12), while black lines with different hatching represent the approximated theory where g−1 (t s ) is obtained from eq. (19).
As a final remark, it is worth to stress that analytical formulas for statistical distributions allow straightforward 120
analysis of write-error rates for memory cells as function of material and geometrical parameters. In fact, by using 11
0
10
-2
10
-4
Write-Error-Rate
10
-6
10
-8
10
-10
10
-12
10
-14
10
0
Jsw = 2*Jsw,crit Jsw = 3*Jsw,crit Jsw = 4*Jsw,crit Jsw = 2*Jsw,crit Jsw = 3*Jsw,crit Jsw = 4*Jsw,crit Jsw = 2*Jsw,crit Jsw = 3*Jsw,crit Jsw = 4*Jsw,crit
5
(mag) (mag) (mag) (analytic) (analytic) (analytic) (asymptotic) (asymptotic) (asymptotic)
10 15 switching time [ns]
20
Figure 8: Comparison between analytical and numerical (µmagnetic) calculations of the write-error rate as function of the switching time, parametrized to the threshold switching current.
eq.(25), one can determine the current pulse duration required in order to have reliable magnetization switching up to a prescribed level. For instance, by observing the results reported in figure 8, one can readily see that a current pulse with amplitude J sw = 4 J sw,crit and duration of about 20 ns is required for a WER of 10−12 , which is a realistic value for the design of magnetic random access memories (MRAM) bit cells. Moreover, by using eq.(26) the slope of 125
the WER curve in logarithmic scale, namely a measure of the sensitivity of the WER to the current pulse amplitude can be estimated. According to the previous example, for a WER of 10−12 the sensitivity is S WER ∼ −2 ns−1 . This means that an increase of 1 ns in the pulse duration affects the WER by a factor of e−2 , which yields a reduction by almost one order of magnitude. In addition, by using eqs.(26)-(29), one can draw the asymptotic WER diagram (also reported in fig.8) in a very straightforward way. It is apparent that there is remarkable agreement between the
130
analytical predictions and the results of micromagnetic simulations. Finally, we observe that a purely statistical computation of the WER curves depicted in fig.8 would require billions of Langevin simulations for a single value of the current amplitude.
Conclusions In this work, theoretical analysis of the magnetization switching process in PMTJs has been performed. Analytical 135
theory for statistical description of the switching process has been developed under the assumption of macrospin approximation and ballistic magnetization dynamics. The validity of macrospin approximation during the switching dynamics has been checked by means of micromagnetic simulations, where the different contributions to the free energy and the amplitude of the average magnetization have been used as indicators. Then, the LLGS equation has 12
been integrated for an ensemble of N = 1000 replicas of the analyzed device, with initial conditions dispersed around 140
mz = +1, both with full micromagnetic and macrospin solvers. From the dynamics, the statistics of switching times have been inferred for different switching current values. Very good agreement has been found between the switching time distributions obtained from micromagnetic and macrospin simulations. The comparison between analytical predictions and numerically computed switching time distributions has been also performed. It is found that the agreement is very good for switching current larger than twice the critical value. Finally, it is demonstrated that the
145
analytical theory provides ready-to-use tools for analyzing Write-Error-Rates, which can favor immediate comparison with experiments for the optimal design of spin-transfer torque MRAMs.
References [1] J. Akerman, Toward a universal memory, Science 308 (2005) 508–510. doi:10.1126/science.1110549. [2] H. W. Schumacher, C. Chappert, P. Crozat, R. C. Sousa, P. P. Freitas, J. Miltat, J. Ferre, Precessional magnetization reversal in microscopic 150
spin valve cells, IEEE Transactions on Magnetics 38 (5) (2002) 2480–2483. doi:10.1109/TMAG.2002.801903. [3] G. Bertotti, I. D. Mayergoyz, C. Serpico, M. d’Aquino, Geometrical analysis of precessional switching and relaxation in uniformly magnetized bodies, IEEE Transactions on Magnetics 39 (5) (2003) 2501–2503. doi:10.1109/TMAG.2003.816453. [4] S. Kaka, S. E. Russek, Precessional switching of submicrometer spin valves, Applied Physics Letters 80 (16) (2002) 2958–2960. doi: 10.1063/1.1470704.
155
[5] M. d’Aquino, W. Scholz, T. Schrefl, C. Serpico, J. Fidler, Numerical and analytical study of fast precessional switching, Journal of Applied Physics 95 (11) (2004) 7055–7057. doi:10.1063/1.1689910. [6] S. Greaves, Y. Kanai, H. Muraoka, Magnetization switching in energy assisted recording, IEEE Transactions on Magnetics 48 (5) (2012) 1794–1800. doi:10.1109/TMAG.2012.2187776. [7] M. H. Kryder, E. C. Gage, T. W. McDaniel, W. A. Challener, R. E. Rottmayer, G. Ju, Y. Hsia, M. F. Erden, Heat assisted magnetic recording,
160
Proceedings of the IEEE 96 (11) (2008) 1810–1835. doi:10.1109/JPROC.2008.2004315. [8] C. Serpico, S. Perna, G. Bertotti, M. d’Aquino, A. Quercia, I. D. Mayergoyz, Noise-induced bifurcations in magnetization dynamics of uniaxial nanomagnets, Journal of Applied Physics 117 (17) (2015) 17A709. doi:10.1063/1.4906961. [9] J. Zhu, X. Zhu, Y. Tang, Microwave assisted magnetic recording, IEEE Transactions on Magnetics 44 (1) (2008) 125–131. doi:10.1109/ TMAG.2007.911031.
165
[10] E. A. Montoya, S. Perna, Y. J. Chen, J. A. Katine, M. d’Aquino, C. Serpico, I. N. Krivorotov, Magnetization reversal driven by low dimensional chaos in a nanoscale ferromagnet, Nature Communications 10 (1). doi:10.1038/s41467-019-08444-2. [11] K. Watanabe, B. Jinnai, S. Fukami, H. Sato, H. Ohno, Shape anisotropy revisited in single-digit nanometer magnetic tunnel junctions, Nature Communications 9 (1) (2018) 663. doi:10.1038/s41467-018-03003-7. [12] N. Perrissin, S. Lequeux, N. Strelkov, A. Chavent, L. Vila, L. Buda-Prejbeanu, S. Auffret, R. Sousa, I. Prejbeanu, B. Dieny, A highly
170
thermally stable sub-20 nm magnetic random-access memory based on perpendicular shape anisotropy, Nanoscale 10 (25) (2018) 12187– 12195. doi:10.1039/C8NR01365A. [13] D. Apalkov, B. Dieny, J. M. Slaughter, Magnetoresistive random access memory, Proceedings of the IEEE 104 (2016) 1796–1830. doi: 10.1109/JPROC.2016.2590142. [14] G. Bertotti, I. D. Mayergoyz, C. Serpico, Nonlinear Magnetization Dynamics in Nanosystems, Elsevier Science, San Diego, (USA), 2009.
175
[15] C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the natural sciences, Springer–Verlag, Berlin, 1985. [16] C. Serpico, G. Bertotti, I. Mayergoyz, M. d’Aquino, R. Bonin, Thermal stability in spin torque driven magnetization dynamics, Journal of applied physics 99 (8). doi:10.1063/1.2158388.
13
[17] M. d’Aquino, V. Scalera, C. Serpico, Analysis of switching times statistical distributions for perpendicular magnetic memories, Journal of Magnetism and Magnetic Materials 475 (2019) 652–661. doi:10.1016/j.jmmm.2018.11.031. 180
[18] G. Bertotti, I. D. Mayergoyz, M. d’Aquino, S. Perna, C. Serpico, Phase-flow interpretation of magnetization relaxation in nanomagnets, IEEE Transactions on Magnetics 50 (11) (2014) 1–4. doi:10.1109/TMAG.2014.2329035. [19] M. d’Aquino, C. Serpico, G. Coppola, I. D. Mayergoyz, G. Bertotti, Midpoint numerical technique for stochastic landau-lifshitz-gilbert dynamics, Journal of Applied Physics 99 (8) (2006) 08B905. doi:10.1063/1.2169472. [20] C. Ragusa, M. d’Aquino, C. Serpico, B. Xie, M. Repetto, G. Bertotti, P. Ansalone, Full micromagnetic numerical simulations of thermal
185
fluctuations, IEEE Transactions on Magnetics 45 (10) (2009) 3919–3922. doi:10.1109/tmag.2009.2021856. [21] M. d’Aquino, C. Serpico, G. Miano, I. D. Mayergoyz, G. Bertotti, Numerical integration of landau-lifshitz-gilbert equation based on the midpoint rule, Journal of Applied Physics 97 (10) (2005) 10E319. doi:10.1063/1.1858784. [22] V. Tsiantos, W. Scholz, D. Suess, T. Schrefl, J. Fidler, The effect of the cell size in langevin micromagnetic simulations, Journal of magnetism and magnetic materials 242 (2002) 999–1001. doi:10.1016/S0304-8853(01)01365-8.
14
S0921-4526(19)30646-5
DOI:
https://doi.org/10.1016/j.physb.2019.411744
Reference:
PHYSB 411744
To appear in:
Physica B: Physics of Condensed Matter
Received Date: 30 June 2019 Revised Date:
30 September 2019
Accepted Date: 2 October 2019
Please cite this article as: M. d'Aquino, S. Perna, C. Serpico, Micromagnetic study of statistical switching in magnetic tunnel junctions stabilized by perpendicular shape anisotropy, Physica B: Physics of Condensed Matter (2019), doi: https://doi.org/10.1016/j.physb.2019.411744. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier B.V. All rights reserved.
Micromagnetic study of statistical switching in magnetic tunnel junctions stabilized by perpendicular shape anisotropy M. d’Aquino1 , S. Perna2 , C. Serpico2 1 Engineering
Department, University of Naples ”Parthenope”, 80143 Naples, Italy. University of Naples Federico II, 80125 Naples, Italy.
2 DIETI,
Abstract Magnetic tunnel junctions with free layer thermally stabilized by the combined action of perpendicular magnetocrystalline and shape anisotropy are considered. The spin-transfer torque driven switching dynamics is studied. Analytical formulas for switching time statistical distributions and write-error rate are derived under the assumption of ballistic dynamics. The analytical results are compared with full micromagnetic and macrospin simulations based on stochastic Landau-Lifshitz-Gilbert-Slonczewski equation. This comparison reveals quantitative agreement with the proposed analytical theory. Keywords: Magnetic tunnel junction, perpendicular shape anisotropy, switching times statistics, write-error rates, micromagnetic simulations
1. Introduction The continual demand in magnetic recording of memory devices with increasingly higher information density and faster operation is pushing through the minimization of the area associated with the single bit cell, while the read/write process is realized via spintronic technology. One of the most promising single bit cell device in this route 5
is the magnetic tunnel junction (MTJ), where the information is written due to the spin-transfer torque mechanism and is read due to the tunnel magneto-resistance effect (TMR). The non-volatile nature of the bit stored, the high TMR which allows large enough signal to noise ratio (SNR) in the reading process for low current values, small dissipation and high speed during the writing process, are strong indicators of the possibility to achieve a universal memory[1]. However, the scaling process has limits well summarized by the trilemma of the magnetic recording. Indeed, reducing
10
the volume of the bit cell, the thermal stability can be achieved only if the anisotropy is increased. This in turn implies the increase of the threshold value of dc field or spin-polarized current for the switching process (critical switching current), which is limited by the technology. To overcome this issue, intense research activity has been performed in the last decades and several strategies have been proposed. For instance, precessional switching[2, 3] relies on magnetization reversal driven by a strong torque produced by tilting the external field with respect to the equilibrium 1 Corresponding
author tel. +390815476767, fax +390815476777, email. [email protected] (M. d’Aquino)
Preprint submitted to Physica B
October 4, 2019
15
magnetization, but suffers the need of extremely precise design of the field pulse duration[4, 5]. Another possibility is energy assisted switching[6], where the magnetization switching is achieved by the spin-polarized dc currents or dc magnetic fields combined with a supplementary external action. The latter action can be a heat pulse (heat assisted magnetic recording, HAMR[7, 8]) or a microwave magnetic field (microwave assisted magnetic recording MAMR[9, 10]). In this way, despite the volume reduction and the consequent anisotropy increment, the switching
20
process can be achieved with sustainable current (field) values. A more recent approach for enhancing the information density is the realization of perpendicular MTJs (PMTJs) where the free layer has thickness size of the order or even bigger than the diameter[11],[12]. In this situation, the shape anisotropy helps the crystalline anisotropy to enforce the thermal stabilization of the magnetic state. With this technology, it is possible to satisfy the industrial requirement of a thermal stability factor ∆ at the temperature T = 300 K such that ∆300 > 60 for sub-20 nm diameter cells (the thermal
25
stability factor is defined as ∆T = E/(kB T ), with E being the energy barrier which separates the two metastable states and kB the Boltzmann constant). Moreover, the increase of the thickness reduces the damping α due to the weakness of the spin-pumping effect. This is reflected in small dissipation and lower switching current. In this work, the free layer of the PMTJ[11] is modeled as a magnetic cylinder with diameter d = 10 nm and thickness t = 15 nm, as sketched in fig.1, with saturation magnetization such that µ0 M s = 1.5 T (µ0 is the vacuum
30
permeability), magneto-crystalline uniaxial anisotropy constant Ku = 1.10 × 105 J/m3 which favours magnetization orientations along the cylinder symmetry axis z, and damping α = 0.00425 (the values of parameters are taken from Ref.[11]). The device is assumed to be in contact with a thermal bath at fixed temperature T = 300 K. For such a device, the magnetic metastable states are spatially-uniform and aligned along the symmetry axis z. In this case, the combined action of the magneto-crystalline anisotropy and the shape anisotropy produces the thermal stability. It has
35
been shown[11] that, for T = 300 K, the thermal stability factor ∆300 ranges between 60 and 80, which yields typical values required for magnetic memories[13]. For such values of ∆T , in the absence of any external driving force, thermal activation acts on a time scale which is much larger than that associated with the driven switching dynamics of interest. This also occurs during the switching process, if the driving force (field or current) is well above the threshold value to achieve the switching.
40
Thus, in this condition, the magnetization dynamics can be assumed to be purely ballistic. On the other hand, the effect of thermal fluctuations at room temperature, in the absence of external actions, determines the fact that the initial magnetization orientation follows suitable statistical distribution. As a consequence, the resulting spin-torque switching time will obey appropriate statistical distribution arising from the combination of the above arguments. In this paper, the statistical description of switching times will be given in terms of their probability distribution function
45
(pdf), cumulative distribution function (cdf) and write-error-rate (WER). The latter quantities will be analytically derived as function of material and geometrical parameters and compared with the results of macrospin and full micromagnetic simulations.
2
𝑑 𝑡
𝑧 𝑦 𝑥
Figure 1: Sketch of the MTJ device: the metallic nano-contacts are depicted in grey, the spin current polarizer in green, the non conductive layer in purple and the free layer in yellow.
2. Analytical theory of switching times statistical distributions Magnetization dynamics is described by the stochastic Landau-Lifshitz-Gilbert-Slonczewski (LLGS) equation[14], which in normalized form reads as: dm = m × ∇Σ gL − α∇Σ φ − νm × hN (t) , dt
(1)
where m = M/M s (T ) is the magnetization unit vector, gL (m, ha ) the normalized Gibbs-Landau free energy, ∇Σ the nabla operator on the unit sphere Σ : m2 = 1, α the damping coefficient and φ = gL − [β/(αc p )] ln(1 + c p m · p) the potential describing damping effects and Slonczewski spin-transfer torque effects, where β = 2η J/J p is the dimensionless spin-polarized current with η being the polarization factor, J p = (eγM s2 t)/(g µB ) a normalization current density (γ is the absolute value of the gyromagnetic ratio, e is the absolute value of the electron charge, t is the free layer thickness, g is the Land´e factor, µB is the Bohr magneton) and c p quantifies the spin-torque asymmetry. The term hN (t) is a three-dimensional Gaussian white noise where ν is a parameter which measures its intensity. Equation (1) is interpreted in the Stratonovic sense[15], such that the stochastic magnetization dynamics occurs on the unit sphere Σ. This equation is associated with a partial differential equation, governing the evolution of the transition probability density w(m, t|m0 , t0 ) [14], which is given by: ∂t w = −∇Σ · J w ,
3
(2)
where w satisfies the following normalization condition: Z w(m, t|m0 , t0 ) dm = 1 , Σ
(3)
with the initial condition w(m, t0 |m0 , t0 ) = δ(m − m0 ), while J w is the probability density current, expressed by: J w = (m × ∇Σ gL − α∇Σ φ) w −
ν2 ∇Σ w . 2
(4)
The stationary probability distribution weq is responsible for a divergenceless probability current. It can be seen [14], that weq is given by the following expression: weq (m) =
1 exp(−µ φ(m)) , Z(µ)
(5)
where Z is the partition function and µ = 2α/ν2 . In order to guarantee that, before the injection of current, the stationary distribution is consistent with the Boltzmann statistics and the normalization used, the following Einstein relation must hold[16]:
1 ν2 kB T . = = µ 2α µ0 M s2 V
(6)
Let us assume that the applied field and the polarizer of the spin current are oriented along the symmetry axis z, namely: ha = haz ez and e p = ez . This produces the following expression for the free energy: gL =
Dk keff 1 m · D · m − ha · m = (1 − m2z ) − haz mz + , 2 2 2
(7)
where keff = D⊥ − Dk (Dk , D⊥ are effective demagnetizing factors along symmetry and transverse axes, respectively), and for the potential associated with non-conservative effects: φ = gL −
β ln(1 + c p mz ) . α cp
(8)
In the absence of injected current (β = 0) and large energy barrier ∆T = µkeff /2, the magnetization orientation is statistically distributed in one energy well (e.g. surrounding the equilibrium mz = 1) according to the following distribution expressed in terms of the tilting angle θ (such that mz = cos θ) [17]: ! keff 2 θ . weq (θ) = µkeff θ exp −µ 2
(9)
The dynamics is forced by the spin polarized current and is assumed to occur in a purely ballistic regime, namely in a phase flow fashion[18]. The more the switching current exceeds the threshold value, the more this approximation is accurate. This means that the dynamics can be studied with eq.(1) by setting ν = 0. Moreover, due to the cylindrical symmetry, only the dynamics of magnetization component mz is relevant for the switching process, which is governed by the following equation:
dmz ∂φ = −α(1 − m2z ) . dt ∂mz 4
(10)
It is shown in ref.[17] that exact integration of eq. (10) can be performed in both cases of symmetric (c p = 0) and asymmetric (c p 6= 0) spin transfer torques. In general, the obtained expression can be written as: ts = 0 where keff = keff + αβ c p , h˜ = haz −
β α
1 (1±c p )
1 0 q(H, mz0 , mz f ) , α keff
(11)
± c p (the ± sign when c p 6= 0 depends on whether parallel to anti-parallel
˜ 0 and such that |H|> 1, t s is the switching time, mz0 (mz f ) is the initial switching is considered or viceversa), H = h/k eff (final) condition and the function q(H, mz0 , mz f ) is: ! ! ! H + mz f 1 − mz f 1 + mz f 1 1 1 q(H, mz0 , mz f ) = − ln ln + ln − 2 . 2(H + 1) 1 − mz0 2(H − 1) 1 + mz0 H + mz0 H −1
(12)
With this notation, assuming the initial conditions distributed around mz = 1 according to eq.(9), the switching time t s is defined as the time needed to reach mz = mz f . For a given value mz f (e.g. mz f = −0.9), eq.(11) can be written as function of the sole mz0 : ts =
1 q(mz0 , mz f ) = g(mz0 ) . α keff
(13)
It is important, for what will be discussed in the sequel, to perform the inversion of the function g(mz0 ). In general, such an inversion has to be performed numerically. Nevertheless, it turns out that an approximate expression of the function g−1 (t s ) can be obtained with algebraic manipulation of the function q defined in eq.(12). In fact, one can be easily convinced that the first term is, in general, dominant over the remaining two. Thus, the idea is to derive an approximation q˜ of q in the following form: ! 1 − mz f 1 q˜ = (H, mz0 , mz f ) = − ln + Q2 + Q3 , 2(H + 1) 1 − mz0
(14)
where Q2 , Q3 are appropriate constants the values of which are determined in order to minimize the difference between q and q. ˜ We propose two possible choices for Q2 , Q3 . The former is based on computing the following averages of the second and third term in eq.(12):
Q2 =
1 1 + mz f
Z
1 Q3 = − 1 + mz f
1
−mz f
Z
1 + mz f 1 ln dmz0 = 2(H − 1) 1 + mz0 !
1
−mz f
2 ln
mz f +1 2
m +1 + ln − mzzff −1 + 1 (mz f − 1) + 2
, 2(H − 1)(mz f + 1) H+mz f H+mz f ! ln − + 1 (H − mz f ) ln + 1 H−mz f H + mz f 1 H+1 ln dm = − + . z0 H + mz0 (H − 1)(1 + mz f ) H2 − 1 (H 2 − 1)(1 + mz f )
(15)
(16)
A simpler possibility is given by taking Q2 , Q3 equal to the values of second and third term in eq.(12) at mz0 = −mz f : ! 1 + mz f 1 ln , 2(H − 1) 1 − mz f ! H + mz f 1 ln Q3 = − 2 . H − mz f (H − 1) Q2 =
5
(17) (18)
mzf=-0.9
0.0 01
0.0 0.0 1 05
0.001
0.001
0.001 -3 0.001
mz0
0.001
-4.5
0.001 0.95
-5 0.8
1
0.01 0.005
1 0.00
01 0.0
0.9
00.1 .05
0.001
05 0.0
-4
1 0.00
05 0.0 0.85
0. 0.01 00 5
01 0.0
-3.5
-4
-5 0.8
1 0.0 5 0 0.0
0.01 5 -2.5 0.00
-3.5 0.0 05
-4.5
.1 00.05
0.1 0.05
-2
H
H
-1.5
01 0.0
0.01 0.005
-2.5 0.005 0.001 -3 0.001
00.1 .05
01 0.0
00..015 1 0.0 05 0 0.
0.1 0.05
-1.5 0.1 0.05 -2
mzf=-0.9
0.85
0.9
mz0
0.95
1
Figure 2: Error analysis for formula (14). Diagram of Relative error of |q(H, mz0 , mz f )− q(H, ˜ mz0 , mz f )| for: (left) Q2 , Q3 values computed according to (15)-(16); (right) Q2 , Q3 values computed according to (17)-(18).
By using the function q˜ in place of q, it is possible to perform exact inversion of the function g(mz0 ), which ends up with the expression: 0 mz0 = g−1 (t s ) = 1 − (1 − mz f ) exp[2(α keff t s − Q2 − Q3 )(H + 1)] .
(19)
It is worth to remark that the functions q, q˜ depend only on the dimensionless quantities H, mz0 , mz f , therefore they 50
describe a dynamical behaviour that is somehow universal in this kind of switching process regardless of the particular choice of material and geometrical parameters of the considered system. The accuracy of such approximations has been tested and some results of error analysis as function of mz0 , H for a given value mz f = −0.9 are reported in figure 2. It is apparent that the determination of Q2 , Q3 according to eqs.(15)-(16) is generally more accurate than that coming from eqs.(17)-(18). Indeed, in the large energy barrier limit,
55
the initial conditions mz0 are very close to mz = 1, with statistical distribution given by eq.(9) (see ref. [17]). For these values of mz0 , the approximation of Q2 , Q3 according to eqs.(15)-(16) works better, although both of them for |H|> 1.5 give a relative error which is in the order of percent or even less. We proceed now to the determination of statistical distributions for switching times. The assumption of ballistic switching is crucial in deriving the switching times probability distribution function (pdf). Indeed, one has: Z w(mz = mz f , t s ) = w(mz = mz f , t s |mz0 ) weq (mz0 ) dmz0 ,
(20)
X0
where X0 is the ensemble of initial conditions (it is understood that we assumed t0 = 0). Given t s , there exists one and only one mz0 = g−1 (t s ) such that after the time t s from the current pulse, one has mz (t s ) = mz f . For this reason, one has w(mz = mz f , t s |mz0 ) = δ(mz0 − g−1 (ts)) and then eq.(20) can be written as: w(mz = mz f , t s ) = weq (g−1 (t s )) . 6
(21)
Thus, denoting the switching times pdf with fT s (t s ), one has: fT s (t s ) = weq (g−1 (ts))
dg−1 (t s ) . dt s
By using eqs.(9), (10), one ends up with[17]: " # keff 2 −1 0 ˜ − (g−1 (t s ))2 ]| , fT s (t s ) = µkeff exp −µ arccos (g (t s )) |α(keff g−1 (t s ) + h)[1 2
(22)
(23)
where g−1 (t s ) can be obtained either from the numerical inversion of eq.(12) or can be directly approximated by eq.(19). Analogously, one can obtain the cumulative distribution function (cdf) of switching times, namely: # " keff arccos2 (g−1 (t s )) . F(t s ) = exp −µ 2
(24)
Once the analytical expression of F(t s ) is known, it is immediate to get the formula for the WER. Indeed, it is given by the following equation, which expresses the probability of not switching after time t s : " # keff 2 −1 WER = 1 − FT s (t s ) = 1 − exp −µ arccos (g (t s )) . 2
(25)
It can be also useful deriving the asymptotic behavior of the WER for large t s . From eq.(25), one can see that, in log scale and for t s → +∞, WER(t s ) approaches a straight line with negative slope S WER and starting from WER0 for t s = 0. By using the fact that g−1 (t s )|ts →+∞ → ±1 (+ sign applies for parallel to anti-parallel switching P → AP and − sign for AP → P), it is readily seen that the asymptotic slope and the offset of the asymptotic straight line (WER0 ), are given by:
0 2αkeff d log WER = ± (±1 + H) , t s →+∞ dt s γ Ms 0 WER0 = lim (log WER − S WER t s ) = log ±µ keff (±1 − mz f ) exp(−2(1 ± H) q0 ) ,
S WER = lim
t s →+∞
(26) (27)
where H ≤ −1 during the P → AP switching process (H ≥ 1 for AP → P switching), and the function q0 is given by: ! ! 1 ± mz f H + mz f 1 1 q0 (H, mz f ) = ± − 2 ln . (28) ln 2(±H − 1) 2 H± 1 H −1 Equations (26),(27),(28) provide a very simple analytical tool for quickly drawing asymptotic write-error rate diagrams, namely: 0 log(WER) ≈ WER0 + S WER t s
0 t s < − WER S WER 0 t s ≥ − WER S WER
(29)
We remark that formulas (19),(23)-(28) are closed-form expressions which depend on material and geometrical parameters of the considered magnetic particle and contain no adjustable parameters. 60
In the next section, the numerical study of the considered device by means of macrospin[19] and full micromagnetic simulations[20] is addressed. Then, the analytical formulas are compared with the numerical results. Finally, the validity and the predictive power of the proposed analytical theory are discussed. 7
1
z
s
0.5
0
-0.5 magnetic macrospin -1 0
2
4
6
8 10 time [ns]
12
14
16
Figure 3: Comparison between magnetization switching dynamics obtained from macrospin simulation (red curve) and full micromagnetic simulation (blue curve).
3. Numerical results Validity of macrospin approximation. Here, the switching dynamics at room temperature T = 300K for the considered 65
system is simulated. We adopt geometrical and material parameters recently used for experimental demonstration of a PMTJ stabilized by shape anisotropy[11]. To this end, the stochastic LLGS equation (1) is integrated starting with initial magnetization mz = +1 and, after 4 ns, the spin-polarized current pulse is applied. The time interval before the current pulse application is chosen to let the magnetization relax in the energy well around mz = 1. Then, magnetization switching driven by symmetric spin-transfer torque (c p = 0) occurs. The spin-polarized current
70
amplitude is chosen equal to twice the critical switching current. In all the simulations, the switching process is considered complete when mz has reached the final value mz f = −0.9. Both macrospin and micromagnetic simulations are performed by using geometrical integration of LLGS based on the midpoint rule[21]. The adopted mesh cell size for micromagnetic simulations is 2 × 2 × 3 nm3 . In figure 3, a comparison of the switching dynamics of the magnetization z-component obtained integrating the macrospin LLGS and the full micromagnetic LLGS is reported
75
for a single realization. Apart from a small shift in the switching time of about 2 ns, the dynamics appears quite similar. Moreover, stronger indicators of the validity of the macrospin approximation can be obtained by looking at the comparison of the respective free energy quantities (anisotropy, demagnetizing and exchange) and at the average magnetization amplitude normalized by the saturation magnetization. These information are shown in figures 4 and 5, respectively.
80
The range of values for the anisotropy and magnetostatic free energies (G AN and G D in figure 4) is almost identical 8
0.3
GD (magnetic) GEX (magnetic) GAN (magnetic)
0.2
GD (macrospin)
0
s
Free Energy [ M2V]
0.25
GEX (macrospin)
0.15
GAN (macrospin)
0.1 0.05 0 0
2
4
6
8 10 time [ns]
12
14
16
Figure 4: Comparison among different free energy terms (see legend for details) during magnetization switching obtained from macrospin and full micromagnetic simulations.
1
s
|
0.99 0.98 0.97 0.96 0.95 0
2
4
6
8 10 time [ns]
12
14
16
Figure 5: Time evolution of spatially-averaged magnetization amplitude, normalized by the saturation magnetization at T = 300 K such that µ0 M s = 1.5 T .
for macrospin and micromagnetic simulations. A substantial energy difference is present on the exchange free energy value which, in the macrospin case is zero, while it is obviously nonzero in the micromagnetic framework. This is due to the effect of the thermal field, which tries to create disuniformities against the exchange field. In this context, the cell size plays an important role. By choosing a smaller cell size, for example a 1 × 1 × 1 nm3 cell, we checked 9
85
that, while the anisotropy and magnetostatic energy are almost identical to those computed for larger mesh size, the exchange energy increases by a large amount (10 times larger). This fact can be reasonably explained from the expression of ν used for micromagnetic integration of the LLGS. The only difference with respect to eq.(6), is the substitution ∆V → V, where ∆V is the computational cell volume[20]. The smaller the cell size, the larger the intensity of thermal fluctuations, and therefore the larger the exchange energy produced by spatial inhomogeneities.
90
Moreover, the increase of ν produces non-uniformities which disappear for bigger cell size. This problem is wellknown in computational micromagnetics[22] and is still relatively open. In figure 5, it is also shown that the average magnetization amplitude computed from the micromagnetic simulation is such that |< M > |≈ 1 throughout the whole switching process. This is a strong numerical evidence for the validity of the macrospin approximation. Similar results have been obtained for higher current values. In addition, it has been
95
checked that, up to J sw ∼ 10 MA/cm2 , the Oersted magnetic field produced by the electric current is µ0 HOe ∼ 10−3 T and, therefore, is negligible with respect to anisotropy and demagnetizing fields which are µ0 HAN , µ0 H M ∼ 10−1 T. Switching times statistical distributions. In order to determine the statistical description of the switching process, the stochastic LLGS equation (1) is repeatedly integrated for a large number N = 1000 of PMTJ replicas, as prescribed by the approach referred to as Langevin dynamics[19]. Each realization of the simulated switching process starts with
100
initial magnetization mz = +1 and, after 4 ns, the spin-polarized current pulse is applied. As mentioned before, the time interval before the current pulse application is chosen to let the magnetization relax in the energy well around mz = 1. Then, magnetization switching driven by symmetric spin-transfer torque (c p = 0) occurs. The spin-polarized current values are chosen J sw = 2 J sw,crit , J sw = 3 J sw,crit , J sw = 4 J sw,crit , where J sw,crit is the critical switching current provided by the analytical prediction βcrit = αkeff obtained by imposing the condition
105
H = −1 = −βcrit /(αkeff ) for P → AP switching (see section 2). In dimensional units, assuming a spin-polarization factor η ∼ 0.1 it follows that J sw,crit ∼ 25 MA/cm2 , which is of the same order of magnitude as the switching current measured in ref.[11]. Figure 6 shows the comparison of the switching time probability distributions obtained from micromagnetic and macrospin simulations. As expected from the aforementioned discussion, the agreement is quantitative. Next, we
110
compare the statistical distributions provided by the analytical theory with the results of macrospin simulations. The function g−1 (t s ) in eq.(23) has been computed both from numerical inversion of eq.(12) and from the approximate expression (19). The comparison is reported in figure 7. One can see that for J sw = 2 J sw,crit the magnetization switching is slightly faster than that predicted by the analytical model. This discrepancy is due to the thermal fluctuations neglected during the switching dynamics, which become relatively more important when the switching current value
115
approaches its threshold value. As expected, the agreement dramatically improves for higher current values where magnetization dynamics is fully compliant with the assumption of ballistic switching. We also observe that the approximate closed-form formula (19) for g−1 (t s ) leads to statistical distributions that are almost indistinguishable from those computed by performing numerical inversion of eq.(12). 10
0.7
Jsw = Jsw = Jsw = Jsw = Jsw = Jsw =
Switching Times pdf
0.6 0.5 0.4
2*Jsw,crit 3*Jsw,crit 4*Jsw,crit 2*Jsw,crit 3*Jsw,crit 4*Jsw,crit
(magnetic) (magnetic) (magnetic) (macrospin) (macrospin) (macrospin)
0.3 0.2 0.1 0 0
5
10 15 switching time [ns]
20
Figure 6: Switching times statistical distributions obtained from micromagnetic simulations (colored dots), and macrospin simulations (dashed lines).
0.7
Jsw = Jsw = Jsw = Jsw = Jsw = Jsw = Jsw = Jsw = Jsw =
Switching Times pdf
0.6 0.5 0.4 0.3
2*Jsw,crit 3*Jsw,crit 4*Jsw,crit 2*Jsw,crit 3*Jsw,crit 4*Jsw,crit 2*Jsw,crit 3*Jsw,crit 4*Jsw,crit
(macrospin) (macrospin) (macrospin) (analytic) (analytic) (analytic) (approx.) (approx.) (approx.)
0.2 0.1 0 0
5
10 switching time [ns]
15
20
Figure 7: Comparison between analytical and numerical (macrospin) calculations of the switching times statistical distributions. Colored dots represent macrospin simulations, colored dashed lines represent the result of the analytic theory where the function g−1 (t s ) is obtained from the numerical inversion of eq. (12), while black lines with different hatching represent the approximated theory where g−1 (t s ) is obtained from eq. (19).
As a final remark, it is worth to stress that analytical formulas for statistical distributions allow straightforward 120
analysis of write-error rates for memory cells as function of material and geometrical parameters. In fact, by using 11
0
10
-2
10
-4
Write-Error-Rate
10
-6
10
-8
10
-10
10
-12
10
-14
10
0
Jsw = 2*Jsw,crit Jsw = 3*Jsw,crit Jsw = 4*Jsw,crit Jsw = 2*Jsw,crit Jsw = 3*Jsw,crit Jsw = 4*Jsw,crit Jsw = 2*Jsw,crit Jsw = 3*Jsw,crit Jsw = 4*Jsw,crit
5
(mag) (mag) (mag) (analytic) (analytic) (analytic) (asymptotic) (asymptotic) (asymptotic)
10 15 switching time [ns]
20
Figure 8: Comparison between analytical and numerical (µmagnetic) calculations of the write-error rate as function of the switching time, parametrized to the threshold switching current.
eq.(25), one can determine the current pulse duration required in order to have reliable magnetization switching up to a prescribed level. For instance, by observing the results reported in figure 8, one can readily see that a current pulse with amplitude J sw = 4 J sw,crit and duration of about 20 ns is required for a WER of 10−12 , which is a realistic value for the design of magnetic random access memories (MRAM) bit cells. Moreover, by using eq.(26) the slope of 125
the WER curve in logarithmic scale, namely a measure of the sensitivity of the WER to the current pulse amplitude can be estimated. According to the previous example, for a WER of 10−12 the sensitivity is S WER ∼ −2 ns−1 . This means that an increase of 1 ns in the pulse duration affects the WER by a factor of e−2 , which yields a reduction by almost one order of magnitude. In addition, by using eqs.(26)-(29), one can draw the asymptotic WER diagram (also reported in fig.8) in a very straightforward way. It is apparent that there is remarkable agreement between the
130
analytical predictions and the results of micromagnetic simulations. Finally, we observe that a purely statistical computation of the WER curves depicted in fig.8 would require billions of Langevin simulations for a single value of the current amplitude.
Conclusions In this work, theoretical analysis of the magnetization switching process in PMTJs has been performed. Analytical 135
theory for statistical description of the switching process has been developed under the assumption of macrospin approximation and ballistic magnetization dynamics. The validity of macrospin approximation during the switching dynamics has been checked by means of micromagnetic simulations, where the different contributions to the free energy and the amplitude of the average magnetization have been used as indicators. Then, the LLGS equation has 12
been integrated for an ensemble of N = 1000 replicas of the analyzed device, with initial conditions dispersed around 140
mz = +1, both with full micromagnetic and macrospin solvers. From the dynamics, the statistics of switching times have been inferred for different switching current values. Very good agreement has been found between the switching time distributions obtained from micromagnetic and macrospin simulations. The comparison between analytical predictions and numerically computed switching time distributions has been also performed. It is found that the agreement is very good for switching current larger than twice the critical value. Finally, it is demonstrated that the
145
analytical theory provides ready-to-use tools for analyzing Write-Error-Rates, which can favor immediate comparison with experiments for the optimal design of spin-transfer torque MRAMs.
References [1] J. Akerman, Toward a universal memory, Science 308 (2005) 508–510. doi:10.1126/science.1110549. [2] H. W. Schumacher, C. Chappert, P. Crozat, R. C. Sousa, P. P. Freitas, J. Miltat, J. Ferre, Precessional magnetization reversal in microscopic 150
spin valve cells, IEEE Transactions on Magnetics 38 (5) (2002) 2480–2483. doi:10.1109/TMAG.2002.801903. [3] G. Bertotti, I. D. Mayergoyz, C. Serpico, M. d’Aquino, Geometrical analysis of precessional switching and relaxation in uniformly magnetized bodies, IEEE Transactions on Magnetics 39 (5) (2003) 2501–2503. doi:10.1109/TMAG.2003.816453. [4] S. Kaka, S. E. Russek, Precessional switching of submicrometer spin valves, Applied Physics Letters 80 (16) (2002) 2958–2960. doi: 10.1063/1.1470704.
155
[5] M. d’Aquino, W. Scholz, T. Schrefl, C. Serpico, J. Fidler, Numerical and analytical study of fast precessional switching, Journal of Applied Physics 95 (11) (2004) 7055–7057. doi:10.1063/1.1689910. [6] S. Greaves, Y. Kanai, H. Muraoka, Magnetization switching in energy assisted recording, IEEE Transactions on Magnetics 48 (5) (2012) 1794–1800. doi:10.1109/TMAG.2012.2187776. [7] M. H. Kryder, E. C. Gage, T. W. McDaniel, W. A. Challener, R. E. Rottmayer, G. Ju, Y. Hsia, M. F. Erden, Heat assisted magnetic recording,
160
Proceedings of the IEEE 96 (11) (2008) 1810–1835. doi:10.1109/JPROC.2008.2004315. [8] C. Serpico, S. Perna, G. Bertotti, M. d’Aquino, A. Quercia, I. D. Mayergoyz, Noise-induced bifurcations in magnetization dynamics of uniaxial nanomagnets, Journal of Applied Physics 117 (17) (2015) 17A709. doi:10.1063/1.4906961. [9] J. Zhu, X. Zhu, Y. Tang, Microwave assisted magnetic recording, IEEE Transactions on Magnetics 44 (1) (2008) 125–131. doi:10.1109/ TMAG.2007.911031.
165
[10] E. A. Montoya, S. Perna, Y. J. Chen, J. A. Katine, M. d’Aquino, C. Serpico, I. N. Krivorotov, Magnetization reversal driven by low dimensional chaos in a nanoscale ferromagnet, Nature Communications 10 (1). doi:10.1038/s41467-019-08444-2. [11] K. Watanabe, B. Jinnai, S. Fukami, H. Sato, H. Ohno, Shape anisotropy revisited in single-digit nanometer magnetic tunnel junctions, Nature Communications 9 (1) (2018) 663. doi:10.1038/s41467-018-03003-7. [12] N. Perrissin, S. Lequeux, N. Strelkov, A. Chavent, L. Vila, L. Buda-Prejbeanu, S. Auffret, R. Sousa, I. Prejbeanu, B. Dieny, A highly
170
thermally stable sub-20 nm magnetic random-access memory based on perpendicular shape anisotropy, Nanoscale 10 (25) (2018) 12187– 12195. doi:10.1039/C8NR01365A. [13] D. Apalkov, B. Dieny, J. M. Slaughter, Magnetoresistive random access memory, Proceedings of the IEEE 104 (2016) 1796–1830. doi: 10.1109/JPROC.2016.2590142. [14] G. Bertotti, I. D. Mayergoyz, C. Serpico, Nonlinear Magnetization Dynamics in Nanosystems, Elsevier Science, San Diego, (USA), 2009.
175
[15] C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the natural sciences, Springer–Verlag, Berlin, 1985. [16] C. Serpico, G. Bertotti, I. Mayergoyz, M. d’Aquino, R. Bonin, Thermal stability in spin torque driven magnetization dynamics, Journal of applied physics 99 (8). doi:10.1063/1.2158388.
13
[17] M. d’Aquino, V. Scalera, C. Serpico, Analysis of switching times statistical distributions for perpendicular magnetic memories, Journal of Magnetism and Magnetic Materials 475 (2019) 652–661. doi:10.1016/j.jmmm.2018.11.031. 180
[18] G. Bertotti, I. D. Mayergoyz, M. d’Aquino, S. Perna, C. Serpico, Phase-flow interpretation of magnetization relaxation in nanomagnets, IEEE Transactions on Magnetics 50 (11) (2014) 1–4. doi:10.1109/TMAG.2014.2329035. [19] M. d’Aquino, C. Serpico, G. Coppola, I. D. Mayergoyz, G. Bertotti, Midpoint numerical technique for stochastic landau-lifshitz-gilbert dynamics, Journal of Applied Physics 99 (8) (2006) 08B905. doi:10.1063/1.2169472. [20] C. Ragusa, M. d’Aquino, C. Serpico, B. Xie, M. Repetto, G. Bertotti, P. Ansalone, Full micromagnetic numerical simulations of thermal
185
fluctuations, IEEE Transactions on Magnetics 45 (10) (2009) 3919–3922. doi:10.1109/tmag.2009.2021856. [21] M. d’Aquino, C. Serpico, G. Miano, I. D. Mayergoyz, G. Bertotti, Numerical integration of landau-lifshitz-gilbert equation based on the midpoint rule, Journal of Applied Physics 97 (10) (2005) 10E319. doi:10.1063/1.1858784. [22] V. Tsiantos, W. Scholz, D. Suess, T. Schrefl, J. Fidler, The effect of the cell size in langevin micromagnetic simulations, Journal of magnetism and magnetic materials 242 (2002) 999–1001. doi:10.1016/S0304-8853(01)01365-8.
14