Article pubs.acs.org/JPCC
Cite This: J. Phys. Chem. C 2018, 122, 20550−20554
Electron Transport from Quantum Kinetic Monte Carlo Simulations Fei Lin,*,†,‡ Jianqiu Huang,† and Celine Hin†,‡ †
Department of Mechanical Engineering and ‡Department of Materials Science and Engineering, Virginia Tech, Blacksburg, Virginia 24061, United States
J. Phys. Chem. C 2018.122:20550-20554. Downloaded from pubs.acs.org by KAOHSIUNG MEDICAL UNIV on 09/09/18. For personal use only.
S Supporting Information *
ABSTRACT: The exact quantum kinetic Monte Carlo method is proposed to calculate electron transport for onedimensional (1D) nanostructure described by the Hubbard model. The method is directly formulated in real time and can be applied to extract time-dependent dynamics of general interacting electronic models in 1D. When coupled with density functional theory calculation and maximally localized Wannier functions analysis, our method can be used to predict electron transport in nanostructures with interfaces. We first applied our method to case study of α-quartz dielectric breakdown.
■
INTRODUCTION Exact calculation of electron transport is one of the most difficult problems in nanoscience. Because of the important role it plays in electric device research, understanding electron transport from microscopic simulations is a subject of intensive ongoing researches.1,2 Interest in such calculations has grown when a large electric field is applied to semiconducting nanodevices, which eventually leads to dielectric breakdown of semiconducting materials.3−6 Two important issues arise in such researches: how can one accurately calculate the electrical conductivity in nanoscale devices, and how can one model the dielectric breakdown in these simulations. One popular approach adopts Boltzmann equation description of electron transport.3 However, exact solution of the Boltzmann equation remains a great challenge, and the method is limited to mesoscopic scale devices. Methods based on nonequilibrium Green’s function formulation have been proposed.2,7−10 But Green’s function methods are perturbative in nature and are most suitable for ballistic transport, as discussed by Luisier.2 The density matrix renormalization group (DMRG) method can also be applied to describe electron transport in a many-body system,11,12 but DMRG evolves a prepared state at zero temperature. A great progress in simulating dynamics of physical systems has been made when kinetic Monte Carlo (KMC) was applied to study electron transport,13 but the method uses phenomenological transition probabilities. A broad research interest also exists in quantum Monte Carlo (QMC) community in developing a feasible method for electric transport. However, the electric transport is essentially a real-time dynamics, which is in sharp contrast to the imaginary-time dynamics inherent in pathintegral formulation of the usual QMC methods. In this article, we describe a new modeling approach to reconcile all of the advantages of the different approaches © 2018 American Chemical Society
above. We propose a quantum kinetic Monte Carlo (QKMC) simulation based on the one-dimensional (1D) Hubbard model in real-time dynamics as opposed to the usual imaginary-time QMC simulation14−17 of the quantum Hamiltonian. These time-dependent simulations enable us to extract the electrical conductivity directly from the real-time evolution of the Hubbard model. We then relate, using density functional theory (DFT)18 coupled with maximally localized Wannier function (MLWF) approach,19,20 the electrical conductivity results to α-quartz (SiO2) in an effort to investigate the dielectric breakdown from purely electronic origin in the microscopic lattice scale. Thus, our QKMC method includes many-body effects, such as electron−electron inelastic scatterings, and performs simulations directly in real time at finite temperatures. We will first describe the QKMC model we developed and then we will apply it to a case study of α-quartz. More details are given below.
■
MODEL We study the following 1D Hubbard model M
∑
/ = −t
(ci†, σci + 1, σ + h. c .) + U ∑ ni ↑ni ↓
i = 1, σ =↑↓
−
∑ μi ni i
i
(1)
c†i,σ(ci,σ)
where is the electron creation (annihilation) operator on site i with spin σ, ni,σ = c†i,σci,σ is the electron density operator with spin σ, t is the hopping integral between nearest neighbors, and U is the on-site Coulomb interaction. The Received: June 4, 2018 Revised: July 24, 2018 Published: August 14, 2018 20550
DOI: 10.1021/acs.jpcc.8b05347 J. Phys. Chem. C 2018, 122, 20550−20554
Article
The Journal of Physical Chemistry C system has M sites with open boundary condition (OBC). A schematic for system setup is shown in Figure 1. On the left side, ML lattice sites have a chemical potential μ = μL and on the right side, MR lattice sites have a chemical potential μ = μR.
x=
Wi p(Wi − 1→Wi ) = Wi − 1 p(Wi →Wi − 1)
■
W′ = W
(8)
Nloop
∏ i=1
p(Wi − 1 → Wi ) p(Wi → Wi − 1)
(9)
However, the new configuration we get in the end of the loop N construction is already generated with probability ∏i=1loopp(Wi−1 → Wi). Thus, the configuration weight ratio we should use in eq 4 for the off-diagonal transition probability is
(2)
where β = 1/T is the inverse temperature, Nb is the number of bonds in the system, and the Hamiltonian eq 1 has been written as a summation over the bond operators /b , which is given by l ̃ n −Ù o Un μi ni if b is diagonal o o i↑ i↓ /b = m o o o−t(ci†σci + 1σ + h. c . ) if b is off‐diagonal (3) n here, Ũ = U/z and μ̃ i = μ/z with z = 2 the coordination number in 1D. A series expansion of the partition function to order Nc, which is large enough to contain all expansion terms at temperature T, is carried out. The resulting Hilbert space represented by a set of particle states and operator lists (called a configuration) with weight W is then sampled by the QKMC method. The probability for all possible transitions from the present configuration with expansion order n is tabulated. One typical probability table includes both diagonal and offdiagonal update probabilities together with their degeneracies. QKMC transition probabilities are given by x P= (4) 1+x
x=
1 N ∏i =loop p(Wi 1
→ Wi − 1)
(10)
Here, the selection probability for a specific state change at the entrance leg has been set to be equally probable. Note that the degeneracy of the off-diagonal transition is 4Nop, where Nop is the number of nonzero operators in the old configuration. Also note that infinite numbers of transition probability solutions p(Wi−1 → Wi) exist for local detailed balance eq 823 However, any set of solutions will give the same simulation result (see Supporting Information). A transition is picked from the probability table using the standard KMC procedure24 by uniformly drawing a random number R in the range (0, QK). The ith transition (either diagonal or off-diagonal) is chosen if Qi−1 ≤ R < Qi, where Q0 = 0 and i
Qi =
∑ djPj ;
i = 1, 2, ..., K
j=1
(11)
here, dj and Pj are the degeneracy and one of the diagonal or off-diagonal transition probabilities discussed above. K is the total number of all possible transitions in the system. Once the ith transition is chosen, a second random number R′ is drawn uniformly in the range (0, 1). The time interval between the (i − 1)th and ith transition is then given by
where W ′Pselect(W ′ → W ) WPselect(W → W ′)
(7)
where i = 1, ..., Nloop with Nloop being the length of the loop constructed. In the end, the ratio of the new configuration weight W′ = WNloop to the old one W = W0 is given by
METHOD In this section, we first describe the QKMC method that we developed for the Hubbard model. We further show that the method can be used to calculate the electric conductivity by applying a potential difference across the system. Following the stochastic series expansion (SSE) method,21,22 we write the partition function of the above Hamiltonian as
x=
Nc − n + 1 βNb⟨α|/b|α⟩
for decreasing expansion order n → n − 1. Here, /b is a diagonal operator picked for one of the Nb bonds in the lattice; J |α⟩ is the propagated state |α⟩ = ∏i = 1 /i|α0⟩, where J is the insertion or deletion point of /b and |α0⟩ is the initial particle state in the SSE expansion. The off-diagonal transition is constructed along a directed loop. All old particle states and operators on the loop will be changed to new particle states and operators. By enforcing detailed balance at the ith transition in the loop, we can calculate the new (Wi) and old (Wi−1) configurations weight ratio by
Figure 1. Schematic picture for 1D system setup. ML and MR are the number of lattice sites (solid dots) on the left and right sides of the system. The dashed line represents interface between these two regions. Hopping integral t is between nearest neighbors. Each site can be empty, singly occupied by a spin-up (up arrow) or -down (down arrow) electron, or doubly occupied by two electrons, in which case an on-site Coulomb interaction U is introduced to the system.
Nb
(6)
for increasing expansion order n → n + 1 and x=
A = Tre−β ∑b=1 /b
βNb⟨α|/b|α⟩ Nc − n
(5)
Δτi = −(τ0/Q K )ln R′
W(W′) are weights of the original(proposed) configuration, Pselect is the selection probability for transitions between W and W′ configurations. Specifically, for the diagonal transition
(12)
where τ0 is a characteristic time constant determined by the physical model.24 The cumulative time τ is then proportional 20551
DOI: 10.1021/acs.jpcc.8b05347 J. Phys. Chem. C 2018, 122, 20550−20554
Article
The Journal of Physical Chemistry C
transient behavior of NR. A linear fit to the transient data yields the current I/e = 0.271(1).
to real time in the quantum system evolution. We have applied our QKMC algorithm to a two-site system in equilibrium. The comparison among QKMC, conventional QMC, and exact results is given in Supporting Information. Our calculations demonstrate the real-time dependence implemented in our QKMC code. To access the electrical conductivity, we first equilibrate our system using the above QKMC procedure. A potential difference (μL < μR) is then created across the interface to mimic the effects of an applied electric field, which drives the electron transport from the left to the right side of the system. We can, therefore, measure the time evolution of total electron numbers on the right side NR(τ). The electric current I across the interface is then given by I=e
dNR dτ
Figure 3. Transient evolution (black circle) of the average electron number NR for t/U = 0.5, V/U = 0.2, and t/U = 0.1. The straight line is a fit to the data. Number inside bracket denotes error bar of the fitted coefficient.
(13)
where e denotes the electron charge. To guarantee that the increase of NR comes from the decrease of NL and not from the environment, we perform QKMC simulations in a canonical system, i.e., NR + NL is a constant. We emphasize that μL and μR in our model eq 1 do not refer to a grand canonical system in QKMC simulations. They merely create a potential difference across the interface.
Similar analysis is performed to obtain current I for different V’s. The results can be plotted in the I−V curve, as shown in Figure 4. From the Figure 4, we observe a linear I−V curve for
■
RESULTS We apply the above QKMC algorithm to the 1D Hubbard model to illustrate the electrical conductivity calculation. We use on-site Coulomb interaction U = 1 as energy unit. We set ML = MR = 10. Each simulation starts with a random configuration and is driven to thermal equilibrium through a typical 107 QKMC steps. The equilibrium state is then used to generate 104 equilibrium configurations, each being separated by 2000 QKMC steps. We then apply a series of electric potential differences V = μR − μL across the interface and measure particle number evolution NR(τ). In Figure 2, we show such an evolution, where the result is averaged over the above 104 equilibrium configurations.
Figure 4. Transient electric current I/e (black circle) across the interface as a function of V across the interface for t/U = 0.5 and 0.1. The straight line is a linear fit to the I−V curve.
small to intermediate V, suggesting an Ohm’s law behavior. A linear fit of the I−V curve yields the slope of the curve as the electric conductivity value σ = 1.44(1) for t/U = 0.5. We further vary the hopping integral t/U in the range (0.05, 0.5) and repeat previous analysis. The resulting electric conductivity σ as a function of t/U is plotted in Figure 5. At temperature T/U = 0.1, we identify the onset of nonzero electric current around t/U = 0.1, signaling the closing of energy gap in the insulating system.28−30 We now proceed to relate the model calculation shown in Figure 5 to the material we are studying. To apply our model Hamiltonian to a real material, we develop the following scheme to extract the effective hopping integral. The energy bands of α quartz are first calculated by the Vienna ab initio simulation software package,25,26 which shows an energy gap of ∼9 eV (see Supporting Information). An electric field is then applied to the DFT results to excite electrons out of the valence bands. Depending on the field strength, different numbers of bands will be involved in the wannierization process,19,20 which is performed with Wannier90 software package.27 This is reflected in different values of hopping integrals from different electric field values. We will follow the usual practice of evaluating hopping integrals using MLWFs as developed by Marzari et al.19,20 Assuming an on-site Coulomb
Figure 2. Average electron number NR as a function of cumulative time τ for different V’s across the interface. Here t/U = 0.5 and 0.1.
To extract the electric current, we focus on the transient and nonequilibrium behavior of NR, i.e., in a short cumulative time after the potential difference V is applied. Our transient measurement is similar to experimental procedure of obtaining I−V curve. However, in our method, we shall not look at the long evolution time limit, since such a limit will produce a new equilibrated state of the system, leading to an absence of electron transport across the interface. Figure 3 shows the 20552
DOI: 10.1021/acs.jpcc.8b05347 J. Phys. Chem. C 2018, 122, 20550−20554
Article
The Journal of Physical Chemistry C
modeling of electric transport is still a challenging problem in nanoscience. However, our current results point out that the purely electronic Hubbard model captures most of the physics of dielectric breakdown, since the critical electric field agrees well with experiments. Consequently, we can predict accurately the dielectric breakdown in nanodevices within a 1D Hubbard model approximation.
■
SUMMARY In this article, we proposed a new and powerful QKMC simulation method for predicting dynamical properties of Hubbard model. The method can calculate electric conductivity based on the real-time evolution of the quantum system. The onset of nonzero electric conductivity at low temperatures as hopping integral is increased also manifests the dielectric breakdown of insulators. When coupled with DFT and MLWF methods, our QKMC algorithm can be employed to predict electron-transport property of real materials. Electric conductivity behavior in strong electric field for α quartz, as an application example, seems to be consistent with the dielectric breakdown found in such materials.31
Figure 5. Electric conductivity σ as a function of hopping integral t for a 20-site system at temperature t/U = 0.1. The red squares are interpolations of electric conductivity for α quartz, whose hopping t is obtained by MLWF analysis in the Supporting Information.
interaction U ∼ 1 eV, we can interpolate the electric conductivity values for α quartz from Figure 5 and overlay the data in the same graph with the red squares. As one can see, the electric conductivity follows a nonlinear curve as hopping integral t increases. For t/U ∼ 0.1, which corresponds to V = 2.5 eV across a unit cell (∼6 Å) (see Supporting Information), the onset of nonzero conductivity suggests the dielectric breakdown of α quartz from an insulating to a conducting state. The critical electric field strength is estimated to be 4 × 109 V/m comparable to 3 × 109 V/m found in experiment.31
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b05347. Technical details in Monte Carlo simulations; comparison between an exact calculation and the Monte Carlo simulations; and the hopping integral calculation from density functional theory (PDF)
■
DISCUSSION QKMC simulation in this article is applied to the 1D Hubbard model with OBC, where the notorious sign problem is absent (see Supporting Information). Extension to higher dimension Fermion models is manageable if the simulation is performed at not too low temperatures, and the average sign is of order 0.1. However, we believe that in electric-transport problem, a 1D model is quite accurate in describing the electron motion. Assuming that an electric field is applied in the z-axis direction, the hopping integral t along z-axis direction becomes dominant in determining the electrical conductivity of the material. A two- or three-dimensional model can open up more parallel channels for electron transport, which can lead to a smaller t value for dielectric breakdown. Since a smaller t value means a smaller applied electric field, our previous critical field strength 4 × 109 V/m will get smaller and closer to the experimental value 3 × 109 V/m. Our simulations are performed on a 20-site system, which leads to the finite-size effect. Since a finite-size system usually overestimates the energy gap, we expect that the critical electric field strength is also overestimated. The actual value of critical field strength could, therefore, be further reduced and get closer to the experimental value. Our result of dielectric breakdown for insulators in this article is general and can be applied to various materials, including metal/oxide interface,32,33 defective interface,34,35 and amorphous36 systems. The increase of hopping integral as electric field increases is also general, but the exact hopping integral value will differ for various atomic configurations of materials. Of course, the mechanism for dielectric breakdown of semiconducting devices might not be just electronic.5,13 Effects of other factors, e.g., phonon37 or defects,38 on electrictransport property remain to be determined. Multiscale
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Fei Lin: 0000-0002-2372-2458 Jianqiu Huang: 0000-0002-8892-6291 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was funded by the Air Force with program name: Aerospace Materials for Extreme Environment and Grant No. FA9550-14-1-0157. We acknowledge Advanced Research Computing at Virginia Tech for providing the necessary computational resources and technical support that have contributed to this work (http://www.arc.vt.edu).
■
REFERENCES
(1) Datta, S. Quantum Transport: Atom to Transistor; Cambridge University Press: Cambridge, U.K., 2005. (2) Luisier, M. Atomistic Simulation of Transport Phenomena in Nanoelectronic Devices. Chem. Soc. Rev. 2014, 43, 4357−4367. (3) Jacoboni, C.; Reggiani, L. The Monte Carlo Method for the Solution of Charge Transport in Semiconductors with Applications to Covalent Materials. Rev. Mod. Phys. 1983, 55, 645−705. (4) Ray, S. An Introduction to High Voltage Engineering; Prentice-Hall of India Private Limited: New Delhi, India, 2004. (5) Chiu, F.-C. A Review on Conduction Mechanisms in Dielectric Films. Adv. Mater. Sci. Eng. 2014, 2014, No. 578168. (6) Luisier, M.; Schenk, A. Two-Dimensional Tunneling Effects on the Leakage Current of MOSFETs with Single Dielectric and High20553
DOI: 10.1021/acs.jpcc.8b05347 J. Phys. Chem. C 2018, 122, 20550−20554
Article
The Journal of Physical Chemistry C kappa Gate Stacks. IEEE Trans. Electron Devices 2008, 55, 1494− 1501. (7) Han, J. E.; Heary, R. J. Imaginary-Time Formulation of SteadyState Nonequilibrium: Application to Strongly Correlated Transport. Phys. Rev. Lett. 2007, 99, No. 236808. (8) Schmidt, T. L.; Werner, P.; Muhlbacher, L.; Komnik, A. Transient Dynamics of the Anderson Impurity Model out of Equilibrium. Phys. Rev. B 2008, 78, No. 235110. (9) Mühlbacher, L.; Rabani, E. Real-Time Path Integral Approach to Nonequilibrium Many-Body Quantum Systems. Phys. Rev. Lett. 2008, 100, No. 176403. (10) Schiró, M.; Fabrizio, M. Real-Time Diagrammatic Monte Carlo for Nonequilibrium Quantum Transport. Phys. Rev. B 2009, 79, No. 153302. (11) Schollwöck, U. The Density-Matrix Renormalization Group. Rev. Mod. Phys. 2005, 77, 259−315. (12) Kirino, S.; Ueda, K. Nonequilibrium Current in the One Dimensional Hubbard Model at Half-Filling. J. Phys. Soc. Jpn. 2010, 79, No. 093710. (13) Jegert, G. C. Modeling of Leakage Currents in High-κ Dielectrics. PhD Thesis, Walter Schottky Institute of Munich Technical University: Munich, Germany, 2011. (14) Ceperley, D. M. Path Integrals in the Theory of Condensed Helium. Rev. Mod. Phys. 1995, 67, 279−355. (15) Martin, R. M.; Reining, L.; Ceperley, D. M. Interacting Electrons, Theory and Computational Approaches; Cambridge University Press: Cambridge, U.K., 2016. (16) Hirsch, J. E. Monte Carlo Study of the Two-Dimensional Hubbard Model. Phys. Rev. Lett. 1983, 51, No. 1900. (17) Sandvik, A. W. A Generalization of Handscomba’s Quantum Monte Carlo Scheme-Application to the 1D Hubbard Model. J. Phys. A: Math. Gen. 1992, 25, 3667−3682. (18) Martin, R. M. Electronic Structure: Basic Theory and Practical Methods; Cambridge University Press: Cambridge, U.K., 2008. (19) Marzari, N.; Vanderbilt, D. Maximally Localized Generalized Wannier Functions for Composite Energy Bands. Phys. Rev. B 1997, 56, No. 12847. (20) Souza, I.; Marzari, N.; Vanderbilt, D. Maximally Localized Wannier Functions for Entangled Energy Bands. Phys. Rev. B 2001, 65, No. 035109. (21) Sandvik, A. W. Stochastic Series Expansion Method with Operator-Loop Update. Phys. Rev. B 1999, 59, No. R14157. (22) Syljuåsen, O. F.; Sandvik, A. W. Quantum Monte Carlo with Directed Loops. Phys. Rev. E 2002, 66, No. 046701. (23) Alet, F.; Wessel, S.; Troyer, M. Generalized Directed Loop Method for Quantum Monte Carlo Simulations. Phys. Rev. E 2005, 71, No. 036706. (24) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. A New Algorithm for Monte Carlo Simulation of Using Spin Systems. J. Comput. Phys. 1975, 17, 10−18. (25) Kresse, G.; Furthmüller, J. Efficiency of Ab-initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (26) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, No. 11169. (27) Mostofi, A. A.; Yates, J. R.; Pizzi, G.; Lee, Y. S.; Souza, I.; Vanderbilt, D.; Marzari, N. An Updated Version of Wannier90: A Tool for Obtaining Maximally-Localised Wannier Functions. Comput. Phys. Commun. 2014, 185, 2309−2310. (28) Kotliar, G.; Savrasov, S. Y.; Palsson, G.; Biroli, G. Cellular Dynamical Mean Field Approach to Strongly Correlated Systems. Phys. Rev. Lett. 2001, 87, No. 186401. (29) Schäfer, T.; Geles, F.; Rost, D.; Rohringer, G.; Arrigoni, E.; Held, K.; Blumer, N.; Aichhorn, M.; Toschi, A. Fate of the False Mott-Hubbard Transition in Two Dimensions. Phys. Rev. B 2015, 91, No. 125109.
(30) Neumayer, J.; Arrigoni, E.; Aichhorn, M.; von der Linden, W. Current Characteristics of a One-Dimensional Hubbard Chain: Role of Correlation and Dissipation. Phys. Rev. B 2015, 92, No. 125149. (31) Harari, E. Dielectric Breakdown in Electrically Stressed Thin Films of Thermal SiO2. J. Appl. Phys. 1978, 49, 2478−2489. (32) Huang, J.; Tea, E.; Li, G.; Hin, C. Hydrogen Release at MetalOxide Interfaces: A First Principle Study of Hydrogenated Al/SiO2 Interfaces. Appl. Surf. Sci. 2017, 406, 128−135. (33) Tea, E.; Huang, J.; Li, G.; Hin, C. Atomic Bonding and Electrical Potential at Metal/Oxide Interfaces, A First Principle Study. J. Chem. Phys. 2017, 146, No. 124706. (34) Tea, E.; Huang, J.; Hin, C. First Principles Study of Band Line up at Defective Metal-Oxide Interface: Oxygen Point Defects at Al/ SiO2 Interface. J. Phys. D: Appl. Phys. 2016, 49, No. 095304. (35) Li, G.; Tea, E.; Huang, J.; Hin, C. First Principle Thermodynamic Study of Oxygen Vacancy at Metal/Oxide Interface, arXiv:1611.03411. arXiv.org e-Print archive. https://arxiv.org/abs/ 1611.03411 (submitted Nov 10, 2016). (36) Huang, J.; Lin, F.; Hin, C. A DFT Based LCAO Approach to Determine Band Offsets and Dielectric Breakdown Properties Across Metal/Crystal Oxide and Metal/Amorphous Oxide Interfaces: an Al/ SiO2 Case Study. J. Mater. Chem. C 2017. (37) O’Dwyer, J. J. The High Temperature Dielectric Breakdown of Alkali Halides. Proc. Phys. Soc., Sect. B 1957, 70, 761. (38) Lombardo, S.; et al. Dielectric Breakdown Mechanisms in Gate Oxides. J. Appl. Phys. 2005, 98, No. 121301.
20554
DOI: 10.1021/acs.jpcc.8b05347 J. Phys. Chem. C 2018, 122, 20550−20554
