Fabian
Ducry
*,
Jan
Aeschlimann
and
Mathieu
Luisier
Integrated Systems Laboratory, ETH Zurich, Gloriastrasse 35, CH-8092, Zurich, Switzerland. E-mail: ducryf@iis.ee.ethz.ch
First published on 19th May 2020
Following the emergence of novel classes of atomic systems with amorphous active regions, device simulations had to rapidly evolve to devise strategies to account for the influence of disordered phases, defects, and interfaces into its core physical models. We review here how molecular dynamics and quantum transport can be combined to shed light on the performance of, for example, conductive bridging random access memories (CBRAM), a type of non-volatile memory. In particular, we show that electro-thermal effects play a critical role in such devices and therefore present a method based on density functional theory and the non-equilibrium Green’s function formalism to accurately describe them. Three CBRAM configurations are investigated to illustrate the functionality of the proposed modeling approach.
The challenges encountered in the scaling of MOSFETs have paved the way for the development of alternative technologies that could complement or even replace these systems in selected cases. This trend is most notable in the data storage segment where emerging non-volatile memories (NVM) are receiving a lot of attention, both from academia and industry.5 Two of the most promising NVM candidates intimately rely on amorphous and disordered materials for their operation, specifically valence change (VCM) and conductive bridging random access memory (CBRAM) cells.6,7 Such structures are sometimes collectively named resistive random access memories (RRAM). A large number of material combinations and device stacks have been demonstrated to exhibit the desired memory behavior. Each individual configuration comes with its own characteristics such as switching speed, retention time, or turn-on voltage. A huge design space must be explored to give rise to a cell that fulfils specific criteria. Simulation and modeling efforts can support this process and assist in the design of optimal RRAMs, thus reducing the required effort. Moreover, in the case of NVM not all operating mechanisms are quantitatively understood, and the origin of certain effects remains in debate.7 An in-depth comprehension of the underlying physics is, however, crucial to enhance the reliability of RRAMs, which is one of the key challenges to address in future research.8 Simulations can provide valuable insight into the processes governing the device operation.
In all of the situations mentioned above any dedicated modeling effort should account for both the atomic configuration of the structures and quantum mechanical effects. Therefore, a quantum transport simulator is essential to reveal the device characteristics and evaluate the performance of nanoscale components incorporating amorphous layers and interfaces between different materials. It could be used to support on-going experimental activities, provided that it satisfies specific requirements: any material combination should be possible, spatial resolution of the structure down to single atoms is necessary, and atomic disorder should be properly treated, not within a so-called virtual crystal approximation.9 Although ultra-scaled devices often operate close to their ballistic limit, many experimental features can only be explained by simultaneously taking electronic as well as thermal and coupled electro-thermal properties into account. This is the case for self-heating, which may be responsible for RRAM failures at high current densities.10 Density functional theory (DFT)11 meets these requirements and has found widespread use in electronic structure calculations. To drive devices out-of-equilibrium and observe electrical or thermal current flows the non-equilibrium Green’s function (NEGF)12 formalism must be coupled to DFT. Together, DFT and NEGF build a powerful and versatile ab initio quantum transport simulation framework.13–15
Memristive metal–insulator–metal cells that implement VCM and CBRAM rely on the random relocation of atoms to change their active switching layer from insulating to conductive and vice versa.7 The switching layer and the mechanism responsible for the insulating to conductive transition are different in both technologies, but they critically depend on atomic disorder in the core of the device. In addition to the inherently non-periodic nature of amorphous phases, the presence of oxygen vacancies (VCM) or metal interstitials (CBRAM) increases the disorder and profoundly affects the functionality of the memory cell. As such, VCM and CBRAM are particularly attractive as test beds to demonstrate the strength of any proposed modeling approach: they take advantage of disordered switching layers, where electrons must flow through complex material interfaces, electro-thermal effects are important, and simulation domains composed of thousands of atoms must be generated to capture the physics. Hence, this review presents the current state-of-the-art in advanced device modeling, starting from the available simulation techniques up to concrete applications where ab initio quantum transport is needed. The associated challenges are addressed in the following sections of this paper.
The assembly of atomistic models for memristive devices is discussed in Section 2.1, highlighting the need for ab initio (transport) techniques, which are briefly summarised in Section 2.2. Ballistic electron and thermal transport equations are described in Section 2.3, while Section 2.4 introduces the coupling between electrons and phonons to account for the self-heating effects arising at large current magnitudes. Section 3 illustrates the power of the reviewed techniques with concrete examples, starting from the simulation of the operation of a single CBRAM cell in Section 3.1. Section 3.2 presents the impact of the oxide thickness on the performance of ultra-scaled CBRAM cells. Finally, Section 3.3 examines the impact of different material combinations on the characteristics and spatial distribution of the current. Selected challenges regarding the accuracy of the physical modeling and computational complexity are highlighted in Section 4 before conclusions are drawn in Section 5.
While CBRAM and VCM cells operate in a similar fashion, they differ with respect to their material stack, filament properties, and underlying physical switching mechanism, as depicted in Fig. 1. The insulating layer in VCM devices usually consists of SrTiOx or binary transition metal oxides such as HfOx, TiOx, or TaOx.16 While it is widely accepted that the switching mechanism of VCM cells is based on the creation and displacement of oxygen vacancies, its exact nature is still in debate.6 It was, for example, proposed that oxygen is extracted from the oxide at the interface to the active electrode through the application of an external voltage, leaving oxygen vacancies behind. The oxygen atoms are either incorporated into an oxygen reservoir at the electrode,17 or recombine with each other forming gaseous O2.18 Therefore, metals with a high oxygen acceptance capability like Ti or Ta are frequently used as active electrodes.7 Meanwhile, the doubly positively charged oxygen vacancies migrate through the oxide layer in the direction of the electric field until they reach the passive, typically electro-inactive electrode where they are reduced. Based on the clustering of the neutral vacancies, a conductive filament is generated within the oxide layer.7 By reversing the voltage polarity, the oxygen vacancies are oxidized again and migrate back towards the active electrode where they recombine with oxygen atoms. This process ruptures the filament and resets the system into its OFF-state.16
For HfOx- and SiOx-based VCM cells, an alternative model has been proposed in which oxygen vacancies are formed in the bulk rather than at the interface, creating anti-Frenkel defect pairs consisting of an oxygen interstitial and an oxygen vacancy.19–21 The oxygen vacancies are assumed to be immobile forming an oxygen-deficient filament. The interstitial oxygen ions, on the other hand, would be mobile and migrate under the applied electric field towards the active electrode forming an interfacial oxide layer. However, it was later shown that anti-Frenkel defect pairs are not stable in bulk hafnium oxide.17 Furthermore, it is energetically more favorable to form oxygen vacancies at HfOx interfaces rather than in the bulk.22
In addition to filamentary type switching, VCM is also proposed to operate in a volume-type mode where the total amount of oxygen vacancies stays constant.16 It relies on the redistribution of oxygen vacancies next to the active electrode during the switching process. If the oxygen vacancy concentration at the electrode–oxide interface is modified, the electrostatic barrier is altered according to the Schottky effect.23,24
The filament in CBRAM cells consists of metal atoms dissolved in the insulating layer. A broad range of materials have been reported as possible candidates for the insulating layer.6 Apart from chalcogenides, amorphous oxide films of SiO2, HfO2, ZrO2, or Al2O3 have attracted a lot of attention due to their compatibility with existing CMOS technologies. Typical materials for the electro-active electrode are Ag and Cu, whereas Au, Pt, or W are preferred for the electro-inactive side. In CBRAM, the application of an external potential triggers an oxidation of the metal atoms at the interface between the active electrode and the insulating layer.6 This oxidation generates metal cations, which migrate along the applied electric field towards the passive electrode where they are reduced. Eventually, a filament of neutral metal atoms forms and grows towards the active electrode.25 Once it bridges the oxide and connects the two electrodes, the electrical resistance drops by several orders of magnitude. The filament can be disrupted by applying a voltage of opposite polarity, leaving metallic ions that only partially span the switching layer. The remaining atoms can then be reused to create a new filament in subsequent SET processes.6,7
Although recent experimental studies26 shed light on the switching principle of RRAMs, the precise mechanisms that control the transition from the OFF- to the ON-state, as well as the nature of the conducting path, are still under intense investigation. Continuum models27 (Fig. 2(a)), in which partial differential equations describe the atomic motions (drift and diffusion), can very accurately reproduce and explain experimental data such as the I–V characteristics during a switching cycle28,29 or the conductive filament life time,30 at low computational cost. However, their efficiency depends on the availability of a large set of material parameters, that must be determined one way, e.g. from higher-order simulations, or another, e.g. through fitting. In addition, any information about the actual atomic configuration is lost, which might become an issue when the stochastic relocation of a few atoms can change the electronic current by several orders of magnitude.31 For instance, doping atoms reduce the formation energy of oxygen vacancies in VCM cells. The oxygen vacancies tend to then accumulate around the dopants, affecting the formation of the RRAM filament and impacting the overall device performance; both the electronic structure, which governs the conductivity of the switching layer, and the switching kinetics are strongly influenced.32 Therefore, atomistic models are needed to highlight the mechanisms underlying the switching behavior of RRAM cells.
One such example is kinetic Monte Carlo33 (KMC) (Fig. 2(b)), a simulation approach that allows us to generate atomistic filament structures and to link them to continuum methods.34 The KMC simulation box is typically discretized into a grid with quadratic tiles representing the atomic positions. The edge length of a square (2-D) or cube (3-D) corresponds to the hopping distance of the filament forming species. In a KMC model, all relevant processes occurring in a RRAM cell, e.g. oxidation and reduction, adsorption and desorption, nucleation as well as ionic migration within the insulating layer or along interfaces, are described by rate equations obeying Arrhenius-type behavior.27 Each rate equation depends on the energy barrier that the specific reaction has to overcome, for instance the activation energy for ionic diffusion or the one for oxidation. Since the activation energy can be lowered by an applied voltage, the processes can be exponentially accelerated by increasing either the applied voltage or the temperature. The rate of each individual process is first calculated and stored in a table. At each step of the KMC algorithm, the event to be executed is randomly chosen based on the occurrence probability of the various processes. After each event, the atomic configuration is potentially modified until a stationary state is reached.
KMC has been successfully applied to grow and dissolve filaments with an atomic resolution in both CBRAM35,36 and VCM37 cells, with excellent agreement with experimental data. Despite valuable insights into the filament dynamics, structures generated with this method suffer from multiple limitations. First of all, most KMC models are two-dimensional, although three-dimensional implementations have been recently demonstrated.36,37 Second of all, cubic grids have difficulty treating amorphous structures, and materials with a non-cubic lattice are only approximately represented. Lastly, the switching layer, contrary to the filament, is described as a continuum, rather than at an atomic level. Thus, more advanced models are needed that can enhance the spatial resolution of KMC and better account for the broad range of material properties encountered in RRAMs.
Classical molecular dynamics (MD) based on force-field (FF) approaches (Fig. 2(c)) meet these requirements and can capture the detailed atomic structure of both the filament and the insulating layer as well as their dynamics. In such simulations, a parameter set describes the different types of atoms and their interactions. Dihedral angles, torsion, and chemical bonds are used to calculate the potential energy of a given system from which the forces acting on each atom can be derived. The parameters are fitted to reproduce reference data from experiments, quantum mechanical calculations, or both.38 The obtained forces are then used to determine the trajectories of the atoms based on Newton’s laws of motion. To model the growth and dissolution of a filament through the switching layer of an RRAM cell, simulation domains containing thousands of atoms must be constructed.39 Additionally, time spans of several nanoseconds must be considered to model a full switching cycle.10 FF-based molecular dynamics achieves that at reasonable computational cost.
Elaborate schemes are needed to construct suitable amorphous structures and interface them with metallic electrodes. For example, a melt-and-quench approach40 can be used for that purpose. Starting with a chunk of crystalline oxide or randomly placed atoms, MD is performed for several hundreds of picoseconds at a temperature above the melting temperature of the oxide. Then, the melt is quenched to 300 K with cooling rates in the order of −15 K ps−1.39 Post-quench annealing at room, or slightly elevated, temperature can be beneficial to eliminate coordination defects and reduce the stress inside the amorphous structure.
The first atomistic simulation of a complete CBRAM switching cycle was demonstrated by Onofrio et al.39 using a so-called reactive FF method. In contrast to traditional FFs, reactive force-fields such as ReaxFF41 are able to describe the formation and breaking of bonds and therefore to model chemical reactions. They rely on the charge equilibration formalism,42 which provides environment-dependent partial atomic charges. Onofrio et al.39 extended ReaxFF MD simulations with electrochemical dynamics including implicit degrees of freedom (EChemDID). This allows for the description of electrochemical reactions driven by the application of an external voltage. The effect of an electrochemical potential difference, ΔΦ, applied between two electrodes was modeled by altering the electronegativity, χ, of the atoms on the left and right electrode to χ + ΔΦ/2 and χ − ΔΦ/2, respectively.
A much simpler model of a conductive filament can be obtained by manually inserting metal atoms into the amorphous insulating layer instead of explicitly growing a structure. A shape must be defined and all atoms within it are replaced by metal. Subsequently, the obtained system is annealed using MD. It can be used as the starting point for reactive MD under an electric field. However, models relying on continuous, rather than localized, electric fields have not lead to realistic filament morphologies so far, at least not for complete ON–OFF switching cycles.43
A similar approach can be applied to VCM cells to create an oxygen vacancy filament in a layer of transition metal oxide. Starting from pristine amorphous oxide, oxygen atoms are removed by visual inspection resulting in a transition metal-rich region.44 Depending on whether the low- or high-resistance state should be modeled, a continuous or disrupted filament can be created.
The parameterization of force-fields is often tailored such that the processes of interest are accurately described, whereas less relevant phenomena are not well accounted for. Therefore, the usage of force-fields to perform MD in complex systems such as RRAM cells, where many different subprocesses are encountered, can result in misleading behavior. A higher level of accuracy can be achieved using ab initio molecular dynamics (AIMD) where the forces acting on each atom are derived from DFT. The latter is a quantum mechanical modeling method that can describe the electronic structure of any given atomic configuration without the need for fitting parameters (Fig. 2(d)). However, the high computational demand limits the time range accessible by AIMD to a few picoseconds and the system size to a few thousand atoms.45
To benefit from the advantages of FFs and DFT, both methods can be combined.40 First, atomistic filamentary type RRAM structures are created using FF approaches as outlined above. Then, the structures are relaxed and optimized using AIMD before a variety of physical properties such as the evolution of the electronic density of states (DOS),45 the activation energy of ion diffusion46 as well as the nucleus formation energy47 in CBRAM cells, or the formation energy of oxygen vacancies in VCM cells32 are extracted with the help of DFT. Due to the disordered nature of the structures, calculating meaningful physical properties can only be achieved by averaging over an ensemble of independent measurements.45 For each of them, the corresponding Hamiltonian can be used in quantum transport simulations to compute electrical and thermal currents as well as self-heating effects, which is the topic of the next section.
The NEGF formalism offers a powerful framework to calculate the non-equilibrium properties of quantum mechanical systems.50 It is widely used to perform quantum transport (QT) simulations.51,52 This approach to non-equilibrium statistical mechanics is based on the work of Kadanoff and Baym,12 and Keldysh.53 It requires a description of the electronic structure of the system under study in the form of a Hamiltonian matrix. A number of strategies have been proposed and tested to construct this quantity. Notable examples are the tight-binding (TB)54 and effective mass approximation (EMA),55 as well as approaches relying on first-principles concepts, where the Hamiltonian matrix is obtained from DFT calculations.13,14,56 Such a coupling of NEGF and DFT to perform ab initio QT calculations was introduced by Lang,56 where the representation of the device was based on DFT and the electrodes modeled using a jellium approximation. Fully atomistic simulations were proposed by Taylor et al.13 featuring an atomistic representation based on DFT for both the device region and the contacts. Since these pioneering examples, several packages capable of treating quantum transport from first-principles have been developed.57–62 Some of them are freely available, others commercially available.
The majority of DFT + NEGF calculations are performed in the ballistic limit of transport, where the energy of each particle is conserved throughout the simulation domain. The effect of inelastic interactions can naturally be incorporated in NEGF through the use of scattering self-energies.63 Besides pure electrical or thermal transport, coupled electro-thermal simulations can be perturbatively carried out through the self-consistent Born-approximation (SCBA).15,64 Owing to the large computational burden induced by the SCBA, such simulations are typically restricted to small systems or require large computational resources. Furthermore, calculating the phonon properties and the electron–phonon coupling from first-principles is a challenging task.65,66 Nevertheless, such calculations have been applied to a wide range of nanoscale devices going from 2-D field-effect transistors (FETs),65 FinFETs,66 or CBRAM cells31 to the modeling of inelastic electron-tunneling spectroscopy (IETS) in molecular junctions.67,68
It should be noted that the SCBA is not the only possibility to account for electro-thermal effects. Lowest order-expansion techniques have been used as well. They have a lower computational burden, but at the cost of additional approximations.67,69 Due to its perturbative nature, SCBA may fail to converge in the presence of strong electron–phonon coupling. An exact, but also computationally more expensive, technique capable of treating such systems is the hierarchical equations of motion.70
The DFT + NEGF framework is not limited to steady-state simulations. It is capable of delivering insight into time- and frequency-resolved quantum transport phenomena as well.71 Due to the heavy computational burden of these simulations they have only recently seen a rise in popularity. While NEGF is frequently used for electrical and thermal QT simulations the formalism can be applied to any quasi-particle that follows the laws of quantum mechanics. It has, among others, also been employed to study the adsorption on metal substrates72 and solid-plasma surfaces.73
The remainder of this paper is dedicated to the description of ab initio electro-thermal QT calculations within the SCBA. The following sections describe the coupling of DFT with NEGF for the case of electrical and thermal transport. Subsequently, the coupling of electrons and phonons via scattering self-energies is explained.
ĤKS|Ψ〉 = E|Ψ〉 | (1) |
(2) |
It contains two terms, first the kinetic operator −∇2/2 and then the effective potential
(3) |
The three terms above correspond to the Hartree, the exchange–correlation, and the external potential with ρ(r) denoting the charge density at position r. The Hartree potential includes the electron–electron Coulomb repulsions. The second term in eqn (3), VXC, accounts for all electron many-body effects, namely electron exchanges and correlations, whereas the influence of external voltage sources are cast into the last one, Vext. DFT takes advantage of the Born–Oppenheimer approximation, which enables a separate treatment of the valence electrons and atomic cores. Electrons freely move in the potential induced by these static cores, which are also described by the external potential.
Multiplying eqn (1) with 〈Ψ| gives rise to the Hamiltonian and overlap matrices
H = 〈Ψ|ĤKS|Ψ) and S = 〈Ψ|Ψ〉 | (4) |
H·ψ(E) = ES·ψ(E), | (5) |
The computational efficiency of solving matrix equations can be directly related to the number of non-zero elements and the sparsity pattern. To maximize the sparsity of H and S a suitable basis must be selected to expand |Ψ〉. Localized basis sets such as Gaussian-type orbitals (GTO)75 are ideal for that as they produce sparsely populated, banded matrices. In contrast, plane-waves, which are popular for electronic structure calculations, lead to dense matrices. While not impossible,76,77 the use of plane-waves is relatively scarce in quantum transport calculations, in particular for large systems. If necessary, plane-waves can still be localized, with the help of e.g. Wannier functions.78 When employing a localized basis set, |Ψ〉 becomes
(6) |
The valence electrons of atom n are expanded in l(n) basis functions ϕln centered at position rn. The number of basis functions per atom l may vary for different chemical species. The cln(E) are the occupation coefficients of the respective basis function. The wave function is then the sum over all individual basis functions of each atom weighted by the occupation factor. With this choice of basis expansion ψ(E) becomes a vector containing all coefficients cln(E), and the size of the matrices H and S is the sum of l(n) over all atoms. If the ϕln are orthogonal to each other, S is the identity matrix and eqn (5) becomes a regular eigenvalue problem (EVP). In most cases the overlap between localized basis functions is non-zero and a generalized EVP must be solved.
It should be noted that the H and S matrices are not unique to DFT. They can also be created based on other methods such as tight-binding (TB),54 where Löwdin orbitals are parameterized to reproduce experimentally measured, or DFT, band structures. The Hamiltonian is not necessarily an atomistic quantity either, it could be expressed in the effective mass approximation (EMA) on a discretized grid.55 The transport equations presented in the next paragraphs apply to Hamiltonian matrices obtained with these methods as well. The continuum nature of EMA and the required TB parameterization, however, make both methods ill-suited to deal with disordered structures. The major difficulty in TB models lies in the derivation of a parameter set that accurately captures amorphous phases or defects such as vacancies and interstitials.
In contrast to electronic structure calculations, which are typically restricted to the ground state of a system, electrons in device simulations must be able to enter and leave an open domain so that a non-equilibrium current can flow. External potentials are applied to contact regions to drive a device out of equilibrium. As a consequence, the boundary conditions applied to eqn (5) require special attention. Whereas periodic (PBCs) or closed boundary conditions (CBCs) are typically used in DFT, open boundary conditions (OBCs)79 are at the core of quantum transport investigations. In OBCs the contacts and the device region are first treated separately, as illustrated in Fig. 3(a). Each contact is modeled as a semi-infinite lead in thermal equilibrium, with a flat electrostatic potential, as shown schematically in Fig. 3(b). It should be represented by at least two identical blocks of atoms. These blocks correspond to the first and last unit cell of the central (device) region. The leads serve as launching pads for electrons or as collectors. They are connected to the device through so-called retarded boundary self-energies ΣR,B(E) that must be introduced into eqn (5). Additionally, an injection vector Inj(E) acts as a source term to model the incoming electrons. OBCs are not limited to two-terminal systems, but can be readily generalized to structures with multiple leads.79
In the presence of OBCs eqn (5) takes the following form
(ES − H − ΣR,B(E))·ψ(E) = Inj(E). | (7) |
This equation must be solved for all discrete energies belonging to the interval of interest, which extends over the Fermi energy of both contacts. In ballistic and coherent simulations, i.e. in the absence of scattering with energy relaxation, all E’s are independent of each other and eqn (7) directly yields ψ(E). Such an approach is known as the quantum transmitting boundary method (QTBM).80 It is computationally attractive as it involves the solution of sparse linear systems of equations with multiple right-handed-sides, but it does not lend itself naturally to the simulation of dissipative transport.
Alternatively, eqn (7) can be rewritten in terms of Green’s functions (GFs) as
(ES − H − ΣR,B(E))·GR(E) = I, | (8) |
G≶(E) = GR(E)·Σ≶,B(E)·GA(E), | (9) |
ψ(E) = GR(E)·Inj(E). | (10) |
Eqn (8) and (9) are known as the non-equilibrium Green’s function (NEGF) formalism.12 Intuitively, the lesser (greater) boundary self-energy, ΣB,≶(E), indicates the probability that a state gets filled (ΣB,<(E)) or emptied (ΣB,>(E)) through interactions with the contact. The off-diagonal elements of the lesser and greater GF, G≶(E), describe the correlation between the involved basis functions, whereas the diagonal entries contain the probability that a state n is occupied (G<nn(E)) or unoccupied (G>nn(E)). Therefore, it is not necessary to convert back the results of eqn (8) and (9) to a wave function ψ(E). All observable quantities such as the charge density ρ(r) at position r and the electrical current Id can be directly derived from selected entries of G≶(E). The latter can be computed efficiently with an iterative algorithm called the recursive GFs (RGF) algorithm.81 Nevertheless, for ballistic simulations this approach is computationally more expensive than QTBM. The strength of NEGF comes from its natural integration of scattering mechanisms through self-energies, as will be introduced in Section 2.4.
Regardless of the transport type, the charge density ρ(r) can be computed as
(11) |
(12) |
Injecting electrons into the device domain drives it out of equilibrium, which changes the distribution of electrons and modifies the charge density ρ(r). This in turn affects the Hamiltonian H through the first two terms in the effective potential Veff(r) in eqn (3), giving rise to a mutual dependence of eqn (2), (8) and (9). It must be resolved in a self-consistent manner until ρ(r) is converged.13
Fully coupled NEGF + DFT simulations come with a heavy computational burden, even in the ballistic case. This can be somewhat alleviated by assuming that the charge density only affects the electron–electron repulsions and by neglecting the change of exchange and correlation. By solving Poisson’s equation
(13) |
(14) |
After convergence of the selected self-consistent loop, physical observables can be extracted. In ballistic simulations the function T(E) describes the transmission probability of an electron from the left to the right side of an open system (or vice versa) at energy E.63 It is calculated according to
T(E) = ΓL(E)GA(E)ΓR(E)GR(E). | (15) |
The broadening function ΓC(E) of contact C depends on the boundary self-energy ΣBC(E) and is defined as
ΓC(E) = i(ΣB,RC(E) − ΣB,R†C(E)), with C = R or L. | (16) |
(17) |
(18) |
(19) |
This formulation of the electrical and energy currents, Id and IdE,el is more general and holds even when no transmission function can be defined, i.e. in the presence of a dissipative scattering mechanism.
Phonons are quasi-particles that arise from the coupled motion of atoms around their equilibrium position.55 An illustration of such a wave of coupled atomic motions is given in Fig. 5(a) for a 1-D wire. What is represented is an excited state of the lattice whose amplitude is related to the crystal temperature. To mathematically describe phonons, the total energy, Etot, of a perturbed atomic system with equilibrium energy E0 should be considered. The displacement of atom m along the cartesian direction i is labeled dim. In the harmonic approximation, Etot is expanded in a Taylor series up to the second order of the displacement
(20) |
In an equilibrium configuration the system resides in a (local) minimum so that dEtot/dim, the force acting on atom m in direction i, vanishes from eqn (20). The second-order term
(21) |
Φ·μ(ω) = ω2μ(ω). | (22) |
Eqn (22) is the phonons equivalent of eqn (5) for electrons. To derive the thermal NEGF equations, OBC are introduced into eqn (22). They can be constructed following the same prescriptions as for electrons illustrated in Fig. 3(a). Their calculation takes either the form of an eigenvalue problem86 or of a complex contour integral.87 By incorporating the resulting phonon boundary self-energies ΠB into eqn (22) and by transforming the wave function expressions into GFs, we obtain the following system of equations to solve:
(ω2I − Φ − ΠR,B(ω))·DR(ω) = I | (23) |
D≶(ω) = DR(ω)·Π≶,B(ω)·DA(ω). | (24) |
The quantities in the equations above are the phonon GFs (DR, DA and D≶) and the boundary self-energies (ΠR,B and Π≶,B). The labeling conventions for retarded, advanced, lesser, and advanced remain the same as in the previous section. Eqn (23) and (24) must be solved for all phonon frequencies of interest. The phonon density and current are then derived from D< using similar expressions as for the electrons. Besides, the energy current carried by phonons can be computed as
(25) |
While eqn (23) and (24) have the same form as the equations for electrons, there is an important difference. Namely, eqn (22) does not depend on the phonon population, therefore eliminating the need for an iterative solution process and reducing the computational cost as compared to electrons. The procedure involved in atomistic thermal transport simulations is depicted in the right-hand side of Fig. 4.
Whereas DFT has become the most widely used method for electronic structure calculations, even in large systems, competing approaches exist to generate the dynamical matrix Φ in eqn (22).89 It can be directly produced from the total energy of a given system with density functional perturbation theory (DFPT).90 Alternatively, it can be observed that Fijmn in eqn (21) corresponds to the first derivative of the forces acting on each atom. It can therefore be calculated from finite differences through the frozen phonon scheme.85 Both DFPT and frozen phonons induce a large computational burden that make them impractical for large atomic systems, if performed at the ab initio level.
The frozen phonon approach is the method of choice for large, disordered systems because of its relatively low computational complexity. It relies on the evaluation of the first derivative of the force fim = dEtot/drim using forward or central differences67,85
(26) |
To couple the electron and phonon populations, the energy of the fermionic and bosonic system must be considered. It can be described by the total Hamiltonian
Htot = H + Hph-kinetic + Hph-harmonic + He–ph. | (27) |
The first term corresponds to the electron Hamiltonian from eqn (5), while the second and third ones are captured by the dynamical matrix Φ. The last Hamiltonian contains the interaction between electrons and phonons. It is treated perturbatively and cast into the scattering self-energies ΣTS and ΠTS of type T ∈{R, A, <, >} for electrons and phonons, respectively. The lesser and greater components can be written as84
(28) |
(29) |
In eqn (28) and (29), all G≶, Σ≶, and ∇H blocks are matrices of size Norb × Norb with the summation over neighbor atoms omitted for brevity. The superscripts i and j denote the entries in ∇H, the phonon GFs, and self-energies corresponding to the cartesian coordinates i and j ∈{x, y, z}. The strength of the electron–phonon coupling is determined by ∇iH, which represents the derivative of the electron Hamiltonian with respect to the displacement of the atoms along the direction i. It thus couples the lattice dynamics created by the phonons to its electronic response. The retarded scattering self-energies ΣR,S and ΠR,S can be derived from Σ≶,S and Π≶,S. Very often, their real part is neglected for simplicity.63
To give an intuitive interpretation of eqn (28) and (29), we first recall that the diagonal elements of the lesser GFs, G<(E) and D<(ω), indicate whether a state at energy E is occupied by an electron and the number of phonons that occupy a state at energy ℏω. The same elements of the greater GFs, G>(E) and D>(ω), determine whether a state is unoccupied or the number of free states at the energy E and ℏω. A specific transition is only possible if an (un-)occupied electron state is available at an energy ℏω above or below the state of interest, as illustrated in Fig. 6(a). If a scattering event is allowed, the likelihood of in-scattering, i.e. an empty state at energy E, G>(E), gets filled is proportional to the lesser scattering self-energy Σ<(E). Such a process can happen through either phonon emission or absorption. An electron at energy E ± ℏω (G<(E ± ℏω)) emits (+) or absorbs (−) a phonon with energy ℏω (D>(ω) for emission, D<(ω) for absorption) and changes its energy to E. The out-scattering probability is given by Σ>(E). An occupied state at energy E, G<(E), gets emptied to E ∓ ℏω (G>(E ∓ ℏω)) by emission (−) or absorption (+) of a phonon (D>(ω) and D<(ω)). A similar interpretation can be made for the phonon scattering self-energies Π≶(ω). They refer to the probabilities that an unoccupied (D>(ω)) or free (D<(ω)) state gets filled (Π<(ω)) or emptied (Π>(ω)) when an electron transitions from one state to the other through phonon emission or absorption.
The diagonal entries of the scattering self-energies describe local interactions, that is, the electron remains on the same atom during the process. The off-diagonal elements, on the other hand, account for non-local transitions where the electron does not only change its energy, but also its position in real space. Non-local scattering events are numerically difficult to handle. In our approach they are neglected for electrons and only close-neighbor interactions are taken into account for phonons. This is necessary to ensure energy conservation in our NEGF calculations. By scaling the strength of the local entries of Σ≶, the influence of the non-local events can be indirectly modeled.93
The derivative of the Hamiltonian matrix from eqn (4) with respect to a perturbation in real space, ∂r, is given by
(30) |
The last two terms appear because the basis changes with perturbation.67 The expression in eqn (30) can be computed at various levels of accuracy. The simplest picture considers the derivative of the Hamiltonian with respect to bond stretching between two neighboring atoms.94 This scheme can be expanded to include more harmonic terms, e.g. bond angles. An alternative approach, similar to the calculation of the dynamical matrix, helps to determine ∇H with respect to the displacement of each individual atom. This is the most accurate method as it does not rely on any assumption regarding the nature of the bonds or angles connecting two atoms, but it comes at the price of increased computational cost. While the bond stretching technique only requires us to perform one additional ground-state DFT calculation under hydrostatic strain, typically 0.01% to 0.1%, the second method is even more expensive than computing the dynamical matrix from first-principles. Because of this, the bond stretching approach is the preferred one for large systems.
The chosen hydrostatic stress is applied to the atomic configuration by increasing the simulation box, thus stretching the bond rmn between atoms m and n without affecting the angle between them. If the change in basis functions induced by the stress is small, the last two terms in eqn (30) vanish. The strained Hamiltonian Hs is then computed and the derivative, with respect to a change in the bond length, evaluated as
(31) |
(32) |
When the scattering self-energies are included in NEGF the equations for electrons become84
(ES − H − [ΣR,B + ΣR,S])·GR = I | (33) |
G≶ = GR·[Σ≶,B + Σ≶,S(G≶,D≶)]·GA, | (34) |
(35) |
For phonons we have
(ω2I − Φ − [ΠR,B + ΠR,S])·DR = I | (36) |
D≶ = DR(ω)·[Π≶,B + Π≶,S(G≶,G≷)]·DA, | (37) |
(38) |
The dependence of the GFs and self-energies on the energy E and frequency ω has been left out in the above equations and substituted by D≶ and G≶ for the scattering self-energies, to emphasize the interplay between the electron and phonon populations. Eqn (33) and (34) now depend on eqn (36) and (37) and vice versa. These two sets of equations must be solved iteratively until convergence is reached. The fulfillment of this property can be verified by looking at the electrical current, eqn (18), and the sum of the electronic and thermal energy currents, eqn (19) and (25). Both quantities have to be conserved along the transport axis of the investigated device, when the GFs do not vary anymore. This process is known as the self-consistent Born-approximation (SCBA). The system of equations to be tackled is depicted in Fig. 6(b).
After converging the electron and phonon densities, physical quantities can be extracted. In addition to the currents that are given by eqn (18) and (25), the lattice temperature is of particular interest to quantify the effect of self-heating. Different possibilities exist to assign a local temperature to individual atoms.84 In the so-called population approach the effective temperature Teffn of atom n is adjusted such that the Bose–Einstein distribution reproduces the phonon population of each individual atom, Neffn. The phonon density is derived from the GF
(39) |
(40) |
The considered atomic CBRAM structures are assembled with the melt-and-quench approach detailed in Section 2.1, and depicted in Fig. 7(a)–(c). The switching layer is made of amorphous SiO2 (a-SiO2) with a density of 2.20 g cm−3. Melting happens at 3000 K for 600 ps. It is followed by a cooling phase with a rate of −30 K ps−1 to reach 300 K. The MD simulations are performed with a ReaxFF force-field,39 as implemented in the QuantumATK 2017.1.62,96 Next, metal electrodes, Cu or Ag, are attached to the a-SiO2 layer and a conical filament is inserted by replacing the Si and O atoms with Cu or Ag ones. The obtained configurations are annealed with AIMD at 800 K for 3–4 ps to relax the stress induced by the insertion of the metallic filament.
Fig. 7 (a) The melt-and-quench procedure to generate samples of a-SiO2. A box with crystalline SiO2 or randomly placed Si and O atoms is melted at 3000 K. Subsequently, it is cooled at a rate of −30 K ps−1 and annealed at 300 K. The melt-and-quench is performed classically, the post-annealing and optimization with DFT. (b) The metal–insulator–metal structure of a pristine CBRAM cell. Metal electrodes are attached to a slab of a-SiO2 obtained with the procedure from (a). (c) Filamentary configuration of a CBRAM cell in the ON-state. A metal filament is inserted into the a-SiO2 from (b) by replacing all Si and O atoms within a cone and annealing the result with DFT. (d) Energy-resolved transmission function around the Fermi energy for the structure shown in (c) using two basis sets, DZVP (solid blue line) and SZV + 3SP (dashed red line). Subplots (b)–(d) are adapted from ref. 31. |
The DFT simulations are executed with the CP2K package.97 This tool employs a Gaussian-type orbital (GTO) basis set with the linear combinations of atomic orbitals (LCAO)74 method. Such a localized basis is suitable for subsequent NEGF quantum transport simulations. All metal atoms are represented by a double zeta-valence polarized (DZVP)98 basis set to construct the atomic structure, while a single zeta-valence (SZV) basis is used to calculate the H and S matrices required for eqn (5). The larger DZVP basis set is required to produce meaningful device geometries, but the small SZV set is accurate enough for transport simulations31 as illustrated in Fig. 7(d). The 3SP parameterization of Zijlstra et al.99 is employed for the Si and O atoms, both in MD and quantum transport. In conjunction with the GTO basis set the atomic cores and inner electrons are described by the pseudopotentials of Goedecker–Teter–Hutter (GTH).100 The exchange–correlation energy is approximated by the PBE functional.101 No k-point grid is used; all calculations are performed at the Γ-point only. After obtaining the final atomic structure, the Hamiltonian and overlap matrices in eqn (5) are created in CP2K. The electron–phonon coupling elements are computed with the proposed bond stretching scheme introduced in Section 2.4 and rely on hydrostatically strained Hamiltonians. This approach is computationally less intense than more elaborate ones, especially for large systems.
The dynamical matrix accounting for the thermal properties is determined with the frozen-phonon approach discussed in Section 2.3.2. Computing the forces required to construct Φ from first-principles is prohibitive for our systems that comprise more than 3000 atoms. Therefore, these calculations were performed using a ReaxFF41 force-field specifically parameterized39 to model the switching behavior of CBRAM cells.
The Hamiltonian, overlap, and dynamical matrices, as well as the derivatives of the Hamiltonian, are imported into the OMEN quantum transport simulator.86,88 The latter implements the NEGF-Poisson scheme for electrons, phonons, and their coupling as detailed in Sections 2.3 and 2.4. The NEGF simulations for electrons are performed in the low field limit, which assumes that the applied bias does not significantly alter the electron density within the considered domain. Thus, the Hamiltonian and overlap matrices from CP2K are not updated based on the non-equilibrium charge and Poisson’s equation is not solved, thereby considerably reducing the computational burden.
The current flowing through the filament configuration of Fig. 8(a) was computed both in the ballistic limit and under the influence of electron–phonon interactions. The resulting I–V characteristics are shown in Fig. 8(b). Ohmic behavior is revealed at low biases, i.e. a linear increase of the current vs. voltage. The electron–phonon limited current reaches about 69% of the value of the ballistic one at room temperature. The impact of scattering on the magnitude of the electrical current is rather low, which is to be expected from such a short device. The resistance, extracted as a linear fit to the I–V curve, is 48.1 kΩ in the ballistic case and 57.6 kΩ in the dissipative one. These values lie in the range of typical CBRAM ON-state resistances.6 The current field lines corresponding to an applied voltage of 0.2 V are rendered in Fig. 8(a). They confirm that the current follows the metallic filament and exhibits its largest density at the tip, i.e. at the thinnest part of the filament. The influence of electron–phonon scattering on the current can be best visualized in the form of a spectral plot representing the energy- and spatial location of the current given in Fig. 8(c). Electrons lose part of their energy when propagating through the device. Such energy dissipation is not possible in ballistic electron transport. Most of the energy loss occurs close to the tip of the filament or in the metallic contact attached to it. This fact agrees well with the observation that the device operates close to its ballistic limit. The electrons cross the oxide layer too fast to relax their energy therein. It should also be noted that a large portion of the electron population enters (leaves) the simulation domain with an energy below (above) the Fermi energy of the respective contact. This indicates that a lot of the power dissipation happens outside the simulated device, namely in the leads. The fraction of power that is consumed within the simulation domain is called the internal dissipation fraction and is labeled α. It can be computed as
(41) |
Fig. 8 (a) Atomic structure of a CBRAM cell containing a metallic filament. The copper atoms are represented by the gray spheres, silicon and oxygen by the orange ones. The current flowing through this configuration upon an external bias of 0.2 V is plotted with blue-red lines. The red color at the tip of the filament (right-hand-side of the oxide layer) indicates a large current density which decreases (dark blue) towards the left. (b) Ballistic (blue triangles) and electron–phonon limited (red diamonds) I–V characteristics of the device in (a). A linear fit to the data points is given with the blue dotted (ballistic) and red dashed (e–ph) lines. (c) Energy- and position-resolved current flowing through the cell in (a) in the presence of electron–phonon scattering. The largest current density is observed just below (above) the equilibrium potential in the left (right) contact, which is indicated by a white line. (d) Atomically-resolved temperature in the device in sub-plot (a) at an applied bias V = 0.2 V. The heating is more prominent in the copper atoms in the right half of the filament, marked by the red color. (e) Atomic temperature projected onto the x-axis (transport). The filament atoms are marked with the red diamonds, the oxide by the red dots. The vertical dashed lines indicate the electrode-device transitions. The black line refers to the average temperature over the cell cross-section (yz-plane). The maximum temperature is taken as the top of this black curve. (f) Maximum temperature and dissipated power as a function of the CBRAM device current. The simulated values are given as symbols. A quadratic fit of the data is superimposed as dashed (power dissipation) and dotted (maximum temperature) lines. Subplots (a)–(f) are adapted from ref. 31. |
Due to energy conservation, a reduction in electron energy must be compensated for by the creation of additional phonons corresponding to an increase of the lattice temperature. Generally, the temperature is an average quantity characterizing an ensemble of particles embedded within a reservoir. It is therefore a macroscopic property. Nevertheless, by relating the excess phonon population on each atom to an equilibrium temperature through the Bose–Einstein distribution as in eqn (39) and (40), self-heating effects can be mapped to the familiar scale of temperature. The heat map of the device from Fig. 8(a) at 0.2 V is rendered in Fig. 8(d)–(e) with a 3-D atomic resolution, and after projection onto the transport axis. The highest temperature of the device lies in the middle of the switching layer, as indicated by the red atoms. This does not correlate with the highest electron–phonon scattering rate, which can be found in the electrodes or at the oxide–metal interface. Evidently, there has to be a second mechanism involved in the heating process. This is identified as the heat extraction rate from the oxide layer, which can be cast into the thermal resistance Rth. We define the latter as
Tmax = T0 + RthαRI2 | (42) |
Fig. 9 (a) CBRAM cell with an oxide length of 3.5 nm. (b) Same as (a), but with an oxide shortened from the left to 2.4 nm. The atomic configuration of the right-hand extremity of the filament is preserved. (c) CBRAM cell with an oxide layer shortened to 1.6 nm in the same fashion as (b). (d) Filament resistance as a function of the SiO2 layer thickness for an external bias of 0.2 V. The diamonds mark the simulated values, the dashed line serves as a guide to the eye. (e) Average temperature along the transport axis of the three CBRAM cells with different SiO2 layer thickness. For convenience, the filament tips are aligned at 7.6 nm. The vertical dashed line indicates the right device-electrode boundary. (f) Internal dissipation fraction α (left axis) and thermal resistance Rth (right axis) as a function of the SiO2 layer thickness. Subplots (a)–(c) and (e) are adapted from ref. 10, the data for (d) and (f) are presented in the ESI of ref. 10. |
As expected, the electrical resistance only slightly depends on the oxide thickness with a 10% reduction when going from the 3.5 nm to the 1.6 nm of a-SiO2 layer, as can be seen in Fig. 9(d). This confirms the assumption that the current is mostly limited by the filament extremity, not by its length. Consequently, self-heating can be compared between the three CBRAM cells using eqn (42), which depends on the current and resistance. The temperature averaged over the device cross-section is displayed in Fig. 9(e) for these three structures at an applied bias of 0.2 V. It is apparent that self-heating diminishes with the oxide thickness. Apart from a higher temperature profile in the thicker oxide, it can be noticed that the peak temperature moves further away from the contact towards the middle of the a-SiO2 layer.
The magnitude of the self-heating depends on two factors according to eqn (42): (1) the amount of power that is converted to heat (αRI2) and (2) the efficiency at which the excess phonon population can be extracted from the active region into the leads (Rth). To determine the origin of the increased temperature in thicker devices, both effects are examined separately. The internal power dissipation factor α versus the oxide thickness is plotted in Fig. 9(f). It can be observed that this quantity decreases with the oxide thickness. Hence, fewer phonons are emitted within the shorter cell. This behavior can be related to the time electrons spend propagating through the filament. The shorter the time, the lower the probability that an electron can interact with the surrounding lattice and dissipate energy.
To examine the thermal resistance of the device, Rth is drawn as a function of the oxide thickness, also in Fig. 9(f). We find that Rth exhibits the same characteristics as the power dissipation, i.e. Rth is smaller in thinner oxides: the phonons emitted by electrons remain longer in thicker oxides because they need more time to escape. The fact that copper is an excellent heat conductor whereas SiO2 is not, explains this behavior. As a result, the heat distribution is different in the three structures. The maximum temperature, which is located at the tip extremity in the shortest device, moves to the middle of the structure when the oxide thickness increases.
As both α and Rth increase with the oxide thickness, a strong dependence of self-heating on the dimensions of the switching layer is observed in the CBRAM cells. If we assume that excessive heating is the main reason for device failure, the experimental findings of ref. 10 can be readily explained: devices with a shorter oxide layer can endure larger currents before reaching the breakdown temperature. Stated differently, at a given current magnitude, self-heating is more pronounced in CBRAM cells with longer filaments because more phonons are emitted and they have more difficulties leaving the amorphous oxide region and attain the metallic contacts.
Constructing a CBRAM cell where one electrode is different from the other requires special treatment due to the breaking of the structure symmetry along the transport axis. Two metals have different lattice constants in general, which creates a mismatch between the end of the right electrode (Ag in Fig. 10(a)) and the beginning of the left electrode (Ag, Pt, or W in our case). The electrodes can be strained to eliminate the lattice difference. To minimize its impact on the results, the metal supercells that build the left and right contacts should be chosen such that their dimensions match closely. To do that, a cross-section of 2.5 × 2.4 nm2 was selected. It ensures that the strain level does not exceed 1% for each electrode. The amorphous SiO2 does not have a lattice constant and hence does not suffer from strain. After setting up the Ag/a-SiO2/Ag structures, a conical Ag filament is inserted manually into the 2 nm thick oxide layer. The whole system is annealed with AIMD so that the filament displays an hourglass shape with a large cone on the counter electrode and a small one on the active electrode side, as shown in Fig. 10(a). The Pt and W electrodes are inserted after the atomic structure was optimized and annealed. One of the Ag electrodes is stripped away and replaced by either metal. Afterwards, the interface between the new metal and the oxide is optimized, while the bulk of the SiO2 and the right Ag electrode remain unchanged. In this way three material stacks with the same atomic filament are obtained. Next, the Hamiltonian and overlap matrices of each configuration are computed to perform ballistic electron transport calculations. Phonons and electron–phonon scattering are not considered here.
Fig. 10 (a) Atomic configuration of a CBRAM cell with two different electrodes and an hourglass filament. White spheres represent silver atoms, the black ones tungsten. The red dashed line delimits the shape of the filament with a large cone on the left and a small one on the silver electrode. (b) Energy-resolved transmission as a function of the energy for the structures illustrated in (a) with Ag (red dotted), Pt (blue dashed), and W (green dash-dotted) as counter electrodes. (c) I–V characteristics of the three devices in (a) with Ag (red diamond), Pt (blue circles), and W (green squares) as counter electrodes. Linear fits are provided for convenience (dashed lines). (d) Current magnitude on the individual silver atoms in the filament and on the surface of the active electrode (black squares). The error bars give a measure of the variance of the current between the three structures. The atoms belonging to the active electrode and the small cone attached to it are indicated by the red circles. Subplots (a)–(d) are adapted from ref. 95. |
The energy-resolved transmission function around the Fermi energy is very similar in all CBRAM cells. It is reported in Fig. 10(b). This means that the probability of electrons of being transmitted through the oxide layer is independent of the contact material. In ballistic simulations the current is directly related to the transmission through eqn (17). Therefore, the I–V characteristics also resemble each other, Fig. 10(c), and confirm that the device current is mostly unaffected by the choice of the electrode metal.
Other than the current magnitude, its distribution can also play an important role in determining the device operation. To asses this quantity, the current vector starting from each atom is computed. The expression for that can be inferred from eqn (18): instead of summing all flows, each individual contribution is recorded and stored in the form of a current vector field. The current magnitude passing through each individual filamentary Ag and electrode surface atom is plotted in Fig. 10(d). The dots refer to the mean value between the three simulations and the error bar measures the standard deviation. Two distinct regimes can be discerned: (1) in the small cone situated on the active electrode, there is no variation in the current between the different structures; (2) some Ag atoms situated inside the large cone attached to the counter electrode (Ag, Pt, W) show a large spread in current. From these results, it can be deduced that the nature of the {Ag, Pt, W}-Ag filament influences the current distribution, but not its magnitude. The latter is determined by the thinnest region of the filament. In other words, the filament resistance depends on its atomic morphology. The current density, on the other hand, is affected by the interface connecting the counter electrode and the filament. Therefore, simulations considering two identical metal electrodes correctly predict the device resistance, but not the spatial distribution of the current. Electro-thermal effects, for instance, strongly depend on the current density, not only on its magnitude. In future work, the atomically resolved temperature of the three investigated CBRAM cells should be calculated to identify possible failure mechanisms.
Approximations to the laws of quantum mechanics are the most challenging ones to address. They permeate the ab initio modeling and its most common form, DFT. The underestimation of the band gap is a known issue. Possible improvements include the use of hybrid functionals that can nowadays be applied even to large atomic systems.102 The HSE06 functional,103 for example, improves the description of the electronic structures, but it does not eliminate all problems. Further research towards more accurate and general exchange–correlation functionals is on-going.104 A related approximation mentioned in this review is the reduction of the calculations to the Γ-point. Because of the large size of the systems, a sparse k-point sampling is sufficient, but not optimal. With increasing available computational power and growing system sizes, the Γ-point limitation will become less prominent as the sampling of the momentum space implicitly increases with the simulation domain.
Improvements can also be brought on the phonon side where the dynamical matrix is often computed from force-fields for large systems. The main challenge to apply DFT is the huge number of force evaluations that are necessary. It scales linearly with the number of atoms, Natom, in addition to the typical O(Natom3) dependence of DFT itself. The fact that the response to the displacement of a single atom is localized could be exploited to develop a scheme where multiple, “far” away from each other, displacements are accounted for in a single force evaluation. This could considerably reduce the computational burden. The same considerations hold true for the calculation of the electron–phonon scattering elements. Instead of calculating them with a very basic model that only includes bond stretching and might neglect important effects the same methodology as for phonons could be used, with the same numerical challenges to overcome.
Within NEGF, generally, only local scattering events are accounted for. The off-diagonal entries of the scattering self-energies are discarded to reduce the computational complexity. To go one step further, the widely used RGF algorithm should be extended to produce the off-diagonal elements of the Green’s function, which are required to calculate other than diagonal components of the scattering self-energies. The computation of these entries becomes more and more difficult as one moves away from the diagonal. Improving the description of the scattering mechanism is definitively possible, but computationally very demanding. If experimental or reference data are available, it is preferable to scale the diagonal entries of the scattering self-energies with a constant factor to reproduce the desired targets.105 Note that additional effects such as anharmonic phonon–phonon interactions might be important. Including them would further increase the computational cost by slowing down the convergence of the SCBA iterations. Moreover, the matrix elements for anharmonic phonon–phonon scattering are related to the second derivatives of the forces acting on the atoms. Obtaining these derivatives is computationally even more challenging than the calculation of the dynamical matrix. Similar considerations as above could render these computations more affordable as long as small systems are investigated.
In terms of parallelization of the workload on high performance computing (HPC) resources solving the NEGF equations in the ballistic limit is an optimally suited problem as all energies (and momentum points in periodic structures) are independent from each other. As such, they can be naturally distributed to different computing units with little to no communication overhead in between. This picture completely changes when considering electron–phonon scattering where all energies (and momenta) are coupled. Recent advances in data-centric programming have shown that great speed-ups can be achieved if the parallelization is done based on the nature of the problem instead of physical intuitions.66 Another active research area is the development of so-called low-dimensional approximations, which drastically reduce the size of the Hamiltonian and overlap matrices in eqn (5).106 Such approximations can achieve cost reductions of multiple orders of magnitude with very little errors in the results. Presently, these approximations suffer from the requirement of devices made of the repetitions of identical cells. This is inherently incompatible with disordered structures as encountered in the oxide layer of RRAMs. Moreover, fully coupled electron–phonon simulations have yet to be demonstrated, when a low-dimensional basis is used. Nevertheless, if successful, such schemes could drastically change the size of the devices that can be simulated with NEGF at the ab initio level.
Limitations of the current methodology have been presented, together with opportunities for improvements. Many challenges remain to be tackled to enable fully ab initio dissipative quantum transport simulations of large-scale nanostructures. Most of them revolve around the computational cost of the DFT calculations that provide input parameters to the NEGF framework. Nevertheless, the demonstrated modeling platform combining first-principle calculations with force-field methods is already very powerful and can provide valuable insights into the physics and operation of nanodevices.
This journal is © The Royal Society of Chemistry 2020 |