A orthogonality relation for the nonorthogonal basis

# First-principles investigation of dynamical properties of molecular devices under a steplike pulse

## Abstract

We report a computationally tractable approach to first principles investigation of time-dependent current of molecular devices under a step-like pulse. For molecular devices, all the resonant states below Fermi level contribute to the time-dependent current. Hence calculation beyond wideband limit must be carried out for a quantitative analysis of transient dynamics of molecules devices. Based on the exact non-equilibrium Green’s function (NEGF) formalism of calculating the transient current in Ref.(6), we develop two approximate schemes going beyond the wideband limit, they are all suitable for first principles calculation using the NEGF combined with density functional theory. Benchmark test has been done by comparing with the exact solution of a single level quantum dot system. Good agreement has been reached for two approximate schemes. As an application, we calculate the transient current using the first approximated formula with opposite voltage in two molecular structures: Al--Al and Al--Al. As illustrated in these examples, our formalism can be easily implemented for real molecular devices. Importantly, our new formula has captured the essential physics of dynamical properties of molecular devices and gives the correct steady state current at and .

###### pacs:
71.15.Mb, 72.30.+q, 85.35.-p 73.23.-b

## I introduction

With the rapid progress in molecular electronics,(2) quantum transport in molecular device has received increasing attention. In particular, the dynamic response of molecular devices to external parameters(3); (4); (5); (6); (7); (8); (9), in which the external time-dependent fields or internal parametric pump potentials drive the electrons to tunnel through the molecular device, is one of the most important issues in molecular electronics. The simplest molecular device structure is the two-probe lead-device-lead (LDL) configuration, where “device” is the molecular device connected to the external probes by the “leads”. In such a device, all the atomic details of the device material can be treated using density functional theory (DFT) and the non-equilibrium physics can be taken into account using non-equilibrium Green’s function (NEGF). Up to now, from an atom point of view, one of the most popular theoretical approaches used to study the quantum transport properties of molecular device is Keldysh nonequilibrium Green’s functions coupled with density-functional theory (NEGF-DFT).(10) Using this approach, the steady state quantum transport properties in molecular devices have been widely studied.

For time dependent response of molecular devices, there have been many different theoretical approaches, such as evolution of time-dependent Schrodinger equation,(11), time development operator approach,(12) and the NEGF technique.(13) These approaches are convenient to deal with dynamic response of time-dependent external field that is sinusoidal (e.g., microwave radiation). Under such an external field an electron can tunnel through the system by emitting or absorbing photons, giving rise to the photon-assisted tunneling (PAT). Concerning the steady state ac response to harmonic external field, the Floquet approach is convenient.(14) For the transient transport, however, the pulse like ac signal is the optimal driven force since they can provide a less ambiguous measure of time scales.(15) In this case, besides PAT, one of the most interesting questions to ask is how fast a device can turn on or turn off a current. With the development of molecular electronics, providing a particular viable switching device has become a key technical issue. Concerning the transient dynamics, different approaches such as path-integral techniques,(16) the solution of Wigner distribution function,(17) the time-dependent numerical renormalization group,(18) time-dependent DFT (TDDFT),(19); (3) and Keldish Green’s funciton(4); (6); (20) have also been developed and applied to different systems. Up to now, most of these approaches can only be implemented in simple systems such as quantum dots(6); (20) or one-dimension tight-binding chains.(3) Numerical calculation of transient current for molecular devices is very difficult at present stage due to the huge computational cost. This is because if we calculate the current as a function of time , the amount of calculation scales as if the time-evolution method is used. This scaling can be reduced to a linear scaling in if the wideband limit is used.(21) As we have demonstrated,(22) the wideband limit is not a good approximation for molecular devices. If one uses the exact solution from NEGF,(6) one can calculate the transient current at a particular time. However, the calculation involves a triple integration over energy which is extremely time consuming. Clearly an approximate scheme that is suitable for numerical calculation of transient properties for real molecular devices while still captured essential physics is needed.

It is the purpose of this paper to provide such a practical scheme. To study transient dynamics, in this paper, we consider a system that consists of a scattering region coupled to two leads with the external time dependent pulse bias potential . For this case, the time-dependent current for a step-like pulse has been derived exactly without using the wide-band limit by Maciejko et al(6). Since the general expression for the current involves triple integrations, it is extremely difficult to perform them in a real systems like molecules devices. So, approximation has to be made. The simplest approximation is the so called wide-band approximation where self-energies are assumed to be constants independent of energy.(23) Unfortunately, this approximation can not give the correct result since in general there are several resonant levels that significantly contribute to the transient current in molecules devices. To go beyond the wideband limit, we propose an approximate scheme of calculating the transient current that is suitable for numerical calculation of real molecules devices.(24) Our scheme is an approximation of the exact solution of Maciejko et al(6). It is very fast computationally and gives the correct limits at and . Since the exact solution of transient current is available for a single level quantum dot system, we have compared our result with the exact solution on the quantum dot system to test our approximate schemes. Good agreement is obtained. Therefore, our approximated scheme maintains essential physics of transient dynamics. Using our scheme, we calculate the transient current for the upward pulse (turn-on) and downward pulse (turn-off) in two molecular structures: Al--Al and Al--Al. We find that different from the single level quantum dot system, upon switching on the current oscillates rapidly in the first a few or tens fs with several characteristic time scales. Furthermore, due to the resonant states in molecular devices, transient currents have a much longer decay time , especially for the molecule device having a complex electronic structure such as Al--Al.

The rest of paper is organized as follows: In Sec.II, starting from the typical molecular device Hamiltonian which is expressed in an non-orthogonal basis, we shall derive a general DC and AC current expressions for a non-orthogonal basis set. It is found that for DC bias, the expressions of current for orthogonal and non-orthogonal basis sets are the same. For ac current, however, the expressions are different as will be demonstrated in Sec.II. The reason that we study the difference between orthogonal and non-orthogonal basis sets is the following. For the NEGF formalism, it is assumed that the basis set is orthogonal. It turns out that for ac transport, the current expression becomes extremely complicated if non-orthogonal basis is used. For DFT calculation, however, most people work in molecular orbitals that are non-orthogonal. Our results show that we must orthogonalizing the nonorthogonal molecular Hamiltonian, so that the present approach in Ref.(6) can be used. In Sec.III, based on the exact solution of Maciejko et al, we derive two approximate expressions for transient current with different levels of approximation. They are all suitable for numerical calculation for real molecular devices. In addition, the initial current and its asymptotic long time limit are shown to be correct. In Sec.IV, in order to appreciate our approximate formulas, we compare our result with the exact result obtained in Ref.(6) for a single-level quantum dot connected to external leads with a Lorentzian linewidth. In Sec.V, we apply our formalism to several molecular devices. Finally, a conclusion is presented in Sec.VI. Two appendices are given at the end of the paper. In Appendix A, we give a detailed derivation of orthogonalization relation for an non-orthogonal basis. This relation is used to derive the effective Green’s function which is the key to approximate exact current expression of Maciejko et al. In Appendix B, we show how to orthogonalize an nonorthogonal Hamiltonian so that the general AC current for real molecules device can be derived.

## Ii general AC current

### ii.1 Hamiltonian

The transport properties of a molecular device can be described by the following general Hamiltonian:

 H=Hc+HT+∑α=L,RHα (1)

where and describe the left and right macroscopic reservoir, respectively; is Hamiltonian of the central molecular device; couples the reservoirs to the molecular device. For a particular basis set, the above Hamiltonian can be written in the following matrix form:

 Hα = ∑μαναc†μα[H0μανα+eVα(t)δμανα]cνα Hc = ∑μcνcd†μc[H0μcνc+Uμcνc(t)]dνc HT = ∑να,νcc†ναT0νανcdνc+h.c. (2)

where is the electron charge, () and () are Fermionic annihilation (creation) operators at the state in the lead- and the state in central molecular device. , are the indices of the given basis set. The Hamiltonian of lead- are divided into two parts: the time independent part and time dependent part due to external bias on the lead-. Here we consider two kinds of step-like bias: upwards pulse (turn-on case) and downwards pulse (turn-off case) , where

 VDα(t)={Vα,   t<00,     t>0,   VUα(t)={0,     t<0Vα,   t>0 (3)

In the adiabatic approximation it is assumed that the single particle energies acquire a rigid time-dependent shift as . The energy shift in the leads is assumed to be uniform throughout. This assumption is reasonable since the pulse rising time is slower than the usual metallic plasma oscillation time, which ensures that the external electric field is effectively screened.(25)

Since Green’s function is obtained by solving Dyson equation from the known history, it is better to set time dependent external bias so that the uncertainty of future can be eliminated.(6) From Eq.(3), this is satisfied only in the downward case. In the following, we will discuss how to eliminate this uncertainty for the upward pulse. To use the Dyson equation, we will separate the Hamiltonian into two pieces: the unperturbed Hamiltonian that can be exactly resolved and the interacting term which contributes to the self energy in Dyson equations. For the downward pulse, we define the non-biased open system as the unperturbed system. It is described by the Hamiltonian . For the upward pulse, however, the situation is different, in which we will set the DC biased open system as the unperturbed Hamiltonian and set as the new time dependent part. Here denotes the coupling between scattering region and biased leads and is the induced coulomb potential due to the external bias. Now, the time dependent bias satisfies , and the uncertainty of the future in the upward case is eliminated. Then, for the downward case, we have and while for the upward case we have and . From now on we will use superscript to denote the unperturbed system that is exactly resolvable.

When the system is biased, the incoming electron will polarize the system. The induced Coulomb potential in the central scattering region consists of two parts: DC and AC parts. The DC part can be put into the exactly resolvable Hamiltonian . The induced time dependent coulomb potential due to the external bias is included as part of the non-equilibrium Hamiltonian. Because the electric field is not screened in the small scattering region where the potential drop occurs, the coulomb potential landscape in the central region is not uniform, which is different from the semi-infinite leads. Note that it is rather difficult to treat the time-dependent coulomb potential and no close formed solution exists if one does not assume wide band limit. In the small bias limit, we can expand the time-dependent coulomb potential to linear order in bias so that the analytic expression for current can be obtained. Here is the characteristic potential.(27) From the gauge invariance, (26) , and is determined from a poisson like equation.(28) In this paper, we consider the symmetric coupling so that for the external bias it is a good approximation to assume that the time dependent coulomb potential is roughly zero in the the molecular device regime.

In the following, we will derive an exact solution of transient current using a non-orthogonal basis set.(29) To facilitate the derivation, we take a unitary transformation to the Hamiltonian (2) with

 ^O(t) = exp{ie∑να∫t0dτ [~Vα(τ)c†ναcνα]}

where for the downward pulse and for the upward pulse. Note that the time in can be negative or positive, and only when . The new Hamiltonian , in which and are different from original ones:

 Hα = ∑μανα¯c†μαH0μανα¯cνα HT = ∑να,νc¯c†ναTνανc(t)dνc+h.c. (4)

where

 ¯cνα=cναexp[ie∑μα∫t0dτ ~Vα(τ)c†μαcμα], Tνανc(t)=T0νανcWα(t) Wα(t)=exp[ie∫t0~Vα(τ)dτ] (5)

For the original Hamiltonian with nonorthogonal basis, the overlap between nonorthogonal basis is expressed as the matrix form . After the unitary transform, annihilation (creation) operators () and consequently the orbital basis in the leads are changed, then overlap matrices between the leads and the scattering region become

 Sνανc(t)=S0νανcWα(t) Sνcνα(t)=W†α(t)S0νcνα. (6)

In the following, we will use the transformed Hamiltonian [Eq.(4,5), in which , are used] to derive the time dependent current expression.

### ii.2 The current

The current operator from a particular lead- to the molecular junction can be calculated from the evolution of the number operator of the electron in the semi-infinite lead-. Assuming there is no direct coupling between the left and right leads, the current operator can be expressed as:(30)

 ^Jα(t) = −e∑ναddt^Nνα(t) (7) = −e∑να[¯c†να(t)ddt¯cνα(t)+(ddt¯c†να(t))¯cνα(t)] = e∑να,νc¯c†να(t)(iTνανc(t)+Sνανc(t)ddt)dνc(t)+H.c.

where ‘H.c.’ denotes the Hermitian conjugate. The current is obtained by taking average over the nonequilibrium quantum state ’’,

 Jα(t)=e ∑να,νc [G<νcνα(t,t′)(Tνα,νc(t′)−Sνα,νc(t′)i∂∂t) − (Tνc,να(t′)−Sνc,να(t′)i´∂∂t)G<νανc(t′,t)]t=t′, (8)

where “” and “” denotes the left and right derivation respectively, and

 G<νc,να(t,t′)=i⟨¯c†να(t′)dνc(t)⟩,  G<να,νc(t′,t)=i⟨d†νc(t)¯cνα(t′)⟩.

Using the Keldysh equation and the theorem of analytic continuation, we have

 G

where

 Bcα(t1) = Tcα(t1)−Scα(t1)i∂∂t (10)

For simplicity, we have dropped the subscript , and keep only the symbol and to indicate the central scattering region and lead-, respectively. In the above expression and in the following, the summation convention on repeated sub-indices is assumed. Substituting Eq.(9) into Eq.(8), we have the general expression for the current:

 Jα(t) = −2eRe∫dt1 Tr (11) [Grcc(t,t1)Bcα(t1)g<αα(t1,t′)Bαc(t′)− G

When the system reaches a stationary state, becomes time independent, from definition Eq.(5), (6) and (10), we can find

 Bcα(t1)XBαc(t)=e−ieVα(t1−t)B0cαXB0αc,

with , where “0” denotes the zero bias system.In addition, all the propagators and depend only on the time difference . Taking the Fourier transformation, from Eq.(8) or Eq.(11), we can easily obtain DC current expressed in the energy representation:

 Jα = ∫dϵ Jα(ϵ) (12) = Re 2e∫dϵ Tr[Gr(ϵ)Σ<α(ϵ)+G<(ϵ)Σaα(ϵ)]

where and are the Green’s function and the self-energy. They have the same matrix dimension as that of the Hamiltonian . The Green’s function and self-energy is defined as

 Gr/a(ϵ)=[ϵI−Hc−Σr/a(ϵ)]−1 Σγα(ϵ)=[T0cα−ϵαS0cα]gγαα(ϵα)[T0αc−ϵαS0αc] (13)

where , is the unitary matrix with same dimension as , , and

 gr/aαα(ϵ) = [((ϵ±i0+)S0αα−H0αα)−1]να∈sur,μα∈sur g<αα(ϵ) = f(ϵ)[gaαα(ϵ)−grαα(ϵ)] (14)

is the surface Green’s function of the semi-infinite periodic lead which can be calculated numerically using a transfer matrix method.(31) Here, is the Fermi distribution. Eq.(12) shows that the dc current expressions are the same for both orthogonal and non-orthogonal basis sets.

When the time dependent field is present, however, the current expressed in energy representation will be very complicated for nonorthogonal basis due to the term in Eq.(8), since can’t be expressed as a function of time difference . One thing is clear, the transient current expressions are different for orthogonal and non-orthogonal basis sets. Instead of deriving a complicated transient current expression using a non-orthogonal basis set, we will eliminate in Eq.(8) and work on an orthogonal basis set. In Appendix B, from the overlap matrix , we derive the orthogonal basis set and new Hamiltonian expressed in this orthogonal basis. With the new orthogonal Hamiltonian, the overlap matrix will be eliminated since the overlap matrix of orthogonal basis . Then, replacing Hamiltonian in Eq.(2) with and go through the derivation leading to Eqs.(2-11) again, we arrive at a new AC current expression:

 Jα(t)=2eRe∫dt1Tr{Grcc(t,t1)[Tcα(t1)g<,exαα(t1−t)Tαc(t)]} +2eRe∫dt1Tr{G

Defining the self-energy on the orthogonal basis

 Σγ=r,a,<α(t,t′)=Tcα(t)gγ,exαα(t−t′)Tαc(t′) (16)

where is the surface Green’s function of semi-infinite lead- in the unperturbed state as defined in the Sec.II.1. For the downward pulse we have set the unperturbed system as the open system at zero bias, in which . For the upward pulse, the unperturbed system means biased open system, in which . From Eq.(15),(16), we have the general current formula

 Jα(t) = 2eRe∫dt1Tr[Gr(t,t1)Σ<α(t1,t)+G<(t,t1)Σaα(t1,t)]

At , AC external bias or time dependent part in Hamiltonian is a constant and the system is in a steady state. Consequently, the total current is known from DC transport theory that is expressed in the form of Eq.(12) but with the Green’s function and self-energy obtained from the orthogonal Hamiltonian defined above. Hence in the following we shall derive only the Ac current when . First, we shall look at the self-energy. From Eq.(5) and (16),

 Σγα(t,t′) = W†α(t)[T0cαgγαα(t,t′)T0αc]Wα(t′) = W†α(t)[∫dϵ2π eiϵ(t−t′)Σγ,exα(ϵ)]Wα(t′) = W†α(t)V†α(t)[∫dϵ2π eiϵ(t−t′)Σγ,0α(ϵ)]Vα(t′)Wα(t′)

where for the downward pulse and for the upward pulse. Here is the self-energy at zero bias, is the self-energy at the unperturbed state defined above. In the downward case ; In the upward case . Setting , and are defined in Eq.(13) with zero and nonzero , respectively. We have . From Eq.(LABEL:generalCur) and (LABEL:self2), we find

 Jα(t) = 2eRe∫dϵ2π∫t−∞dt1  eiϵ(t−t1) (19) [Gr(t,t1)~Σ<α(ϵ,t1,t)+G<(t,t1)~Σaα(ϵ,t1,t)]

where the first term is the current flowing into the molecular device while the second one is the current flowing from the molecular device, and

 ~Σγα(ϵ,t1,t) = W†α(t1)Σγ,0α(ϵ)Wα(t) (20)

where . Here is the self-energy of lead- at zero bias. The lesser Green’s function is given by

 G<(t,t′)=∫dt1∫dt2 Gr(t,t1)⎡⎣∑βΣ<β(t1,t2)⎤⎦Ga(t2,t′) (21) = ∑β∫dϵ2π e−iϵ(t−t′)[∫t−∞dt1 eiϵ(t−t1)Wβ(t)Gr(t,t1)W†β(t1)] Σ<,0β(ϵ)[∫t′−∞dt2 e−iϵ(t′−t2)Wβ(t2)Ga(t′,t2)W†β(t)]

Substitute Eq.(20) and (21) into Eq.(19) and introducing a spectrum function

 Aα(t,ϵ)=∫t−∞dt1 eiϵ(t−t1)Wα(t)Gr(t,t1)W†α(t1) (22)

we have

 Jinα(t) = 2eRe∫dϵ2π Aα(t,ϵ)Σ<,0α(ϵ) (23) Joutα(t) = 2eRe∫dϵ2π ∑βAβ(t,ϵ)Σ<,0β(ϵ)~Fβα(t,ϵ) (24)

where

 ~Fβα(t,ϵ) = ∫t−∞dt′ e−iϵ(t−t′)∫dE2π eiE(t−t′) (25) A†β(t′,ϵ)W†α(t′)Σa,0α(E)Wα(t)

Very often, is singular at , such as the quantum dot system with the wide-band limit , or the superconducting-quantum dot-normal metal system, and so on. In these cases, we should be careful with Eq.(25),

 ~Fβα(t,ϵ) = Fβα(t,ϵ)+¯Fβα(t,ϵ) (26) = (∫t−−∞+12∫t+t−)dt′ e−iϵ(t−t′)∫dE2π eiE(t−t′) A†β(t′,ϵ)W†α(t′)Σa,0α(E)Wα(t)

The first integral is the same as Eq.(25), the second integral now becomes , where we have defined

 Δr/aα = 12∫t+t−dt′ [∫dE2π Σr/a,0α(E)] (27) = 12∫t+t−dt′ Σr/a,0α(0)

Then, Eq.(24) becomes

 Joutα(t) = 2eRe∫dϵ2π ∑βAβ(t,ϵ)Σ<,0β(ϵ)Fβα(t,ϵ) (28) + 2eRe∫dϵ2π ∑βAβ(t,ϵ)Σ<,0β(ϵ)A†β(t,ϵ)Δaα

We note that Eq.(28) is the same as that derived in Ref.(6). Different from Ref.(6), we have split the expression into two terms. The first term corresponds to the non-wideband limit, i.e., when the linewidth function goes to zero at large energy. The second term of Eq.(28) is related to the wideband limit. Hence, for a quantum dot with a Lorentzian linewidth function(6), only the first term is nonzero while for the system in contact with a superconducting lead both terms are nonzero.

So far, we have discussed the ac conduction current under the time dependent bias derived from the evolution of the number operator of the electron in the semi-infinite lead-. Now we wish to address the issue of charge accumulation in the scattering region. In principle, this can be done by including the self-consistent Coulomb potential due to ac bias.(28) However, at finite voltages, there is no close form expression for ac current if Coulomb potential is included. Alternatively, one can treat Coulomb potential phenomenologically as follows. From the continuity equation, , we see that the conduction current is not a conserved quantity. In the presence of ac bias, the displacement current due to the charge pileup inside the scattering region becomes important and must be considered. Since we have neglected the Coulomb interaction in our calculation, we can use the method of current partition(32); (33) to include the displacement current. This can be done by partitioning the total displacement current into each leads giving rise to a conserving total current . For symmetric systems like what we shall study below, it is reasonable to assume that from which we find . Hence the total current is given by (25) which satisfies the current conservation .

## Iii transient AC current

Up to now, we have derived the general expression for time dependent current, Eq.(22,23,25,28) which can be used for orthogonal as well as nonorthogonal basis set. To calculate the transient current we have to solve the retarded Green’s function and integrate it over time to find and . For the pulse-like voltage , we can obtain the Green’ function by solving Dyson equation from the known history in the time domain. Depending on what is the chosen unperturbed system that can be solved exactly, the Dyson equation can be written in a different but equivalent form. In the study of time-dependent transport, it is better to treat the time-independent, open steady state system as the unperturbed system as described in Sec.II.1, and treat the time dependent part and as a perturbation. As a result, the effective self-energy , which is due to the ac bias, would have two sources: the perturbation in leads and the induced Coulomb interaction in molecular device . Then,

 Gr(t,t′) = Gr,ex(t,t′)+∫0−∞dt1 Gr,ex(t,t1)U(t1)Gr(t1,t′) + ∫dt1 dt2 Gr,ex(t,t1)[∑α¯Σrα(t1,t2)]Gr(t2,t′)

where is the response of the molecular device that is due to the Coulomb interaction when the time-dependent voltage is turned on. Here we have assumed an adiabatic response since most of time the variance of the applied electric field is much slower than the particles’ intrinsic lifetime inside the scattering region. Then we have for downward case and upward case with .

 ∫dt1 dt2=(∫0−∞dt1∫t1−∞dt2+∫t0dt1∫0−∞dt2)
 ¯Σrα(t,t′)=Σrα(t,t′)−Σr,exα(t−t′) Σr,exα(t−t′)=V†α(t)Σr,0α(t−t′)Vα(t′)

### iii.1 Exact expression of Aβ(t,ϵ) and Fβα(t,ϵ)

Following the derivations in Ref.(6), we can get the exact expression for and with the aid of the expressions and :

 ADβ(t,ϵ) = Gr,0(ϵ)+∫dE2π ei(ϵ−E)t (29) × Gr,0(E)[Z(ϵβ)−Z(ϵ)+PDGr,V(ϵβ)] FDβα(t,ϵ) = ∫dE2π Z∗(ϵ)Ga,0(ϵ)Σa,0α(E)+∫dE2π e−i(ϵ−E)t (30) × {[Z∗(ϵβ)−Z∗(ϵ)+Ga,V(ϵβ)P†D]Ga,0(E)QD(E) + [Z∗(ϵβα)Ga,V(ϵβ)−Z∗(ϵ)Ga,0(ϵ)]Σa,0α(E)}
 AUβ(t,ϵ) = Gr,V(ϵβ)+∫dE2π ei(ϵβ−E)t (31) × Gr,V(E)[Z(ϵ)−Z(ϵβ)+PUGr,0(ϵ)] FUβα(t,ϵ)