A fast model based on muffin-tin approximation to study charge transfer effects in time-dependent quantum transport simulations: doped Si-SiO2 quantum-dot systems
aa r X i v : . [ phy s i c s . pop - ph ] D ec APS/
A fast model based on muffin-tin approximation to study charge transfer effects intime-dependent quantum transport simulations: doped Si-SiO quantum-dot systems I-Lin Ho a) ChiMei Visual Technology Corporation, Tainan 741, Taiwan,R.O.C. (Dated: 12 December 2018)
In order to quickly study quantum devices in transient problems, this work demon-strates an analytical algorithm to solve the Hartree potential associated withcharge fluctuations in the time-dependent non-equilibrium green function (TDNEGF)method. We implement the calculations in the heterojunction system of gold met-als and silicon quantum dots for applications of photoelectric semiconductors in thefuture. Numerical results for the transient solutions are shown to be valid by compar-ing with the steady solutions calculated by the standard time-independent densityfunctional method.PACS numbers: 73.63.Kv, 05.60.Gg, 85.65.+hKeywords: quantum dot, time-dependent non-equilibrium green function, poissonequation a) Electronic mail: [email protected] . INTRODUCTION Photoelectric bioengineering - the use of photoelectric semiconductors as functional enti-ties in biological systems - is heralded as an alternative option for signaling communicationsbetween organisms and physical devices in future biomedicines. In particular, researchon quantum dots has already revealed a variety of biologically-oriented applications, e.g.drug discovery , disease detection , protein tracking , and intracellular reporting . Aqualitative understanding of these complex processes has been accessed by perturbativeelectron-photon interactions associated with strong electron correlations , but the quanti-tative agreement between the first-principles theory and experiments is still unsatisfactoryfrom the perspective of the ground-state density functional theory (DFT) .In recent years, the majority of studies for quantum-dot electronics have focused on thetime-dependent density functional theory (TDDFT) that provides a more rigorous theo-retical foundation . The formalism may also be easily extended to cover the interaction ofelectrons with light or under environments in open quantum systems by the time-dependentnon-equilibrium green function (TDNEGF) technique , e.g. for the photon-assisted trans-port and fluorescence of contacted atomic devices.Several scenarios of open quantum systems implemented with TDDFT have so far beensuggested, including the ring-topology of the electronic circuit and the approximately-isolated atomic device . This present work adopts more general set-ups to study systemscomposed of functional atomic devices and environmental clusters as shown in Figure 1.Here, the device in the central region focuses on the Si-SiO core-shell quantum dot due toits wide application spectrum. Moreover, the energetically-favored phosphorus impurities areconsidered for low-voltage operations. Nuria reported relevant analyses of doped Si-SiO quantum dots in detail. The effect of the neighboring Au(111) electrodes upon devices isalso exactly accounted for through properly defined self-energies. For numerical treatments,the initial Kohn-Sham (KS) single-particle Hamiltonians and the overlap matrices for theSi-SiO quantum dots and Au electrodes are obtained by DFT calculations in SIESTAprograms . On the basis of the holographic electron density theorem and Runge-Grosstheorem, the time-dependent electron dynamics are then determined by solving equationsof motion for the devices using the TDNEGF technique in own-implemented fortranprograms . 2omputations by TDNEGF of realistic devices having a large amount of atoms are nu-merically demanding, because all electron motions have to be fully resolved, leading toconsiderable degrees of freedom in the sub-fs time resolution. To arrive at a computation-ally efficient but still predictive method, this study follows the work of Chen for openquantum systems driven by time-dependent bias potentials. Furthermore, to consider theeffects of the charge variations through devices, this work demonstrates a capacitive networkmodel as an analytical Poisson solution in the muffin-tin approximation. Rather than thenumerically-demanding iterative Poisson-equation solution using discretized spatial grids ,the model can quickly solve the Hartree potential associated with the charge fluctuations inthe time-dependent non-equilibrium green function (TDNEGF) formulae. Numerical resultsare shown to be valid by comparing with the steady solutions calculated by the standardtime-independent density functional method using SIESTA.This paper is organized as follows. Section II describes the employed computationalalgorithms. Section III demonstrates the density of state of the gold electrode for validatingthe wide-band limit approximation, and depicts the characteristics of energy spectrumsfor quantum-dot devices. We then compute the electronic transmission functions of thequantum device coupled to the semi-infinite electrodes by the SIESTA program so as tostudy the steady transport dynamics. The properties of the transient electronic transportfor the open quantum system are simulated by using the a fortran program, and are shownto be good versus the steady solutions in a long time limit. Finally, Section IV presentsconcluding remarks. The mathematical derivations and relevant physical approximationsfor the time-dependent non-equilibrium green function (TDNEGF) formulae are stated inthe appendix. II. TIME-DEPENDENT NON-EQUILIBRIUM GREEN FUNCTION FORQUASI-ONE-DIMENSIONAL OPEN QUANTUM SYSTEMS
Figure 1 shows a regular open quantum system, including semi-infinite electrodes andatomistic devices. The system is partitioned by several electronically-functional areas,named as L-electrode (L), device (D), and R-electrode (R). The equation of motion (EOM)for electrons can be described by quantum dynamics: i ˙ σ ( t ) = [ h ( t ) , σ ( t )] (1)3 IG. 1. Schematic representation of general simulation setups for open quantum systems includingAu(111) electrodes and atomic devices. where h ( t ) is the Kohn-Sham Hamiltonian matrix, and the square bracket on the right-handside (RHS) denotes a commutator. The matrix element of the single-electron density σ isdefined by σ ij ( t ) = D a † j ( t ) a i ( t ) E , where a † j ( t ) and a i ( t ) are the creation and annihilationoperators for atomic orbitals j and i at time t , respectively. On the basis of the atomicorbital sets for electrons, the matrix representation of σ and h can be written as h = h L h LD h DL h D h DR h RD h R , σ = σ L σ LD σ LR σ DL σ D σ DR σ RL σ RD σ R (2)Here, m L , m D , and m R ( m ∈ { h , σ } ) represent the matrix blocks corresponding to left-electrode L , device D , and right-electrode R partitions, respectively. Moreover, h LR and h RL are ignored due to the distant separation between L and R electrodes in common appli-cations. We note that the holographic electron density theorem and Runge-Gross theoremare adopted for time-dependent electron dynamics , stating that the initial ground-statedensity of the subsystem σ D ( t ) can determine all physical properties of systems at anytime t . Hence h and σ can be approximately expressed as functions of σ D ( t ) for a formallyclosed-form equation of motion as described below.Placing Eq. (2) into Eq. (1), we can write the equation of motion for σ D as i ˙ σ D,mn = X ℓ ∈ D ( h D,mℓ σ D,ℓn − σ D,mℓ h D,ℓn ) − i X α = L,R,N Q α,mn (3) Q α,mn ≡ i X k α ∈ α ( h Dα,mk α σ αD,k α n − σ Dα,mk α h αD,k α n ) (4)Here, m and n denote the atomic orbital in partition D, k α denotes the state of α ( α =L, R)electrode, and Q α is the dissipation term due to the contacts of the device with electrodes4 and R. The transient current through an electrode’s interfaces can be calculated by: I α ∈{ L,R } ( t ) = − Z α d r ∂ t ρ ( r , t ) = − X k α ∈ α ∂ t σ k α k α ( t )= i X k α ∈ α X ℓ ∈ D ( h Dα,k α ℓ σ αD,ℓk α − σ Dα,k α ℓ h αD,ℓk α )= − tr [ Q α ( t )] (5) A. Expressions of the dissipation function Q α using the Green functionformalism To calculate the dissipation term Q α in EOM and transient current equation, we usethe time-dependent non-equilibrium Green function (TDNEGF) formalism. We note thatthe overlap matrix s is treated as an identity matrix I when deriving the Green functionformalism in the appendices. This replacement of the overlap-matrix by the identity matrixhas been verified to be valid if the hamiltonian matrix is modified according to mathematicaltechniques in the textbook (Ch. 8.1.2), e.g. h D − E s D = h D − E ( s D − I ) − E I = h ′ D − E I .Appendix A demonstrates relevant approximations and derivations, and gives the formulaeof Q α : Q α,mn ( t ) = − X ℓ ∈ D Z ∞−∞ dτ (cid:2) G For the open quantum system, the device Hamiltonian can take the perturbative form of: h D = h D ( q ) + δ h D ( δq ) (17)Here, the change of electron density can be calculated by the density matrix σ D in Eq. (3) δn ( ~r ) = X µν Re [ ρ µν χ µ ( ~r ) χ ∗ ν ( ~r )] − n ( r ) (18)as a spatial distribution function, or, alternatively, by δq i = X µ ∈{ i } X ν Re [ ρ µν s D,νµ ] − q ,i (19)for atom i . Here, n ( r ) and q ,i are the reference atomic charges chosen for neutrality, s D is the device overlap matrix, and χ i ( ~r ) is a set of local basis functions used in thetight-binding formulation. According to a Taylor expansion of the total energy around the7eference density, this change of charge density can result in corrections to the Hartreeand the exchange-correlation potentials for the device Hamiltonian as in Eq. (17), andit is continuously renewed with the density matrix in Eq. (3). In this work we simplifythe correction of the device Hamiltonian δ h D by retaining only the Hartree potential δV H (assuming the exchange-correlation term is insignificant in the mean-field scope), whichobeys the three-dimensional Poisson equation ∇ δV H ( r ) = − δn ( ~r ) (20)with the boundary conditions imposed by the lead potentials. The conventional Poissonsolution is based on spatially-discretized grids for numerically iterative processes, and canbe time-consuming for large systems. Thus, it is convenient to develop an approximatelyanalytical model.Extending the idea of the muffin-tin (MT) approximation, each atom i can define aspherical region (MT-sphere) that bounds the total excess charges δq i from Eq. (19). Thepaired parts of charges inside different neighboring MT-spheres are considered as capacitanceeffects . All these spheres of atoms in the system then further construct a capacitancenetwork architecture that supplies an analytical solution for the Poisson equation . Inprinciple, the capacitances are treated as a combination of the electrostatic capacitance c e and quantum capacitance c Q . Herein, the quantum capacitance is assumed to be lessdominant than the electrostatic capacitance for δq i and is ignored in our work.Replacing the spatial solution ( ∇ r ) of Poisson equation by the atomic-site solution ( ∇ i )of the capacitive model , we can rewrite Eq. (20) by a matrix-form equation ˆC ˜V = ˜Q c ij = 4 πǫ ¯ a ij | r ij | (cid:18) a ij | r ij | − a ij + ... (cid:19) (21)ˆ C ij = X k ∈{ NN } i,con δ i,j c ik + X k ∈{ NN } i δ i,j c ik − δ j,k c ij (22) ~Q i = e · δq i + e · δq d,i + X k ∈{ NN } i,con. δ j,k c ij V con.j (23)Here, the matrix elements of ˆC are calculated in a two-center approximation as proposed inthe tight-binding formulation , obeying the formal condition eδq i + eδq d,i = P j c ij ( δV i − δV j ).The notation { N N } i is the group of the first nearest-neighbor (NN) atoms in the deviceregion for atom i , and { N N } i,con. is the group of the first nearest-neighbor atoms in the leadregion for device atom i . Extending the solution by including more capacitively-coupling8erms { nN N } ( n ∈ , , ... ) is reasonable, because the additional capacitance terms c ij ( n ≥ 2) diminish with the increasing separation | r ij | as indicated in Eq. (21). Herein, c ij defines the capacitance between two ideal metal spheres, | r ij | is the spatial distance betweenatoms i and j , and ¯ a ij is the effective muffin-tin radius for atoms i and j and is definedby ¯ a ij = ( r MT,i + r MT,j ) / ˜V ≡ ( δV , δV , ..., δV N ) is the potential vectorwith the components being the electrostatic potential for device atoms i ∈ { , ..., N } . V con.j is the potential of lead atom j imposed by boundary conditions. δq i is the charge densityobtained by Eq. (19), and δq d,i represents the defect charge for atom i . By linear algebra thepotential vector ˜V can be easily solved using ˜V = ˆC − ˜Q . For instance, in a 1-dimensionalhomogeneous system having 4 atoms L-A-A-R, the capacitance between nearby atoms isdenoted as c , and the biases are denoted as v L and v R for lead atoms L and R, respectively.There are no excess charges ( δq = 0) inside the MT-sphere of device atoms A. In this way,the 2x2 capacitance matrix has components ˆ C = ˆ C = 2 c and ˆ C = ˆ C = − c , andthe charge vector is ˜Q t = [ cv L cv R ]. One then can obtain the electrostatic potentials forthe two atoms A as ˜V t = [ 2 v L + v R v L + 2 v R ] / 3, which agree with the free-space Poissonsolution.Figures (2-3) illustrate 2-dimensional examples with a comparison between the numeri-cally iterative solution and the analytical capacitance model. In order to solve the spatialPoisson equation in Eq. (20), the density function δn i ( r ) for two-dimensional systems isassumed to be: δn i ( r ) = δq i πη e − | r − Ri | η (24)This is according to the symmetry assumption , where R i is the position for atom i , and η is associated with the effective radius of the MT-sphere by η ∝ r MT,i (use η = r MT,i forexamples in Figs. (2-3)). The obtained potential δV H ( r ) is projected on the atomic sitesthrough δV i = R d r δV H ( r ) e − | r − Ri | η R d r e − | r − Ri | η (25)for a comparison with the analytical solution ˜V in this work. We note that the analyticalmodel associated with orientatingly capacitive couplings implies a spatial density functionbeyond the symmetry assumption. As indicated in Figs. (2-3), the analytical solutionshows comparable results with that from the numerically-iterative method. Otherwise, itis emphasized that the analytical model turns inefficient at large biases or strong density9 a) numerical Poisson solution (b) analytical solution V r =-1V l =1 q=0 V=0.33 V=0.39 V=-2.32 V=0.81 q=0 V=0.59 q=0 V=0.46 q=0.05 V=-2.07 V=0.40 FIG. 2. Profile of Hartree potential δV i for the first structure with area 9 × A , solved by (a)the numerical Poisson solution and (b) the analytical solution. Four atoms with specified charges δq (see Figure (a)) are placed between two leads and have r MT = 1˚ A . In (a) the additionalspatial function δV ( r ) is illustrated by the contour curves. In (b) each boundary condition of theelectrostatic potential is represented by four lead atoms. variations, because the MT sphere cannot accurately account for the distorted and displacedelectron density distribution from the nucleus.The relevant parameters of the MT radius used herein are r MT ( Si ) = 1 . A , r MT ( O ) = 0 . A , r MT ( Au ) = 1 . A , and r MT ( P ) = 1 . A . All computations areoperated on a workstation having 2xCPU(E5-2690 v2) and 128G of DRAM. Fortran sourcecodes can be downloaded online . III. NUMERICAL RESULTSA. Atomic Electrodes: Au(111) Nanotubes This research uses Au(111) nanotubes as atomic electrodes. The length ℓ of the Au-Aubond is determined with geometry relaxations of the Au bulk in the SIESTA program ,and the obtained value is ℓ =2.8785 ˚ A (lattice constant a= √ ℓ =4.0708 ˚ A , which is simi-lar to the experimental value of 4.0782 ˚ A ). The effects of core electrons are evaluatedwith norm-conserving pseudopotentials in the local density approximation (Ceperley-Alderexchange-correlation potential ), which are generated by the ATOM program . Thevalence electrons of Au are calculated in the s-d hybridized configuration . All the calcu-10 a) numerical Poisson solution (b) analytical solution V r =-1V l =1 q=0.04 V=-1.03 q=-0.07 V=1.91 q=-0.03 V=0.21 q=0.05 V=-0.28 V=2.19 V=-1.58 V=-0.43 V=0.34 FIG. 3. Profile of Hartree potential δV i for the second structure with area 9 × A , solved by (a)the numerical Poisson solution and (b) the analytical solution. Relevant setups are the same withthat in Fig. 2, except for atom charge δq . (a) (b) (c) R PL PL PL PLq- q- FIG. 4. Ball-stick representation of the Au(111) nanotube in (a) the longitudinal perspective andin (b-c) two lateral perspectives. The radius of the cross section in (a) is set as R=2a. Two ofthe principle layers (PL) arranged along the longitudinal (transport) direction in (b-c) represent asegment of the nanotube. q − indicates the quasi one dimensional (red arrow) charge transport. lations for nanotubes are performed on 8 × × E F is the Fermi level corre-11 10 -8 -6 -4 -2 0 2 E-E F (eV) no r m a li ze d DO S bulktube(R=0.5a)tube(R=2.0a)tube(R=4.0a) FIG. 5. Normalized density of states (DOS) for Au bulk and infinite Au(111) nanotubes, in whichthe radiuses of the nanotubes are set as R=0.5a, R=2.0a, and R=4.0a. sponding to the mentioned system. In Fig. 5, DOS of the Au bulk shows metallic propertiesas the literature reports. When increasing the cross-section radius R, the DOS functionsof Au(111) nanotubes at energies near E F change from discrete to uniform distributions,depicting the transfer of systems from 1D-nanotube to 3D-bulk structures. In this work, weuse Au(111) nanotubes with R=2a for semi-infinite electrodes in transport problems. Thisadoption (setting R=2a) meets the requirement of slowly-varying DOS for the wide-bandlimit (WBL) condition , and demands computation resources that are affordable. B. Doped Si-SiO quantum dots This research investigates the silicon quantum dots with diameters around 1.0 nm thatare embedded in a β -cristobalite SiO matrix. The dopant phosphorus (P) atoms are placedinside quantum dots due to their energetically-favored formation of structures (see Fig.6). Lattice constants are determined with geometry relaxations in the SIESTA program(set orbital bases s and p for species Si, O, and P). The obtained values are 5.5001 ˚ A (5.4306 by experiment ) for the Si diamond structure and 7.46831 ˚ A (7.160-7.403 ˚ A intextbooks ) for β -cristobalite silica. We investigate the energy band diagram of Si-SiO -slabs heterojunctions by using Anderson’s rule through Fig. 7, in which the vacuum levels(green dotted lines) of Si and SiO slabs are aligned at the same energy. Here, the vacuum12 X FIG. 6. Schematics of Si Quantum Dot (red atoms) embedded in SiO matrix (light cyan-yellowatoms). The phosphorus atom P (blue atom) is doped inside the quantum dot for the 1P-dopingcondition. Two red atoms with X marks denote the doping locations inside the quantum dot andat the interface, respectively, for the 2P-doping condition. level is defined as the effective potential φ (adding local pseudopotential, Hartree potential,and exchange-correlation potential) at zero-density points near the surface of slabs having35 atomic layers. All calculations are performed at Γ-point of the reciprocal space. Asindicated in Fig. 7, the vacuum levels are 1.064 eV and 1.626 eV for Si-slab and SiO -slab, respectively, corresponding to working functions W Si =4.46 eV and W SiO =4.52 eV.The experimental value is 4 . ≤ W Si ≤ . 91 eV. The computed energy gaps are 1.17eV for bulk silicon and 7.7 eV for β -cristobalite silica, which can be compared with theexperimental values of 1.1 eV and 9.0 eV, respectively. The valence band offset (VBO)and conduction band offset (CBO) for Si-SiO heterojunctions are estimated to be 3.18eV and 3.31 eV, respectively. The obtained VBO values are smaller than experimentalmeasurements with VBO=4.6 eV and CBO=3.1 eV. Several theoretical works usinghopping mechanisms give VBO ≈ ≈ quantum-dot device in Fig. 6 is con-structed from a 3 × × β -cristobalite silica by removing O atoms in a cut-offbox. Figure 8 reports the eigenvalue spectra for the undoped, 1P-doping, and 2P-doping13 e n e r gy ( e V ) DOS of Si DOS of SiO E g,Si =1.17 VBO=3.18CBO=3.31 E g,SiO =7.7 FIG. 7. Band diagrams of Si-SiO -slabs heterojunction by Anderson’s rule. The density of state atequilibrium is arranged according to a hypothetical flat vacuum level. The computed energy gapsare E g,Si =1.17 eV and E g,SiO =7.7 eV. The valence band offset VBO is 3.18 eV and the conductionband offset CBO is 3.31 eV. structures after relaxation processes, using the corresponding initial geometries in Fig. 6.The spectrum energies are aligned using the SiO states (deep valence states), and the en-ergy axes show the common origin according to the fermi level of the undoped structure.Black and gray circles mark the highest occupied molecular orbital (HOMO) and lowestunoccupied molecular orbital (LUMO) states, respectively. The green dotted line representsthe fermi level of the corresponding structure. In Fig. 8(a), the undoped structure exhibitsa distinguished energy spectrum from that of the slab-heterojunction in Fig. 7, revealing theessential mechanism for strain-induced electron levels . For the 1P-doping system, the oddnumber of electrons leads to the spin-dependent energy spectrum in Fig. 8(b), which depictsa clear donor behavior and agrees well with previous works . This study adopts the 2P-dopping structure in Fig. 8(c) due to the following considerations: (i) high conductivity ata low bias V owing to the rising fermi level and the decreasing energy gap, compared to theundoped structure; and (ii) having spin independence for better computational efficiencyand the negligible spin-flip mechanism, compared to the 1P-doping case.14 IG. 8. Spin-up and spin-down spectra of (a) undoped, (b) 1P-doping,and (c) 2P-doping systems.Energies are aligned using the embedding SiO states, and are shifted with the reference of thefermi level of the undoped structure. Black and gray circles mark HOMO and LUMO states,respectively. The green dotted line represents the fermi level of the corresponding system. C. Steady (time-independent) electron transport in open quantum-dotsystems Fig. 9 constructs the open transport system of quantum dots. The device region, i.e.Si-based quantum dot (red spheres) and SiO matrix (small light cyan-yellow spheres), isenclosed by two semi-infinitely long Au wires. Two dopant atoms (phosphorus; blue spheres)are placed inside the quantum dot and at the Si-SiO interface, respectively, according totheir energetically-favored formation energy . It is assumed that the positions of atoms ofAu electrodes are under constraint by the experimental set-ups, while the atoms of the dopedSi-SiO quantum dot are in equilibrium according to geometry relaxations by the SIESTAprogram. The distance between the nearest cross sections of silica and gold boundariesbefore geometry relaxations is initially set to be 1.8 ˚ A in this work.By applying non-equilibrium green functions for steady transport problems , the trans-mission function of the quantum-dot system is obtained as shown in Fig. 10. In Fig. 10(a),the transmission function T (blue curve) is calculated by SIESTA::Transiesta programs, andis compared with the projected density of state (PDOS; gray curve) of the Si-SiO quantum15 IG. 9. Schematics of a Si-based (red spheres) quantum dot embedded in β -cristobalite SiO matrix (small light cyan-yellow spheres), where two dopant atoms (blue phosphorus spheres) areplaced inside the quantum dot and at the Si-SiO interface, respectively. The device is enclosedbetween two semi-infinitely-long Au(111) wires (larger yellow spheres) with an applied voltage. dot. The red-curve transmission function is calculated by the fortran program which extractsthe relevant hamiltonian and overlap matrices from SIESTA for modeling the tight-bindingformulation, and the numerical framework is employed in time-dependent non-equilibriumgreen functions for transient problems below. Numerical results demonstrate that (i) theconductance channels in the transmission function T are associated with the density of statesof the Si-SiO quantum dots as expected, and (ii) the tight-binding formulation works wellsince the comparable transmission functions by SIESTA and fortran programs. Figure 10(b)depicts transmission functions for several systems calculated by SIESTA programs with dif-ferent bias setups. It is found that the transmission profiles of the devices non-linearly varywith biases, and reveal considerable effects of charge fluctuations inside the device. It is em-phasized that SIESTA is the standard time-independent density functional program withoutthe mentioned approximations for this work. We also note that the quantum-dot systempresents non-zero conductance at near zero bias, which is similar to the analysis in Nuria’swork . This conductance associated with the finite transmission function at the fermi level,however, decreases with an increasing bias as shown in Fig. 10(b).16 (a) E-E f (eV) T r a n s m i ss i on T for V=0.0 eV by fortranT for V=0.0 eV by siestaPDOS for V=0.0 eV by siesta -2 0 2 (b) E-E f (eV) T r a n s m i ss i on T for V=0.0 eV by siestaT for V=0.1 eV by siestaT for V=0.5 eV by siesta FIG. 10. (a) Transmission functions of the quantum-dot system (V=0.0 eV) using SIESTA andthe fortran programs, to be compared with the projected density of state (PDOS) of Si-SiO QDs.(b) Transmission functions of the quantum-dot system calculated by SIESTA programs at differentvoltage biases. D. Transient (time-dependent) electron transport in open quantum-dotsystems This section studies the time-dependent electron transport for the quantum-dot systemin Fig. 9. Additional parameters and numerical methods are as follows: time step δt = 5 as ,voltage function V f ( t ) = V dc (cid:2) − exp − t/τ (cid:3) + V ac cos ( ωt ) with τ = 2 f s , globally-adaptiveintegrator treating singularities in the energy domain, and the fourth-order Runge Kuttamethods (RK4) for solving the time-differential equation. Here, we adopt the linear extrap-olation of the density matrix σ D (Eq. 3) during the iterative process by RK4.Figure 11(a) shows transient currents by our TDNEGF codes with and without cor-rections for charge transfer effects (CTE). The voltage functions are symmetrically set by V L = V f,V dc =0 . V,V ac =0 V and V R = V f,V dc = − . V,V ac =0 V for the left and right electrodes, re-spectively. Moreover, the time-independent solutions for steady currents are derived viaLandauer Buttiker formula , an integral of the transmission functions in Fig. 10(b). Wenote that the integrals using the transmission function for V = 0 . V = 0 . (a) time/ps -3-2-1012 c u rr e n t I ( A ) -4 I L,V=0.5 (t) wo CTEI R,V=0.5 (t) wo CTEI L,V=0.5 (t) w CTEI R,V=0.5 (t) w CTE -6 t-indep. I L wo CTEt-indep. I R wo CTEt-indep. I L w CTEt-indep. I R w CTE (b) time/ps t r ace () / e ii (t) wo CTE ii (t)- I L dt- I R dt wo CTE ii (t) w CTE ii (t)- I L dt- I R dt w CTE FIG. 11. (a) Transient current I of quantum-dot devices at a symmetry DC bias V dc =0.5Vwith/without including charge transfer effect (CTE). The inset diagram shows that the currentsasymptotically approach the value of the steady solution by SIESTA (green curves). (b) Tran-sient charge numbers of the quantum-dot device with/without including CTE. After applying bias,a part of the electrons (Σ ρ ii , solid curves) participate in the inter-orbital transferring process( R I L dt + R I R dt ). The dot curves show the conservation of total charges. observe that the transient currents asymptotically approach the values of the correspondingsteady solutions by SIESTA (see the inset diagram), no matter the charge transfer effectsare considered or not. This concludes the validation of the TDNEGF program as well asthe proposed analytical model for treating CTE. Moreover, the inclusion of charge transfereffects presents considerable corrections for the convergence of the transient current, sug-gesting non-trivial variations/excitations of charge in the device. The curves also depictthat the calculation including CTE requires a much longer time to bring the system intothe steady sate, inferring a slow redistribution process of the charge density. Figure 11(b)shows the transient electron number of the device and the integrals of boundary currents,obeying the continuity equation for the device region.In Fig. 12, we study the transient currents for the quantum-dot devices under AC biasvoltages. To observe the diffusion of CTE through the device, the voltage functions areasymmetrically set by V L = V f,V dc =0 . V,V ac =0 V and V R = V f,V dc = − . V,V ac =0 . V for the leftand right electrodes, respectively. Here the AC frequency is ω = 0 . × . In Fig.12(a), for the quantum-dot device with charge transfer effects, the interfacial current I L (a) time/ps -50510 c u rr e n t I ( A ) -4 I L,V=0.1 wo CTEI R,V=-0.1+0.4cos( t) wo CTEI L,V=0.1 w CTEI R,V=-0.1+0.4cos( t) w CTE (b) time/ps -505 a v e r a g e c u rr e n t I a vg . ( A ) -5 I L,V=0.1 wo CTEI R,V=-0.1+0.4cos( t) wo CTEI L,V=0.1 w CTEI R,V=-0.1+0.4cos( t) w CTE0.03 0.04 0.05-4-202 10 -6 (c) time/ps t r ace () / e ii (t) wo CTE ii (t)- I L dt- I R dt wo CTE ii (t) w CTE ii (t)- I L dt- I R dt w CTE FIG. 12. (a) Transient current I of quantum-dot devices at AC bias V L =0.1V and V R =-0.1+0.4 cos ( ωt ) V, with/without including a charge transfer effect (CTE). Here the AC frequencyis ω = 0 . × . By including CTE, I L exhibits continuous oscillations of interfacial currents.(b) Effective currents I avg. by averaging I (t) in (a) over one period T = 2 π/ω . The inset diagramshows that the effective currents, with/without including CTE, asymptotically converge into sim-ilar values. (c) Transient charge numbers of the quantum-dot device with and without includingCTE. The dot curves show the conservation of total charges. of the left electrode exhibits continuous oscillations; while I L quickly declines to a steadyvalue for the case without charge transfer effects. The interfacial current I R of the rightelectrode, however, is always oscillatory due to the driving of the AC potentials at the local(right) electrode. This observation identifies the AC-induced oscillation of charge densitiesinside the quantum-dot device via CTE algorithms. To correlate with the alternative DCmeasurements, we calculate the average current I avg. ( t ) by averaging I ( t ) over one period T = 2 π/ω . Numerical results of I avg. ( t ) are shown in Fig. 12(b). In its inset diagram,the effective currents, no matter with and without charge transfer effects, asymptoticallyconverge into similar values. We attribute the similarity of the asymptotical currents to thelow DC bias, which contributes insignificant charge fluctuations on average, as suggested bythe alike curves for V=0.0 and 0.1 eV in Fig. 10(b). Figure 12(c) validates the algorithmsby the continuity equation in AC cases. 19 V. CONCLUSIONS This research proposes an approximate analytical model to efficiently calculate the tran-sient properties of quantum-dot systems under time-dependent external potentials. Numer-ical results in the low DC bias and long-time limits present good agreements with the cor-responding steady solutions, no matter the charge transfer effects are included or excluded.For the cases using asymmetric AC biases, numerical calculations for transient currents alsoshow distinct characteristics between the systems with and without charge transfer effects,revealing the essential oscillations/excitations of charge densities inside the device. V. ACKNOWLEDGEMENT This work was supported by ChiMei Visual Technology Corporation. Appendix A: Equation of motion for Green functions The Hamiltonian operator ˆh for open transport systems without spin notations can begiven by ˆh = X k α h α,k α n k α + X m,n h D,mn a † m a n + X m,k α h Dα,mk α a † m a k α + h αD,k α m a † k α a m (A1)The first term describes the α th electrode with state k α , the second term is for the devicein geometry region D, and the third term is for the coupling between the device and theelectrode α .In this appendix, the algorithm of the time-dependent non-equilibrium green functionis addressed for systems under the following conditions: during t < t , the system is inthermal equilibrium at an inverse temperature β and chemical potential µ . At times t ≥ t ,the system departs from the equilibrium conditions after applying external voltages. Here,the response of the spin to external fields is ignored, resulting in diagonal green-functionand self-energy matrices with respect to the spin parameter. The initial condition is definedby the ground state of the system.For time-dependent electron transport problems, one begins with the one-particle greenfunction on the Keldysh contour γ K (see Fig. 13). This green function is defined as the20 K ∞ − + − iβ t R t I FIG. 13. The Keldysh contour γ K is an oriented contour having endpoints at 0 − and − iβ . β isthe inverse temperature. The contour is composed of a forward branch going from t = 0 − to t= ∞ ,a backward branch coming back from t = ∞ to t = 0+, and a vertical (thermic) track on theimaginary times axis between 0 + and − iβ . z and z ′ define variables along γ K . ensemble average of the contour-ordered product of electron creation and annihilation op-erators in the Heisenberg picture, G rs ( z, z ′ ) = − i (cid:10) T C (cid:2) a r ( z ) a † s ( z ′ ) (cid:3)(cid:11) (A2)Here, r, s present states of L(left-electrode), R(right-electrode), and D(device). z and z ′ define complex variables along the contour γ K . T C is the time-ordering operator. Thecreation operator a † and annihilation operator a obey the equation of motion i ddz a ( z ) = [ a ( z ) , h ( z )] = h ( z ) a ( z ) (A3) i ddz a † ( z ) = (cid:2) a † ( z ) , h ( z ) (cid:3) = − h ( z ) a † ( z ) (A4)Here, the anticommutator relation for fermions have been used, i.e. { a r , a s } = (cid:8) a † r , a † s (cid:9) = 0and (cid:8) a r , a † s (cid:9) = h r | s i with orthonormal state bases r and s . Combining Eqs. (A3-A4) and(A2), the equation of motion for the green function can be given as: (cid:20) i ddz − h ( z ) (cid:21) G ( z, z ′ ) = δ ( z, z ′ ) 1 (A5) G ( z, z ′ ) (cid:20) − i ddz ′ − h ( z ′ ) (cid:21) = δ ( z, z ′ ) 1 (A6)Here, the green function follows Kubo-Martin-Schwinger (KMS) boundary conditions on theimaginary axis in γ K . h ( z ) is the single-particle Hamiltonian.21pplying the definition in Eq. (2), the matrix structure of G has block matrix form: G = G L G LD G LR G DL G D G DR G RL G RD G R (A7)Equation (A5) in matrix form is hence given by i ddz G ( z, z ′ ) − h ( z ) G ( z, z ′ ) = δ ( z, z ′ ) (A8)Here, the equations for components G αD and G D are (cid:20) i ddz − h D ( z ) (cid:21) G D ( z, z ′ ) = δ ( z, z ′ ) + X α ∈ L,R h Dα ( z ) G αD ( z, z ′ ) (A9) (cid:20) i ddz − h α ( z ) (cid:21) G αD ( z, z ′ ) = h αD ( z ) G D ( z, z ′ ) (A10)By multiplying Eq. (A10) with the green function G α , i.e. (cid:2) − i ddz ′ − h α ( z ′ ) (cid:3) G α ( z, z ′ ) = δ ( z, z ′ ) in Eq. (A5), we can obtain G αD ( z, z ′ ) as Z γ K d ¯ z G α ( z, ¯ z ) (cid:20) i dd ¯ z − h α (¯ z ) (cid:21) G αD (¯ z, z ′ ) = Z γ K d ¯ z (cid:20)(cid:18) − i dd ¯ z − h α (¯ z ) (cid:19) G α ( z, ¯ z ) (cid:21) G αD (¯ z, z ′ )= Z γ K d ¯ zδ ( z, ¯ z ) G αD (¯ z, z ′ )= Z γ K d ¯ z G α ( z, ¯ z ) h αD (¯ z ) G D (¯ z, z ′ ) G αD ( z, z ′ ) = Z γ K d ¯ z G α ( z, ¯ z ) h αD (¯ z ) G D (¯ z, z ′ ) (A11)We apply integration by parts and assume the disappearance of electrons at infinite distance.Inserting equation (A11) into equation (A9), the equation of motion for G D ( z, z ′ ) can beobtained as (cid:20) i ddz − h D ( z ) (cid:21) G D ( z, z ′ ) = δ ( z, z ′ ) + Z γ K d ¯ z "X α h Dα ( z ) G α ( z, ¯ z ) h αD (¯ z ) G D (¯ z, z ′ )The term P α h Dα ( z ) G α ( z, ¯ z ) h αD (¯ z ) = P α Σ α is defined as the coupling self-energy Σ ( z, ¯ z ) and the equation is reformulated as (cid:20) i ddz − h D ( z ) (cid:21) G D ( z, z ′ ) = δ ( z, z ′ ) + Z γ K d ¯ z Σ ( z, ¯ z ) G MM (¯ z, z ′ ) (A12)22 . Kadanoff-Baym equations The equations for the device’s green function G D ( z, z ′ ) are summarized as: (cid:20) i ddz − h ( z ) (cid:21) G ( z, z ′ ) = δ ( z, z ′ ) + Z γ K d ¯ z Σ ( z, ¯ z ) G (¯ z, z ′ ) (A13) G ( z, z ′ ) (cid:20) − i ddz ′ − h ( z ′ ) (cid:21) = δ ( z, z ′ ) + Z γ K d ¯ z G ( z, ¯ z ) Σ (¯ z, z ′ ) (A14)Here, the subscript D is dropped in this subsection for simplicity. Because the lesser greenfunction G < is directly related to observable physical quantities, i.e. electron densitiesand currents, its integro-differential equation is described first. Using the definition of G < ( t − , t + ) = G ( z = t − , z ′ = t + ) with t − < t + and separating the Keldysh contour byreal and imaginary segments in Eq. (A13), one gets: (cid:20) i ddt − (cid:21) G < ( t − , t + ) − h ( t − ) G < ( t − , t + )= Z Reγ K d ¯ t [ Σ ( t − , ¯ t ) G (¯ t, t + )] − i Z Imγ K dτ [ Σ ( t − , t − iτ ) G ( t − iτ, t + )] (A15)Adopting common notations f for green functions G and self-energy Σ in the Keldysh space,we arrive at: f ( t, t ′ ) | f ∈ G , Σ = f δ ( t ) δ ( t − t ′ ) + Θ ( t − t ′ ) f > ( t, t ′ ) + Θ ( t ′ − t ) f < ( t, t ′ ) (A16) f R ( t, t ′ ) | f ∈ G , Σ = f R,δ ( t ) δ ( t − t ′ ) + Θ ( t − t ′ ) [ f > ( t, t ′ ) − f < ( t, t ′ )] (A17) f A ( t, t ′ ) | f ∈ G , Σ = f A,δ ( t ) δ ( t − t ′ ) − Θ ( t ′ − t ) [ f > ( t, t ′ ) − f < ( t, t ′ )] (A18) f ⌉ ( t, τ ) | f ∈ G , Σ = f < ( t, t − iτ ) (A19) f ⌈ ( τ, t ) | f ∈ G , Σ = f > ( t − iτ, t ) (A20)Equation (A15) can be rewritten as i ddt − G < ( t − , t + ) − h ( t − ) G < ( t − , t + )= Z ∞ t d ¯ t Σ R ( t − , ¯ t ) G < (¯ t, t + ) + Z ∞ t d ¯ t Σ < ( t − , ¯ t ) G A (¯ t, t + ) − i Z β dτ Σ ⌉ ( t − , τ ) G ⌈ ( τ, t + )Alternatively, we have i ddt − G < ( t − , t + ) − h ( t − ) G < ( t − , t + ) = (cid:2) Σ R · G < + Σ < · G A + Σ ⌉ ⋆ G ⌈ (cid:3) ( t + , t − ) (A21)23ith notations [ f · g ] ( t, t ′ ) = R ∞ t d ¯ tf ( t, ¯ t ) g (¯ t, t ′ ) and [ f ⋆ g ] ( t, t ′ ) = − i R β d ¯ tf ( t, τ ) g ( τ, t ′ ).Equations for the greater green function G > can be obtained by similar processes: i ddt − G < ( t − , t + ) − h ( t − ) G < ( t − , t + ) = (cid:2) Σ R · G < + Σ < · G A + Σ ⌉ ⋆ G ⌈ (cid:3) ( t − , t + ) (A22) i ddt + G > ( t + , t − ) − h ( t + ) G > ( t + , t − ) = (cid:2) Σ R · G > + Σ > · G A + Σ ⌉ ⋆ G ⌈ (cid:3) ( t + , t − ) (A23)and, − i ddt + G < ( t − , t + ) − h ( t + ) G < ( t − , t + ) = (cid:2) G R · Σ < + G < · Σ A + G ⌉ ⋆ Σ ⌈ (cid:3) ( t − , t + )(A24) − i ddt − G > ( t + , t − ) − h ( t − ) G > ( t + , t − ) = (cid:2) G R · Σ > + G > · Σ A + G ⌉ ⋆ Σ ⌈ (cid:3) ( t + , t − )(A25)Equations (A22)-(A25) are the Kadanoff-Baym equations with symmetry relations of func-tions f ∈ G , Σ : f ≷ ( t, t ′ ) | f ∈ G , Σ = − (cid:2) f ≷ ( t ′ , t ) (cid:3) † (A26) f ⌉⌈ ( t, t ′ ) | f ∈ G , Σ = − (cid:2) f ⌉⌈ ( t ′ , t ) (cid:3) † (A27) G > ( t, t ) = − i + G < ( t, t ) , at equal time (A28) G A ( t, t ′ ) = (cid:2) G R ( t ′ , t ) (cid:3) † (A29) 2. Approximate equations for fast numerical implementation by neglectingthe complex-axis integral For the equation of motion for the retarded green function G R ( t, t ′ ), one can differentiateEq. (A17) with respect to t , ignoring the G δ δ ( t − t ′ ) term and the complex path in theKeldysh contour, i ddt G R ( t, t ′ ) = iδ ( t − t ′ ) [ G > ( t, t ′ ) − G < ( t, t ′ )]+Θ ( t − t ′ ) (cid:20) i ddt G > ( t, t ′ ) − i ddt G < ( t, t ′ ) (cid:21) (A30)Together with Eqs. (A22) and (A23), Eq. (A30) can be rewritten as i ddt G R ( t, t ′ ) = δ ( t − t ′ ) + h ( t ) Θ ( t − t ′ ) [ G > ( t, t ′ ) − G < ( t, t ′ )]+Θ ( t − t ′ ) (cid:2) Σ R · G > + Σ > · G A − Σ R · G < − Σ < · G A (cid:3) ( t, t ′ )= δ ( t − t ′ ) + h ( t ) G R + Σ R · G R (A31)24or the equation of motion for the lesser green function G < ( t, t ′ ), Eq. (A22), by ignoringcomplex integration, gives i ddt G < ( t, t ′ ) − h ( t ) G < ( t, t ′ ) = (cid:2) Σ R · G < + Σ < · G A (cid:3) ( t, t ′ ) (A32)For the evaluation of the dissipation term Q α in Eq. (4), we calculate the equation ofmotion for the lesser green function σ ( t ) = − i G < ( t, t ) in the equal-time limit. Substitut-ing Eq. (A14) from Eq. (A13) and applying the limit condition t − ≈ t + during similarderivations of Eq. (A21), one obtains i ddt G < ( t, t ) − [ h ( t ) , G < ( t, t )] = (cid:2) Σ R · G < + Σ < · G A − G R · Σ < − G < · Σ A (cid:3) ( t, t )= − (cid:2) G R · Σ < + G < · Σ A (cid:3) ( t, t ) + h.c. (A33)Here, the complex integration has been ignored and the relations in Eq. (A26) are used. Bycomparing Eq. (A33) with Eq. (3) and using σ ( t ) = − i G < ( t, t ), the dissipation term canbe given by Q = − (cid:2) G R · Σ < + G < · Σ A (cid:3) ( t, t ) + h.c. (A34) Appendix B: Wide-Band Limit approximation for the dissipation term Q α By applying the assumptions of the wide-band limit approximation, the advanced self-energy for L and R in Eq. (9) becomesΣ Aα,mn ( t, t ′ ) = i Θ ( t ′ − t ) X k α h Dα,mk α ( t ) exp ( i Z t ′ t ǫ k α + V α (¯ t ) d ¯ t ) h αD,k α n ( t ′ ) ≃ Θ ( t ′ − t ) Z ∞−∞ dǫe iǫ ( t ′ − t ) " i · h Dα,mk α ( t ) e exp ( i Z t ′ t V α (¯ t ) d ¯ t ) h αD,k α n ( t ′ ) = δ ( t ′ − t ) [Λ α,mn + i Γ α,mn ] (B1)where the matrix in the square bracket in the last line is approximated by the initial Σ Aα ( ǫ F ) at fermi level of the unbiased system . V α ( t ) is the external potential that isturned on at t > t , resulting in time-dependent level shifts of α ∈ { L, R } . The summa-tion over all single-electron states in the electrodes is replaced by an integration over theentire energy, i.e. P k α → R ∞−∞ dǫ . The retarded/advanced self-energies are Σ R,Aα ( t, t ′ ) =[Λ α,mn ∓ i Γ α,mn ] δ ( t ′ − t ). 25he lesser self-energy in Eq. (10) isΣ <α,mn ( t, t ′ ) = X k α h Dα,mk α ( t ) G <α,k α h αD,k α m ( t ′ )= X k α h Dα,mk α ( t ) h αD,k α m ( t ′ ) h i · f α ( ǫ k α ) e iǫ ( t ′ − t ) e i R t ′ t V α (¯ t ) d ¯ t i = 2 iπ Γ α,mn e i R t ′ t V α (¯ t ) d ¯ t Z ∞−∞ f α ( ǫ ) e iǫ ( t ′ − t ) dǫ (B2)From Eq. (7), the lesser green function can be solved as G RD ( t, t ′ ) = − i Θ ( t − t ′ ) e − i R t [ h D (¯ t )+ P α ( Λ α − i Γ α ) ] d ¯ t e − i R t ′ [ h D (¯ t )+ P α ( Λ α − i Γ α ) ] d ¯ t (B3)Inserting Eqs. (B1)-(B3) into Eq. (6), the dissipation term for electrodes L and R can begiven by Q α ( t ) = − Z ∞−∞ dτ (cid:2) G 2Θ ( t − τ ) π Z ∞ dτ e − i R tτ [ h D (¯ t )+ P α ( Λ α − i Γ α ) − V α (¯ t ) ] d ¯ t Z ∞−∞ f α ( ǫ ) e iǫ ( t − τ ) dǫ Γ α = − iπ U α ( t ) Z ∞−∞ f α ( ǫ ) e iǫt ǫ − h D (0) − P α ( Λ α − i Γ α ) dǫ Γ α − iπ Z ∞−∞ (cid:2) I − U α ( t ) e iǫt (cid:3) f α ( ǫ ) ǫ − h D ( t ) − P α ( Λ α − i Γ α ) + V α ( t ) I dǫ Γ α (B5)with U α ( t ) = e − i R t [ h D (¯ t )+ P α ( Λ α − i Γ α ) − V α (¯ t ) I ] d ¯ t (B6)Conclusively, the dissipation term now is Q α ( t ) = K α ( t ) + K † α ( t ) + { Γ α , σ ( t ) } + i [ Λ α , σ ( t )] (B7)with the definition of K α ( t ) in Eq. (B5). 26 EFERENCES S. J. Rosenthal, J. C. Chang, O. Kovtun, J. R. McBride, I. D. Tomlinson, Chemistry &Biology 18, 10 (2011). S. Jin, and K. Ye, Biotechnol. Prog. 23, 32 (2007). I. L. Medintz, H. Mattoussi, and A. R. Clapp, Int. J. Nanomed. 3, 151 (2008). X. Michalet, F. F. Pinaud, L. A. Bentolila, J. M. Tsay, S. Doose, J. J. Li, G. Sundaresan,A. M. Wu, S. S. Gambhir, and S. Weiss, Science 307, 538 (2005). X. Gao, Y. Cui, R. M. Levenson, L. W. K. Chung, and S. Nie, Nat. Biotech. 22, 969(2004). T. Pons, and H. Mattoussi, Ann. Biomed. Eng. 37, 1934 (2009). F. Pinaud, S. Clarke, A. Sittner, and M. Dahan, Nat. Meth. 7, 275 (2010). A. M. Derfus, W. C. W. Chan, and S. N. Bhatia, Adv. Mat. 16, 961 (2004). G. Ruan, A. Agrawal, A. I. Marcus, and S. Nie, J. Am. Chem. Soc. 129, 14759 (2007). H. Dong, T. Hou, X. Sun, Y. Li, and S. T. Lee, Appl. Phys. Lett. 103, 123115 (2013). Y. Matsumoto, A. Dutt, G. S. Rodrguez, J. S. Salazar, and M. A. Mijares, Appl. Phys.Lett. 106, 171912 (2015) S. M. Lindsay, and M. A. Ratner, Adv. Mater. 19, 23 (2007). S. Nazemi, M. Pourfath, E. A. Soleimani, and H. Kosina, J. App. Phys. 119, 144302 (2016). G. Stefanucci, C.O. Almbladh, Phys. Rev. B 69(19), 195318(2004). Y. Wang, C. Y. Yam, Th. Frauenheim, G.H. Chen, T.A. Niehaus, Chemical Physics 391,69 (2011). X. Zheng, F. Wang, C. Y. Yam, Y. Mo, and G. H. Chen, Phys. Rev. B 75, 195127 (2007). K. Burke, R. Car, and R. Gebauer, Phys. Rev. Lett. 94(14), 146803 (2005). M. Koentopp, C. Chang, K. Burke, and R. Car, J. Phys. Condens. Matter 20, 083203(2008). C. L. Cheng, J. S. Evans, and T. V. Voorhis, Phys. Rev. B 74, 155112 (2006). J. Evans, O. Vydrov, and T. Van Voorhis, J. Chem. Phys. 131, 034106 (2009). N. Garcia-Castello, S. Illera, J. D. Prades, S. Ossicini, A. Cirera, and R. Guerra, Nanoscale7, 12564 (2015). J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon and D. Sanchez-Portal, J. Phys.: Condens. Matter 14, 2745 (2002).27 P. Ordejon, E. Artacho and J. M. Soler, Phys. Rev. B: Condens. Matter 53, R10441 (1996). E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984). preparing. I. L. Ho, T. H. Chou, and Y. C. Chang, Comput. Phys. Commun. 185, 1383 (2014). J. C. Cuevas, and E. Scheer, Molecular Electronics: An Introduction to Theory and Ex-periment , World Scientific, 2010 R. Tuovinen, E. Perfetto, G. Stefanucci, and R. van Leeuwen, Phys. Rev. B 89, 085131(2014). A. P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 50, 5528 (1994). Christian Oppenlander, Time-dependent density functional tight binding combined with theLiouville-von Neumann equation applied to AC transport in molecular electronics , Disser-tation, Faculty of Physics, University of Regensburg, (2014). S. Yokojima, G. Chen, R. Xu, and Y. Yan, Chem. Phys. Lett., 369(3), 495 (2003). M. P. Lopez-Sancho, J. M. Lopez-Sancho, and J. Rubio, J. Phys. F: Met. Phys. 14, 1205(1984); 15, 851 (1985). A. Pecchia, G. Penazzi, L. Salvucci, and A. D. Carlo, New J. Phys. 10, 065022 (2008). W. M. C. Foulkes and R. Haydock, Phys. Rev. B 39 12520 (1989). M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai,and G. Seifert, Phys. Rev. B 58 7260 (1998). S. Datta, Quantum Transport: Atom to Transistor , New York, Cambridge University Press(2005). H. Grabert, and M. H. Devoret, Single Charge Tunneling: Coulomb Blockade PhenomenaIn Nanostructures , New York, Springer Science & Business Media (1992). J. A. Melsen, U. Hanke, H. O. Muller, and K. A. Chao, Phys. Rev. B 55 (1997),10638-10642. I. L. Ho, D. S. Chung, M. T. Lee, C. S. Wu, Y. C. Chang, and C. D. Chen, J. Appl. Phys.111 (2012), 064501. H. Dreysse (editor), Electronic Structure and Physical Properties of Solids: The Uses ofthe LMTO Method , p. 122, Springer-Verlag, Berlin Heidelberg (2000). G. Nazir, A. Ahmad, M. F. Khan, and S. Tariq, Comp. Cond. Mat. 4, 32 (2015). M. F. Thomas, J. M. Williams, and T. C. Gibb (editors), Hyperfine interactions (C) , p.186, Springer Science+Business Media Dordrecht (2002).28 Z. Dai, W. Jin, J. X. Yu, M. Grady, J. T. Sadowski, Y. D. Kim, J. Hone, J. I. Dadap, J.Zang, R. M. Osgood, Jr., and K. Pohl, Phys. Rev. Materials 1, 074003 (2017). W. M. Haynes, CRC Handbook of Chemistry and Physics , 96th Edition (2015). D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). J. P. Perdew and A. Zunger, Phys. Rev. B: Condens. Matter 23, 5048 (1981). M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos, Rev. Mod.Phys. 64, 1045 (1992). L. S. Wang, Phys. Chem. Chem. Phys. 12, 8694 (2010). C. G. Sanchez, E. P. M. Leiva, and W. Schmickler, Electrochem. Commun. 5(7), 584(2003). C. Sevik, and C. Bulutay, J Mater. Sci. 42, 6555 (2007). G. Conibeer, M. A. Green, D. Konig, I. Perez-Wurfl, S. Huang, X. Hao, D. Di, L. Shi, S.Shrestha, B. PuthenVeetil, Y. So, B. Zhang and Z. Wan, Prog. Photovoltaics: Res. Appl.19, 813 (2011). G. Seguini, S. Schamm-Chardon, P. Pellegrino and M. Perego, Appl. Phys. Lett. 99, 082107(2011). K. Seino, F. Bechstedt and P. Kroll, Phys. Rev. B: Condens. Matter 82, 085320 (2010). K. Seino, F. Bechstedt and P. Kroll, Phys. Rev. B: Condens. Matter 86, 075312 (2012). R. Guerra, E. Degoli and S. Ossicini, Phys. Rev. B: Condens. Matter 80, 155332 (2009). M. Mavros, D. A. Micha and D. S. Kilin, J. Phys. Chem. C 115, 19529 (2011).57