Non-equilibrium Green's Function and First Principle Approach to Modeling of Multiferroic Tunnel Junctions

TL;DR

This paper presents a comprehensive multi-physics modeling framework for multiferroic tunnel junctions using NEGF, Landau-Khalatnikov, and LLG equations, validated by experimental data and analyzing device parameter effects.

Contribution

It introduces a novel integrated simulation approach combining quantum transport, ferroelectric, and magnetic dynamics for MFTJs, filling a gap in multi-physics studies.

Findings

Simulation results agree with experimental data.

Analysis of TER and TMR dependence on device parameters.

First-principles estimates of electrode screening lengths.

Abstract

Recently, multiferroic tunnel junctions (MFTJs) have gained significant spotlight in the literature due to its high tunneling electro-resistance together with its non-volatility. In order to analyze such devices and to have insightful understanding of its characteristics, there is a need for developing a multi-physics modeling and simulation framework. The simulation framework discussed in this paper is motivated by the scarcity of such multi-physics studies in the literature. In this study, a theoretical analysis of MFTJs is demonstrated using self-consistent analysis of spin-based non-equilibrium Green's function (NEGF) method to estimate the tunneling current, Landau-Khalatnikov (LK) equation to model the ferroelectric polarization dynamics, together with landau-Lifshitz-Gilbert's (LLG) equations to capture the magnetization dynamics. The spin-based NEGF method is equipped with a magnetization dependent Hamiltonian that eases the modeling of the tunneling electro-resistance (TER), tunneling magneto-resistance (TMR), and the magnetoelectric effect (ME) in MFTJs. Moreover, we apply the first principle calculations to estimate the screening lengths of the MFTJ electrodes that are necessary for estimation of tunneling current. The simulation results of the proposed framework are in good agreement with the experimental results. Finally, a comprehensive analysis of TER and TMR of MFTJs and their dependence on various device parameters is illustrated.

Equations

$$\displaystyle\phi(z)=\left\{\begin{array}[]{cc}\frac{\sigma_{s}\delta_{1}}{\epsilon_{0}\epsilon_{r1}}e^{-\frac{|z|}{\delta_{1}}}&\quad\quad z\leqslant 0\\ \\ -\frac{\sigma_{s}\delta_{2}}{\epsilon_{0}\epsilon_{r2}}e^{-\frac{|z-d|}{\delta_{2}}}&\quad\quad z>0,\\ \end{array}\right.$$

$$E_{FE}=\frac{({\sigma }_s-P)}{{ϵ}_0},$$

$$\frac{{\sigma }_s{\delta }_1}{{ϵ}_0{ϵ}_{r1}}+\frac{{\sigma }_s{\delta }_1}{{ϵ}_0{ϵ}_{r2}}+E_{FE}{t}_{FE}=0.$$

$${\sigma }_s=\frac{P{t}_{FE}}{\frac{{\delta }_1}{{ϵ}_{r1}}+\frac{{\delta }_2}{{ϵ}_{r2}}+t_{FE}}.$$

$$\frac{\partial m_i}{\partial \tau }$$

$$H_{eff,i}=H_{ext}+H_{therm}+H_{anis}+H_{exch,i},$$

$$F=$$

$$\lambda \frac{\partial P}{\partial t}=-\frac{\partial F}{\partial P},$$

$$H_{L(R)}$$

$$H_{ch}$$

$$G={\left[EI-H-{\Sigma }_L-{\Sigma }_R\right]}^{-1},$$

$${\Sigma }_L=\left[\begin{array}{ccccccc}-t_L{e}^{i{K}_L^{\uparrow }a}& 0& 0& .& .& 0& \\ 0& -t_L{e}^{i{K}_L^{\downarrow }a}& 0& .& .& .& \\ 0& 0& .& .& .& .& \\ .& .& .& .& .& .& \\ .& .& .& .& .& .& \\ 0& .& .& .& .& 0& \\ & & & & & & \end{array}\right],$$

$${\Sigma }_R=\left[\begin{array}{ccccccc}0& .& .& .& .& 0& \\ .& .& .& .& .& .& \\ .& .& .& .& .& .& \\ .& .& .& .& 0& .& \\ .& .& .& 0& -t_R{e}^{i{K}_R^{\uparrow }a}& 0& \\ 0& .& .& .& 0& -t_R{e}^{i{K}_R^{\downarrow }a}& \end{array}\right],$$

$$K_{L(R)}^{\uparrow }=co{s}^{-1}\left(1-\frac{E\pm \frac{q{V}_a}{2}-\frac{h^2{K}_t^2}{2m_{L(R)}^{\ast }}+\frac{{\Delta }_{L(R)}}{2}}{2t_{L(R)}}\right),$$

$$K_{L(R)}^{\downarrow }=co{s}^{-1}\left(1-\frac{E\pm \frac{q{V}_a}{2}-\left(\frac{h^2{K}_t^2}{2m_{L(R)}^{\ast }}\right)-\frac{{\Delta }_{L(R)}}{2}}{2t_{L(R)}}\right).$$

$$J=\frac{-e}{2{\pi }^{2}h}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }d{k}_{x}d{k}_y\int dEt\left(f_L-f_R\right),$$

$$t=Trace\left({\Gamma }_{L}G{\Gamma }_R{G}^{†}\right),$$

$${\Gamma }_{L(R)}=i\left({\Sigma }_{L(R)}-{\Sigma }_{L(R)}^{†}\right).$$

$$f_{L(R)}=\frac{1}{e^{(E-{\mu }_{L(R)})/K_{B}T}+1}.$$

$$J^{spin}=$$

$$G^n=G({\Gamma }_L{f}_L+{\Gamma }_R{f}_R)G^{†}.$$

$$V(z,t)$$

$$V_0$$

$$a(t)$$

$$J_{exch,i,j}^{\prime }=|J_{exch,i,j}|\left(1-|a(t){|}^2\right)-...$$

$$|J_{exch,i,j}||a(t){|}^2,$$

Full text

††preprint: APS/123-QED

I Introduction

Over the last few decades, the complementary metal-oxide-semiconductor (CMOS) technology has been continuously downscaled following Moore’s law Moore (1998). However, the static power dissipation and the threshold voltage variations of downscaled short channel transistor have become dominating factors that limit the static random access memory (SRAM) performance Stt-ram et al. (2013); Ye et al. (2010); Chun et al. (2013). Consequently, the high static power dissipation of SRAM inspired the exploration of alternative memory technologies like spin transfer torque magnetic memory (STT-MRAM). However, the limited tunneling magneto-resistance (TMR) of magnetic tunnel junction (MTJ) together with the threshold voltage fluctuations of the short channel access transistor affect the STT-MRAM read error rate. Therefore, the read performance of STT-MRAM has become a fundamental limiting factor in its applicability. Consequently, a new family of tunnel junctions, called ferroelectric tunnel junctions (FTJs), have emerged in literature Chanthbouala et al. (2012a); Garcia and Bibes (2014); Yin and Li (2017).

An FTJ consists of a ferroelectric insulator sandwiched between two different metal electrodes, as illustrated in Fig. 1. The information is stored in the electric polarization of the insulator. The FTJ resistance is a function of the electric polarization of the insulator. The electric polarization of the ferroelectric insulator modulates the FTJ resistance, and hence the information can be extracted by sensing the FTJ resistance. The tunneling electro-resistance (TER) of FTJ is defined as $TER=\frac{|R_{\rightarrow}-R_{\leftarrow}|}{min(R_{\rightarrow},R_{\leftarrow})}$, where $R_{\rightarrow}$ and $R_{\leftarrow}$ are the resistance of positive and negative electric polarization states, respectively Esaki et al. (1971). The physical origin of TER is discussed in detail in section II. The charge current of FTJ consists of three main components: Fowler–Nordheim tunneling, direct tunneling, and thermionic emission Pantel and Alexe (2010).

On the other hand, the multiferroic tunnel junctions (MFTJ) is a nonvolatile tunnel junction that consists of two ferromagnetic layers separated by a ferroelectric insulator, as illustrated in Fig. 1. Intuitively, from the structure of an MFTJ, we can predict that an MFTJ combines the resistive switching mechanism of FTJ and MTJ to constitute a four-state device. However, it turns out that the MFTJ has more advantages over its constituent devices due to the magnetoelectric effect (ME). The ME effect at the FM/FE interface is observed in LaSrMnO3(LSMO)/LaCaMnO3(LCMO)/BaTiO3 (BTO)/LSMO MFTJ Pantel and Alexe (2010). It originates from the modulation of the screening charges at the LCMO side by the bound charges at the BTO interface. The change in the electron concentration at the LCMO interface affects the LCMO magnetic configuration. The magnetic alignment of the LCMO layer is switched from the ferromagnetic (FM) to the antiferromagnetic (AFM) alignment due to the change in electron concentration Burton and Tsymbal (2009); Yin et al. (2015). The transition to the AFM alignment shifts the density of states (DOS) of the majority spin carriers to higher energy levels, and hence limits the majority spin current. In brief, the overall influence of the ME effect is to improve the TER ratio, as explained in detail in section II and section III.

A detailed review of the state of the art in FTJs and MFTJs could be found in Garcia and Bibes (2014); Yin and Li (2017). However, a brief review of the progress in FTJs and MFTJs literature is provided in the following discussion. Although the FTJ has been predicted by Esaki et al. Esaki et al. (1971) in 1971 with the name ”polar switch”, the FTJ has not been realized until recently. The lack of the knowledge of fabrication techniques of ferroelectric ultra-thin films had prevented the FTJ realization. However, due to the breakthrough that has been achieved by Zembilgotov et al. Zembilgotov et al. (2002), the ferroelectric ultra-thin film has been realized followed by many other experimental studies Qu et al. (2003); Rodrı́guez Contreras et al. (2003). Zhuravlev et al. Zhuravlev et al. (2005, 2009) have explained the dependence of the barrier height on electric polarization with the help of Thomas-Fermi equation and Landau tunneling current formula Appelbaum (1970). The Wenzel-Kramer-Brillouin (WKB) approximation and one-band model Appelbaum (1970); Kane (1969); Wolf (1985) have been used to calculate the tunneling current through the FE insulator in Kohlstedt et al. (2005). Hinsche et al. have used Landauer-Büttiker formula and the WKB approximation together with ab initio calculation to model the FTJ characteristics Hinsche et al. (2010). Fechner et al. have used the ab initio method to study the electric polarization dependent phase transition in Fe/ATiO3 interface Fechner et al. (2008). On the other hand, the non-equilibrium Green’s function (NEGF) method along with Landau-Khalatnikov (LK) equation have been used to estimate the FTJ I-V characteristics in Chang et al. (2017). However, the study did not consider the magnetization dynamics or the Hamiltonian dependence on the magnetization.

To conclude, the scarcity of multi-physics simulation studies that capture the MFTJ magnetization dynamics, along with TMR, and TER effects motivates the modeling and simulation framework applied in this study. In this study, the spin-based NEGF is applied to model the tunneling current, along with landau-Lifshitz-Gilbert’s (LLG) equation are applied to model the magnetization dynamics, and the LK equation is applied to describe the FE motion. However, the accuracy of these models depends on the parameters used to model various materials. In our study, we use the density functional theory (DFT) to estimate the electrostatic potential, and hence the screening lengths of the electrode that are used in the NEGF transport simulations. The simulation results are compared to experimental results of the MFTJ in Yin et al. (2015); Chanthbouala et al. (2012b) to confirm the validity of the method.

The quantum transport model adopted by this study is based on the mean field approximation. In addition, the proposed model is based on single-band effective mass approximation of the complex band structure of the material. A detailed discussion of the advantages and limitations of the adopted quantum transport model could be found in Datta (2012). However, the effective mass approximation is a computationally efficient method compared to other computationally intensive methods that account for the complex band structure of the material. The self-consistent solution of the magnetization dynamics, electric polarization, and quantum transport requires thousands of evaluations of the quantum transport model.

The paper is organized as follows. Section II is dedicated to explaining the difference between MTJ, FTJ, and MFTJ along with the origin of TER and TMR effects. Section III is devoted to DFT simulations of the ME effect at the FM/FE interface. The magnetization dynamics, ferroelectric dynamics, and quantum transport (NEGF) are illustrated in sections IV, V, and VI, respectively. Section VII is dedicated to the formulation of the magnetic exchange coefficient as a function of the electric polarization based on time-dependent perturbation theory. The simulation procedure is explained in section VIII. Finally, section IX is assigned to the simulation results followed by conclusions in section X.

II Multiferroic Tunnel Junctions

The MFTJ structure combines the FM electrodes of MTJ together with the FE insulator of FTJ to produce a four-state device. We start by explaining the TMR effect in the MTJ along with the TER effect in FTJ before describing the MFTJ characteristics. Furthermore, the ME effect at the FM/FE interface is a unique property of MFTJs that enhances the TER effect, as explained in detail in this section and in section III.

The TMR effect could be explained in the light of spin dependent transport illustrated in Fig. 2 Julliere (1975); Mathon and Umerski (2001); Butler et al. (2001). In such FM materials, the lower band of the density of states (DOS) of the majority and minority spin carriers have energy shift, as illustrated in Fig. 2 (a). The energy splitting is dependent on the magnetization direction. Therefore, in the case of anti-parallel alignment of the electrode magnetization, the majority spin carriers that migrate from the left electrode are restricted by the shortage of matched spin states at the right electrode. Consequently, the overall charge current is reduced in the case of anti-parallel alignment of the electrode magnetization. In the case of parallel alignment of the magnetization, the majority and minority spin carriers migrate from the left electrode and are absorbed by the matched spin states that are sufficiently available at the right electrode. Consequently, the overall charge current is not limited by the availability of the spin states in the case of parallel alignment of magnetization. In other words, the MTJ resistance changes according to the magnetization alignment of the electrodes that control the DOS energy splitting between the majority and the minority spin carriers. Finally, the TMR is defined as $TMR=\frac{R_{AP}-R_{P}}{R_{P}}$, where $R_{AP}$ is the resistance of anti-parallel aligned magnetization state, and $R_{P}$ is the resistance of parallel aligned magnetization state Julliere (1975).

The TER effect of FTJ could be explained by the help of charge screening phenomena in the metal electrode Zhuravlev et al. (2005). The electric polarization of FE insulator induces bound charges at the metal/FE interface. The bound charges are partially screened by the free electron gas in the metal side, as illustrated in Fig. 2 (b). The uncompensated charges at the interface result in a constant electric field inside the insulator, and hence linear potential. Moreover, the polarization direction controls the polarity of the bound charge, and hence the polarity of the potential drop in the FE insulator $(\phi_{1}-\phi_{2})$, where $\phi_{1}$ and $\phi_{2}$ are the potential at the left and right metal/FE interfaces, respectively. Therefore, the barrier height increases by $|\phi_{1}+\phi_{2}|$ in the case of positive electric polarization. In contrast, the barrier height is reduced by $|\phi_{1}+\phi_{2}|$ in the case of negative electric polarization. Finally, the large TER value of FTJ is a natural result of the exponential dependence of tunneling current on the barrier height. The asymmetry of the electrodes screening lengths is an important factor for MFTJs to exhibit a nonzero TER. It is important to mention that the TER also depends on the barrier effective thickness, which can change upon the polarization reversal due to the change from the insulating to metal phase, and both interface terminations Quindeau et al. (2015); Borisov et al. (2015).

The aforementioned qualitative discussion could be formulated quantitatively, as illustrated in Zhuravlev et al. (2005). The screening charges and potential distribution in the metal side is described by the Thomas-Fermi formalism:

[TABLE]

where $\phi$ is the electrostatic potential, $\sigma_{s}$ is the surface charge density of free charges, $\epsilon_{0}$ is permittivity of free space, $\epsilon_{r1}$ ($\epsilon_{r2}$) is the permittivity of the first (second) electrode, and $\delta_{1}$ ($\delta_{2}$) is the screening length of first (second) electrode. According to Thomas-Fermi relation, the charge and the potential of any point in the electrodes decrease as an exponential function of the distance between the point and the interface. Moreover, the potential values at the interface are defined as $\phi_{1}=\frac{\sigma_{s}\delta_{1}}{\epsilon_{0}\epsilon_{r2}}$ and $\phi_{2}=\frac{-\sigma_{s}\delta_{2}}{\epsilon_{0}\epsilon_{r2}}$. By applying the Gauss’s law at the metal/FE interface, we get the expression

[TABLE]

where $P$ is the polarization vector and $E_{FE}$ is the electric field in the FE layer. The potential drop $\phi_{1}-\phi_{2}$ is equal to the constant electric field inside the FE insulator multiplied by $t_{FE}$ as given by

[TABLE]

Finally, from (5) and (6), the $\sigma_{s}$ that satisfies the continuity of the potential at the interface is defined as

[TABLE]

However, for the limiting case of $t_{FE}\gg\frac{\delta_{1}}{\epsilon_{r1}}+\frac{\delta_{2}}{\epsilon_{r2}}$, the free charge density $\sigma_{s}$ is equal to $P$, and hence the potential drop is zero which eliminates the TER effect Zhuravlev et al. (2005). Therefore, a mandatory constraint is imposed on the maximum FE thickness that maintains the TER effect. However, the stability of the ferroelectricity imposes a lower limit on the FE thickness. Therefore, the FE layer should have an optimal thickness that maintains the ferroelectricity and provides high TER ratio at the same time.

An MFTJ combines the TER effect of FTJ along with the TMR effect of MTJ to produce a four-state device as illustrated in Fig. 2 (c). However, it has been experimentally observed that the LCMO electrode of LSMO/LCMO/BTO/LSMO MFTJ Yin et al. (2015) goes through phase transition from the FM state to the AFM phase by the influence of the electric polarization switching. To understand the effect of the FM to the AFM phase transition on the TER, let us assume that both electrodes have the magnetization in the positive $z$ direction. In the case of positive electric polarization (high resistance state), the LCMO left electrode has an AFM configuration. The lower band of the DOS of the spin-up carriers shifts to higher energy levels. Therefore, the spin-up carriers that migrate from the right electrode are restricted by the shortage of spin-up states at the left electrode. Thus, the MFTJ high resistance increases, and hence the TER increases. In case of negative electric polarization (low resistance state), the LCMO has an FM configuration. The spin-up carriers that migrate from the right electrode are absorbed by the matched spin states without any restriction. Consequently, the overall TER of the MFTJ improves. The details and origin of this ME effect are explained in the following section.

III The First Principle Calculations of the FM/FE interface

Recently, many different forms of magnetoelectric effects have been observed in the literature such as electric field manipulated magnetization, electric field induced magnetic phase transition, and voltage controlled magnetic anisotropy Weisheit et al. (2007); Papers (2009); Martin and Ramesh (2012). In this study, we focus on the magnetoelectric effect that happens in the interface between $La_{1-x}A_{x}MnO_{3}/BaTiO_{3}$, where A is a divalent cation, i.e., Ca, Ba, and Sr and $x$ is the chemical doping concentration. $LAMO$’s phase diagram exhibits a phase change between the ferromagnetic state and antiferromagnetic state as a function of hole carrier concentration $x$ Dagotto et al. (2001). The transition between the FM and AFM phases and its dependence on hole concentration could be explained by the existence of two competing interactions that happen between the adjacent Mn sites in LAMO: superexchange interaction and double exchange interaction. In contrast to superexchange interaction that prefers AFM alignments, the double exchange interaction favors FM alignment De Gennes (1960). The doping concentration $x,$ modulates the density of electrons in Mn $e_{g}$ orbitals that mediate the double exchange interaction. Note, the doping concentration supports one of the interactions over the other, and hence favors one of the configurations over the other.

The electrostatic doping created by electric polarization could change the electron concentrations similar to chemical doping Ahn et al. (2003). Since the bound charges induced by electric polarization of BTO at the interface modulates the screening charges at the LAMO side, the electric polarization could control the magnetization phase transition similar to chemical doping. The magnetoelectric effect in LAMO/BTO interface is illustrated in Fig. 3. The second Mn site in LAMO exhibits AFM (FM) alignment in case of positive (negative) polarization state. However, the chemical doping concentration has to be fixed at the magnetic phase transition point ($x=0.5$) between the FM and the AFM phases to facilitate the magnetic phase transition by electrostatic doping Moritomo et al. (1998).

III.1 Simulation Procedure and Parameters

We applied DFT method to extract the electrostatic potential profile of LAMO/BTO and Co/BTO structures. The screening lengths of LAMO and Co electrodes are estimated from the electrostatic potential. The generalized gradient approximation (GGA) method Perdew et al. (1996) implemented in Quantum-ESPRESSO package Giannozzi et al. (2009) is used to perform all of the DFT calculations in this study. The Vanderbilt’s ultrasoft pseudopotential Vanderbilt (1990) is used along with virtual crystal approximation (VCA) Colombo et al. (1991) to handle the La-A doping. The VCA method is used by Burton et al. Burton and Tsymbal (2009) to perform DFT calculations for typical structure with acceptable accuracy. The energy cutoff of 400 eV and Monkhorst-Pack grid of 12x12x1 of k-points are used for all the DFT simulations in this study.

The supercell used to simulate LAMO/BTO interface consists of 4.5 unit cells of LAMO and 5.5 unit cells of BTO, as illustrated in Fig. 3. The structure is stacked along (001) direction of the perovskite cell. The stacking sequence at LAMO/BTO interface is $AO-BO_{2}$ Burton and Tsymbal (2009). The supercell illustrated in Fig. 4 is utilized to model the Co/BTO interface. The structure consists of 4.5 unit cells of Co and 5.5 unit cells BTO along (001) direction and is rotated $45^{\circ}$ in the $x-y$ plane. The most stable interface has TiO2 termination as described in Cao et al. (2011). We did not include any vacuum regions in these structures. As both structures are epitaxial growth on a SrTiO3 substrate that has a bulk in-plane lattice constant of $a=3.94\AA$, the lattice constant in the lateral direction is constrained to $a=3.94\AA$ for all the layers of LAMO, Co and BTO. The lateral strain results in tetragonal distortion in the longitudinal direction ($z$ direction). To estimate the tetragonal distortion, the DFT calculation of a single LAMO unit cell is repeated with different longitudinal lattice constants $c$. The lattice constant that has the minimum total energy is selected for further supercell simulations. The longitudinal lattice constant of LAMO that has minimum total energy is $c/a=0.99$. Similarly, the longitudinal lattice constant of BTO and Co are estimated to be $c/a=1.05$ and $c/a=0.83$, respectively. Next, both supercells of LAMO/BTO and Co/BTO with the in-plane constraint and the corresponding tetragonal distortion are relaxed until the total force on the atoms is less than $10^{-3}$ $Ryd/au$. As the total force on the atoms reaches the limit of $10^{-3}$ $Ryd/au$, the atoms reach their equilibrium positions. Further optimization beyond this limit results in a negligible change in the positions of the atoms.

III.2 DFT Simulation Results

The magnetic configurations that minimize the total energy of the LAMO/BTO interface are illustrated in Fig. 3. The second Mn site in LAMO (left interface) exhibits AFM (FM) alignment in case of positive (negative) polarization state. The magnetic configuration of LAMO (left interface), that minimize the total energy, is AFM in the case of positive polarization and FM in the case of negative polarization Burton and Tsymbal (2009). The DFT simulation results for the magnetization of Mn sites are 2.57, -2.66, 2.76, 2.87, and 3.02 $\mu_{B}$. In case of positive polarization state, the magnitude of the magnetization of the Mn atoms is lower at left interface and increases for the atoms away from that left interface Burton and Tsymbal (2009). In contrast, the magnetic configuration of the Co/BTO interface, that minimizes the total energy, is FM configuration independent of electric polarization of BTO. The Co/BTO interface exhibits a constant magnetic configuration independent of the electric polarization of BTO Cao et al. (2011), as illustrated in Fig. 4. The DFT simulation result for the magnetization of bulk Co is $1.74$ $\mu_{B}$.

The atomistic electrostatic potential of LAMO/BTO and Co/BTO are illustrated in Fig. 5 and Fig. 6, respectively. The macroscopic potential is estimated from the atomistic potential by a moving window integral Colombo et al. (1991) and fitted by a spline function. Finally, the screening lengths $\delta/\epsilon_{0}\epsilon_{r}$ of La0.7A0.3MO, L0.5A0.5MO, and Co are estimated from electrostatic potential to be 1.06, 1.05, and 1.14 $m^{2}/F$, respectively. The estimated screening lengths are used in the spin dependent transport calculations, as explained in section VI.

IV The Magnetization Dynamics

The Landau-Lifshitz-Gilbert (LLG) equation formulates the precessional and damped motion of magnetization induced by the magnetic field and spin current Gilbert (2004, 2004); D’Aquino (2004). The single domain LLG equation is used along with NEGF self-consistently in many studies in literature. The single domain solution of LLG equation is not appropriate for LAMO Because the LAMO material has atoms in the FM order and other atoms in the AFM order at the same time. Similarly, the solution of the magnetization as a continuum fails as well, because it requires a second order derivative of the magnetization with respect to space. The second order derivative appears in the definition of the exchange interaction effective field $H_{exch}=\frac{2A}{\mu_{0}M_{S}}\frac{\partial^{2}m}{\partial x^{2}}$, where $A$ is a constant. However, the abrupt change of the magnetization at that atomistic scale makes the derivative with respect to space is not possible. The usual solution in case of AFM material is to replace the magnetization by the total magnetization $l=m_{j+1}+m_{j}$ and the AFM Neel field $n=m_{j+1}-m_{j}$ that are continuous variables in case of AFM material. However, LAMO has the FM and AFM orders that exist at the same time. The Neel field will be discontinuous at the area between the FM and the AFM phase.

We adopted a discrete multi-domain version of the LLG equation. The main difference between the continuum and discrete LLG equation is the definition of the exchange field. The definition of the exchange field in the discrete LLG equation does not require differentiation with respect to spatial coordinates. The discrete multi-domain LLG equation is similar to the atomistic LLG equation Evans et al. (2014). However, the discrete multi-domain LLG equation models the lateral direction as a single domain to reduce the computational effort. The lateral single domain assumption does not affect the accuracy of the method because the cross-section area of the junction is large enough to neglect the effect of edges. The discrete domains have a thickness equal to a single unit cell in the normal direction.

IV.1 The LLG Equation

The LLG equation Evans et al. (2014) can be expressed as

[TABLE]

where $m$ is a unit vector in direction of magnetization, $\tau$ is defined as $\tau=\frac{|\gamma|}{(1+{\alpha}^{2})}dt$, t is the time, $H_{eff}$ is the effective magnetic field, $\alpha$ is the Gilbert damping constant, $\gamma$ is the gyromagnetic ratio,i is index over the atoms along the $x$ axis, and STT is the spin transfer torque. The first term of (8) is the precessional motion of the magnetization due to the effective magnetic field. The second term models the damped part of magnetization oscillation. The third term is the spin torque exerted by the spin current on the magnetization.

LAMO is modeled as a 1D chain of discrete domains with a magnetization variable $m_{i}$ assigned to each domain. Each mesh cell has a length equal to the lattice constant and cross section area equals the total cross-section area of the MFTJ. In other words, we assumed that the MFTJ cross-section area is large enough. Therefore, we can neglect the effect of the boundary cell on the magnetization dynamics.

IV.2 The Effective Magnetic Field

The effective magnetic field is given by

[TABLE]

where $H_{ext}$ is the external magnetic field, $H_{exch,i}$ is the exchange interaction effective field, $H_{therm}$ models the random thermal variations, and $H_{anis}$ is the magnetic anisotropy. The magnetic anisotropy is defined as $H_{anis}=\frac{2K_{U}}{M_{s}}$, where $K_{U}$ is the anisotropy constant, and $M_{S}$ is the saturation magnetization.

The exchange interaction field is defined as $H_{exch,i}=\frac{1}{\mu_{0}M_{S}}\sum_{j}^{N}J^{\prime}_{exch,i,j}m_{j}$ Evans et al. (2014), where $J^{\prime}_{exch,i,j}$ is the material magnetic exchange coefficient that is averaged by the FM to AFM transition probability as explained in VII, and N is the number of nearest neighbors.

IV.3 The Thermal Fluctuations

The random variation in the magnetization due to thermal excitations is modeled by the effective magnetic field $H_{therm}$ defined as $H_{therm}=\zeta\sqrt[]{\frac{2\alpha KT}{|\gamma|M_{S}V_{cell}dt}}$, where $K$ is Boltzmann’s constant, $T$ is the temperature, $V_{cell}$ is the mesh cell volume that equals to the lattice constant multiplied by the total cross-section area of the MFTJ, $dt$ is the numerical time step, and $\zeta$ is a vector with random components which are selected from standard normal distribution Brown (1963).

The thermal fluctuations term makes the LLG equation stochastic differential equation (SDE). The integration of the thermal fluctuations results in Wiener stochastic process Berkov (2007); Gardiner (1986) that is not differentiable with respect to time. Therefore, the Stieltjes integral is used instead of the Riemann integral to integrate the thermal term Berkov (2007); Gardiner (1986). The Stieltjes integral of the thermal term is defined in terms of the differential increments of Wiener process that has a variance proportional to the integration time step. Therefore, time-step $dt$ appears in the denominator of the thermal field. The details of the integration of the LLG equation as a stochastic differential equation is explained in Scholz et al. (2000); Berkov (2007); Evans et al. (2014).

IV.4 The spin transfer torque

The STT term is defined as $STT=\frac{\hbar}{2\mu_{0}M_{s}a}m\times(m\times J^{spin})$ Datta et al. (2012), where $a$ is lattice constant, and $J^{spin}$ is the spin current that is calculated from the quantum transport, as explained in VI. The definition of STT term adopted in this paper is preferred over the Slonczewski STT term in the context of the quantum transport. The Slonczewski STT term has parameters that depend on the material and geometry of the junction to account for the junction efficiency of producing spin current. These effects are already included in the quantum transport. Therefore, we avoided using the Slonczewski STT term in favour of the quantum transport formulation of the spin current.

V The Ferroelectric Dynamics

The Landau-Devonshire (LD) expression of the free energy of FE material that describes the dependence of free energy of the FE material on the electric polarization and the electric field is defined as

[TABLE]

where $\alpha_{1}$, $\alpha_{11}$, and $\alpha_{111}$ are the free-energy expansion coefficients for bulk material Liu et al. (2013); Tan et al. (2001). The polarization of the material could be determined by minimizing the free energy (F) with respect to the electric polarization in (10). However, the LD expansion describes the static relation between electric field and polarization. The dynamic behavior of FE and its dependence on time is described by Landau-Khalatnikov (LK) equation:

[TABLE]

where $\lambda$ is the viscosity coefficient that represents the resistance of FE polarization motion toward the free energy minimum state.

VI Quantum Transport: Non-equilibrium Green’s Function modeling of MFTJ

The NEGF models the magnetization dependent tunneling current by splitting the device into two independent channels for the spin-up and spin-down carriers. The schematic diagram in Fig. 7 shows the device meshing and the magnetization dependent DOS. The spin-based channel Hamiltonian $H_{ch}$ and the left (right) contact Hamiltonian $H_{L(R)}$ Datta et al. (2012, 2010) of the MFTJ are defined as

[TABLE]

where $a$ is the length of the mesh element, N is the number of the mesh elements, q is the electron charge, $V_{a}$ is the applied voltage, $\phi_{BI}$ is the built-in potential, I is the identity matrix, i is the horizontal index of the Hamiltonian matrix, j is the vertical index of the Hamiltonian matrix, $\sigma$ are the Pauli spin matrices, $\Delta_{L(R)}$ is the splitting energy of the left (right) contact as illustrated in Fig. 7, respectively, and $M_{L(R)}$ is the normalized magnetization of the left (right) contact. The Hamiltonian tight binding parameters are defined as $\alpha_{ch}(K_{t})=2t_{ch}+\frac{\hbar^{2}K_{t}^{2}}{2m_{ch}^{*}}$, $\alpha_{L(R)}(K_{t})=2t_{L(R)}+\frac{\hbar^{2}K_{t}^{2}}{2m_{L(R)}^{*}}$, $t_{ch}=\hbar/(2m_{ch}^{*}a^{2})$, and $t_{L(R)}=\hbar/(2m_{L(R)}^{*}a^{2})$, where $K_{t}$ is the momentum vector in the transverse direction, $m_{L(R)}^{*}$ is the electron effective mass of left (right) electrodes, and $m_{ch}^{*}$ is the electron effective mass of the channel.

The Hamiltonian $H_{L(R)}$ dependence on the magnetization direction is modeled by the term $\left(I-\sigma.M_{L(R)}(i)\right)\frac{\Delta_{L(R)}}{2}$. The term $(qV_{a}+\phi_{BI}+\phi_{1}-\phi_{2})(\frac{N+1-2i}{2N+2})I$ in $H_{ch}$ linearly interpolates the applied field and built-in potential $V_{a}+\phi_{BI}+\phi_{1}-\phi_{2}$ over the channel. $U_{B}$ is the barrier height relative to the conduction band. The term $\phi_{1}-\phi_{2}$, that is estimated by (4-7), represents the dependence of the potential on the electric polarization of the FE insulator. The screening lengths estimated by the DFT method and the electric polarization estimated by LK equation are plugged in Thomas-Fermi relation (4) to determine $|\phi_{1}-\phi_{2}|$. The term $|\phi_{1}-\phi_{2}|$ is necessary for evaluating the NEGF Hamiltonian (17).

The Green’s function $G$ is defined as

[TABLE]

where $H$ is the full device Hamiltonian, $\Sigma_{L(R)}$ is the left (right) contact self-energy that is defined as

[TABLE]

where $K^{\uparrow}_{L(R)}$ is the left (right) contact longitudinal wave vector of spin-up electron given by

[TABLE]

and $K^{\downarrow}_{L(R)}$ is the left (right) contact longitudinal wave vector of spin-down electron given by

[TABLE]

Finally, the Landau’s current formula is defined as

[TABLE]

where t is the transmission coefficient of the channel given by the expression

[TABLE]

and $\Gamma_{L(R)}$ is the left (right) broadening function defined by

[TABLE]

The Fermi-Dirac distribution $f_{L}(R)$ is defined as

[TABLE]

The spin current is defined as Datta et al. (2010)

[TABLE]

where Gn defined as

[TABLE]

VII Time-dependent Formulation of Exchange interaction Coefficient based on Time-dependent Perturbation Theory

In the following discussion, we formulate a time-dependent formulation of the evolution from the FM to AFM phase induced by electric polarization. The sign of exchange interaction constant $J_{exch,i,j}$ is responsible for the magnetic order. The FM to AFM phase transition is controlled by the electrical polarization of the BTO. The change in the potential results from polarization switching can be considered as a small perturbation. Therefore we can use the time-dependent perturbation theory to model the time evolution of the exchange coefficient due to the perturbation potential Griffiths and Schroeter (2018). The time-dependent perturbation potential can be formulated from the Thomas-Fermi relation as

[TABLE]

where $V_{0}$ is the initial value of the potential. Assuming that the FM to AFM phase transition results from the transition from wave function $\psi_{b}$ to $\psi_{a}$. The time dependent wave function can be written as $\Psi(z,t)=a(t)\psi_{a}(z)e^{-i\frac{E_{a}}{\hbar}t}+b(t)\psi_{b}(z)e^{-i\frac{E_{b}}{\hbar}t}$. The time evolution of a(t) can be formulated as Griffiths and Schroeter (2018)

[TABLE]

where $\omega=\frac{E_{b}-E_{a}}{\hbar}$, $E_{a}$ is the energy of $\psi_{a}$, and $E_{b}$ is the energy of $\psi_{b}$. finally, the transition probability from state $\psi_{a}$ to state $\psi_{b}$ is $P_{b\rightarrow a}=|a(t)|^{2}$. The term $\left<\psi_{a}|e^{-\frac{|z|}{\delta_{1}}}|\psi_{b}\right>$ can be determined from the normalization of the probability density. However, the polarization as function of time P(t) is not known analytically. Therefore the integration (37) has to be numerically evaluated. The time evolution of the exchange interaction coefficient can be formulated as

[TABLE]

where $J_{exch,i,j}$ is the material-dependent magnetic exchange constant.

VIII Simulation Procedure and Parameters selection

In the previous discussion, we explained the quantum transport model, the magnetization LLG equation, and the LK equation, separately. In the following discussion, we explain the methodology we used to solve these models together to get the MFTJ characteristics. The steady-state characteristic of the MFTJ is calculated by the following procedure. Given the initial polarization $P$ and the external applied voltage $V_{a}$, the term $\frac{dF}{dP}$ is calculated by differentiating the LD equation analytically. The $\frac{dF}{dP}$ obtained in the previous step is substituted in the LK equation to get $\frac{dP}{dt}$. Then the forward difference formula is used to update the polarization $P(t+dt)=P(t)+\frac{dP}{dt}dt$. The previous steps are repeated iteratively until the electric polarization reaches its steady-state value. After the steady-state electric polarization is obtained, the electrostatic potential is calculated from (1) and (4). Next the current is calculated from the quantum transport model. However, the solution of the transport model is dependent on the magnetization directions of the electrodes which are calculated by the LLG equation. At the same time, the solution of the LLG equation depends on the current obtained from the NEGF equation. This raises the need for a self-consistent solution of the quantum transport and the LLG equation iteratively until the steady state current and magnetization are reached. The LLG equation is solved using Huen’s method Brown (1963); Evans et al. (2014). The existence of the thermal fluctuation in the LLG equation makes it a stochastic differential equation (SDE). The integration of the SDE is explained in Brown (1963); Evans et al. (2014); Berkov (2007).

The time-dependent response of the MFTJ is calculated using the following procedure. Step 1: the term $\frac{dF}{dP}$ is calculated from the LD equation, given the initial polarization $P$ and the external applied voltage $V_{a}$. The $\frac{dF}{dP}$ obtained in the previous step is substituted in the LK equation to get $\frac{dP}{dt}$. Then the forward difference formula is used to update the polarization as following $P(t+dt)=P(t)+\frac{dP}{dt}dt$. Step 2: the electrostatic potential is calculated from (1) and (4). Step 3: The spin current is calculated from the quantum transport model. The magnetization obtained from the solution of LLG and the electrostatic potential obtained in the previous step are used to solve the quantum transport. Step 4: the exchange coefficient $J^{\prime}_{exch,i,j}$ is calculated from (37)-(38) using the electric polarization at time $t$ obtained from the LK equation. Step 5: the LLG equation is solved to get $m(t+dt)$ using the spin current obtained from the quantum transport. Finally, the steps (Step 1) to (Step 5) are repeated iteratively at each time step.

The transport parameters are usually estimated by fitting the parameters on the experimental I-V characteristic Quindeau et al. (2015); Gruverman et al. (2009); Datta et al. (2012). However, the estimation process is not straight forward and more than one solution could produce the same transport properties. In this study, we try to go beyond that and predict some of the parameters from DFT calculations to reduce the complexity of the estimation process. On the other hand, some transport processes cannot be included in the DFT calculation. Note, the DFT by definition describes the system at the ground minimum energy state. In contrast, the quantum transport model exhibits nonequilibrium conditions by definition. Therefore, it is better to estimate certain parameters from experimental data to account for these limitations of DFT. Because of the aforementioned discussion, we adopted a combination of estimating the parameters directly from experimental data and estimating the parameters from DFT to improve the quality of parameter estimation.

The simulation parameters used in this study are selected according to the following criterion. The saturation magnetization, magnetic exchange constant are estimated from DFT calculations. The magnetic anisotropy is selected from the experimental study Yin et al. (2015). The magnetic damping factor is set within the acceptable range of similar structure. The screening lengths and the splitting energy are estimated by the DFT calculation. The effective mass are tuned within the acceptable range in literature to produce the experimental results. This is a very common procedure for selecting effective mass Quindeau et al. (2015); Gruverman et al. (2009); Datta et al. (2012, 2010).

The Landau-Devonshire equation parameters $\alpha_{1}$, $\alpha_{11}$, and $\alpha_{111}$ are calculated from the critical voltage at which the electric polarization is switched and the values of the polarization at zero voltage. The values of the critical voltage and polarization at zero voltage are known from the experimental results in Yin et al. (2015). The two values of the polarization at zero voltage are local minimum points of the free energy. The free energy has a maximum point at zero polarization. The maximum and minimum points are located at the zeros of the first derivative of the free energy and impose constraints on the sign of the second derivative of the free energy. In addition, the coefficient of the highest order term in the free energy has to be positive because the free energy has to reach positive infinity as the polarization reaches $\pm\infty$. We used a numerical grid search to solve for $\alpha_{1}$, $\alpha_{11}$, and $\alpha_{111}$ that considers the aforementioned constraints. The viscosity coefficient is estimated by a grid search to get $5ns$ switching time. Finally, the energy $E_{a}$ and $E_{b}$ in (37) are estimated from DFT calculations.

IX Simulation Results and Analysis

IX.1 Comparison with Experimental MFTJ Characteristics

The main advantage of spin-based NEGF is the ability to model the four resistance states of the MFTJ. The LSMO/LCMO/BTO/LSMO MFTJ in Yin et al. (2015) is the most appropriate device to demonstrate these physical characteristics. The simulation results of LSMO/LCMO/BTO/LSMO are illustrated in Fig. 8. The proposed framework can capture the majority of the MFTJ I-V characteristic for positive (low resistance) and negative (high resistance) polarization states, as illustrated in Fig. 8. The simulation results for LSMO/BTO/Co FTJ are in agreement with the experimental results Chanthbouala et al. (2012a), as illustrated in Fig. 9. The simulation parameters of the LK and the LLG equations are $\alpha_{1}=-1.0654\times 10^{9}$ $m/F$, $\alpha_{11}=-6.0878\times 10^{9}$ $m^{5}/(C^{2}.F)$, $\alpha_{111}=5.0499\times 10^{10}$ $m^{9}/(C^{4}.F)$, $dt=1\times 10^{-14}$ s, $K_{U}=7.8\times 10^{4}$ $erg/cm^{3}$, $\alpha=0.05$, $\lambda=1.8$ $s/F$, $j_{exch,i,j}=4.5\times 10^{5}$ $erg/cm^{3}$, $|E_{a}-E_{b}|=0.029$ $eV$ and $M_{s}=414.15$ $emu/cm^{3}$.

As explained in section III, the LCMO layer goes through a phase transition from the FM to the AFM alignment by the influence of the electric polarization switching. The phase transition is confirmed by the following experimental procedure Yin et al. (2015). Starting by applying an external magnetic field to the MFTJ, the change in the resistance of the MFTJ due to the increase of the external magnetic field is measured. In the case of positive polarization, the device resistance diminishes due to the increase of the external magnetic field. In contrast, the negative polarization state exhibits a constant resistance independent of the magnetic field. This behavior of the MFTJ resistance is explained by the influence of the external magnetic field on the AFM aligned Mn site and the ability of the external magnetic field to switch the AFM aligned Mn site back to FM alignment Yin et al. (2015). Fig. 10 illustrates the effect of the external magnetic field on the MFTJ resistance that is produced by the quantum transport and magnetization dynamics. The simulation mimics the same physical device behavior because the FM electrode Hamiltonian has a magnetization dependent term (23).

IX.2 Analysis of Various MFTJ Parameters

The TER is estimated at different values of the splitting energy, as illustrated in Fig. 11. The TER dependence on splitting energy originates from the ME effect that happens in the LCMO/BTO interface. Moreover, the $TMR_{\rightarrow}$ is consistently lower than $TMR_{\leftarrow}$ as illustrated in Fig. 8. This asymmetric behavior is due to the antiferromagnetic alignment of the Mn sites of the LCMO electrode in the case of $P_{\rightarrow}$ state that reduces the TMR effect at that state. In contrast, the LCMO exhibits FM alignment in the case of $P_{\leftarrow}$ state, and hence the $TMR_{\leftarrow}$ is higher compared to $TMR_{\rightarrow}$.

The asymmetry in the electrodes screening lengths is necessary for an FTJ to exhibit a TER effect, as explained in section II. However, the NEGF transport simulations show a significant TER ratio for a hypothetical device that has electrodes of identical screening lengths. Fig. 12 shows the TER ratio as a function of the built-in potential $\phi_{BI}$ along with the electrostatic potential of the positive and negative polarization states at different values of $\phi_{BI}$. Although the electrodes screening lengths are identical, the TER ratio increases significantly due to the increase of the built-in potential. The origin of the TER effect in the case of symmetric electrodes screening lengths could be explained by observing the potential profile of the positive and negative polarization states. In the case of $\phi_{BI}=0$, the potential profiles of $P_{\rightarrow}$ and $P_{\leftarrow}$ are symmetric, and hence the TER ratio is zero as expected. However, as the built-in potential $\phi_{BI}$ increases, the potential profile of $P_{\rightarrow}$ and $P_{\leftarrow}$ start to deviate from the symmetric shapes to asymmetric potential profiles that have different average barrier height, as illustrated in the inset of Fig. 12. In other words, the built-in potential introduces another source of asymmetry that allows the modulation of the barrier average height by the electric polarization.

The exponential dependence of the TER on the left and right electrode screening lengths ratio is demonstrated in Fig. 13. In order to understand the TER behavior as a function of $\frac{\delta 1}{\delta 2}$, we have reordered the TER definition as $TER=\frac{j_{\leftarrow}}{j_{\rightarrow}}-1$, where $j_{\rightarrow}$ and $j_{\leftarrow}$ are the current at positive (high resistance) and negative (low resistance) polarization states, respectively. The strength of the barrier height modulation, induced by polarization switching, is enhanced by increasing the difference between the screening lengths of the electrodes. Therefore, the current $j_{\rightarrow}$ increases and the current $j_{\leftarrow}$ decreases, as illustrated in Fig. 13. Consequently, the TER improves exponentially as the ratio $\frac{\delta 1}{\delta 2}$ reaches zero. The same conclusion can be quantitatively derived from the Thomas-Fermi relation (4) that formulates the potential at the interface as $\phi_{1}=\frac{\sigma_{s}\delta_{1}}{\epsilon_{0}\epsilon_{r1}}$. Therefore, the potential $\phi_{1}$ decreases as $\delta_{1}$ shrinks, and hence the the potential difference $||\phi_{1}|-|\phi_{2}||$ rises along with the ratio $\frac{j_{\leftarrow}}{j_{\rightarrow}}$. As a result of $\frac{j_{\leftarrow}}{j_{\rightarrow}}$ exponential increase, the TER exponentially improves.

Interestingly, the TER shows exponential dependence on the ratio $\frac{\delta_{1}}{\delta_{2}}$, but weaker dependence on the barrier height $U_{B}$. The rationale behind the difference in the dependence of TER on $\frac{\delta_{1}}{\delta_{2}}$ and $U_{B}$ is explained in the following comparative analysis. The decay of the ratio $\delta_{1}/\delta_{2}$ results in increasing $j_{\rightarrow}$ and decreasing $j_{\leftarrow}$ that exponentially enhance the TER. In contrast, the increase in the barrier height $U_{B}$ reduces both $j_{\rightarrow}$ and $j_{\leftarrow}$ but with different rates. Therefore, the TER changes with a weak rate because it is proportional to $\frac{j_{\leftarrow}}{j_{\rightarrow}}$, as illustrated in Fig. 14. However, the TER curve looks approximately linear because of the narrow range of $U_{B}$ along with the weak exponential dependence of the TER on $U_{B}$. The MFTJ high and low resistances is exponentially augmented as $U_{B}$ elevates, as observed from the currents $j_{\rightarrow}$ and $j_{\leftarrow}$ in Fig. 14. The barrier height is dependent on the insulator and electrodes work functions. Therefore, the electrodes work functions together with the screening lengths have a strong influence on the resistance and the TER of MFTJ, respectively.

IX.3 MFTJ Dynamic Characteristics

The time-dependent response of the MFTJ switching is illustrated in Fig. 15. The MFTJ dynamic response is calculated according to the procedure explained in section VIII. The electric polarization takes around $5ns$ to switch from from negative to positive value and saturate. The transition probability elevates from zero to one and saturates following the electric polarization, as illustrated in Fig. 15. The transition probability of the FM phase to the AFM phase is calculated from (37) that is based on the time-dependent perturbation theory. The exchange coupling coefficient follows the transition probability according to (38). Therefore, the exchange coefficient change from a positive value (FM alignment) to a negative value (AFM alignment). The magnetization of the second Mn atom switches from positive to negative magnetization following a magnetic precession motion because of the change in the exchange coupling. The precession motion of the magnetization causes the oscillation of the MFTJ current, as illustrated in Fig. 16.

Interestingly, the electric current decreases significantly after $4ns$ from the start of the switching process. The large variations of the current are due to the electric polarization switching. In contrast, the change in magnetic configuration lags the electric polarization switching, as illustrated in Fig. 15 and Fig. 16. During the switching process, the electric polarization changes from negative to positive direction passing by $P=0$. The electrostatic potential is modulated by the electric polarization as described by Thomas-fermi relation. The current reaches its minimum value because the tunneling current changes from Nordheim tunneling to direct tunneling at that point. The current changes back to Nordheim tunneling after the minimum point. Therefore, the current starts to increase after the minimum current point. Note, non-ideal contacts are assumed at high switching voltage to allow a constant voltage drop at each contact of $0.2V$ during quantum transport.

In contrast, the oscillations that start after the current minimum point are due to the precession motion of the Mn atoms. The precession motion of the Mn atoms is derived by the magnetic exchange torque of neighbor atoms. Due to the time-dependent perturbation potential caused by electric polarization switching, the probability of magnetization switching to AFM alignment increases from zero to one. The magnetization switching lags the switching probability by hundreds of picoseconds. The reason behind the delay in the magnetization switching is that the magnetization and the effective exchange field have almost an angle of $\pi$ at the initial position. Therefore, the magnetization motion under the effect of magnetic torque is slow at the beginning. The thermal fluctuations could assist the switching process at the slow-starting part of the switching process. However, the thermal fluctuations are small compared to the magnetic anisotropy due to the large area of the MFTJ.

X Conclusion

In this study, we propose a modeling and simulation framework that captures the behavior of MFTJ as a four-state device. Furthermore, the DFT method is used to estimate the screening length of the electrodes that have a strong influence on the TER. The estimated screening lengths are used in the quantum transport calculations to mimic the realistic device behavior at different bias voltages. The TER and TMR estimated by the proposed framework is compared with the experimental results of LSMO/LCMO/BTO/LSMO Yin et al. (2015). The quantum transport and magnetization dynamics could show the dependence of the device resistance on the external applied magnetic field. The dependence of the MFTJ resistance on the external magnetic field is in agreement with the experimental results in Yin et al. (2015) that confirms the transition from the FM to the AFM phase in the LCMO electrode.

Our analysis illustrates that not only the TMR but also the TER of the MFTJ depends on splitting energy because of the magnetoelectric effect at the interface that originates from the LCMO electrode phase transition from the FM to the AFM phase. On the other hand, the contrast between the weak and strong exponential dependence of the TER on the barrier height and electrodes screening length ratio, respectively is analyzed. The barrier height that is dependent on the electrodes and insulator work functions, could increase the MFTJ high and low resistance. Consequently, the power and speed of the MFTJ sensing could be enhanced by choosing the insulator and electrodes that have the appropriate work functions. However, the ratio of the electrodes screening length $\delta_{1}/\delta_{2}$ could exponentially enhance the TER ratio. Finally, our analysis reveals that the TER effect is improved by the asymmetry exhibited by the built-in potential that results in average barrier height modulation by electric polarization.

Based on the time-dependent perturbation theory, we could derive a mathematical formulation that relates the magnetic exchange interaction coefficient to the time evolution of electric polarization (37)-(38). This formulation is an important step toward a consistent model of MFTJs. The formulation emphasizes that the magnetization switching from FM to AFM alignment induced by polarization reversal follows a precessional motion. The transient response of the MFTJ exhibits a transition from Nordheim tunneling to direct tunneling and back to Nordheim tunneling current during the polarization switching. The transient response of the MFTJ demonstrates a set of oscillations due to the magnetic precessional motion. Although the switching from the FM to AFM alignment is induced by the electric polarization switching, the thermal magnetic fluctuations still assist the magnetization motion especially at the start of the switching.

Acknowledgement

The work was supported in part by Semiconductor Research Corporation(SRC), the National Science Foundation, Intel Corporation and by the DoD Vannevar Bush Fellowship.