Insights into Electron Transport in a Ferroelectric Tunnel Junction

The success of a ferroelectric tunnel junction (FTJ) depends on the asymmetry of electron tunneling as given by the tunneling electroresistance (TER) effect. This characteristic is mainly assessed considering three transport mechanisms: direct tunneling, thermionic emission, and Fowler-Nordheim tunneling. Here, by analyzing the effect of temperature on TER, we show that taking into account only these mechanisms may not be enough in order to fully characterize the performance of FTJ devices. We approach the electron tunneling in FTJ with the non-equilibrium Green function (NEGF) method, which is able to overcome the limitations affecting the three mechanisms mentioned above. We bring evidence that the performance of FTJs is also affected by temperature–in a non-trivial way–via resonance (Gamow-Siegert) states, which are present in the electron transmission probability and are usually situated above the barrier. Although the NEGF technique does not provide direct access to the wavefunctions, we show that, for single-band transport, one can find the wavefunction at any given energy and in particular at resonant energies in the system.


Introduction
The asymmetric ferroelectric tunnel junction (FTJ), a thin ferroelectric (FE) layer sandwiched between two dissimilar electrodes-or with different interfaces in the case of the same electrode material-is a promising electron device for many applications such as low power memories [1] or neuromorphic computing [2]. The salient mechanism of FTJ is based on the tunneling electroresistance (TER) effect, i.e., a change in the electrical resistivity when the ferroelectric polarization is reversed under an external electric field [3,4]. In other words, the electrical resistance of FTJ switches from a high conduction (ON) state to a low conduction (OFF) state or vice-versa when a voltage pulse is applied. Ferroelectrics like BaTiO 3 , PbZr 0.2 Ti 0.8 O 3 , BiFeO 3 [5][6][7], as well as high-k dielectrics like HfO 2 and Hf 0.5 Zr 0.5 O 2 [8,9] have been shown to work as FTJs. The latter are particularly attractive since they are compatible with the CMOS technology. The TER effect can be characterized by the TER ratio where J ON/OFF are the current densities in the two different states. The higher the TER ratio and J ON , the better the FTJ performance. The extent of the barrier modulation depends on numerous factors, such as: (i) the thickness and spontaneous polarization of the ferroelectric film, (ii) the bias voltage, (iii) the differences between work functions and screening lengths of the electrodes, and (iv) the built-in field and screening of the polarization charge, as well as on the variation of barrier thickness due to piezoelectricity [10]. In some FTJ implementations, an ultrathin dielectric (DE) layer is also introduced, as shown in Figure 1, in order to ensure an additional way to tune the FTJ performance. In semi-empirical methods based on NEGF, the parameters can be adjusted to match the experimental data [14]. In many cases, in order to characterize the electron transport and interpret the I-V curves, various models of transport mechanisms, involving direct tunneling, thermionic emission, and Fowler-Nordheim tunneling, are used [16]. Direct tunneling presumes low energy for incident electrons such that the tunneling probability does not depend on the profile of the barrier, which is assumed to be rectangular on average in the case of Simmons formula [17,18] or trapezoidal in the case of Brinkman et al. formula [18,19]. The thermionic emission is that part of the current carried over the top of the barrier by thermally activated charge carriers [20]. The last mechanism, the Fowler- Figure 1. Schematic representation of a FTJ with a composite barrier, consisting of a ferroelectric (FE) layer and a dielectric (DE) layer between two metallic electrodes (ML and MR). The dotted line delimitates the device that is considered in calculations, λ 1,2 are the screening lengths in electrodes, ε FE,DE the relative dielectric constants of the ferroelectric and dielectric layers, ε 1,2 the relative dielectric constants of the electrodes, and t FE , DE , ML , MR mark the margins of the layers. In the case of a simple FE barrier, the DE layer is not present.
To understand, describe, and model the TER effect, the ab-initio methods may provide an accurate image of the physics governing FTJ systems including the geometry, polarization state, electronic structure [11], and transport properties [12,13]. Nevertheless, ab-initio methods require large computational resources, which may prevent their intensive use in the early stages of FTJs design, when there is a requirement for a fast scan in parameter spaces to obtain a device with required basic functionalities. In the search for optimized devices, the semi-empirical methods are fast and, despite all simplifications, are extremely useful for the treatment of the charge transport as long as they are fed with appropriate physical parameters [14,15]. These semi-empirical approaches are based on the non-equilibrium Green function (NEGF) method that can calculate exactly, in principle, the tunneling current and the I-V characteristic of a FTJ [15].
In semi-empirical methods based on NEGF, the parameters can be adjusted to match the experimental data [14]. In many cases, in order to characterize the electron transport and interpret the I-V curves, various models of transport mechanisms, involving direct tunneling, thermionic emission, and Fowler-Nordheim tunneling, are used [16]. Direct tunneling presumes low energy for incident electrons such that the tunneling probability does not depend on the profile of the barrier, which is assumed to be rectangular on average in the case of Simmons formula [17,18] or trapezoidal in the case of Brinkman et al. formula [18,19]. The thermionic emission is that part of the current carried over the top of the barrier by thermally activated charge carriers [20]. The last mechanism, the Fowler-Nordheim tunneling, is used for triangular barrier profiles, which are encountered at finite applied bias voltages [21]. The treatment of the above mechanisms is based on several approximations for electron transport. The Simmons and Brinkman et al. as well as the Fowler-Nordheim formulae are based on the semi-classical Wentzel-Kramers-Brillouin (WKB) approximation [17,21], while thermionic emission tacitly assumes a homogeneous transmission probability of one for all electrons with energy above the top of the barrier [20]. In this work, by the exact treatment of the tunneling within the NEGF method, we show that the above models (describing direct tunneling, thermionic emission, and Fowler-Nordheim tunneling) are not accurate enough for performance characterization of FTJs. These models cannot capture resonant features emerging in FTJ structures that might play an important role in carrier transport and ultimately in the performances of the device. The electric current of an impinging electron at a given energy is determined by the electron transmission probability at that energy weighted (multiplied) by its Fermi-Dirac distribution function. Hence, even though they seem unlikely to play a role in the tunneling due to the relatively simple structure of the barriers, the resonances can contribute significantly to the electric current since their contribution to transmission probability may offset the temperature-dependent Fermi-Dirac factor. In the following, we will show that the resonances emerging close or above the top of the barriers play a decisive role for transport at room temperature, a fact that is not captured by the models previously mentioned. Moreover, resonances may appear in FTJs with composite (both a ferroelectric and a dielectric) barriers. The NEGF method can also accurately describe these type of heterostructures. Generally, the NEGF approach cannot provide full access to the wavefunctions, especially for multi-band transport [22]. However, for a single-band transport, we show that this is not the case; the eigenvectors of the spectral functions are just the wavefunctions of the system. At this stage, we are able to identify the Green function of the device as the discrete version of the outgoing Green function, which is used in diverse scattering problems like nuclear reactions or electron transport in nanostructures [23,24]. The advantage of such correspondence is that the NEGF of the device and the transmission probability can be expanded as a sum of resonances or Gamow-Siegert states [23][24][25][26]. Thus, we are able to sort out the resonances that count for the electron transport in FTJs.
The paper is organized as follows: The next section deals with the theoretical background of the calculations of applied voltage induced electrostatic potential across the FTJ, of transport quantities like I-V curve and conductance by NEGF method, and the evaluation of wavefunctions and resonance states from NEGF calculations. Section 3 presents numerical results regarding the effect of temperature on electric conductance and TER ratio for several practical cases of simple and composite BaTiO 3 -based FTJs. In addition, the wavefunctions and resonance states for such FTJ structures are analyzed. Section 4 summarizes the conclusions of this work, and lastly, in Appendix A we present the steps to obtain a discrete tight-binding Hamiltonian from the BenDaniel and Duke Hamiltonian [22].

The Profile of Potential Barrier
In semi-empirical methods, the calculation of tunneling current and electric conductance is performed by simultaneously solving the electrostatic and transport problems. They are intricately interrelated since the charge density in Poisson equation is calculated self-consistently from the NEGFs [22]. However, when one deals with free carriers, only in the metallic contacts they can decouple each other. Thus, let us first deal with the electrostatics and the profile of the potential barrier. We assume that the dielectric and ferroelectric are located at x between −t DE and 0 and x between 0 and t FE , respectively, while the metallic contact ML is at x < −t DE and the metallic contact MR is at x > t FE , as illustrated in Figure 1. A widely used and good approximation of electrostatics in metals is the Thomas-Fermi approximation [14,17]. Within this framework, the electrical fields in ML and MR are given by: In Equations (2) and (3) τ S is the screening charge at ML/dielectric and MR/ferroelectric interfaces, λ 1 (ε 1 ) and λ 2 (ε 2 ) are the Thomas-Fermi screening lengths (dielectric constants) of ML and MR, respectively, and ε 0 is the vacuum permittivity. In the following, we denote by ε FE , E FE , and P the dielectric constant, the electric field, and the intrinsic polarization, respectively, of the ferroelectric, and by E DE and ε DE the electric field and the dielectric constant of the dielectric, respectively. The electrostatic equations are in fact the conti- Additionally, the bias voltage, V, across the structure obeys the equation: where V BI is the built-in voltage due to mismatch of the conduction band discontinuities and of Fermi energies E FL (E FR ) of ML (MR).
In (6), e is the elementary electric charge, ϕ 1 is the band discontinuity at the first interface between ML and the dielectric, ϕ 2 is the band discontinuity at the second interface between the ferroelectric and MR, and ϕ C is the band discontinuity at the interface between the dielectric and ferroelectric when the ferroelectric is unpolarized. Eliminating E FE and E DE from (4) and (5), we obtain the screening charge τ S Knowing that the electric field is homogeneous in both the dielectric and the ferroelectric, now it is straightforward to obtain the electric potential. Then, the tunneling barrier profile for electrons has the following form: In Equation (8), P is considered positive when points from ML to MR, and negative when it points the other way around.

Transport by NEGF
Transport properties like electric current density and conductance can be calculated by solving the Schrödinger equation with scattering boundary conditions, i.e., either an incoming wave from the left or from the right. For a single band, the equation is merely 1D with a BenDaniel and Duke Hamiltonian, where the effective mass of electrons can vary across the structure [22]: In Equation (9), it is assumed that the energy band is parabolic with an isotropic effective mass in each layer of the heterostructure, k is the transverse wavevector (parallel to each interface), m* is the effective mass, and U(x) the potential energy given by (8). Equation (9) can be discretized into a 1D problem in a tight-binding representation [22]. The details are presented in Appendix A. In general, in the tight binding representation, the Hamiltonian has a block matrix form, or more precisely, a block tridiagonal form, as in Equation (10): For a single band model, the block tridiagonal matrix form of H becomes a simple tridiagonal matrix. In this format, one can easily see that the Hamiltonian has three parts: the semi-infinite left (L) and right (R) metallic electrodes and the device (D) that is defined by H D , an n D × n D matrix. Both electrodes act as reservoirs, hence they have well-defined chemical potentials and temperatures. We can define the retarded Green function (G R ) of the system at energy E as the inverse matrix of In the NEGF method, one eliminates the degrees of freedom of contacts by introducing self-energies into the projected Green functions of the device (D).
The self-energy Σ R,A B (E) = Σ R,A L (E) + Σ R,A R (E) has two components due to the coupling to the left and right contacts. It replaces the boundary conditions that otherwise would be fulfilled by the construction of a Green function in the device region. The self-energy due to the left contact has the following expression: where is the Green function of the semi-infinite left contact.
Additionally, Σ R,A R (E) has a similar expression. Here we notice that due to the fact that we deal with a nearest-neighbor tight-binding Hamiltonian, the self-energies Σ R L and Σ R R are n D × n D matrices with just one non-zero element, the (1,1) element for Σ R L and (n D ,n D ) for Σ R R . Thus, the retarded Green function is: The poles of the device Green function are no longer real-an attribute of an open quantum system. Starting from Equation (11), the calculation of σ R L (E) involves the calculation of matrix element (0,0) of g R L (E), which is the surface Green function of the left contact. Bearing in mind that D L and t L are the tight-binding parameters of the left contact and defining a longitudinal wavevector k l , the energy can be parameterized according to the single-band dispersion relation for the left contact E = D L − 2t L cos(k l ∆). One can find that the matrix element (0,0) of [14,22,27]. The expression of g R L (E) is calculated with outgoing boundary condition [27], hence the self-energy is for this kind of boundary condition.
One can further define: (a) the spectral function A = i(G r − G a ), which also has a matrix form-whose diagonal is just the local density of states up to a 2π factor-and (b) the broadening function due to the coupling to the left and right contacts Using the projection operator on the device D, one can show that the projection of the full spectral function A on the device space is [28].
Moreover, one can further show that partial spectral densities are spectral densities due to incident Bloch waves that come from the left (L) and from the right (R), respectively [28]. Thus, G R D (E) contains information about both solutions of the scattering problem. Furthermore, it can be shown that the current flow from an incident Bloch wave that comes from the left electrode into the right electrode is is just the matrix of transmission probability T(E) that has a schematic representation in Figure 2. From Equation (17), the expression of the total current takes the form of the Landauer-Büttiker formula [15,22]: In addition, in the linear regime (small bias voltages V) and at temperature of 0 K, we obtain the Landauer conductance formula [15] ( ) ( ) The behavior of an incident Bloch wave Ψ nk+ coming from the left. One part is reflected with the coefficient r and the other part is transmitted with a coefficient t. In case of a single-band transport, the coefficients r and t are simple scalars.
In Equation (18), the Tr () operation includes both the trace over the T(E) matrix and the integration over the transverse wavevector k. The functions f L,R are the Fermi-Dirac distribution functions of the L and R electrodes. Performing the trace operation only over the T matrix, we obtain the transmission probability coefficient t, and the Landauer-Büttiker formula becomes: In addition, in the linear regime (small bias voltages V) and at temperature of 0 K, we obtain the Landauer conductance formula [15]

Retrieving the Wavefunction from the Spectral Function. Resonance States
Let us consider the spectral function . It signifies the projected spectral function on the device for incident waves coming from the left. As it was pointed out in Ref. [28], in general the eigenvectors of A L cannot be identified with the wavefunction in the device region since there are several eigenvectors of A L with non-zero eigenvalues. This is rather obvious in the tight-binding representation, but for the multiband problem. For a single-band problem, however, this is not the case. A L has just one non-zero eigenvalue. Its corresponding eigenvector is just the function that is proportional to the wavefunction in the device region. This statement can be proven directly. The matrix Γ L (E) has just one element different from 0, like its corresponding self-energies Σ R,A L (E) is an eigenvector of A L (E) with the eigenvalue The fact that λ L = Tr(A L (E)) ensures that λ L is the only non-zero eigenvalue of A L (E). In a bra and ket notation, we thus write A L (E) as since the norm of v L is just λ L /γ L . Equations (22) and (24) of v L and A L (E) guarantee that v L is proportional to the solution of Schrödinger equation for incoming waves from the left projected on the device space: Similar calculations can be performed for A R (E), explicitly, the matrix form is the eigenvalue and the simple bra and ket form Also Equations (26) and (28) of v R and A R (E) guarantee that v R is proportional to the solution of Schrödinger equation for incoming waves from the right, projected on the device space: The total spectral function A D (E) is the sum of A L (E) and A R (E), hence its range is spanned by v L and v R with two eigenvalues λ 1 and λ 2 . They obey the following equation: Also it is easy to find that the squared modulus of the overlap between |Ψ L D and |Ψ R D is From the analysis we are going to perform in the next section on realistic examples, we will see that for energies in the direct tunneling regime there are two distinct solutions with low overlap. In this case, λ 1 and λ 2 are close to λ L and λ R . In the opposite case, at resonance, one of the two eigenvalues λ 1 or λ 2 is much larger than the other, hence the overlap is large. In a similar manner, we can obtain the explicit expression of the transmission probability coefficient in terms of the Green function Equation (31) has previously been deduced using an iterative procedure to calculate the Green function [22,29]. It additionally shows that the transmission probability coefficient has the spectral properties of the Green function. There is a large body of work in which the Green function G R D (E) is set in a meaningful representation. Such a representation is given by the expansion of the Green function in resonance (Siegert-Gamow) states [23,24], which lead to the Breit-Wigner formula for resonances in transmission. Here we will outline a few results about resonant states representation that are connected with our results discussed above. The expansion of the Green function in resonance states is the sum over these states plus a background, and it may take the following form: where u n (x) is a resonant state that satisfies the Schrödinger equation with outgoing boundary conditions at the device boundary. Since these boundary conditions are not hermitian, the eigenvalues are rather complex, i.e., E n = E rn − iΓ n /2. Resonant states come into pairs with a negative and positive imaginary part [23,25,26]. Now, it is rather obvious that for energy E far from any resonant energy E rn , the solution of the Schrödinger equation for incoming waves from the left is different from the solution of the Schrödinger equation for incoming waves from the right due to different mixing of the tails of various resonances. However, for energy E near a resonant energy E rn , both solutions from the left and from the right are similar since the dominant term is that given by u n (x). Finally, it is quite straightforward to notice by using Equation (29) that the transmission probability coefficient has a multi-resonance Breit-Wigner-like form [24] t The term is the Breit-Wigner expression, T nm (E) of the interference term between resonances, and R(E) is the term resulting from the background.

Temperature Influence on Conductance and TER Ratio
Numerically, we calculated all quantities defined in the previous section using the NEGF formalism. The calculation of I-V characteristics is performed with Equation (19), which is able to reproduce the non-linear regime for large bias voltages and finite temperatures. Equation (20) describes the linear regime at 0 K and is often used (especially in ab-initio calculations, where a full I-V curve is costly) in the calculation of the TER ratio defined by Equation (1). In the following we shall analyze a few practical FTJ structures.
3.1.1. Pt/BaTiO 3 /SrRuO 3 FTJ The first system to deal with is that of a BaTiO 3 ferroelectric barrier on top of a SrRuO 3 substrate which acts as the right contact layer in Figure 1, and Pt is the left contact to the barrier. The physical parameters are: λ 1 = 0.45Ǻ, λ 2 = 0.75Ǻ, ε 1 = 2, ε 2 = 8.45, ε FE = 125, E FL = E FR = 3 eV, ϕ 1 = ϕ 2 = 3.6 eV, P = 16 µC/cm 2 . The effective masses are: 5 m 0 for SrRuO 3 , 2 m 0 for BaTiO 3 , and m 0 for Pt, where m 0 is the free electron mass [30]. The Fermi energy is set to E F = E FL = E FR = 3 eV. One should note that the electron effective masses are different in the three layers of the structure, hence the calculations need full numerical integration over transverse k wavevector in Equations (19) and (20). In Figure 3 we show the potential profiles of a 2 nm thick BaTiO 3 barrier for both directions of polarization. The transmission probability coefficients, also depicted in the right of Figure 3, exhibit resonances above the very top of the barriers. The resonances are associated with the peaks in the transmission spectra and if they are well resolved then they obey Equation (35). In Figure 4a, we plot the conductance of the system, and in Figure 4b the TER ratio at 0 K and 300 K.
is the Breit-Wigner expression, ( ) nm T E of the interference term between resonances, and ( ) R E is the term resulting from the background.

Temperature Influence on Conductance and TER Ratio
Numerically, we calculated all quantities defined in the previous section using the NEGF formalism. The calculation of I-V characteristics is performed with Equation (19), which is able to reproduce the non-linear regime for large bias voltages and finite temperatures. Equation (20) describes the linear regime at 0 K and is often used (especially in abinitio calculations, where a full I-V curve is costly) in the calculation of the TER ratio defined by Equation (1). In the following we shall analyze a few practical FTJ structures.
3.1.1. Pt/BaTiO3/SrRuO3 FTJ The first system to deal with is that of a BaTiO3 ferroelectric barrier on top of a SrRuO3 substrate which acts as the right contact layer in Figure 1, and Pt is the left contact to the barrier. The physical parameters are: λ1 = 0.45 Ǻ, λ2 = 0.75 Ǻ, ε1 = 2, ε2 = 8.45, εFE = 125, EFL = EFR = 3 eV, φ1 = φ2 = 3.6 eV, P = 16 μC/cm 2 . The effective masses are: 5 m0 for SrRuO3, 2 m0 for BaTiO3, and m0 for Pt, where m0 is the free electron mass [30]. The Fermi energy is set to EF = EFL = EFR = 3 eV. One should note that the electron effective masses are different in the three layers of the structure, hence the calculations need full numerical integration over transverse k wavevector in Equations (19) and (20). In Figure 3 we show the potential profiles of a 2 nm thick BaTiO3 barrier for both directions of polarization. The transmission probability coefficients, also depicted in the right of Figure 3, exhibit resonances above the very top of the barriers. The resonances are associated with the peaks in the transmission spectra and if they are well resolved then they obey Equation (35). In Figure 4a, we plot the conductance of the system, and in Figure 4b the TER ratio at 0 K and 300 K.  The calculation of the conductance at 300 K, see Figure 4a, was performed with an applied bias of 0.0001 V, which is low enough to satisfy the linear regime conditions. Both the conductance and TER ratio behave differently at 300 K with respect to 0 K. At 300 K we observe two regimes depending on the ferroelectric thickness-one regime up to 2.5 nm and another one beyond this value. The difference between the two regimes is explained in Figure 5. At 0 K, the channels open to electron flow are those up to the Fermi energy, represented by dashed line. These channels are those of direct tunneling such that the Simmons and Brinkman formulae are appropriate [17][18][19]. For the barrier of 1.6 nm thickness (Figure 5a), the transport occurs around Fermi energy even at room temperature (300 K), hence the similar behavior of the conductance at 300 K with respect to 0 K, see Figure 4a. There is a small contribution to the transmitted current from the resonance states, but this contribution is orders of magnitude smaller. However, for thicker barriers like the one shown in Figure 5b, things may change. First, the direct tunneling current is much smaller since it has an exponential dependence on the barrier thickness. Second, the resonance levels are spaced much closer, hence their contribution to transmission increases considerably. Even if their occupancy-given by the Fermi-Dirac function-might be small, it is offset by the large transmission probability. In Figure 5b it can easily be seen that the contribution from the resonance states is overwhelmingly larger than the contribution from the states near the Fermi energy. Moreover, electrons with energies closer to the barrier top encounter a triangular barrier profile and so they may still reach a resonance state; this mechanism cannot be described within the semiclassical WKB theory of Fowler-Nordheim [21]. A similar behavior of the TER ratio at room temperature was also observed in experimental data [31]; a degradation of TER was even found when increasing the ferroelectric thickness. However, performing their analysis based on Brinkman model at finite bias voltage, the authors attributed the TER degradation to the rather high levels of noise in the measurements.  The calculation of the conductance at 300 K, see Figure 4a, was performed with an applied bias of 0.0001 V, which is low enough to satisfy the linear regime conditions. Both the conductance and TER ratio behave differently at 300 K with respect to 0 K. At 300 K we observe two regimes depending on the ferroelectric thickness-one regime up to 2.5 nm and another one beyond this value. The difference between the two regimes is explained in Figure 5. At 0 K, the channels open to electron flow are those up to the Fermi energy, represented by dashed line. These channels are those of direct tunneling such that the Simmons and Brinkman formulae are appropriate [17][18][19]. For the barrier of 1.6 nm thickness (Figure 5a), the transport occurs around Fermi energy even at room temperature (300 K), hence the similar behavior of the conductance at 300 K with respect to 0 K, see Figure 4a. There is a small contribution to the transmitted current from the resonance states, but this contribution is orders of magnitude smaller. However, for thicker barriers like the one shown in Figure 5b, things may change. First, the direct tunneling current is much smaller since it has an exponential dependence on the barrier thickness. Second, the resonance levels are spaced much closer, hence their contribution to transmission increases considerably. Even if their occupancy-given by the Fermi-Dirac function-might be small, it is offset by the large transmission probability. In Figure 5b it can easily be seen that the contribution from the resonance states is overwhelmingly larger than the contribution from the states near the Fermi energy. Moreover, electrons with energies closer to the barrier top encounter a triangular barrier profile and so they may still reach a resonance state; this mechanism cannot be described within the semiclassical WKB theory of Fowler-Nordheim [21]. A similar behavior of the TER ratio at room temperature was also observed in experimental data [31]; a degradation of TER was even found when increasing the ferroelectric thickness. However, performing their analysis based on Brinkman model  3.1.2. Pt/ SrTiO3/BaTiO3/SrRuO3 Composite Barrier FTJ The second system we analyze is a composite FTJ, where a dielectric barrier (SrTiO3) is added besides the BaTiO3 ferroelectric barrier. The electrodes are Pt (left) and SrRuO3 (right). The role of the dielectric layer is to increase the asymmetry of the system, hence it

Pt/ SrTiO3/BaTiO3/SrRuO3 Composite Barrier FTJ
The second system we analyze is a composite FTJ, where a dielectric barrier (SrTiO3) is added besides the BaTiO3 ferroelectric barrier. The electrodes are Pt (left) and SrRuO3 (right). The role of the dielectric layer is to increase the asymmetry of the system, hence it is expected a higher TER ratio [32]. Physical parameters are slightly changed with respect to the previous case [32,33]: λ1 = 0.45 Ǻ, λ2 = 0.8 Ǻ, ε1 = 2, ε2 = 8.45, εFE = 90, εDE = 90, EFL = EFR = 3 eV, φ1 = φ2 = 3.6 eV, φC = 0 eV, P = 20 μC/cm 2 . The effective masses are: 5 m0 for SrRuO3, 2 m0 for BaTiO3 and SrTiO3, and m0 for Pt. In the following, BaTiO3 thickness is fixed to 2.4 nm and we varied the SrTiO3 thickness between 0.5 and 3 nm. The results of the calculations are presented in Figure 6. At 0 K, the conductance for both polarization directions shows an exponential dependence on the dielectric thickness. At room temperature, on the other hand, it appears that this exponential dependence is no longer valid, as revealed by the TER ratio. In this case, a decrease in TER takes place when the dielectric thickness increases (Figure 6b). The At 0 K, the conductance for both polarization directions shows an exponential dependence on the dielectric thickness. At room temperature, on the other hand, it appears that this exponential dependence is no longer valid, as revealed by the TER ratio. In this case, a decrease in TER takes place when the dielectric thickness increases (Figure 6b). The explanation of TER degradation is provided by the results presented in Figure 7. One may observe that the electric charges are mainly transported through quantum states that are close or above the barrier top, the contribution of the states near the Fermi energy being negligible. The total barrier thickness is 3 nm at least, a value at which the resonance states start to play a significant role. As the barrier thickness increases, the contribution of the resonance states is enhanced. Nevertheless, those resonances located above the barrier (see the F-D weighted curves in Figure 7 right panel) become less sensitive to the barrier profile, and hence the decline in the TER ratio with the barrier thickness decrease.

Metal/CaO/BaTiO 3 /Metal Composite Barrier FTJ
The third system studied herein is also a composite FTJ, with a CaO dielectric barrier added to the BaTiO 3 ferroelectric barrier. The electrodes are of the same generic metal, Me. In this case the asymmetry is ensured just by the presence of the dielectric. We have set the physical parameters to [32]: λ 1 = λ 2 = 1.0Ǻ, ε 1 = ε 2 = 1, ε FE = 90, ε DE = 10, E FL = E FR = 3 eV, ϕ 1 = 5.5 eV, ϕ 2 = 3.6 eV, ϕ C = 1.9 eV, P = 40 µC/cm 2 . The effective masses are equal to m 0 for all materials. We have also kept the thickness of BaTiO 3 to 2.4 nm and we have varied the thickness of CaO from 0.5 to 3 nm. In comparison to the previous system, in this particular case the barrier is much higher on the dielectric. The calculated conductance and TER ratio are presented in Figure 8.
observe that the electric charges are mainly transported through quantum states that are close or above the barrier top, the contribution of the states near the Fermi energy being negligible. The total barrier thickness is 3 nm at least, a value at which the resonance states start to play a significant role. As the barrier thickness increases, the contribution of the resonance states is enhanced. Nevertheless, those resonances located above the barrier (see the F-D weighted curves in Figure 7 right panel) become less sensitive to the barrier profile, and hence the decline in the TER ratio with the barrier thickness decrease.

Metal/CaO/BaTiO3/Metal Composite Barrier FTJ
The third system studied herein is also a composite FTJ, with a CaO dielectric barrier added to the BaTiO3 ferroelectric barrier. The electrodes are of the same generic metal, Me. In this case the asymmetry is ensured just by the presence of the dielectric. We have set the physical parameters to [32]: λ1 = λ2 = 1.0 Ǻ, ε1 = ε2 = 1, εFE = 90, εDE = 10, EFL = EFR = 3 eV, φ1 = 5.5 eV, φ2 = 3.6 eV, φC = 1.9 eV, P = 40 μC/cm 2 . The effective masses are equal to m0 for all materials. We have also kept the thickness of BaTiO3 to 2.4 nm and we have varied the thickness of CaO from 0.5 to 3 nm. In comparison to the previous system, in this particular case the barrier is much higher on the dielectric. The calculated conductance and TER ratio are presented in Figure 8. In contrast to Pt/SrTiO3/BaTiO3/SrRuO3, the conductance and the TER ratio of Me/CaO/BaTiO3/Me composite FTJ exhibit an exponential dependence on the dielectric thickness at both 0 K and 300 K. Therefore, qualitatively at least, a 0 K analysis remains valid also at room temperature. In Figure 9, we show the data for a quantitative explanation of conductance and TER ratio behavior with temperature. One can see that the resonance states have a minor contribution to the current; the vast majority of carriers are transported through states around Fermi energy, although there are many resonances below the top of the barrier. Due to the higher dielectric barrier, these resonant states are just weakly coupled to the contacts, hence they show low transmission probability coefficients and they have a modest contribution to conductance. In contrast to Pt/SrTiO 3 /BaTiO 3 /SrRuO 3 , the conductance and the TER ratio of Me/CaO/BaTiO 3 /Me composite FTJ exhibit an exponential dependence on the dielectric thickness at both 0 K and 300 K. Therefore, qualitatively at least, a 0 K analysis remains valid also at room temperature. In Figure 9, we show the data for a quantitative explanation of conductance and TER ratio behavior with temperature. One can see that the resonance states have a minor contribution to the current; the vast majority of carriers are transported through states around Fermi energy, although there are many resonances below the top of the barrier. Due to the higher dielectric barrier, these resonant states are just weakly coupled to the contacts, hence they show low transmission probability coefficients and they have a modest contribution to conductance. valid also at room temperature. In Figure 9, we show the data for a quantitative explana-tion of conductance and TER ratio behavior with temperature. One can see that the resonance states have a minor contribution to the current; the vast majority of carriers are transported through states around Fermi energy, although there are many resonances below the top of the barrier. Due to the higher dielectric barrier, these resonant states are just weakly coupled to the contacts, hence they show low transmission probability coefficients and they have a modest contribution to conductance.

The Tunneling Wavefunctions. The Wavefunctions of Resonances
In the previous section, it was discussed that the energy-dependent transmission in FTJs has resonant features besides an exponential background. In this section we show that the wavefunctions also exhibit general features, some of which-e.g., those corresponding to resonances-are also observed in other nanostructures like quantum wells, multi-barrier structures, etc. Usually, as a scattering problem, the electron transport in FTJs has two solutions at a given energy: one solution for an incident electron wave coming from the left and the other for the electron wave coming from the right. In the following we will illustrate the wavefunctions at some representative energy values for two FTJs: Pt/BaTiO 3 /SrRuO 3 and Pt/SrTiO 3 /BaTiO 3 /SrRuO 3 .

The Wavefunctions of Pt/BaTiO 3 /SrRuO 3 FTJ
We illustrate the wave functions at the Fermi energy (3 eV) and at some resonant energies that can be taken from Figure 10.

The Tunneling Wavefunctions. The Wavefunctions of Resonances
In the previous section, it was discussed that the energy-dependent transmission in FTJs has resonant features besides an exponential background. In this section we show that the wavefunctions also exhibit general features, some of which-e.g., those corresponding to resonances-are also observed in other nanostructures like quantum wells, multi-barrier structures, etc. Usually, as a scattering problem, the electron transport in FTJs has two solutions at a given energy: one solution for an incident electron wave coming from the left and the other for the electron wave coming from the right. In the following we will illustrate the wavefunctions at some representative energy values for two FTJs: Pt/BaTiO3/SrRuO3 and Pt/SrTiO3/BaTiO3/SrRuO3.

The Wavefunctions of Pt/BaTiO3/SrRuO3 FTJ
We illustrate the wave functions at the Fermi energy (3 eV) and at some resonant energies that can be taken from Figure 10. It is obvious that up to about 3.5 eV, the transmission probability has an exponential dependence with respect to energy. In Figure 11, we plotted the wavefunctions of a 1.6 nm thick BaTiO3 barrier at Fermi energy, as well as the first few resonance energies for both directions of polarizations. In order to compare these wavefunctions we plotted the "normalized" eigenvectors of AL, AR, and AD. In other words, for instance, It is obvious that up to about 3.5 eV, the transmission probability has an exponential dependence with respect to energy. In Figure 11, we plotted the wavefunctions of a 1.6 nm thick BaTiO 3 barrier at Fermi energy, as well as the first few resonance energies for both directions of polarizations. In order to compare these wavefunctions we plotted the "normalized" eigenvectors of AL, AR, and AD. In other words, for instance, |Ψ L,R D are divided by λ L,R /2π. The device region is comprised of the barrier and several layers of contact regions, such that the region outside the device should exhibit a flat electrostatic potential. In practice this is achieved when a few nanometers of contacts are added to the device region. At Fermi energy, the wavefunctions are almost real and their overlap is almost zero. They decay exponentially in the barrier, hence the Simmons or Brinkman formula applies well for states around Fermi energy.   Figure 12. The transmission probability for polarizations is shown in Figure 12a, from which we can extract the resonances that play a role in the electron transport at room temperature. The Since the barrier is thin, the resonances are well separated. From Figure 11a,b,d,e we see that the "normalized" |Ψ L D is almost identical with the complex conjugate of "normalized" |Ψ R D , however the values of λ L,R /2π provide the levels of electron density inside the device for the corresponding wavefunction. These states exhibit a confining character within the barrier, in contrast to the states shown at Fermi energy. These solutions belong to resonance and anti-resonance states in the complex wavevector plane [23,25,26]. Moreover, the real and the imaginary parts of |Ψ L,R D can be found in the normalized eigenvectors of A D (Figure 11c,f). Finally, by analyzing Figures 5 and 11, we notice that just the first resonance would participate to the electron transport at room temperature.

The Wavefunctions of Pt/SrTiO 3 /BaTiO 3 /SrTiO 3 FTJ
The plots of the wave functions for this composite barrier FTJ are shown in Figure 12. The transmission probability for polarizations is shown in Figure 12a, from which we can extract the resonances that play a role in the electron transport at room temperature. The physical parameters are those that have already been used in the calculations. The effective thickness of the barrier is larger than in simple FTJ presented in the previous subsection; hence the resonances are closely spaced. Still, at Fermi energy the wavefunctions decay exponentially in the barrier and are almost real, while their overlap is almost zero. The wavefunctions of the first four resonances in transmission have a much smaller imaginary part, which is not shown here. As in the previous section, these wavefunctions manifest confining character in the barrier. Moreover, the dielectric induces a much smaller coupling to one of the two contacts, such that the first two resonances in transmission contain a significant background contribution as a decaying wavefunction in the barrier seen in the solution from the left (Figure 12b, positive polarization) and in the solution from the right (Figure 12d, negative polarization). Analyzing Figure 12c,e, one can notice that the overlap between the solution from the left and that from the right is almost zero for the first two resonances in transmission. The next two resonances in transmission, however, can be distinguished from the background, the normalized eigenvectors of A L and A R being quite similar, yet the corresponding wavefunctions |Ψ L D and |Ψ R D having quite different amplitudes. The difference in amplitudes starts to close in as we move to higher energies, since at sufficiently higher energy the asymmetry of the barrier will not play any major role in the electron transmission. As a final comment, all four resonances shown here participate to temperature-activated transport ( Figure 7). However, the main contribution to temperature-dependent transport is given by the first three resonances, even though the couplings to the contacts of the first two resonances are not as strong as the couplings of the third one, which is further apart in energy.

Conclusions
In conclusion, the semi-empirical model of electrostatic and NEGF calculations can provide a detailed picture of electron transport in systems with FTJs. It treats on equal footing several transport mechanisms that are usually invoked when studying these devices, such as direct tunneling, thermionic emission, and Fowler-Nordheim tunneling. This feature of NEGF allows us to assess in detailed form the role of temperature in

Conclusions
In conclusion, the semi-empirical model of electrostatic and NEGF calculations can provide a detailed picture of electron transport in systems with FTJs. It treats on equal footing several transport mechanisms that are usually invoked when studying these devices, such as direct tunneling, thermionic emission, and Fowler-Nordheim tunneling. This feature of NEGF allows us to assess in detailed form the role of temperature in electron transport, and how temperature affects the TER ratio. We have found that for simple or composite barrier BaTiO 3 -based FTJs, the transport through temperature activated resonant (Gamow-Siegert) states may become dominant, affecting both the conductance and the TER ratio; thus, the more resonances are activated and the stronger their coupling to the contacts, the more powerful the effect of temperature. The effect of temperature is obvious in thicker FTJs, since the resonant states-especially those above the barrier-are closer one another, and hence more of them may participate in the transport. More insight into this phenomenology may be acquired by calculating the wavefunctions at any given energy. We show that this is possible for NEGF calculations of single-band transport. Thus, the transport by direct tunneling takes place through states whose wavefunctions have a decaying shape in the barrier. These states belong to the background generated by all resonant states. Furthermore, where the resonances are strong and well-separated, the wavefunctions show confinement character in the barrier. In the intermediate regime, where the resonances are weak, the confining character of the resonance competes with the decaying character of the background. This analysis is valid not only in the linear regime, but also for larger bias voltages characteristic to non-linear regime and relevant to applications; the only change is a rigid shift in the screening charge and in the potential energy of the right contact, as can be seen from Equations (7) and (8). Lastly, we suggest that these results may be usefull in the optimization process of the FTJ design for various applications, particularly at high temperature.
where m * L is the effective mass in the left contact and Equation (A2) is discretized and the discrete form of (A2) can be mapped into a tightbinding Hamiltonian. Thus the continuous variable x is transformed in the discrete version ∆·n, where ∆ is the discretization step and n is an integer running from −∞ to ∞. We define a localized orbital n, R L with ∆·n and R L as longitudinal and transverse positions. Moreover, we construct transverse Bloch orbitals as a sum over N localized orbitals in the transverse plane In this Bloch basis the Hamiltonian has the following expression n, k H n ,k = D n (k)δ n,n − t n,n δ n ,n±1 ,