{'fraction_non_alphanumeric': 0.04471565827057907, 'fraction_numerical': 0.017242274733835366, 'mean_word_length': 4.176344086021506, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear bias dependence of spin-transfer torque from atomic first principles

We report first-principles analysis on the bias dependence of spin-transfer torque (STT) in Fe/MgO/Fe magnetic tunnel junctions. The in-plane STT changes from linear to nonlinear dependence as the bias voltage is increased from zero. The angle dependence of STT is symmetric at low bias but asymmetric at high bias. The nonlinear behavior is marked by a threshold point in the STT versus bias curve. The high-bias nonlinear STT is found to be controlled by a resonant transmission channel in the anti-parallel configuration of the magnetic moments. Disorder scattering due to oxygen vacancies in MgO significantly changes the STT threshold bias.


I. INTRODUCTION
Magnetic tunnel junction (MTJ) has attracted great attention due to its importance in magnetic randomaccess memory (MRAM) and read-sensor technology. A basic MTJ is made of two ferromagnetic (FM) layers sandwiching a thin insulating material. Digital information is coded by magnetic moments of the FM layer being in parallel or anti-parallel configurations (PC or APC). In commercial MRAMs, switching between PC and APC is achieved by external magnetic fields. An emerging trend is to switch by spin-transfer torque (STT) predicted theoretically 1 and demonstrated experimentally 2-4 where a spin-polarized current transfers angular momentum to a FM layer causing magnetic moment reversal. STT holds great promise to simplify the MRAM structure, down-scale device size, and reduce power consumption.
One of the most important issues concerning STT is its dependence on external bias voltage V b that drives the spin-polarized current in the first place. Understanding this issue is rather difficult not only because torque is a vector [5][6][7][8][9] but also because STT operates at nonequilibrium due to the flow of spin currents. A physical picture concerning the bias dependence of STT, from the linear to nonlinear bias regime, is yet to be established from atomic first principles. Theoretically, V b dependence of STT was investigated by model analysis such as the tight-binding model 10,11 and free-electron model. [12][13][14] More recently, Heiliger and Stiles 15 analyzed STT for an MgO based MTJ based on equilibrium electronic structure obtained from density-functional theory (DFT). The in-plane STT was calculated 15 up to V b = ±0.5 V and showed essentially a linear dependence on V b , and the out-of-plane STT was found to be quadratic in V b , which is consistent with tight-binding 10 and experimental results. 5,16 Existing calculations provided valuable understanding of STT at the small V b regime. To establish a complete picture, it is important to investigate STT at higher bias from atomic first principles. Indeed, recent experimental data 8 already showed a strong nonlinear V b dependence of in-plane STT, and observable non-quadratic out-ofplane STT beyond V b ≈ 0.2V. It is the purpose of this work to report a first-principles analysis of STT in the most technologically important Fe/MgO/Fe MTJ from low-to high-bias regimes, focusing on its nonlinear V b dependence. We found that the in-plane STT changes from linear to nonlinear dependence as V b is increased from zero. The angle dependence of STT is symmetric at low bias but asymmetric at high bias. The nonlinear behavior is marked by a threshold point in the STT versus bias curve which is controlled by a resonant transmission channel in the APC current. We also found that the nonlinear STT is affected by atomic defects−oxygen vacancies in MgO−hence, in principle, can be tuned by interface engineering. Results are compared to the corresponding experimental data. This paper is organized as follows. In Sec. II, we give the details of our calculation based on Keldysh nonequilibrium Green's function (NEGF) and wave-functionmatching function method. In Sec. III, we present our results on Fe/MgO/Fe(001) MTJs with clean and disordered interfaces. Section IV is our summary.

II. ELECTRONIC STRUCTURE AND TRANSPORT CALCULATION
The MTJ we consider consists of an MgO barrier sandwiched by two semi-infinite Fe leads shown in Fig.1. The system is periodic in the x-z plane and current flows along y corresponding to the (001) material-growth direction. A very small lattice mismatch between Fe and MgO is neglected by fixing the interfacial atoms at their bulk bcc positions. Our calculation is from first principles where DFT is carried out within the Keldysh nonequilibrium Green's function (NEGF) method. 17,18 Here, DFT determines the Hamiltonian of the open device structure, and NEGF determines the nonequilibrium statistics of the device operation and the nonequilibrium density matrix. The NEGF-DFT is solved self-consistently under finite bias during current flow. For devices hav- ing oxygen vacancies (OV, the most energetic favorable defects), we further apply the nonequilibrium-vertexcorrection (NVC) theory 18 for disorder averaging at both the nonequilibrium-density-matrix level and transportcalculation level. The NEGF-DFT-NVC formalism allows one to self-consistently calculate quantum transport properties from atomic first principles without phenomenological parameters. We refer interested readers to the original work. 18 Moreover, The self-consistent nonequilibrium potentials were used as input to a TB-MTO wave-functionmatching calculation. The scattering wave functions of the whole system were obtained explicitly. Because we cannot embody the nonequilibrium NVC into the outof-plane STT calculations at present, a supercell is used to assess out-of-plane STT in the presence of interfacial disorder. In the calculations, an 800 × 800 k -mesh is used to sample the two-dimensional Brillouin zone (BZ) in the x-z plane (see Fig. 1) to ensure accurate convergence. Other numerical details are similar to those of Ref. 19 and 20. To check our electronic structure of Fe/MgO/Fe, we calculate the zero-bias TMR for 7 monolayers (ML) MgO barrier, which is 5100% with the perfect Fe/MgO interface while decreases dramatically to 350% (78%) with 3% (7%) OV at both interfaces. The TMR of a clean junction is consistent with the published first-principles calculation, 21 while the TMR of junctions with vacancies are in the range of experimental measurement. [22][23][24] The results indicate that we have got the right electronic structure.

III. BIAS DEPENDENCE OF STT
Before presenting results, let's consider some general trends. Neglecting any spin dependent scattering in the tunnel barrier, the in-plane STT can be obtained from conservation of spin current, 10 where θ is the relative angle between magnetization M 1 and M 2 (see Fig. 1). I (s) is that of APC. By structural symmetry, I where γ 1 and γ 2 are constants. From these considerations, Eq. (1) suggests I (s) z (π) to dominate at higher V b . Hence, T || should start from a linear dependence and turn into a nonlinear dependence as V b is increased; the turning point defines a "threshold" bias whose value is determined roughly by the material-specific constants γ 1,2 . At low bias, the outof-plane STT demonstrates a simple quadratic relation 10 that we can write as T ⊥ = γ ⊥ V 2 b . After the NEGF-DFT self-consistent calculation is converged, 18 spin currents I ↑ and I ↓ can be calculated by NEGF 18 and the in-plane STT T || obtained from Eq.
(1). The out-of-plane STT T ⊥ (also T || ) can be obtained from a scattering-wave-function approach following Ref. 19. Figure 2 shows plots of the calculated T || and T ⊥ versus V b for a perfect Fe/MgO/Fe MTJ (no disorder) having seven layers of MgO, where magnetization M 1 ⊥ M 2 , i.e., θ = 90 o (see Fig. 1). Seven-layer MgO is used because it is the thickness of experimental device 16 with which we shall compare results (see below). The blue (up-triangle, T || ) and black (solid square, T ⊥ ) curves were calculated from the scattering wave functions. 19 The red curve (solid circle, T || ) was calculated from Eq. (1), which agrees with the blue curve very well at least for lower V b . 25 At lower bias up to ±0.5V, T || is linear in V b and it is in agreement with Ref. 15. As V b increases, T || becomes gradually nonlinear and a clear threshold at V b ∼ 0.65V (e.g., turning point of the curves) is seen, as expected from the general considerations discussed above.
The inset of Fig. 2(a) plots the spin-current density versus V b . Spin-polarized current density I (s) is, as expected, an even function of V b and, by fitting to a quadratic form γ 2 V 2 b , we found that γ 2 ≈ 5µJV −2 m −2 for V b ≤ 0.5 V. The small value of γ 2 makes I (s) z (π) to have a very weak V b dependence, until its sudden increase at V b ∼ 0.6 V, resulting in the threshold of T || at ∼0.65 V.
The out-of-plane STT shows quadratic behavior at the lower bias with γ ⊥ ≈ −13µJV −2 m −2 , and it is very close to the published calculated value of −14µJV −2 m −2 with a six-ML MgO barrier. 15 Quantitatively, T ⊥ < T || for V b < 0.5 V but reaches a similar scale when V b > 0.5 V.
Not only the bias dependence of STT [see Fig. 2(a)] is changed by the sudden increase of APC current I (s) z (π), the angular dependence of STT is also modified. As shown in Fig. 2(b), the angular dependence behaves as sinθ at a small bias V b = −0.4 V. It however deviates from the sine dependence at a higher bias V b = −0.75 V [see Fig. 2(c)]. Experimental evidence of asymmetric angular-dependent STT at higher bias had been already observed before. 16 Can the angular dependence of in-plane STT be understood by the angular dependence of the charge current or TMR? We give the angular dependence of TMR at a small bias V b = −0.4 V [see Fig. 2 Fig.2(e)]. The angular dependence of TMR at low bias follows the cosine function 26 and deviates from the cosine function at high bias. At high bias, the in-plane STT shows an asymmetry and the peak is close to the parallel side, while the angular-resolved TMR decreases always faster than cosine function as the magnetization goes from parallel to antiparallel. Notice here, we have negative TMR at high bias. So, the angular dependence of TMR cannot explain our calculated angular dependence of in-plane STT.

(d)], and a large bias
It is the sudden increase of APC current I (s) z (π) at a higher bias V b = 0.6 V [inset of Fig. 2(a)] that led to the nonlinear bias dependence of STT. To understand the microscopic origin of this sudden increase, we plot spinresolved transmission coefficient versus electron energy E (relative to the Fermi energy E F ) in Fig. 3, calculated at V b = 0.75V. The ↑ or ↓ denotes the majority or minority spin channel in the left Fe lead. For APC, the ↑ channel of the left lead transmits to the ↓ channel of the right lead through the MgO. Figure .3 shows that there is an abrupt increase in ↑ channel transmission at about E = 0.3 eV. For a symmetric MTJ, roughly half V b is dropped at each Fe/MgO interface. Hence the abrupt increase of transmission at E ≈ 0.3 eV should contribute to APC current when the bias V b = 0.6 V. This is indeed what is found [inset of Fig .2(a)]. We conclude that the sudden increase of I (s) z (π) is due to the abrupt increase of the ↑ channel transmission, which, importantly, is also due to bias-dependent transport features as we explain now.
The inset (a) of Fig.3 plots the transmission hot spot of APC(↑) in the two-dimentional (2D) BZ, which is the transverse momentum (k || ) resolved transmission coefficient. There is a clear bright ring surrounding the Γ point. The abrupt increase of the APC ↑ channel transmission is due to this bright ring. Focusing on one particular k point on this ring and scanning the electron energy, a single resonant peak is discovered in the transmission as shown in the inset (b) of Fig. 3.
Importantly, this sharp resonance is different from the well-known resonance states on Fe surfaces. 27 In particular, we found that this sharp APC resonance is only established at high bias and is located between the ∆ 1 band in the left Fe and the ∆ 5 , ∆ 2′ bands in the right Fe. The contribution of ∆ 5 and ∆ 2′ bands to the total transmission is k || dependent, giving rise to the ring in the hot spot. A higher bias significantly changes the shape of the tunnel barrier: we checked that by artificially flattening the barrier potential from the biased slop shape, this sharp APC resonant channel disappears.
Actually, the existence of resonance is not only responsible for the sudden increase of STT but also responsible for the deviation of angular-dependent STT as shown in the inset (c) of Fig. 3. When at a resonant ring such as A point (k x = k y = 0.278Å −1 ), the angular-dependent in-plane STT deviates from a sine function and follows similar behavior as that in Fig. 2(c). When far from the resonant ring such as B point (k x = k y = −1.624 A −1 ), the angular dependent in-plane STT follows a sine function as well. Under most situations, at high bias, the resonance dominates the in-plane STT. Without resonant states, the Bloch wave should decay very fast in the tunnel barrier, the multiple scattering should be negligible, and Eq. 1 should hold in this case. However, in the presence of resonant states, the wave function of these states will decay little in the barrier and Eq. 1 is not valid anymore.
Having understood the microscopic physics behind the nonlinear bias dependence of STT, we note that the ob-tained threshold bias at V b ∼ 0.6 V for the ideal junction is higher than that reported experimentally, 6 which was about 0.3 V. Since experimental devices are never ideal, an investigation of defect scattering and its effect on STT is warranted. There are several possible defects such as oxygen vacancy in MgO, interfacial roughness, 15 material imperfections, etc. It was reported experimentally that due to compressive strain during crystal growth, oxygen vacancies (OV) are accommodated in MgO. 29 Both theory 20 and experiment 30 have shown that existence of OV can drastically influence spin-polarized transport and its effect on STT is therefore expected.
We calculated the in-plane STT from Eq. (1) by putting OVs in the interfacial MgO layer immediately adjacent to Fe. An alloy model O 1−x Va x is used, where Va stands for vacancy and x is the OV concentration. The disorder average is carried out by the NVC theory 18 mentioned above. Results are plotted in Fig. 4 versus V b , and several observations are in order.
First, at low V b , T || is linear for both disordered (x = 0) and clean samples (x = 0)-hence, the linear regime is not qualitatively altered by OVs. Second, both symmetric junctions where the same OV concentration is on both Fe/MgO interfaces (black curve with solid squares), and asymmetric junctions where OVs only exist at one Fe/MgO interface (red curve with solid circles), yield almost the same T || . In particular, its slope τ || ≡ dT || /dV ≈ −5.5 × 10 13 τ o . This value is only slightly larger than that of the ideal junction, which has τ || = −4.7 × 10 13 τ o . Hence, a small amount of OV (e.g. 3%) does not alter T || significantly. These results should be compared to the experimental value of the slope τ || = −3.2 × 10 13 τ o measured at V b = ±0.3 V. 16 Giving the uncertainties in comparing our atomic structure to the experimental one, the quantitative consistency of our results is quite satisfactory. Third, in Ref. 20, it was shown that filling the OVs with nitrogen atoms can significantly reduce disorder-induced diffusive scattering and drastically increase the TMR ratio. Here, filling the OVs with N using the alloy model O 1−x N x , we found that the T || is somewhat decreased, namely, its slope τ || = −1.9 × 10 13 τ o (blue curve with up-triangles).
Similar to the ideal junction, systems with OV also show abrupt increase of T || at higher bias, which, we found, is also due to the quantum resonance in the 2D BZ discussed above. Different from the ideal case, existence of OVs shifts the threshold bias of T || to lower values thus bringing it closer to the experimental result. 6 As shown in the inset of Fig. 4, the effect of OV on T || can also be understood by investigating transmission of the ↑ channel in APC.
The OVs induce two new effects. First, interfacial OVs broaden the area having high transmission values in the hot spots and also suppress the transmission peak from 1.0 e 2 /h to about 0.35 e 2 /h (at 3% OV on both interfaces). The overall effect is to enhance T || slightly. Second, interfacial OVs shift the resonant peak in the ↑ channel of APC to lower energy and, as a result, the threshold bias is now shifted to lower values.
Remarkably, from Fig. 4 one can determine which Fe/MgO interface dominates STT. When OVs exist at both interfaces (black curve with solid squares), the threshold bias shifts to lower values from those of the clean sample for both V b polarity. When OVs exist only at the right Fe/MgO interface (red curve with solid circles), the threshold bias shifts to lower value only in the negative V b . When OVs are filled by N atoms (blue curve with up-triangles), the threshold bias shifts to higher V b for both polarities.
Quantitatively, the threshold bias is ∼0.45V in the presence of 3% OV at interface. It is consistent with the experimental value ∼0.4 V (see Ref. 16 ) with a 1.25-nm (6-ML) MgO barrier. A more recent experiment measured a threshold bias ∼0.2 V (see Ref. 8 ) with a 1.0-nm (5-ML) MgO barrier. The experimentally measured resistance is 1.5 Ωµm 2 also close to our calculated value for a junction with a 5-ML MgO barrier.
The effects of OVs on the out-of-plane STT are also very important. We have examined a junction using the wave-function-matching method. The calculation was done in a 7 × 7 lateral supercell containing five randomly distributed OVs at each Fe-MgO interface (thus x ∼ 10%). For V b = 0.5 V and E = E F +0.25 eV, the calculated in-plane and out-of-plane STT are 20.30 × 10 14 τ o and 7.32 × 10 14 τ o , respectively. These values are larger than that of the ideal junction which are 1.06 × 10 14 τ o and 0.94 × 10 14 τ o for in-plane and out-of-plane STT, respectively. The larger in-plane and out-of-plane STT in the 10% OV sample is attributed to the resonant tun-neling discussed above. These results indicate that the interfacial OV disorder is not detrimental to the value of the out-of-plane STT. In comparison, interface roughness disorder quenches the out-of-plane STT as reported in Ref. 15.

IV. SUMMARY
In summary, the bias dependence of STT was investigated from atomic first principles, covering both lowand high-bias regimes. STT shows weak and linear bias dependence at small V b but strong and nonlinear dependence at higher V b . The angle dependence of STT is symmetric at low bias but asymmetric at high bias. The nonlinear behavior is marked by a threshold bias in the STT versus bias curve and is controlled by a resonant transmission channel in APC. Very importantly, disorder scattering by small amount of oxygen vacancies in MgO lowers the STT threshold bias but not the STT value, suggesting that the nonlinear STT can be tuned by clever interfacial engineering as also reported in Ref. 30 where crystal-growth temperature changes the OV contents.