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

Yanxia Xing1,2, Bin Wang1 and Jian Wang1,∗ 1Department of Physics and the
Center of Theoretical and Computational Physics, The University of Hong Kong,
Pokfulam Road, Hong Kong, China.
2Department of Physics, Beijing Institute of Technology, Beijing 100081, China

###### 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.Maciejko, , 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
$V_{L}(t)=-V_{R}(t)$ in two molecular structures: Al-${\rm C}_{5}$-Al and
Al-${\rm C}_{60}$-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 $t=0$ and
$t\rightarrow\infty$.

###### pacs:

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

##  I introduction

With the rapid progress in molecular electronics,Book  quantum transport in
molecular device has received increasing attention. In particular, the dynamic
response of molecular devices to external parametersgross3 ; Zhu ; yip ;
Maciejko ; diventra ; chen1 ; gross2 , 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).mcdcal  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,Schro , time development operator approach,evolution  and
the NEGF technique.NEGF  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.Floquet  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.time  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,pathintegral  the solution of Wigner distribution
function,Wigner  the time-dependent numerical renormalization group,NRG  time-
dependent DFT (TDDFT),TDDFT ; gross3  and Keldish Green’s funcitonZhu ;
Maciejko ; pulse  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 dotsMaciejko ; pulse  or one-dimension tight-binding
chains.gross3  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 $t$, the
amount of calculation scales as $t^{3}$ if the time-evolution method is used.
This scaling can be reduced to a linear scaling in $t$ if the wideband limit
is used.zheng  As we have demonstrated,WangBin  the wideband limit is not a
good approximation for molecular devices. If one uses the exact solution from
NEGF,Maciejko  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 $V_{\alpha}(t)=\theta(\pm t)V_{\alpha}$. 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 alMaciejko . 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 $\Sigma^{r,a}$ are assumed to be
constants independent of energy.Jauho1  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.foot1  Our scheme is an approximation of the exact
solution of Maciejko et alMaciejko . It is very fast computationally and gives
the correct limits at $t=0$ and $t=\infty$. 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-${\rm C}_{5}$-Al
and Al-${\rm C}_{60}$-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 $\tau$, especially for the molecule device having a complex
electronic structure such as Al-${\rm C}_{60}$-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.Maciejko,  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.Maciejko,  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:

$\displaystyle H=H_{c}+H_{T}+\sum_{\alpha=L,R}H_{\alpha}$ (1)

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

$\displaystyle H_{\alpha}$ $\displaystyle=$
$\displaystyle\sum\limits_{\mu_{\alpha}\nu_{\alpha}}c^{\dagger}_{\mu_{\alpha}}\left[{\mathbf{H}}^{0}_{\mu_{\alpha}\nu_{\alpha}}+eV_{\alpha}(t)\delta_{\mu_{\alpha}\nu_{\alpha}}\right]c_{\nu_{\alpha}}$
$\displaystyle H_{c}$ $\displaystyle=$
$\displaystyle\sum\limits_{\mu_{c}\nu_{c}}d^{\dagger}_{\mu_{c}}\left[{\mathbf{H}}^{0}_{\mu_{c}\nu_{c}}+{\mathbf{U}}_{\mu_{c}\nu_{c}}(t)\right]d_{\nu_{c}}$
$\displaystyle H_{T}$ $\displaystyle=$
$\displaystyle\sum\limits_{\nu_{\alpha},\nu_{c}}c^{\dagger}_{\nu_{\alpha}}{\mathbf{T}}^{0}_{\nu_{\alpha}\nu_{c}}d_{\nu_{c}}+h.c.$
(2)

where $e$ is the electron charge, $c_{\nu_{\alpha}}$
($c_{\nu^{\dagger}_{\alpha}}$) and $d_{\nu_{c}}$ ($d^{\dagger}_{\nu_{c}}$) are
Fermionic annihilation (creation) operators at the state $\nu$ in the
lead-$\alpha$ and the state $\nu$ in central molecular device. $\nu_{\alpha}$,
$\nu_{c}$ are the indices of the given basis set. The Hamiltonian of
lead-$\alpha$ are divided into two parts: the time independent part ${\bf
H}_{\alpha}^{0}$ and time dependent part due to external bias $V_{\alpha}(t)$
on the lead-$\alpha$. Here we consider two kinds of step-like bias: upwards
pulse (turn-on case) $V_{\alpha}^{U}(t)$ and downwards pulse (turn-off case)
$V_{\alpha}^{D}(t)$, where

$\displaystyle
V^{D}_{\alpha}(t)=\left\\{\begin{array}[]{cc}V_{\alpha},~{}~{}~{}t<0\\\
0,~{}~{}~{}~{}~{}t>0\end{array}\right.,~{}~{}~{}V^{U}_{\alpha}(t)=\left\\{\begin{array}[]{cc}0,~{}~{}~{}~{}~{}t<0\\\
V_{\alpha},~{}~{}~{}t>0\end{array}\right.$ (7)

In the adiabatic approximation it is assumed that the single particle energies
acquire a rigid time-dependent shift as ${\bf H}_{\alpha}^{0}+{\bf
I}V_{\alpha}(t)$. 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.Jauho2

Since Green’s function ${\bf G}^{r}(t,t^{\prime})$ is obtained by solving
Dyson equation from the known history, it is better to set time dependent
external bias $V_{\alpha}(t>0)=0$ so that the uncertainty of future can be
eliminated.Maciejko  From Eq.(7), 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 ${\bf H}^{0}={\bf
H}^{0}_{\alpha}+{\bf H}^{0}_{c}+{\bf H}^{0}_{T}$. For the upward pulse,
however, the situation is different, in which we will set the DC biased open
system ${\bf H}^{V}=[{\bf H}^{0}_{\alpha}+V_{\alpha}{\bf I}]+[{\bf
H}^{0}_{c}+{\bf U}^{V}]+{\bf H}^{V}_{T}$ as the unperturbed Hamiltonian and
set ${\tilde{V}}^{U}_{\alpha}(t)=V^{U}_{\alpha}(t)-V_{\alpha}$ as the new time
dependent part. Here ${\bf H}^{V}_{T}$ denotes the coupling between scattering
region and biased leads and ${\bf U}^{V}$ is the induced coulomb potential due
to the external bias. Now, the time dependent bias ${\tilde{V}}^{U}$ satisfies
${\tilde{V}}^{U}(t>0)=0$, and the uncertainty of the future in the upward case
is eliminated. Then, for the downward case, we have
${\tilde{V}}^{D}_{\alpha}(t)=V^{D}_{\alpha}(t)$ and ${\bf H}^{er}={\bf H}^{0}$
while for the upward case we have
${\tilde{V}}^{U}_{\alpha}(t)=V^{U}_{\alpha}(t)-V_{\alpha}$ and ${\bf
H}^{ex}={\bf H}^{V}$. From now on we will use superscript $``ex"$ 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 ${\bf H}^{ex}$. The induced time dependent coulomb potential ${\bf
U}(t)$ due to the external bias ${\tilde{V}}_{\alpha}(t)$ 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 ${\bf U(t)}$ 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 ${\bf
U}(t)=e\sum_{\alpha}{\bf u}_{\alpha}{\tilde{V}}_{\alpha}(t)$ so that the
analytic expression for current can be obtained. Here $u_{\alpha}$ is the
characteristic potential.buttiker9  From the gauge invariance, c2gauge
$\sum_{\alpha}{\bf u}_{\alpha}={\bf I}$, and ${\bf u}_{\alpha}$ is determined
from a poisson like equation.yadong  In this paper, we consider the symmetric
coupling so that for the external bias
${\tilde{V}}_{L}(t)=-{\tilde{V}}_{R}(t)$ it is a good approximation to assume
that the time dependent coulomb potential $U(t)$ 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.foot2  To facilitate the derivation, we take a
unitary transformation $\hat{O}(t)$ to the Hamiltonian (2) with

$\displaystyle\hat{O}(t)$ $\displaystyle=$ $\displaystyle{\rm
exp}\left\\{ie\sum_{\nu_{\alpha}}\int_{0}^{t}d\tau~{}\left[{\tilde{V}}_{\alpha}(\tau)c^{\dagger}_{\nu_{\alpha}}c_{\nu_{\alpha}}\right]\right\\}$

where ${\tilde{V}}_{\alpha}(\tau)=\theta(-\tau)V_{\alpha}$ for the downward
pulse and ${\tilde{V}}_{\alpha}(\tau)=-\theta(-\tau)V_{\alpha}$ for the upward
pulse. Note that the time $t$ in $\hat{O}(t)$ can be negative or positive, and
$\hat{O}(t)=1$ only when $t>0$. The new Hamiltonian
${\mathcal{H}}=\hat{O}H\hat{O}^{\dagger}(t)+i(\frac{\partial}{\partial
t}\hat{O}(t))\hat{O}^{\dagger}(t)$, in which ${\mathcal{H}}_{\alpha}$ and
${\mathcal{H}}_{T}$ are different from original ones:

$\displaystyle{\mathcal{H}}_{\alpha}$ $\displaystyle=$
$\displaystyle\sum\limits_{\mu_{\alpha}\nu_{\alpha}}{\bar{c}}^{\dagger}_{\mu_{\alpha}}{\mathbf{H}}^{0}_{\mu_{\alpha}\nu_{\alpha}}{\bar{c}}_{\nu_{\alpha}}$
$\displaystyle{\mathcal{H}}_{T}$ $\displaystyle=$
$\displaystyle\sum\limits_{\nu_{\alpha},\nu_{c}}{\bar{c}}^{\dagger}_{\nu_{\alpha}}{\mathbf{T}}_{\nu_{\alpha}\nu_{c}}(t){d}_{\nu_{c}}+h.c.$
(8)

where

$\displaystyle{\bar{c}}_{\nu_{\alpha}}=c_{\nu_{\alpha}}\exp[{ie\sum_{\mu_{\alpha}}\int_{0}^{t}d\tau~{}{\tilde{V}}_{\alpha}(\tau){c}^{\dagger}_{\mu_{\alpha}}{c}_{\mu_{\alpha}}}],$
$\displaystyle{\mathbf{T}}_{\nu_{\alpha}\nu_{c}}(t)={\mathbf{T}}^{0}_{\nu_{\alpha}\nu_{c}}{\mathfrak{W}}_{\alpha}(t)$
$\displaystyle{\mathfrak{W}}_{\alpha}(t)=\exp[{ie\int_{0}^{t}{\tilde{V}}_{\alpha}(\tau)d\tau}]$
(9)

For the original Hamiltonian with nonorthogonal basis, the overlap between
nonorthogonal basis is expressed as the matrix form ${\bf
S}^{0}_{\mu\nu}=\langle\mu|\nu\rangle$. After the unitary transform,
annihilation (creation) operators $c_{\alpha}$ ($c^{\dagger}_{\alpha}$) and
consequently the orbital basis $\mu_{\alpha}$ in the leads are changed, then
overlap matrices between the leads and the scattering region become

$\displaystyle{\bf S}_{\nu_{\alpha}\nu_{c}}(t)={\bf
S}^{0}_{\nu_{\alpha}\nu_{c}}{\mathfrak{W}}_{\alpha}(t)$ $\displaystyle{\bf
S}_{\nu_{c}\nu_{\alpha}}(t)={\mathfrak{W}}^{\dagger}_{\alpha}(t){\bf
S}^{0}_{\nu_{c}\nu_{\alpha}}.$ (10)

In the following, we will use the transformed Hamiltonian [Eq.(8,9), in which
${\bar{c}}_{\nu_{\alpha}}$, ${d}_{\nu_{c}}$ are used] to derive the time
dependent current expression.

###  II.2 The current

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

$\displaystyle\hat{J}_{\alpha}(t)$ $\displaystyle=$
$\displaystyle-e\sum_{\nu_{\alpha}}\frac{d}{dt}\hat{N}_{\nu_{\alpha}}(t)$ (11)
$\displaystyle=$
$\displaystyle-e\sum_{\nu_{\alpha}}\left[{\bar{c}}^{\dagger}_{\nu_{\alpha}}(t)\frac{d}{dt}{\bar{c}}_{\nu_{\alpha}}(t)+\left(\frac{d}{dt}{\bar{c}}^{\dagger}_{\nu_{\alpha}}(t)\right){\bar{c}}_{\nu_{\alpha}}(t)\right]$
$\displaystyle=$ $\displaystyle
e\sum_{\nu_{\alpha},\nu_{c}}{\bar{c}}^{\dagger}_{\nu_{\alpha}}(t)\left(i{\mathbf{T}}_{\nu_{\alpha}\nu_{c}}(t)+{\mathbf{S}}_{\nu_{\alpha}\nu_{c}}(t)\frac{d}{dt}\right)d_{\nu_{c}}(t)+H.c.$

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

$\displaystyle J_{\alpha}(t)=e~{}\sum_{\nu_{\alpha},\nu_{c}}$
$\displaystyle\left[\mathbf{G}^{<}_{\nu_{c}\nu_{\alpha}}(t,t^{\prime})\left(\mathbf{T}_{\nu_{\alpha},\nu_{c}}(t^{\prime})-{\bf
S}_{\nu_{\alpha},\nu_{c}}(t^{\prime})i\grave{\frac{\partial}{\partial
t}}\right)\right.$ $\displaystyle-$
$\displaystyle\left.\left(\mathbf{T}_{\nu_{c},\nu_{\alpha}}(t^{\prime})-{\bf
S}_{\nu_{c},\nu_{\alpha}}(t^{\prime})i\acute{\frac{\partial}{\partial
t}}\right)\mathbf{G}^{<}_{\nu_{\alpha}\nu_{c}}(t^{\prime},t)\right]_{t=t^{\prime}},$
(12)

where “$\grave{\frac{\partial}{\partial t}}$” and
“$\acute{\frac{\partial}{\partial t}}$” denotes the left and right derivation
respectively, and

${\bf
G}^{<}_{\nu_{c},\nu_{\alpha}}(t,t^{\prime})=i\left\langle{\bar{c}}^{\dagger}_{\nu_{\alpha}}(t^{\prime})d_{\nu_{c}}(t)\right\rangle,~{}~{}{\bf
G}^{<}_{\nu_{\alpha},\nu_{c}}(t^{\prime},t)=i\left\langle
d^{\dagger}_{\nu_{c}}(t){\bar{c}}_{\nu_{\alpha}}(t^{\prime})\right\rangle.$

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

$\displaystyle{\bf G}^{<}_{c\alpha}(t,t^{\prime})=\int dt_{1}$
$\displaystyle\left[{\bf G}^{r}_{cc}(t,t_{1}){\bf B}_{c\alpha}(t_{1}){\bf
g}^{<}_{\alpha\alpha}(t_{1},t^{\prime})+\right.$ (13) $\displaystyle\left.{\bf
G}^{<}_{cc}(t,t_{1}){\bf B}_{c\alpha}(t_{1}){\bf
g}^{a}_{\alpha\alpha}(t_{1},t^{\prime})\right]$

where

$\displaystyle{\bf B}_{c\alpha}(t_{1})$ $\displaystyle=$ $\displaystyle{\bf
T}_{c\alpha}(t_{1})-{\bf S}_{c\alpha}(t_{1})i\grave{\frac{\partial}{\partial
t}}$ (14)

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

$\displaystyle J_{\alpha}(t)$ $\displaystyle=$ $\displaystyle-2e{\rm Re}\int
dt_{1}~{}{\rm Tr}$ (15) $\displaystyle\left[{\bf G}^{r}_{cc}(t,t_{1}){\bf
B}_{c\alpha}(t_{1}){\bf g}^{<}_{\alpha\alpha}(t_{1},t^{\prime}){\bf B}_{\alpha
c}(t^{\prime})-\right.$ $\displaystyle\left.{\bf G}^{<}_{cc}(t,t_{1}){\bf
B}_{c\alpha}(t_{1}){\bf g}^{a}_{\alpha\alpha}(t_{1},t^{\prime}){\bf B}_{\alpha
c}(t^{\prime})\right]_{t=t^{\prime}}$

When the system reaches a stationary state, $V_{\alpha}(t)=V_{\alpha}$ becomes
time independent, from definition Eq.(9), (10) and (14), we can find

${\bf B}_{c\alpha}(t_{1})X{\bf B}_{\alpha c}(t)=e^{-ieV_{\alpha}(t_{1}-t)}{\bf
B}^{0}_{c\alpha}X{\bf B}^{0}_{\alpha c},$

with ${\bf B}^{0}_{c\alpha/\alpha c}={\bf T}^{0}_{c\alpha/\alpha
c}-i\grave{\frac{\partial}{\partial t}}{\bf S}^{0}_{c\alpha/\alpha c}$, where
“0” denotes the zero bias system.In addition, all the propagators ${\bf G}$
and ${\bf g}$ depend only on the time difference $t_{1}-t$. Taking the Fourier
transformation, from Eq.(12) or Eq.(15), we can easily obtain DC current
expressed in the energy representation:

$\displaystyle J_{\alpha}$ $\displaystyle=$ $\displaystyle\int
d\epsilon~{}\mathcal{J}_{\alpha}(\epsilon)$ (16) $\displaystyle=$
$\displaystyle\rm{Re}~{}2e\int d\epsilon~{}{\rm
Tr}\left[\mathbf{G}^{r}(\epsilon)\mathbf{\Sigma}^{<}_{\alpha}(\epsilon)+\mathbf{G}^{<}(\epsilon)\mathbf{\Sigma}^{a}_{\alpha}(\epsilon)\right]$

where ${\mathbf{G}}$ and ${\mathbf{\Sigma}}$ are the Green’s function and the
self-energy. They have the same matrix dimension as that of the Hamiltonian
${\bf H}_{c}$. The Green’s function ${\bf G}^{r/a}$ and self-energy
${\bf\Sigma}^{r/a}$ is defined as

$\displaystyle{\bf G}^{r/a}(\epsilon)=\left[\epsilon{\bf I}-{\bf
H}_{c}-{\bf\Sigma}^{r/a}(\epsilon)\right]^{-1}$
$\displaystyle\mathbf{\Sigma}^{\gamma}_{\alpha}(\epsilon)=\left[\mathbf{T}^{0}_{c\alpha}-\epsilon^{\alpha}{\bf
S}^{0}_{c\alpha}\right]\mathbf{g}^{\gamma}_{\alpha\alpha}(\epsilon^{\alpha})\left[\mathbf{T}^{0}_{\alpha
c}-\epsilon^{\alpha}{\bf S}^{0}_{\alpha c}\right]$ (17)

where $\epsilon^{\alpha}=\epsilon-eV_{\alpha}$, ${\bf I}$ is the unitary
matrix with same dimension as ${\bf H}_{c}$, $\gamma=r,a,<$, and

$\displaystyle{\bf g}^{r/a}_{\alpha\alpha}(\epsilon)$ $\displaystyle=$
$\displaystyle\left[\left((\epsilon\pm i0^{+}){\bf S}^{0}_{\alpha\alpha}-{\bf
H}^{0}_{\alpha\alpha}\right)^{-1}\right]_{\nu_{\alpha}\in{\rm
sur},\mu_{\alpha}\in{\rm sur}}$ $\displaystyle{\bf
g}^{<}_{\alpha\alpha}(\epsilon)$ $\displaystyle=$ $\displaystyle
f(\epsilon)\left[{\bf g}^{a}_{\alpha\alpha}(\epsilon)-{\bf
g}^{r}_{\alpha\alpha}(\epsilon)\right]$ (18)

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

When the time dependent field $V_{\alpha}(t)$ is present, however, the current
expressed in energy representation will be very complicated for nonorthogonal
basis due to the term ${\bf S}(t^{\prime})i\frac{\partial}{\partial t}$ in
Eq.(12), since ${\bf B}(t_{1})X{\bf B}(t)$ can’t be expressed as a function of
time difference $t_{1}-t$. 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 ${\bf S}_{c\alpha/\alpha
c}(t^{\prime})i\frac{\partial}{\partial t}$ in Eq.(12) and work on an
orthogonal basis set. In Appendix B, from the overlap matrix ${\bf S}$, we
derive the orthogonal basis set and new Hamiltonian ${\tilde{H}}$ expressed in
this orthogonal basis. With the new orthogonal Hamiltonian, the overlap matrix
${\bf S}_{c\alpha/\alpha c}(t^{\prime})$ will be eliminated since the overlap
matrix of orthogonal basis ${\bf S}^{orth}={\bf I}$. Then, replacing
Hamiltonian ${\bf H}$ in Eq.(2) with $\tilde{\bf H}$ and go through the
derivation leading to Eqs.(2-15) again, we arrive at a new AC current
expression:

$\displaystyle J_{\alpha}(t)=2e{\rm Re}\int dt_{1}{\rm
Tr}\left\\{\mathbf{G}_{cc}^{r}(t,t_{1})\left[{\bf T}_{c\alpha}(t_{1}){\bf
g}^{<,ex}_{\alpha\alpha}(t_{1}-t){\bf T}_{\alpha c}(t)\right]\right\\}$
$\displaystyle+2e{\rm Re}\int dt_{1}{\rm
Tr}\left\\{\mathbf{G}^{<}_{cc}(t,t_{1})\left[{\bf T}_{c\alpha}(t_{1}){\bf
g}^{a,ex}_{\alpha\alpha}(t_{1}-t){\bf T}_{\alpha c}(t)\right]\right\\}$ (19)

Defining the self-energy on the orthogonal basis

$\displaystyle{\bf\Sigma}^{\gamma=r,a,<}_{\alpha}(t,t^{\prime})={\bf
T}_{c\alpha}(t){\bf g}^{\gamma,ex}_{\alpha\alpha}(t-t^{\prime}){\bf T}_{\alpha
c}(t^{\prime})$ (20)

where ${\bf
g}^{\gamma,ex}_{\alpha\alpha}(t-t^{\prime})=\int\frac{d\epsilon}{2\pi}~{}e^{-i\epsilon(t-t^{\prime})}{\bf
g}^{\gamma,ex}_{\alpha\alpha}(\epsilon)$ is the surface Green’s function of
semi-infinite lead-$\alpha$ 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 ${\bf
g}^{\gamma,ex}_{\alpha\alpha}(\epsilon)=\left[\epsilon-H^{0}_{\alpha}+i0^{+}\right]^{-1}_{\alpha\in{\rm
sur}}$. For the upward pulse, the unperturbed system means $V_{\alpha}$ biased
open system, in which ${\bf
g}^{\gamma,eq}_{\alpha\alpha}(\epsilon)=\left[\epsilon-
eV_{\alpha}-H^{0}_{\alpha}+i0^{+}\right]^{-1}_{\alpha\in{\rm sur}}$. From
Eq.(19),(20), we have the general current formula

$\displaystyle J_{\alpha}(t)$ $\displaystyle=$ $\displaystyle 2e{\rm Re}\int
dt_{1}{\rm
Tr}\left[\mathbf{G}^{r}(t,t_{1})\Sigma^{<}_{\alpha}(t_{1},t)+\mathbf{G}^{<}(t,t_{1})\Sigma^{a}_{\alpha}(t_{1},t)\right]$

At $t<0$, AC external bias $V_{\alpha}(t)$ or time dependent part in
Hamiltonian ${\tilde{V}}_{\alpha}(t)$ 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.(16) 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 $t>0$. First, we
shall look at the self-energy. From Eq.(9) and (20),

$\displaystyle\mathbf{\Sigma}^{\gamma}_{\alpha}(t,t^{\prime})$
$\displaystyle=$
$\displaystyle\mathfrak{W}^{\dagger}_{\alpha}(t)\left[\mathbf{T}^{0}_{c\alpha}\mathbf{g}_{\alpha\alpha}^{\gamma}(t,t^{\prime})T^{0}_{\alpha
c}\right]\mathfrak{W}_{\alpha}(t^{\prime})$ $\displaystyle=$
$\displaystyle\mathfrak{W}^{\dagger}_{\alpha}(t)\left[\int\frac{d\epsilon}{2\pi}~{}e^{i\epsilon(t-t^{\prime})}\mathbf{\Sigma}_{\alpha}^{\gamma,ex}(\epsilon)\right]\mathfrak{W}_{\alpha}(t^{\prime})$
$\displaystyle=$
$\displaystyle\mathfrak{W}^{\dagger}_{\alpha}(t)\mathfrak{V}^{\dagger}_{\alpha}(t)\left[\int\frac{d\epsilon}{2\pi}~{}e^{i\epsilon(t-t^{\prime})}\mathbf{\Sigma}_{\alpha}^{\gamma,0}(\epsilon)\right]\mathfrak{V}_{\alpha}(t^{\prime})\mathfrak{W}_{\alpha}(t^{\prime})$

where $\mathfrak{V}_{\alpha}(t)=1$ for the downward pulse and
$\mathfrak{V}_{\alpha}(t)=e^{ieV_{\alpha}t}$ for the upward pulse. Here
$\mathbf{\Sigma}_{\alpha}^{\gamma,0}(\epsilon)$ is the self-energy at zero
bias,
$\mathbf{\Sigma}_{\alpha}^{\gamma,ex}(\epsilon)=\mathbf{T}^{0}_{c\alpha}\mathbf{g}_{\alpha\alpha}^{\gamma,ex}(\epsilon)\mathbf{T}^{0}_{\alpha
c}$ is the self-energy at the unperturbed state defined above. In the downward
case ${\bf\Sigma}^{\gamma,ex}_{\alpha}={\bf\Sigma}^{\gamma,0}_{\alpha}$; In
the upward case
${\bf\Sigma}^{\gamma,ex}_{\alpha}={\bf\Sigma}^{\gamma,V}_{\alpha}$. Setting
${\bf S}^{0}_{\alpha c}={\bf S}^{0}_{c\alpha}=0$,
${\bf\Sigma}^{\gamma,0}_{\alpha}$ and ${\bf\Sigma}^{\gamma,V}_{\alpha}$ are
defined in Eq.(17) with zero and nonzero $V_{\alpha}$, respectively. We have
${\bf\Sigma}^{r/a,V}_{\alpha}(\epsilon)={\bf\Sigma}^{r/a,0}_{\alpha}(\epsilon-
eV_{\alpha})$. From Eq.(LABEL:generalCur) and (LABEL:self2), we find

$\displaystyle J_{\alpha}(t)$ $\displaystyle=$ $\displaystyle
2e\rm{Re}\int\frac{d\epsilon}{2\pi}\int^{t}_{-\infty}dt_{1}~{}~{}e^{i\epsilon(t-t_{1})}$
(23)
$\displaystyle\left[\mathbf{G}^{r}(t,t_{1})\tilde{\mathbf{\Sigma}}^{<}_{\alpha}(\epsilon,t_{1},t)+\mathbf{G}^{<}(t,t_{1})\tilde{\mathbf{\Sigma}}^{a}_{\alpha}(\epsilon,t_{1},t)\right]$

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

$\displaystyle\tilde{\mathbf{\Sigma}}^{\gamma}_{\alpha}(\epsilon,t_{1},t)$
$\displaystyle=$
$\displaystyle\mathcal{W}_{\alpha}^{\dagger}(t_{1})\mathbf{\Sigma}^{\gamma,0}_{\alpha}(\epsilon)\mathcal{W}_{\alpha}(t)$
(24)

where
$\mathcal{W}_{\alpha}(t)=\mathfrak{V}_{\alpha}(t)\mathfrak{W}_{\alpha}(t)$.
Here ${\bf\Sigma}^{\gamma,0}_{\alpha\alpha}$ is the self-energy of
lead-$\alpha$ at zero bias. The lesser Green’s function is given by

$\displaystyle\mathbf{G}^{<}(t,t^{\prime})=\int dt_{1}\int
dt_{2}~{}\mathbf{G}^{r}(t,t_{1})\left[\sum_{\beta}\mathbf{\bf{\Sigma}}^{<}_{\beta}(t_{1},t_{2})\right]\mathbf{G}^{a}(t_{2},t^{\prime})$
(25) $\displaystyle=$
$\displaystyle\sum_{\beta}\int\frac{d\epsilon}{2\pi}~{}e^{-i\epsilon(t-t^{\prime})}\left[\int^{t}_{-\infty}dt_{1}~{}e^{i\epsilon(t-t_{1})}\mathcal{W}_{\beta}(t)\mathbf{G}^{r}(t,t_{1})\mathcal{W}^{\dagger}_{\beta}(t_{1})\right]$
$\displaystyle\mathbf{\Sigma}^{<,0}_{\beta}(\epsilon)\left[\int^{t^{\prime}}_{-\infty}dt_{2}~{}e^{-i\epsilon(t^{\prime}-t_{2})}\mathcal{W}_{\beta}(t_{2})\mathbf{G}^{a}(t^{\prime},t_{2})\mathcal{W}^{\dagger}_{\beta}(t)\right]$

Substitute Eq.(24) and (25) into Eq.(23) and introducing a spectrum function

$\displaystyle{\bf
A}_{\alpha}(t,\epsilon)=\int_{-\infty}^{t}dt_{1}~{}e^{i\epsilon(t-t_{1})}\mathcal{W}_{\alpha}(t){\bf
G}^{r}(t,t_{1})\mathcal{W}^{\dagger}_{\alpha}(t_{1})$ (26)

we have

$\displaystyle J^{in}_{\alpha}(t)$ $\displaystyle=$ $\displaystyle
2e\rm{Re}\int\frac{d\epsilon}{2\pi}~{}{\bf
A}_{\alpha}(t,\epsilon){\bf\Sigma}_{\alpha}^{<,0}(\epsilon)$ (27)
$\displaystyle J^{out}_{\alpha}(t)$ $\displaystyle=$ $\displaystyle
2e\rm{Re}\int\frac{d\epsilon}{2\pi}~{}\sum_{\beta}{\bf
A}_{\beta}(t,\epsilon){\bf\Sigma}^{<,0}_{\beta}(\epsilon)\tilde{\bf
F}_{\beta\alpha}(t,\epsilon)$ (28)

where

$\displaystyle\tilde{\bf F}_{\beta\alpha}(t,\epsilon)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{t}dt^{\prime}~{}e^{-i\epsilon(t-t^{\prime})}\int\frac{dE}{2\pi}~{}e^{iE(t-t^{\prime})}$
(29) $\displaystyle{\bf
A}_{\beta}^{\dagger}(t^{\prime},\epsilon){\mathcal{W}}_{\alpha}^{\dagger}(t^{\prime}){\bf\Sigma}_{\alpha}^{a,0}(E){\mathcal{W}}_{\alpha}(t)$

Very often, ${\bf\Sigma}^{r/a}(t-t^{\prime})$ is singular at $t=t^{\prime}$,
such as the quantum dot system with the wide-band limit
${\bf\Sigma}^{r/a}(0)=\int\frac{dE}{2\pi}~{}{\bf\Sigma}^{r/a}(E)=\delta(0)(\mp\Gamma/2)$,
or the superconducting-quantum dot-normal metal system, and so on. In these
cases, we should be careful with Eq.(29),

$\displaystyle\tilde{\bf F}_{\beta\alpha}(t,\epsilon)$ $\displaystyle=$
$\displaystyle{\bf F}_{\beta\alpha}(t,\epsilon)+\bar{\bf
F}_{\beta\alpha}(t,\epsilon)$ (30) $\displaystyle=$
$\displaystyle\left(\int_{-\infty}^{t^{-}}+\frac{1}{2}\int_{t^{-}}^{t^{+}}\right)dt^{\prime}~{}e^{-i\epsilon(t-t^{\prime})}\int\frac{dE}{2\pi}~{}e^{iE(t-t^{\prime})}$
$\displaystyle{\bf
A}_{\beta}^{\dagger}(t^{\prime},\epsilon){\mathcal{W}}_{\alpha}^{\dagger}(t^{\prime}){\bf\Sigma}_{\alpha}^{a,0}(E){\mathcal{W}}_{\alpha}(t)$

The first integral $\int_{-\infty}^{t^{-}}$ is the same as Eq.(29), the second
integral $\frac{1}{2}\int_{t^{-}}^{t^{+}}$ now becomes $\bar{\bf
F}_{\beta\alpha}(t,\epsilon)={\bf
A}_{\beta}^{\dagger}(t,\epsilon){\bf\Delta}^{a}_{\alpha}$, where we have
defined

$\displaystyle{\bf\Delta}^{r/a}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{t^{-}}^{t^{+}}dt^{\prime}~{}\left[\int\frac{dE}{2\pi}~{}{\bf\Sigma}_{\alpha}^{r/a,0}(E)\right]$
(31) $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{t^{-}}^{t^{+}}dt^{\prime}~{}{\bf\Sigma}_{\alpha}^{r/a,0}(0)$

Then, Eq.(28) becomes

$\displaystyle J^{out}_{\alpha}(t)$ $\displaystyle=$ $\displaystyle
2e\rm{Re}\int\frac{d\epsilon}{2\pi}~{}\sum_{\beta}{\bf
A}_{\beta}(t,\epsilon){\bf\Sigma}^{<,0}_{\beta}(\epsilon){\bf
F}_{\beta\alpha}(t,\epsilon)$ (32) $\displaystyle+$ $\displaystyle
2e\rm{Re}\int\frac{d\epsilon}{2\pi}~{}\sum_{\beta}{\bf
A}_{\beta}(t,\epsilon){\bf\Sigma}^{<,0}_{\beta}(\epsilon){\bf
A}^{\dagger}_{\beta}(t,\epsilon){\bf\Delta}^{a}_{\alpha}$

We note that Eq.(32) is the same as that derived in Ref.Maciejko, . Different
from Ref.Maciejko, , we have split the expression into two terms. The first
term corresponds to the non-wideband limit, i.e., when the linewidth function
$\bf\Gamma$ goes to zero at large energy. The second term of Eq.(32) is
related to the wideband limit. Hence, for a quantum dot with a Lorentzian
linewidth functionMaciejko , 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 $J_{\alpha}(t)$ under the
time dependent bias derived from the evolution of the number operator of the
electron in the semi-infinite lead-$\alpha$. 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.yadong
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,
$\sum_{\alpha}J_{\alpha}(t)+dQ(t)/dt=0$, we see that the conduction current is
not a conserved quantity. In the presence of ac bias, the displacement current
$J_{\alpha}^{d}$ due to the charge pileup $dQ/dt$ 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
partitionbuttiker4 ; wbg  to include the displacement current. This can be
done by partitioning the total displacement current
$\sum_{\alpha}J^{d}_{\alpha}=dQ/dt$ into each leads giving rise to a
conserving total current $I_{\alpha}=J_{\alpha}+J^{d}_{\alpha}$. For symmetric
systems like what we shall study below, it is reasonable to assume that
$J^{d}_{L}=J^{d}_{R}$ from which we find $J^{d}_{\alpha}=-(J_{L}+J_{R})/2$.
Hence the total current is given by $I_{L}=(J_{L}-J_{R})/2$Jauho2  which
satisfies the current conservation $I_{L}+I_{R}=0$.

##  III transient AC current

Up to now, we have derived the general expression for time dependent current,
Eq.(26,27,29,32) 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 ${\bf G}^{r}(t,t^{\prime})$ and integrate it over time to
find ${\bf A}_{\beta}(t,\epsilon)$ and $\tilde{\bf
F}_{\beta\alpha}(t,\epsilon)$. For the pulse-like voltage
${\tilde{V}}_{\alpha}(t)=\pm\theta(-t)$, we can obtain the Green’ function
${\bf G}^{r}(t,t^{\prime})$ by solving Dyson equation ${\bf G}^{r}={\bf
G}^{r,eq}+{\bf G}^{r,eq}{\bf\Xi}{\bf G}^{r}$ 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
$\tilde{V}_{\alpha}(t)$ and ${\bf U}(t)$ as a perturbation. As a result, the
effective self-energy ${\bf\Xi}$, which is due to the ac bias, would have two
sources: the perturbation in leads $\bar{\bf\Sigma}^{r}_{\alpha}$ and the
induced Coulomb interaction in molecular device ${\bf U}(t)$. Then,

$\displaystyle\mathbf{G}^{r}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\mathbf{G}^{r,ex}(t,t^{\prime})+\int_{-\infty}^{0}dt_{1}~{}\mathbf{G}^{r,ex}(t,t_{1})\mathbf{U}(t_{1})\mathbf{G}^{r}(t_{1},t^{\prime})$
$\displaystyle+$ $\displaystyle~{}\int
dt_{1}~{}dt_{2}~{}\mathbf{G}^{r,ex}(t,t_{1})\left[\sum_{\alpha}\mathbf{\bar{\Sigma}}^{r}_{\alpha}(t_{1},t_{2})\right]\mathbf{G}^{r}(t_{2},t^{\prime})$

where ${\bf U}(t)$ 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 ${\bf U}(t)=\pm{\bf U}\theta(-t)$ for
downward case and upward case with ${\bf U}={\bf H}^{V}_{c}-{\bf H}^{0}_{c}$.

$\int
dt_{1}~{}dt_{2}=\left(\int^{0}_{-\infty}dt_{1}\int^{t_{1}}_{-\infty}dt_{2}+\int^{t}_{0}dt_{1}\int^{0}_{-\infty}dt_{2}\right)$
$\displaystyle\mathfrak{\bar{{\bf\Sigma}}}^{r}_{\alpha}(t,t^{\prime})={\bf\Sigma}^{r}_{\alpha}(t,t^{\prime})-{\bf\Sigma}_{\alpha}^{r,ex}(t-t^{\prime})$
$\displaystyle{\bf\Sigma}_{\alpha}^{r,ex}(t-t^{\prime})={\mathfrak{V}}^{\dagger}_{\alpha}(t){\bf\Sigma}_{\alpha}^{r,0}(t-t^{\prime}){\mathfrak{V}}_{\alpha}(t^{\prime})$

###  III.1 Exact expression of ${\bf A}_{\beta}(t,\epsilon)$ and ${\bf
F}_{\beta\alpha}(t,\epsilon)$

Following the derivations in Ref.Maciejko, , we can get the exact expression
for ${\bf A}_{\beta}(t,\epsilon)$ and ${\bf F}_{\beta\alpha}(t,\epsilon)$ with
the aid of the expressions $\epsilon_{\beta}=\epsilon+eV_{\beta}$ and
$\epsilon_{\beta\alpha}=\epsilon+eV_{\beta}-eV_{\alpha}$:

$\displaystyle{\bf A}^{D}_{\beta}(t,\epsilon)$ $\displaystyle=$
$\displaystyle{\bf
G}^{r,0}(\epsilon)+\int\frac{dE}{2\pi}~{}e^{i(\epsilon-E)t}$ (33)
$\displaystyle\times$ $\displaystyle{\bf
G}^{r,0}(E)\left[Z(\epsilon_{\beta})-Z(\epsilon)+{\bf P}_{D}{\bf
G}^{r,V}(\epsilon_{\beta})\right]$ $\displaystyle{\bf
F}^{D}_{\beta\alpha}(t,\epsilon)$ $\displaystyle=$
$\displaystyle\int\frac{dE}{2\pi}~{}Z^{*}(\epsilon){\bf
G}^{a,0}(\epsilon){\bf\Sigma}^{a,0}_{\alpha}(E)+\int\frac{dE}{2\pi}~{}e^{-i(\epsilon-E)t}$
(34) $\displaystyle\times$
$\displaystyle\left\\{\left[Z^{*}(\epsilon_{\beta})-Z^{*}(\epsilon)+{\bf
G}^{a,V}(\epsilon_{\beta}){\bf P}^{\dagger}_{D}\right]{\bf G}^{a,0}(E){\bf
Q}_{D}(E)\right.$ $\displaystyle+$
$\displaystyle\left.\left[Z^{*}(\epsilon_{\beta\alpha}){\bf
G}^{a,V}(\epsilon_{\beta})-Z^{*}(\epsilon){\bf
G}^{a,0}(\epsilon)\right]{\bf\Sigma}^{a,0}_{\alpha}(E)\right\\}$
$\displaystyle{\bf A}^{U}_{\beta}(t,\epsilon)$ $\displaystyle=$
$\displaystyle{\bf
G}^{r,V}(\epsilon_{\beta})+\int\frac{dE}{2\pi}~{}e^{i(\epsilon_{\beta}-E)t}$
(35) $\displaystyle\times$ $\displaystyle{\bf
G}^{r,V}(E)\left[Z(\epsilon)-Z(\epsilon_{\beta})+{\bf P}_{U}{\bf
G}^{r,0}(\epsilon)\right]$ $\displaystyle{\bf
F}^{U}_{\beta\alpha}(t,\epsilon)$ $\displaystyle=$
$\displaystyle\int\frac{dE}{2\pi}~{}Z^{*}(\epsilon_{\beta\alpha}){\bf
G}^{a,V}(\epsilon_{\beta}){\bf\Sigma}^{a,0}_{\alpha}(E)+\int\frac{dE}{2\pi}~{}e^{-i(\epsilon_{\beta}-E)t}$
$\displaystyle\times$
$\displaystyle\left\\{\left[Z^{*}(\epsilon)-Z^{*}(\epsilon_{\beta})+{\bf
G}^{a,0}(\epsilon){\bf P}^{\dagger}_{U}\right]{\bf G}^{a,V}(E){\bf
Q}_{U}(E)\right.$ $\displaystyle+$
$\displaystyle\left.e^{ieV_{\alpha}t}\left[Z^{*}(\epsilon){\bf
G}^{a,0}(\epsilon)-Z^{*}(\epsilon_{\beta\alpha}){\bf
G}^{a,V}(\epsilon_{\beta})\right]{\bf\Sigma}^{a,0}_{\alpha}(E)\right\\}$

where

$\displaystyle{\bf P}_{D}$ $\displaystyle=$ $\displaystyle
Z(\epsilon_{\beta}){\bf
U}+\sum_{\delta}\left[Z(\epsilon_{\beta})-Z(\epsilon_{\beta\delta})\right]\left[{\bf\Sigma}^{r,0}_{\delta}(\epsilon_{\beta\delta})-{\bf\Sigma}^{r,0}_{\delta}(E)\right]$
$\displaystyle{\bf P}_{U}$ $\displaystyle=$ $\displaystyle-Z(\epsilon){\bf
U}+\sum_{\delta}\left[Z(\epsilon)-Z(\epsilon_{\delta})\right]\left[{\bf\Sigma}^{r,0}_{\delta}(\epsilon)-{\bf\Sigma}^{r,0}_{\delta}(E-V_{\delta})\right]$
(37) $\displaystyle{\bf
Q}_{D}(E)=\int\frac{d\epsilon^{\prime}}{2\pi}~{}\left[1-e^{i(\epsilon^{\prime}-E)t}\right]Z(\epsilon^{\prime}){\bf\Sigma}^{a,0}_{\alpha}(\epsilon^{\prime})$
$\displaystyle{\bf
Q}_{U}(E)=\int\frac{d\epsilon^{\prime}}{2\pi}~{}\left[1-e^{i(\epsilon^{\prime}_{\alpha}-E)t}\right]Z(\epsilon^{\prime}_{\alpha}){\bf\Sigma}^{a,0}_{\alpha}(\epsilon^{\prime})$

with

$\displaystyle Z(\epsilon)=[i(E-\epsilon-i0^{+})]^{-1}$ (38)

In the absence of the ac bias, the quantity $A_{\alpha}$ is the Fourier
transform of the retarded Green’s function while the quantity
$F_{\beta\alpha}$ is related to the Fourier transform of the advanced Green’s
function. They are all expressed in terms of the unperturbed Green’s functions
${\bf G}^{r/a,0/V}$ and self energy ${\bf\Sigma}^{0/V}$ which have been widely
studied in molecular device using the NEGF-DFT formalism. ${\bf G}^{r/a,0/V}$
and self energy ${\bf\Sigma}^{0/V}$ can be expressed as

$\displaystyle{\bf G}^{r/a,0/V}(\epsilon)=\left[\epsilon{\bf I}-{\bf
H}^{0/V}_{c}-{\bf\Sigma}^{r/a,0/V}(\epsilon)\right]^{-1}$
$\displaystyle\mathbf{\Sigma}^{\gamma,0}_{\alpha}(\epsilon)=\left[\mathbf{T}^{0}_{c\alpha}-\epsilon{\bf
S}^{0}_{c\alpha}\right]\mathbf{g}^{\gamma}_{\alpha\alpha}(\epsilon)\left[\mathbf{T}^{0}_{\alpha
c}-\epsilon{\bf S}^{0}_{\alpha c}\right]$
$\displaystyle\mathbf{\Sigma}^{\gamma,V}_{\alpha}(\epsilon)=\left[\mathbf{T}^{0}_{c\alpha}-\epsilon^{\alpha}{\bf
S}^{0}_{c\alpha}\right]\mathbf{g}^{\gamma}_{\alpha\alpha}(\epsilon^{\alpha})\left[\mathbf{T}^{0}_{\alpha
c}-\epsilon^{\alpha}{\bf S}^{0}_{\alpha c}\right]$

where $\gamma=r,a,<$, $\epsilon^{\alpha}=\epsilon-eV_{\alpha}$. Obviously,
${\bf\Sigma}^{\gamma,V}_{\alpha}(\epsilon)={\bf\Sigma}^{\gamma,0}_{\alpha}(\epsilon-
eV_{\alpha})$. In the wideband limit, Eq.(33-LABEL:FU1) will reduce to the
formula first derived by Jauho et al.Jauho2  With ${\bf A}$ and ${\bf F}$
obtained we can, in principle, solve the AC current biased by downwards or
upwards pulse exactly. In practice, however, its computational cost is
expensive for a realistic molecular device. For example, to calculate
$J_{\alpha}^{out}(t)$, we have to do triple integrals over energy and repeat
this procedure to collect data for all time sequence. In the numerical
calculation especially in ab-initio modeling, it is practically very difficult
if not impossible to calculate the transient current for the complex structure
in molecular devices. So approximation must be made so that Eq.(33-LABEL:FU1)
can be simplified.

###  III.2 Approximate scheme of ${\bf A}_{\beta}(t,\epsilon)$ and ${\bf
F}_{\beta\alpha}(t,\epsilon)$

The approximate solution of ${\bf A}_{\beta}(t,\epsilon)$ and ${\bf
F}_{\beta\alpha}(t,\epsilon)$ in Eq.(33-LABEL:FU1) have to satisfy the
following requirements. First, it has to greatly reduce the calculational
cost. Second, it has to keep essential physics of transient dynamics. Third,
it must have the correct initial current at $t=0$ and approach the correct
asymptotic limit at $t\rightarrow\infty$. The first goal is realized by
eliminating double energy integral using a reasonable ansatz, with which the
dynamical properties of molecular device is maintained.

To find such an ansatz, we first assume that ${\bf\Sigma}^{a,0}(E)$ changes
smoothly and slightly with $E$ and is analytic in the upper half plane, so
that the typical integral like $\int d\epsilon
dE~{}\frac{e^{i(\epsilon-E)t}}{-i(E-\epsilon+i0^{+})}{\bf\Sigma}^{a,0}(E)$ is
roughly zero due to the different phase in $e^{\i(\epsilon-E)t}$. Then the
last term of ${\bf F}^{U/D}$ and the second term of ${\bf Q}^{U/D}$ disappear.
Considering the following identity,

$\displaystyle\int\frac{dE}{2\pi}~{}\frac{{\bf\Sigma}^{a}_{\alpha}(E)}{-i(E-\epsilon+i0^{+})}$
$\displaystyle=$
$\displaystyle\left[\int_{-\infty}^{0^{-}}+\frac{1}{2}\int^{0^{+}}_{0^{-}}\right]d\tau~{}{\bf\Sigma}^{a}_{\alpha}(\tau)\int\frac{dE}{2\pi}~{}\frac{e^{iE\tau}}{-i(E-\epsilon+i0^{+})}$
$\displaystyle=$
$\displaystyle\left[\int_{-\infty}^{0^{+}}-\frac{1}{2}\int^{0^{+}}_{0^{-}}\right]d\tau~{}e^{i\epsilon\tau}{\bf\Sigma}^{a}_{\alpha}(\tau)={\bf\Sigma}^{a}_{\alpha}(\epsilon)-{\bf\Delta}^{a}_{\alpha}$

and defining
${\bf\Sigma}^{a}_{\alpha}(E,\Delta)={\bf\Sigma}^{a}_{\alpha}(E)-{\bf\Delta}^{a}_{\alpha}$,
the first term of ${\bf F}_{U/D}$ and ${\bf Q}_{U/D}$ in Eqs.(37) can be
simplified, ${\bf F}_{U/D}$ now becomes

$\displaystyle{\bf F}^{D}_{\beta\alpha}$ $\displaystyle\simeq$
$\displaystyle{\bf
G}^{a,0}(\epsilon){\bf\Sigma}^{a,0}_{\alpha}(\epsilon,\Delta)+\int\frac{dE}{2\pi}~{}e^{-i(\epsilon-E)t}$
(39) $\displaystyle\times$
$\displaystyle\left[Z^{*}(\epsilon_{\beta})-Z^{*}(\epsilon)+{\bf
G}^{a,V}(\epsilon_{\beta}){\bf P}^{\dagger}_{D}\right]{\bf
G}^{a,0}(E){\bf\Sigma}^{a,0}_{\alpha}(E,\Delta)$ $\displaystyle{\bf
F}^{U}_{\beta\alpha}$ $\displaystyle\simeq$ $\displaystyle{\bf
G}^{a,V}(\epsilon_{\beta}){\bf\Sigma}^{a,0}_{\alpha}(\epsilon_{\beta\alpha},\Delta)+\int\frac{dE}{2\pi}~{}e^{-i(\epsilon_{\beta}-E)t}$
$\displaystyle\times$
$\displaystyle\left[Z^{*}(\epsilon)-Z^{*}(\epsilon_{\beta})+{\bf
G}^{a,0}(\epsilon){\bf P}^{\dagger}_{U}\right]{\bf
G}^{a,V}(E){\bf\Sigma}^{a,0}_{\alpha}(E-eV_{\alpha},\Delta)$

We note that, in the wide-band limit, Eq.(39,LABEL:newFU) is exact. With our
approximation we have eliminated one of the energy integrals in $J^{out}$, and
${\bf A}$ and ${\bf F}$ now have similar structures since $\tilde{\bf
F}\sim{\bf A}^{\dagger}{\bf\Sigma}^{a}$.

With the approximation defined in Eq.(39,LABEL:newFU), the current can be
written in a compact form (see section C) if we introduce the effective
Green’s function

$\displaystyle\tilde{\bf G}^{r/a,0}(E,\epsilon)=\left[E{\bf S}-{\bf
H}^{0}_{c}-\sum_{\alpha}{\bf\Sigma}^{r/a,0}_{\alpha}(\epsilon)\right]^{-1}$
(41) $\displaystyle\tilde{\bf G}^{r/a,V}(E,\epsilon)=\left[E{\bf S}-{\bf
H}^{V}_{c}-\sum_{\alpha}{\bf\Sigma}^{r/a,V}_{\alpha}(\epsilon)\right]^{-1}$
(42)

In general we have to consider the overlap matrix ${\bf S}$. However, we
should keep in mind that in the deriving of the time dependent current, we
have to orthogonalize the basis set, which would lead to ${\bf S}={\bf I}$.
Here $\tilde{\bf G}^{r/a}(E,\epsilon)$ can be regarded as the Green’s
functions at energy $E$ and constant parameter $\epsilon$ for open system with
the effective Hamiltonian ${\bf H}^{r/a}_{eff}={\bf
H}_{c}+{\bf\Sigma}^{r}_{\alpha}(\epsilon)$. For a given ${\bf H}_{eff}$,
Eqs.(41,42) are equivalent to

$\displaystyle(E{\bf S}-{\bf H}^{r}_{eff})\tilde{\bf G}^{r}={\bf I}$ (43)

On the other hand, Green’s function can be expanded in terms of the
eigenfunctions of the corresponding Hamiltonian,

$\tilde{\bf G}^{r}=\sum_{n}{\bf\Psi}^{n}C_{n}.$ (44)

where ${\bf H}_{eff}{\bf\Psi}^{n}=E_{n}(\epsilon){\bf\Psi}_{n}$. Substituting
Eq.(44) into Eq.(43), and using the general orthogonality relation
${\bf\Phi}^{n,\dagger}{\bf S}{\bf\Psi}^{m}=C_{m}\delta_{nm}$ [see Appendix A]
and the eigenvalue equation ${\bf
H}_{eff}{\bf\Psi}^{n}=E_{n}(\epsilon){\bf\Psi}^{n}$, we have

$\displaystyle\tilde{\bf
G}^{r}(E,\epsilon)=\sum_{n}\frac{{\bf\Psi}^{n}{\bf\Phi}^{n,\dagger}}{[E-E_{n}(\epsilon)]{\bf\Phi}^{n,\dagger}{\bf
S}{\bf\Psi}^{n}}$ (45)

Obviously, this Green’s function can be calculated by finding the residues
${\rm Res}_{n}={\bf\Psi}^{n}{\bf\Phi}^{n,\dagger}/{{\bf\Phi}^{n,\dagger}{\bf
S}{\bf\Psi}^{n}}$ at various poles $E=E_{n}(\epsilon)$.

Then, we replace $Z(\epsilon){\bf G}^{r/a}(E)$ in Eqs.(33,35,39,LABEL:newFU)
by $Z(\epsilon)\tilde{\bf G}^{r/a}(E,\epsilon)$. Although $\tilde{\bf
G}^{r/a}(E,\epsilon)$ is different from initial Green’s function ${\bf
G}^{r/a}(E)=\left[E-{\bf H}_{c}-{\bf\Sigma}^{r/a}(E)\right]^{-1}$, this
substitution is reasonable since the major contribution of the integration in
Eqs.(33-LABEL:FU1) comes from the pole $\epsilon$ in $Z(\epsilon)$ (see
Eq.(38)). Similarly, considering the major contribution of the pole of
$Z(\epsilon)$, we replace $Z(\epsilon){\bf\Sigma}^{a,0}(E)$ in
Eqs.(33,35,39,LABEL:newFU) by $Z(\epsilon){\bf\Sigma}^{a,0}(\epsilon)$. Since
${\bf\Sigma}(\epsilon)$ in $\tilde{\bf G}^{r}(E,\epsilon)$ is independent of
energy $E$, we can perform contour integration over energy $E$ in Eqs.(33) and
(35) by closing a contour on lower half plane and perform the integration over
energy $E$ in Eqs.(34) and (LABEL:FU1) by closing a contour on upper half
plane. Thus, energy integration over $E$ can be analytically performed. It
should be noted that the self energy ${\bf\Sigma}^{r/a}$ is not independent of
energy in contrast to the wide-band limit, this energy dependence is on
$\epsilon$ but not on $E$. In this way, we can reduce the computational cost
and keep the essential physics of the dynamics as we will show later.

###  III.3 Approximate expression of ${\bf A}_{\beta}(t,\epsilon)$ and ${\bf
F}_{\beta\alpha}(t,\epsilon)$

Now, considering the initial current and the asymptotic long time limit, we
can write the approximate expression of ${\bf A}_{\beta}(t,\epsilon)$ and
${\bf F}_{\beta\alpha}(t,\epsilon)$ from Eqs.(33,35,39,LABEL:newFU):

$\displaystyle{\bf A}^{D/U}_{\beta}(t,\epsilon)={\bf A}^{D/U}_{\beta,1}+{\bf
A}^{D/U}_{\beta,2}$ (46) $\displaystyle{\bf
F}^{D}_{\beta\alpha}(t,\epsilon)={\bf
A}^{D,\dagger}_{\beta,1}{\bf\Sigma}^{a,0}_{\alpha}(\epsilon_{\beta\alpha},\Delta)+{\bf
A}^{D,\dagger}_{\beta,2}{\bf\Sigma}^{a,0}_{\alpha}(\epsilon,\Delta)$ (47)
$\displaystyle{\bf F}^{U}_{\beta\alpha}(t,\epsilon)={\bf
A}^{U,\dagger}_{\beta,1}{\bf\Sigma}^{a,0}_{\alpha}(\epsilon,\Delta)+{\bf
A}^{U,\dagger}_{\beta,2}{\bf\Sigma}^{a,0}_{\alpha}(\epsilon_{\beta\alpha},\Delta)$
(48)

with

$\displaystyle{\bf
A}^{D}_{\beta,1}=\int\frac{dE}{2\pi}~{}e^{i(\epsilon-E)t}\left[Z(\epsilon_{\beta})\tilde{\bf
G}^{r,0}(E,\epsilon_{\beta})\left({\bf I}+{\bf\Xi}^{D}{\bf
G}^{r,V}(\epsilon_{\beta})\right)\right]$ (49) $\displaystyle{\bf
A}^{D}_{\beta,2}={\bf
G}^{r,0}(\epsilon)-\int\frac{dE}{2\pi}~{}e^{i(\epsilon-E)t}\left[Z(\epsilon)\tilde{\bf
G}^{r,0}(E,\epsilon)\right]$ (50) $\displaystyle{\bf
A}^{U}_{\beta,1}=\int\frac{dE}{2\pi}~{}e^{i(\epsilon_{\beta}-E)t}\left[Z(\epsilon)\tilde{\bf
G}^{r,V}(E,\epsilon)\left({\bf I}+{\bf\Xi}^{U}{\bf
G}^{r,0}(\epsilon)\right)\right]$ (51) $\displaystyle{\bf
A}^{U}_{\beta,2}={\bf
G}^{r,V}(\epsilon_{\beta})-\int\frac{dE}{2\pi}~{}e^{i(\epsilon_{\beta}-E)t}\left[Z(\epsilon_{\beta})\tilde{\bf
G}^{r,V}(E,\epsilon_{\beta})\right]$ (52)

where

$\displaystyle{\bf\Xi}^{D}$ $\displaystyle=$ $\displaystyle{\bf
U}+\sum_{\delta}\left[{\bf\Sigma}^{r,0}_{\delta}(\epsilon_{\beta\delta})-{\bf\Sigma}^{r,0}_{\delta}(\epsilon_{\beta})\right]$
$\displaystyle=$ $\displaystyle{\bf
U}+\sum_{\delta}\left[{\bf\Sigma}^{r,V}_{\delta}(\epsilon_{\beta})-{\bf\Sigma}^{r,0}_{\delta}(\epsilon_{\beta})\right]$
$\displaystyle{\bf\Xi}^{U}$ $\displaystyle=$ $\displaystyle-{\bf
U}+\sum_{\delta}\left[{\bf\Sigma}^{r,0}_{\delta}(\epsilon)-{\bf\Sigma}^{r,0}_{\delta}(\epsilon-
eV_{\delta})\right]$ (53) $\displaystyle=$ $\displaystyle-{\bf
U}+\sum_{\delta}\left[{\bf\Sigma}^{r,0}_{\delta}(\epsilon)-{\bf\Sigma}^{r,V}_{\delta}(\epsilon)\right]$

This is the second level of approximation. As we will see later that it is
better than the first level approximation described below. Now we can make
further approximation (the first level). To do this, we note that the Green’s
function ${\bf G}^{r}$ can be obtained using the Dyson equation,

${\bf G}^{r,tot}={\bf G}^{r,ex}+{\bf G}^{r,ex}{\bf\Xi}{\bf G}^{r,tot}$ (54)

where ${\bf G}^{r,tot}$ is the Green’s function of system denoted by ${\bf
H}^{tot}={\bf H}^{ex}+{\bf H}^{\prime}$, ${\bf G}^{r,ex}$ is the unperturbed
Green’s function corresponding to ${\bf H}^{ex}$ that can be exactly solved,
${\bf\Xi}$ is the effective self energy describing ${\bf H}^{\prime}$. If we
set ${\bf H}^{ex}$ and ${\bf H}^{tot}$ as zero biased open system and
$V_{\alpha}$ biased open system respectively, we have

${\bf G}^{r,tot}={\bf G}^{r,V}(\epsilon)={\bf G}^{r,0}(\epsilon)+{\bf
G}^{r,0}(\epsilon){\bf\Xi}^{D}{\bf G}^{r,V}(\epsilon)$ (55)

Similarly, if we treat ${\bf H}^{ex}$ and ${\bf H}^{tot}$ as $V_{\alpha}$
biased open system and zero biased open system, respectively, we obtain
another Dyson equation

${\bf G}^{r,tot}={\bf G}^{r,0}(\epsilon)={\bf G}^{r,V}(\epsilon)+{\bf
G}^{r,V}(\epsilon){\bf\Xi}^{U}{\bf G}^{r,0}(\epsilon)$ (56)

Similar to the derivation of the second level of approximation, we can also
replace ${\bf G}^{ex}(\epsilon)$ by $\tilde{\bf G}^{ex}(E,\epsilon)$ in
Eq.(55,56) which leads to

$\displaystyle\tilde{\bf G}^{r,V}(E,\epsilon)\simeq\tilde{\bf
G}^{r,0}(E,\epsilon)\left[{\bf I}+{\bf\Xi}^{D}{\bf G}^{r,V}(\epsilon)\right]$
$\displaystyle\tilde{\bf G}^{r,0}(E,\epsilon)\simeq\tilde{\bf
G}^{r,V}(E,\epsilon)\left[{\bf I}+{\bf\Xi}^{U}{\bf G}^{r,0}(\epsilon)\right]$
(57)

Then, Eqs.(49) and (51) can be further approximated as

$\displaystyle{\bf
A}^{D}_{\beta,1}=\int\frac{dE}{2\pi}~{}e^{i(\epsilon-E)t}\left[Z(\epsilon_{\beta})\tilde{\bf
G}^{r,V}(E,\epsilon_{\beta})\right]$ (58) $\displaystyle{\bf
A}^{U}_{\beta,1}=\int\frac{dE}{2\pi}~{}e^{i(\epsilon_{\beta}-E)t}\left[Z(\epsilon)\tilde{\bf
G}^{r,0}(E,\epsilon)\right]$ (59)

This is the first level of approximation. It is easy to confirm that when the
self-energy is energy independent these two approximations lead to exactly the
same expression of transient current in the wide-band limit. In the next
section we will numerically compare these two approximations with the exact
solution.

###  III.4 initial and asymptotic currents

We now show that the currents calculated from Eqs.(27,32,46-52) and from
Eqs.(27,32,46-48,50,52,58,59) satisfy the correct current limit at initial
$t=0$ and asymptotic limit $t\rightarrow\infty$ times. Note that the initial
current and asymptotic currents can be calculated from a standard DC transport
nonequilibrium Green’s function analysis. It is expected that the asymptotic
current for the downward pulse $J^{D}_{\alpha}(t\rightarrow\infty)$ and
initial current for the upward pulse $J^{U}_{\alpha}(t=0)$ are zero since
there is no bias in the system. Now we discuss the limiting cases for two
versions of approximations developed in section IIIC.

When $t=0$, $e^{i(\epsilon-E)t}=1$, we can perform integration over energy $E$
in Eqs.(49-52) by closing a contour at upper half plane, where only a single
residual exists at an energy pole of $Z$. At $t=0$, $\tilde{\bf
G}^{r/a,0/V}(E,\epsilon)={\bf G}^{r/a,0/V}(\epsilon)$, therefore Eqs.(49,51)
and Eqs.(58,59) are equivalent. Now we focus on the current obtained from
Eqs.(27,32,46-48,50,52,58,59). After integrating over $\epsilon$, the two
terms in Eqs.(50,52) cancels to each other, then from Eq.(58, 59), ${\bf
A}^{D/U}_{\beta}(t=0)$ becomes

$\displaystyle{\bf A}^{D}_{\beta}(t=0)$ $\displaystyle=$
$\displaystyle\tilde{\bf G}^{r,V}(\epsilon_{\beta},\epsilon_{\beta})={\bf
G}^{r,V}(\epsilon_{\beta})$ (60) $\displaystyle{\bf A}^{U}_{\beta}(t=0)$
$\displaystyle=$ $\displaystyle\tilde{\bf G}^{r,0}(\epsilon,\epsilon)={\bf
G}^{r,0}(\epsilon)$ (61)

For ${\bf F}_{\beta\alpha}$, we can perform integration over energy $E$ by
closing a contour at lower half plane. Similarly, there also exists only a
single residual on energy pole $E_{Z}$ of $Z^{*}$ in the lower half plane, and

$\displaystyle{\bf F}^{D}_{\beta\alpha}(t=0)$ $\displaystyle=$
$\displaystyle\tilde{\bf
G}^{a,V}(\epsilon_{\beta},\epsilon_{\beta}){\bf\Sigma}^{a,0}_{\alpha}(\epsilon_{\beta\alpha},\Delta)={\bf
G}^{a,V}(\epsilon_{\beta}){\bf\Sigma}^{a,V}_{\alpha}(\epsilon_{\beta},\Delta)$
$\displaystyle{\bf F}^{U}_{\beta\alpha}(t=0)$ $\displaystyle=$
$\displaystyle\tilde{\bf
G}^{a,V}(\epsilon,\epsilon){\bf\Sigma}^{a,0}_{\alpha}(\epsilon,\Delta)={\bf
G}^{a,0}(\epsilon){\bf\Sigma}^{a,0}_{\alpha}(\epsilon,\Delta)$ (63)

Substituting Eq.(60-63) into Eq.(27,32), and considering

$\displaystyle{\bf\Sigma}^{\gamma,0}_{\beta}(\epsilon)={\bf\Sigma}^{\gamma,V}_{\beta}(\epsilon_{\beta})$
$\displaystyle{\bf G}^{<,0/V}(\epsilon)={\bf
G}^{r,0/V}(\epsilon)\left[\sum_{\beta}{\bf\Sigma}^{<,0/V}_{\beta}(\epsilon)\right]{\bf
G}^{a,0/V}(\epsilon)$
$\displaystyle{\bf\Sigma}^{<,0}_{\beta}(\epsilon)=f(\epsilon)\left[{\bf\Sigma}^{a,0}_{\beta}(\epsilon)-{\bf\Sigma}^{r,0}_{\beta}(\epsilon)\right]$
$\displaystyle{\bf\Sigma}^{<,V}_{\beta}(\epsilon)=f(\epsilon-
eV_{\beta})\left[{\bf\Sigma}^{a,V}_{\beta}(\epsilon)-{\bf\Sigma}^{r,V}_{\beta}(\epsilon)\right]$
(64)

where $f(\epsilon)$ is Fermi distribution function, we have initial current at
$t=0$

$\displaystyle J^{D}_{\alpha}=2e\rm{Re}\int\frac{d\epsilon}{2\pi}~{}{\bf
G}^{r,V}(\epsilon){\bf\Sigma}_{\alpha}^{<,V}(\epsilon)+{\bf
G}^{<,V}(\epsilon){\bf\Sigma}^{a,V}_{\alpha}(\epsilon)$ (65) $\displaystyle
J^{U}_{\alpha}=2e\rm{Re}\int\frac{d\epsilon}{2\pi}~{}{\bf
G}^{r,0}(\epsilon){\bf\Sigma}_{\alpha}^{<,0}(\epsilon)+{\bf
G}^{<,0}(\epsilon){\bf\Sigma}^{a,0}_{\alpha}(\epsilon)$ (66)

Eqs.(65) and (66) are the same as the formal DC current expression in the case
of nonzero bias and zero bias, respectively. $J^{U}_{\alpha}(t=0)$ in Eq.(66)
is exactly zero since the Fermi distribution in ${\bf\Sigma}^{<}_{\alpha}$ and
${\bf G}^{<}$ are equal for $\alpha=L$ and $\alpha=R$.

When $t\rightarrow\infty$, by virtue of the Riemann-Lebesgue lemma,lemma  the
Fourier integral over $\epsilon$ vanishes, i.e.,
$\int\frac{d\epsilon}{2\pi}~{}e^{-i\epsilon t}{\bf G}^{r}{\bf\Sigma}^{r}...$
equal to zero at $t\rightarrow\infty$ since there always exist poles in lower
half plane. With this in mind, we have

$\displaystyle{\bf A}^{D}_{\beta}(t\rightarrow\infty,\epsilon)$
$\displaystyle=$ $\displaystyle{\bf G}^{r,0}(\epsilon)$ (67)
$\displaystyle{\bf F}^{D}_{\beta\alpha}(t\rightarrow\infty,\epsilon)$
$\displaystyle=$ $\displaystyle{\bf
G}^{a,0}(\epsilon){\bf\Sigma}^{a,0}_{\alpha}(\epsilon,\Delta)$ (68)
$\displaystyle{\bf A}^{U}_{\beta}(t\rightarrow\infty,\epsilon)$
$\displaystyle=$ $\displaystyle{\bf G}^{r,V}(\epsilon_{\beta})$ (69)
$\displaystyle{\bf F}^{U}_{\beta\alpha}(t\rightarrow\infty,\epsilon)$
$\displaystyle=$ $\displaystyle{\bf
G}^{a,V}(\epsilon_{\beta}){\bf\Sigma}^{a,0}_{\alpha}(\epsilon_{\beta\alpha},\Delta)={\bf
G}^{a,V}(\epsilon_{\beta}){\bf\Sigma}^{a,V}_{\alpha}(\epsilon_{\beta},\Delta)$

From Eq.(67-LABEL:FU4) and Eq.(27,32), we have the asymptotic current

$\displaystyle J^{D}_{\alpha}=2e\rm{Re}\int\frac{d\epsilon}{2\pi}~{}{\bf
G}^{r,0}(\epsilon){\bf\Sigma}_{\alpha}^{<,0}(\epsilon)+{\bf
G}^{<,0}(\epsilon){\bf\Sigma}^{a,0}_{\alpha}(\epsilon)$ (71) $\displaystyle
J^{U}_{\alpha}=2e\rm{Re}\int\frac{d\epsilon}{2\pi}~{}{\bf
G}^{r,V}(\epsilon){\bf\Sigma}_{\alpha}^{<,V}(\epsilon)+{\bf
G}^{<,V}(\epsilon){\bf\Sigma}^{a,V}_{\alpha}(\epsilon)$ (72)

It is easy to see, Eqs.(71) and (72) are the formal DC current expression in
the case of zero bias and nonzero bias, respectively, and
$J^{D}_{\alpha}(t\rightarrow\infty)$ in Eq.(71) is exactly zero.

##  IV comparison with the exact result in quantum dot system

Now we consider a system composed of a single-level quantum dot connected to
external leads with a Lorentzian linewidth. This system can be solved exactly
to give a transient current for pulse-like biasMaciejko . We can obtain
transient current using three methods: (i) the exact current expressed by
Eqs.(27,32, 33-LABEL:FU1), (ii) the first level of approximation from
Eqs.(27,32,46-48,50,52,58,59) and (iii) the second level of approximation from
Eqs.(27,32,46-52). We will compare the current obtained from these three
methods. The system is described by the following simple Hamiltonian

$\displaystyle
H=\sum_{k_{\alpha}}\epsilon_{k_{\alpha}}(t)c^{\dagger}_{k_{\alpha}}c_{k_{\alpha}}+\epsilon_{d}(t)d^{\dagger}d+\sum_{k_{\alpha}}(t_{k_{\alpha}}c^{\dagger}_{k_{\alpha}}d+h.c.)$
(73)

where $\epsilon_{d}(t)=\epsilon_{d}^{0}+U(t)$ and
$\epsilon_{k_{\alpha}}(t)=\epsilon^{0}_{k_{\alpha}}+V_{\alpha}(t)$. Because
the scattering region has only one state with energy level $\epsilon^{0}_{d}$,
the Green’s functions $G(\epsilon)$ and self energy $\Sigma(\epsilon)$ thus
become scalars instead of matrices. If we choose linewidth function
$\Gamma_{\alpha}(\omega)\equiv 2\pi\rho_{\alpha}(\omega)|t_{k_{\alpha}}|^{2}$
to be Lorentzian with the linewidth amplitude $\Gamma^{0}_{\alpha}$,

$\Gamma_{\alpha}(\omega)=\frac{W^{2}}{\omega^{2}+W^{2}}\Gamma^{0}_{\alpha}$

then $G^{\gamma}(\epsilon)$ and $\Sigma^{\gamma}(\epsilon)$ can be expressed
as

$\displaystyle
G^{r/a,0}(\epsilon)=\left[\epsilon-\epsilon_{d}^{0}-\sum_{\alpha}\Sigma^{r/a,0}(\epsilon)\right]^{-1}$
$\displaystyle
G^{r/a,V}(\epsilon)=\left[\epsilon-\epsilon_{d}^{0}-U^{V}-\sum_{\alpha}\Sigma^{r/a,V}(\epsilon)\right]^{-1}$
$\displaystyle
G^{<,0/V}(\epsilon)=G^{<,0/V}(\epsilon)\left[\sum_{\alpha}\Sigma^{<,0/V}(\epsilon)\right]G^{<,0/V}(\epsilon)$
$\displaystyle\Sigma^{r/a,0}_{\alpha}(\epsilon)=\int\frac{d\omega}{2\pi}~{}\Gamma_{\alpha}(\omega)/(\epsilon-\omega\pm
i0^{+})$
$\displaystyle\Sigma^{r/a,V}_{\alpha}(\epsilon)=\int\frac{d\omega}{2\pi}~{}\Gamma_{\alpha}(\omega)/(\epsilon-
eV_{\alpha}-\omega\pm i0^{+})$
$\displaystyle\Sigma^{<,0}_{\alpha}(\epsilon)=f(\epsilon)\left[\Sigma^{a,0}_{\alpha}(\epsilon)-\Sigma^{r,0}_{\alpha}(\epsilon)\right]$
$\displaystyle\Sigma^{<,V}_{\alpha}(\epsilon)=f(\epsilon-
eV_{\alpha})\left[\Sigma^{a,V}_{\alpha}(\epsilon)-\Sigma^{r,V}_{\alpha}(\epsilon)\right]$

Using the theorem of residual, we can analytically perform integral in
$A_{\beta}$ and $F_{\beta\alpha}$ for either exact formula or two approximate
formulas. In the calculation, we set $\Gamma=\Gamma^{0}_{L}+\Gamma^{0}_{R}$ as
the energy unit, and set $\Gamma^{0}_{L}=\Gamma^{0}_{R}=0.5$.

We first consider the transient current induced by opposite voltage
$V_{L}(t)=-V_{R}(t)$. In this case, the equilibrium coulomb potential in
quantum dot $U^{0,V}=0$, and the time dependent perturbation coming from
coulomb response $U(t)$ is assumed to be zero. It is a reasonable assumption
since the coulomb potential in scattering region is canceled by the opposite
voltage in left and right lead. In Fig.1, we plot two approximated transient
currents and exact transient current in downward [panel (a), (b), (c)] and
upward [panel (d), (e), (f)] case vs time for different bandwidth $W$. We find
that for all bandwidth $W$, the approximated current and exact current have
the same dynamical behaviors. Fig.2 gives direct comparison where we merge
panels (a), (b) and (c) in Fig.1 as panel (a) in Fig.2, and merge panels (d),
(e) and (f) in Fig.1 as panel (b) in Fig.2. We can see that for the downward
pulse [panel (a)], transient current using three formulas are almost
indistinguishable. This means that in the opposite voltage, our approximation,
the first approximation [Eqs.(50,52,58,59)] and the second approximation
[Eqs.(49-52)] are all very good for studying transient dynamics. For the
upward pulse, although the approximations are not as good as in downward case,
the currents calculated from approximate scheme are still in good agreement
with the exact solution especially for the second approximation. Hence we may
conclude that the two approximations are all reasonable in the opposite
voltage $V_{L}(t)=-V_{R}(t)$. They can be used to study transient dynamics in
the real molecular device to speed up the calculation.

Figure 1:  (Color online) Time dependent current $J(t)$ corresponding to an
opposite downward pulse or upward pulse in three versions: the exact solution
and two approximations. The different black lines are for different bandwidth
$W$. The red line is for $W=\infty$, i.e., the wide-band limit. The current is
in the unit of $e\Gamma$, and the time is in the unit of $2\pi/\Gamma$.
$eV_{L}=-eV_{R}=5.$ Figure 2: (Color online) Merged version of Fig.1 for
$W=1$, $2$, $5$ and $20$. Panel (a) corresponding to the downward pulse
current comes from panel (a), (b) and (c) in Fig.1, panel (b) corresponding to
upward pulse current comes from panel (d), (e) and (f). Along the black arrow,
the bandwidth are $W=1$, $2$, $5$ and $20$, respectively.

Next, we focus on the asymmetric voltage, i.e., $V_{L}(t)\neq V_{R}(t)$. In
this case, the equilibrium coulomb potential in quantum dot $U^{0/V}$, and the
time dependent perturbation coming from coulomb response $U(t)$ can’t be
canceled by the voltage in left and right lead. In principle, perturbation
$U(t)$ should be calculated by solving time dependent Schrödinger equation, it
will be very difficult and computational demanding therefore can’t be
implemented in real molecular device. As an alternative scheme, we have set
$U(t)=[eV_{L}(t)\Gamma^{0}_{L}+eV_{L}(t)\Gamma^{0}_{L}]/\Gamma$. For the
single level quantum dot system, this is exact because the central scattering
region now is expressed in a scalar instead of matrices, which leads to the
same transient current for the opposite voltage $V_{L}(t)=-V_{R}(t)$ and
asymmetric voltage $V_{L}(t)=V(t)$, $V_{R}(t)=0$ or $V_{L}(t)=0$,
$V_{R}(t)=-V(t)$ in the exact solution.

For the first approximation the poles in time dependent term
$e^{i(\epsilon-E)t}$ are different from that in the second level
approximation, i.e., the poles of $\tilde{\bf G}^{r,0}$ in Eq.49 and
$\tilde{\bf G}^{r,V}$ in Eq.51 are replaced by the poles of $\tilde{\bf
G}^{r,V}$ in Eq.58 and $\tilde{\bf G}^{r,0}$ in Eq.59, respectively. Because
of this, the time evolution process are not as accurate in the first
approximation, especially for the large $V_{\alpha}$. So, for the asymmetric
voltage, the second approximation is better. In Fig.3 and Fig.4, we compare
the transient current obtained from the second approximation [panel (b-d)] for
opposite or asymmetric voltage with the exact transient current [panel (a)] in
response to the downward pulse and upward pulse, respectively. We find that
all transient currents from the second approximation in Fig.3 and Fig.4 [panel
(b)] are very close to the exact result [panel (a)]. Moreover, in Fig.3 and
Fig.4, the approximate transient current in panel (b), (c), (d) have almost
the same behavior. It is safe to say that our approximations have kept
essential physics of dynamical transport properties.

Figure 3: (Color online) Panel (a): exact time dependent current $J(t)$
corresponding to downward pulse for $dV=V_{L}-V_{R}=5$. Panel (b-d) are
corresponding to the second approximate transient current corresponding to
downward pulse for opposite voltage $V_{L}=-V_{R}=2.5$, asymmetric voltage
$V_{L}=5$, $V_{R}=0$ and $V_{L}=0$, $V_{R}=-5$, respectively. The different
black lines are for different bandwidths $W$. The red line is wide-band limit
for $W=\infty$. Figure 4: (Color online) Same to Fig.3, transient current
corresponding to upward pulse vs time are plotted.

##  V several examples for real molecular devices

In this section, we implement our approximate formula in two representative
molecular devices including a short carbon chain coupled to aluminum leads and
a $C_{60}$ molecule coupled to aluminum leads. These systems were chosen
because they are typical in first-principles calculation and their practical
importance to nano-electronics. In Fig.5(a) and Fig.5(b), we show the
structure of Al-${\rm C}_{5}$-Al and Al-${\rm C}_{60}$-Al, respectively, where
Al leads are along (100) direction, one unit cell of Al lead consists of 9 Al
atoms and total 40 atoms were included in the simulation box. For the Al-C5-Al
device, the nearest distance between Al leads and the carbon chain is 3.781
a.u. and the distance of C-C bond is 2.5 a.u.(1 a.u.=0.529Å). In Al-C60-Al
device, the distance between the Al atom and the nearest C atom equal to 3.625
a.u..

To calculate the dynamic response of molecular devices, we have used the
first-principles quantum transport package MATDCAL.matdcal  Considering the
complicated coulomb response in scattering region, we set
$V_{L}(t)=-V_{R}(t)$. In this case, the first approximation is simple but as
good as the second one. So, in the following, the first approximate formula
[Eqs.(27,32,46-48,50,52,58,59)] is used. In principle, the calculation
involves the following steps: (1) calculate the device Hamiltonian including
central scattering Hamiltonian and lead Hamiltonian using NEGF-DFT package to
get two potential landscapes $U^{0}$ at zero bias and $U^{V}$ at $V_{\alpha}$
bias, respectively. They are originally expressed in a nonorthogonal fireball
basis. (2) orthogonalize the nonorthogonal device Hamiltonian using the
approachthesis  introduced in Appendix B so that they are finally expressed in
an orthogonal basis. (3) with the orthogonal lead Hamiltonian $H_{\alpha}$,
one calculates zero biased self energy ${\bf\Sigma}^{r/a,0}_{\alpha}$ and
$V_{\alpha}$ biased self energy ${\bf\Sigma}^{r/a,V}_{\alpha}$ from
Eqs.(17,18) using the transfer matrix method.transfer  (4) with orthogonalized
central scattering Hamiltonian ${\bf H}^{0}_{c}$ and ${\bf H}^{V}_{c}$ and
self energy ${\bf\Sigma}^{r/a,0}_{\alpha}$ and ${\bf\Sigma}^{r/a,V}_{\alpha}$
obtained from two potential landscapes $U^{0}$ and $U^{V}$, one solves the
effective Green’s function ${\bf G}^{r/a,0/V}$ using Eqs.(41,42) by
calculating its poles and residuals from Eq.(45). Step (1)-(4) are time
independent processes and easy to perform. (5) calculate time dependent
quantities ${\bf A}^{D/U}_{\beta,1}$ and ${\bf A}^{D/U}_{\beta,2}$ from
Eqs.(58,59) and Eqs.(50,52). Then ${\bf A}_{\beta}$ and ${\bf
F}_{\beta\alpha}$ can be calculated from Eqs.(46-48). (6) integrate over
$\epsilon$ and obtain the final AC current
$J^{D/U}(t)=[J^{D/U}_{L}(t)-J^{D/U}_{R}(t)]/2$ from Eqs.(27,32).

Figure 5: (Color online) Panel (a): Structure of Al-C5-Al. Panel(b): structure
of Al-${\rm C}_{60}$-Al.

First we study the Al-${\rm C}_{5}$-Al structure. In Fig.6, we plot the
transient current $J(t)$ corresponding to the upward pulse [panel (a) and (b)]
and the downward pulse [panel (c) and (d)] for different external voltages
$V_{R}=-V_{L}=0.001a.u.$ [panel (a) and (c)] and $V_{R}=-V_{L}=0.01a.u.$
[panel (b) and (d)] in Al-${\rm C}_{5}$-Al structure. Following observations
are in order: (1) as we have discussed in Sec.III.4, for all bias voltages
$V_{\alpha}$ the transient currents indeed reach the correct limits at $t=0$
and $t\rightarrow\infty$. For the upward pulse, $J(t=0)=0$ and
$J(t\rightarrow\infty)=J_{dc}$ while for the downward pulse we have
$J(t=0)=J_{dc}$ and $J(t\rightarrow\infty)=0$. (2) for both upward pulse
(turn-on voltage) and downward pulse (turn-off voltage), once the bias voltage
is switched, currents oscillate rapidly in the first a few or tens fs and then
gradually approach to the steady-state values, i.e., $J_{dc}$ for turn-on
voltage and zero for turn-off voltage. The larger the voltage $V_{\alpha}$,
the more rapid the current oscillates. (3) concerning the long time behavior,
the time dependent current oscillates with a frequency proportional to
$|V_{\alpha}|$.WangBin  This is because the time dependent term
$e^{i(\epsilon-E)t}$ in Eqs.(50,52,58,59) are $V_{\alpha}$ dependent. For the
upward pulse, $e^{i(\epsilon_{\alpha}-E)t}\propto e^{iV_{\alpha}t}$, which
directly leads to the oscillating frequency proportional to $|V_{\alpha}|$.
For the downward pulse, although $e^{i(\epsilon-E)t}$ is $V_{\alpha}$
independent, in the energy integral on $E$, the pole $E_{n}$ of $\tilde{\bf
G}^{r}(E,\epsilon)$ are determined by the self energy
${\bf\Sigma}^{r,V}_{\alpha}$. Since ${\bf\Sigma}^{r,V}_{\alpha}$ depends on
$V_{\alpha}$, this leads to $V_{\alpha}$ dependent oscillating frequency. In
addition, we notice that although the properties of dc conductance of short
carbon chains are different for the chains with odd and even number
atomsWangBin2  due to the completely different electronic structure near Fermi
level, the ac signals are similar (see Ref.WangBin,  where Al-${\rm C}_{4}$-Al
structure was analyzed). This indicates that in AC transport, all states with
energy from $-\infty$ to the Fermi energy are contributing, which is very
different from dc case where only the states near Fermi level contribute to
transport processes.

Figure 6: (Color online) Time dependent current $J(t)$ corresponding to the
upward pulse [panel (a) and (b)] and the downward pulse [panel (c) and (d)]
for different external voltages $V_{\alpha}$ in Al-${\rm C}_{5}$-Al device.
The inset of panel (a) shows the long time behavior of the time-dependent
current. The red (gray in print) dashed lines in panels indicate asymptotic
current $J(t\rightarrow\infty)$ which the DC current biased by $V_{L/R}$
labeled in corresponding panels for the upward pulse, and arrive at zero for
the downward case.

Next, we study the second sample: the Al-${\rm C}_{60}$-Al structure. In
Fig.7, the transient current $J(t)$ of the structure corresponding to an
upward pulse [panel (a) and (b)] and a downward pulse [panel (c) and (d)] for
different external voltages $V_{R}=-V_{L}=0.001a.u.$ [panel (a) and (c)] and
$V_{R}=-V_{L}=0.01a.u.$ [panel (b) and (d)] are plotted. Similar to the
Al-${\rm C}_{5}$-Al structure, correct initial current $J(t=0)$ and asymptotic
current $J(t\rightarrow\infty)$ are also obtained in Al-${\rm C}_{60}$-Al
structure. In addition, there are also rapidly oscillations at short times
after the switch although the oscillation is not as rapid as that in the
Al-${\rm C}_{5}$-Al structure. Furthermore, similar to Al-${\rm C}_{5}$-Al
structure, in gradually reaching the steady-state values, the current
oscillates with a frequency proportional to $|V_{\alpha}|$ but its decay rate
is much slower than that in Al-${\rm C}_{5}$-Al structure. It indicates that
there are much more quasi-resonant state that contribute to the transient
current in Al-${\rm C}_{60}$-Al structure which is reasonable considering the
complex electronic structure of isolated ${\rm C}_{60}$. In the following, we
will analyze in detail how the current decays for the Al-${\rm C}_{60}$-Al
structure.

Figure 7: (Color online) Time dependent current $J(t)$ corresponding to the
upward pulse [panel (a) and (b)] and the downward pulse [panel (c) and (d)] in
Al-${\rm C}_{60}$-Al device for different $V_{\alpha}$. In panel (a) and (c),
$V_{R}=-V_{L}=0.001a.u.$. In panel (b) and (d), $V_{R}=-V_{L}=0.01a.u.$. Same
to Fig.6, the red (gray in print) dashed lines in panels indicate asymptotic
current $J(t\rightarrow\infty)$. The long time AC current or detailed short
time AC current are shown in inset of panels.

Physically, decay time of current corresponds to the width of the quasi-bound
state. In molecular devices, because the linewidth function
${\bf\Gamma}(\epsilon)$ are complex and energy dependent matrix, we can’t
extract characteristic time scale directly from $1/{\bf\Gamma}$. As such, the
transmission coefficient $T(\epsilon)$ is needed to understand the resonant
state and corresponding characteristic time scale. In Fig.8(a), we plot
transmission coefficient $T(\epsilon)$ in the energy range from the energy
band bottom to the Fermi energy for Al-${\rm C}_{60}$-Al structure at zero
bias. Here, the sharp peaks [some of them, see red crossed signed peaks in
Fig.8(a)] correspond to resonant states with large lifetimes. At a particular
resonant state, the incoming electron can dwell for a long time, which
contributes to a much more slowly decaying current than other non-resonant
states. In Fig.8(b), (c) and (d), we amplify the first, second and forth
labeled quasi-resonant transmission, respectively, where the peaks’ width
$\Gamma_{peak}\sim 10^{-5}a.u.$ are indicated, corresponding to a decay time
$\tau\sim 2400fs$ from the expression $\Gamma_{peak}t=1$. In Fig.8(e)-(g),
corresponding to different $\epsilon$ where the resonant peaks in Fig.8(b)-(d)
are located, we plot long time behavior of current element $J_{L}(\epsilon)$.
Here $J_{L}(\epsilon)$ is the time dependent current for each energy
$\epsilon$, the integration over which gives the final current
$J_{\alpha}(t)$. We can see that for each resonant state the current
$J_{L}(\epsilon)$ keeps oscillating in a long time comparable to the decay
time $\tau\sim 2400fs$. Furthermore the intensity of the oscillation $\Delta
J\sim 0.2\mu A$ is not very small comparing to the DC signal $J_{dc}=5.1\mu
A$.

Figure 8: (Color online) Panel (a): transmission coefficient $T(\epsilon)$ in
the energy range from the energy band bottom to the Fermi energy. In the whole
energy range, there are some resonant states corresponding to the very sharp
transmission coefficient $T(\epsilon)$, as we have indicated (see red cross)
and labeled (by 1, 2, 3 and 4) in panel (a), some of them contribute to the
current at long time. We amplify the first, second and forth labeled resonant
transmission in panel (b), (c) and (d), respectively. In panel (e)-(g), we
plot the long time behavior of current $J_{L}(\epsilon)$ at a fixed $\epsilon$
for the first, second and forth resonant states. The external voltage
$V_{\alpha}=0.001a.u.$.

After integration over energy, these slowly decaying currents
$J_{L}(\epsilon)$ due to the resonant states may cancel to each other
partially due to the difference in their phases. However, we should keep in
mind that it is these resonant peaks that may give rise to convergence
problem. Hence in the calculation, we should first scan the equilibrium and
non-equilibrium transmission coefficient (100,000 energy points for example)
to resolve sharp resonant peaks in the whole energy range from minimum energy
to Fermi energy. Then, for each sharp resonant peak, enough (100 for example)
energy points should be chosen to converge the integration of the current
$J_{L}(\epsilon)$ over $\epsilon$, i.e., $\int d\epsilon J(\epsilon)$. For the
non-resonant state, i.e., the smoothly changed region in $T(\epsilon)$, the
current $J(\epsilon)$ are integrated using less energy points.

As we have discussed that the resonant states are important for the transient
current and they must be carefully treated in calculation. However, in the
calculation of the effective Green’s function $\tilde{\bf G}^{r/a,0/V}$, a
small imaginary part that is usually added to the real energy
$\epsilon\rightarrow\epsilon+i\eta$ to help resolving the retarded or advanced
self-energies. This in turn introduces pseudo resonant states. In order to
eliminate the pseudo resonant state in effective Green’s function $\tilde{\bf
G}^{r/a,0/V}$ [Eqs.(41,42)], one has to calculate the self-energy by setting
$\eta=0$ and resolve the retarded or advanced self-energies with the aid of
the group velocity $v_{k}=(\partial E(k)/\partial k)$.Sanvito

##  VI conclusion

By orthogonalizing the Hamiltonian expressed in the nonorthogonal basis and
considering the singularity of self-energy ${\Sigma}^{r/a}(t,t^{\prime})$ at
$t=t^{\prime}$, we have generalized the solution [ developed in Ref.Maciejko,
] of the transient current driven either by a downward step voltage pulse or
by a upward step pulse. This generalized result can be applied to both the
quantum dot model and real molecular device. Based on the exact solution given
in Ref.Maciejko, , we derived two approximate formulas that are suitable for
numerical calculation of the transient current for molecular devices. We have
tested our approximate formula in a quantum dot system where exact numerical
solution exists. For the quantum dot system, we chose a Lorentzian linewidth
(beyond wideband limit) and compared the time-dependent current calculated
using both exact formula and our approximate formula. We found that for the
opposite voltage $V_{L}(t)=-V_{R}(t)$, the results obtained from the exact
formalism and two approximate scheme agree very well with each other
especially in the downward pulse case. For the nonsymmetric voltage
$V_{L}(t)=V(t)$, $V_{R}(t)=0$ or $V_{L}(t)=0$, $V_{R}(t)=-V(t)$, the second
approximation is better. This shows that our approximate formulas captured the
essential physics of the transient current. In addition, it gives the correct
initial current at $t=0$ and correct asymptotic current at
$t\rightarrow\infty$. Since we have reduced the calculation from triple
integral to single integral over the energy, the approximated approach reduces
the computational cost drastically and it can be easily implemented in first
principles calculation for molecular devices. To demonstrate this, we
calculated the transient current using the first approximated scheme with an
opposite voltage $V_{L}(t)=-V_{R}(t)$ for two molecular structures: Al-${\rm
C}_{5}$-Al and Al-${\rm C}_{60}$-Al. Different from the quantum dot system,
because of the complex electronic structure in molecular devices, transient
currents oscillate rapidly in the first a few or tens fs as the bias voltage
is switched, then gradually approach to the steady-state values. Furthermore,
due to the resonant state in molecular devices, transient currents have a very
long decay time $\tau$.

##  Appendix A orthogonality relation for the nonorthogonal basis

For a system described by $H$, the time independent eigenvalue equation is
written as:

$\displaystyle H|n\rangle$ $\displaystyle=$ $\displaystyle E_{n}|n\rangle$
(74)

the eigenvectors $|n\rangle$ form an orthogonal complete basis set. However,
in many systems such as a molecular device connected to external leads, the
basis set constructed by eigen vectors is not convenient. We usually expand
the eigen vector $|n\rangle$ in other basis $|\mu\rangle$, which is non-
orthogonal complete set (or nearly complete).

$\displaystyle|n\rangle\simeq\sum_{\mu}|\mu\rangle\langle\mu|n\rangle$ (75)

the eigenvalue equation now becomes

$\displaystyle\sum_{\mu}H|\mu\rangle\langle\mu|n\rangle$ $\displaystyle=$
$\displaystyle E_{n}\sum_{\mu}|\mu\rangle\langle\mu|n\rangle$
$\displaystyle\sum_{\mu}\langle\nu|H|\mu\rangle\langle\mu|n\rangle$
$\displaystyle=$ $\displaystyle
E_{n}\sum_{\mu}\langle\nu|\mu\rangle\langle\mu|n\rangle$
$\displaystyle\sum_{\mu}{\bf H}_{\nu\mu}{\bf\Psi}_{\mu}^{n}$ $\displaystyle=$
$\displaystyle E_{n}\sum_{\mu}{\bf S}_{\nu\mu}{\bf\Psi}_{\mu}^{n}$ (76)

where ${\bf S}_{\nu\mu}=\langle\nu|\mu\rangle$. In the matrix form, we have
${\bf H}{\bf\Psi}^{n}=E_{n}{\bf S}{\bf\Psi}^{n}$. If we use the self-energy to
replace the effect of leads the effective Hamiltonian for the open system
becomes ${\bf H}={\bf H}_{0}+{\bf\Sigma}^{r}$. Since the effective Hamiltonian
is not Hermitian, we can define the adjoint operator ${\bf H}^{\dagger}={\bf
H}={\bf H}_{0}+{\bf\Sigma}^{a}$ and corresponding eigen-equation becomes ${\bf
H}^{\dagger}|\phi_{n}\rangle=E^{*}_{n}{\bf S}|\phi_{n}\rangle$. Then

$\displaystyle{\bf\Phi}^{m,\dagger}{\bf H}{\bf\Psi}^{n}$ $\displaystyle=$
$\displaystyle E_{n}{\bf\Phi}^{m,\dagger}{\bf S}{\bf\Psi}^{n},$ (77)
$\displaystyle{\bf\Psi}^{n,\dagger}{\bf H}^{\dagger}{\bf\Phi}^{m}$
$\displaystyle=$ $\displaystyle E^{*}_{m}{\bf\Psi}^{n,\dagger}{\bf
S}^{\dagger}{\bf\Phi}^{m}$ (78)

Taking hermitian conjugate of Eq.(78),

$\displaystyle{\bf\Phi}^{m,\dagger}{\bf H}{\bf\Psi}^{n}$ $\displaystyle=$
$\displaystyle E_{m}{\bf\Phi}^{m,\dagger}{\bf S}{\bf\Psi}^{n}$ (79)

From (77) and (79), we have

$\displaystyle{\bf\Phi}^{n,\dagger}{\bf S}{\bf\Psi}^{m}=C_{m}\delta_{nm}$ (80)

For the normalized wave function $|\psi_{n}\rangle$ and $|\phi_{n}\rangle$,

$\displaystyle{\bf\Phi}^{\dagger}{\bf S\Psi}={\bf I}$ (81)

It is the usual orthogonality relation for eigenvectors expressed in a
nonorthogonal basis set. For an hermitian Hamiltonian ${\bf H}={\bf
H}^{\dagger}$, $|\psi_{n}\rangle=|\phi_{n}\rangle$, we have

${\bf\Psi}^{\dagger}{\bf S\Psi}={\bf I}.$

##  Appendix B Orthogonalize Hamiltonian expressed in nonorthogonal basis

In this appendix, we will show how to construct a new orthogonal basis from
the atomic real-space nonorthogonal basis. We will transform the original
Hamiltonian ${\bf H}$ which is expressed in the nonorthogonal basis into
Hamiltonian $\tilde{\bf H}$ expressed in the new orthogonal basis. Of course,
instead of ${\bf S}$, the overlap matrix in the new basis will be ${\bf I}$.

Denoting nonorthogonal basis $|\mu\rangle$ and orthogonal basis $|j\rangle$,
they are related by unitary transform operator ${\bf\mathcal{U}}$

$\displaystyle|{\mu}\rangle=\sum_{j}|j\rangle\langle
j|{\mu}\rangle=\sum_{j}|j\rangle{\bf\mathcal{U}}_{j\mu}$
$\displaystyle{\bf\mathcal{U}}_{j\mu}=\langle j|{\mu}\rangle$ (82)

where we have used the completeness of orthogonal basis $|j\rangle$. Using the
orthogonality $\langle i|j\rangle=\delta_{ij}$

$\displaystyle\sum_{\mu\nu}\langle
i|\mu\rangle\langle\mu|\nu\rangle\langle\nu|j\rangle$ $\displaystyle=$
$\displaystyle\sum_{\mu\nu}{\bf{\mathcal{U}}}_{i\mu}{\bf
S}_{\mu\nu}{\bf{\mathcal{U}}}^{\dagger}_{\nu j}=\delta_{ij}$

where we have used the completeness of nonorthogonal basis. In the matrix
form, ${\bf{\mathcal{U}}}{\bf S}{\bf{\mathcal{U}}}^{\dagger}={\bf I}$. We can
formally define

${\bf{\mathcal{U}}}={\bf
S}^{-\frac{1}{2}},~{}~{}{\bf{\mathcal{U}}}^{\dagger}=\left[{\bf
S}^{-\frac{1}{2}}\right]^{\dagger}.$

Then new Hamiltonian ${\tilde{\bf H}}$ expressed in basis $|i\rangle$ can be
expressed as:

$\displaystyle{\tilde{\bf H}}_{ij}$ $\displaystyle=$ $\displaystyle\langle
i|H|j\rangle$ (83) $\displaystyle=$ $\displaystyle\sum_{\mu\nu}\langle
i|\mu\rangle\langle\mu|H|\nu\rangle\langle\nu|j\rangle$ $\displaystyle=$
$\displaystyle\sum_{\mu\nu}{\bf{\mathcal{U}}}_{i\mu}{\bf
H}_{\mu\nu}{\bf{\mathcal{U}}}^{\dagger}_{\nu j}$

In the matrix form, ${\tilde{\bf H}}={\bf S}^{-\frac{1}{2}}{\bf H}\left[{\bf
S}^{-\frac{1}{2}}\right]^{\dagger}$.

We now discuss how to find the matrix ${\bf S}^{-\frac{1}{2}}$. Without loss
generality, we assume the real overlap matrix ${\bf S}$ satisfies eigen
function ${\bf SV}={\bf V}{\rm diag}(\lambda_{1},...,\lambda_{n})$ with the
eigenvalues $\lambda_{1},...,\lambda_{n}$ and eigenvectors ${\bf
V}=[v_{1},...,v_{n}]$. Since ${\bf S}$ is real and symmetric, the eigenvectors
are real and orthogonal, and it thus holds that ${\bf V}^{\dagger}{\bf V}={\bf
V}{\bf V}^{\dagger}={\bf I}$. Then

$\displaystyle{\bf S}$ $\displaystyle=$ $\displaystyle{\bf V}{\rm
diag}(\lambda_{1},...,\lambda_{n}){\bf V}^{\dagger}$ $\displaystyle=$
$\displaystyle{\bf V}{\rm diag}(\sqrt{\lambda_{1}},...,\sqrt{\lambda_{n}}){\bf
V}^{\dagger}{\bf V}{\rm diag}(\sqrt{\lambda_{1}},...,\sqrt{\lambda_{n}}){\bf
V}^{\dagger}$

It follows that

$\displaystyle{\bf S}^{\frac{1}{2}}={\bf V}{\rm
diag}(\sqrt{\lambda_{1}},...,\sqrt{\lambda_{n}}){\bf V}^{\dagger}$ (84)

From ${\bf S}^{-\frac{1}{2}}{\bf S}^{\frac{1}{2}}={\bf I}$ and Eq.(84), we
have

$\displaystyle{\bf S}^{-\frac{1}{2}}{\bf V}{\rm
diag}(\sqrt{\lambda_{1}},...,\sqrt{\lambda_{n}}){\bf V}^{\dagger}={\bf I}$
$\displaystyle{\bf S}^{-\frac{1}{2}}{\bf V}{\rm
diag}(\sqrt{\lambda_{1}},...,\sqrt{\lambda_{n}}){\bf V}^{\dagger}{\bf V}{\rm
diag}(\frac{1}{\sqrt{\lambda_{1}}},...,\frac{1}{\sqrt{\lambda_{n}}}){\bf
V}^{\dagger}$ $\displaystyle=$ $\displaystyle{\bf S}^{-\frac{1}{2}}={\bf
V}{\rm
diag}(\frac{1}{\sqrt{\lambda_{1}}},...,\frac{1}{\sqrt{\lambda_{n}}}){\bf
V}^{\dagger}$ (85)

In general, the dimension of matrix ${\bf S}$ is infinity, we can’t calculate
its eigenvalue $\lambda_{i}$ and eigenvector $v_{i}$ by diagonalizing ${\bf
S}$. However, in the tight-binding representation, the state $\mu$ and $\nu$
hardly overlap when their separation is large enough in real space, i.e.,
${\bf S}_{\mu\nu}\approx 0$ for most of off-diagonal elements. Considering the
periodic properties in semi-infinite leads, we can select a block matrix which
is large enough to include all the overlap between leads and central molecular
regions. For the non-orthogonal basis including several unit cell of atomic
leads as a buffer layer into the central scattering region is enough to get a
good screening for dc transport calculation. In transforming the Hamiltonian
to the orthogonal basis needed for ac transport calculation, however, it turns
out that we have to include at least 10 unit cells of atomic leads into the
central scattering region. Partly because the overlap of orthogonal basis has
longer range than that of non-orthogonal basis. With this large simulation box
(finite dimension), we can calculate the overlap matrix ${\bf
S}^{\frac{1}{2}}$ therefore transform ${\bf H}$ into $\tilde{\bf H}$. The
accuracy of transformed Hamiltonian $\tilde{\bf H}$ should be examined by
comparing dc conductance obtained from the original Hamiltonian ${\bf H}$ and
the transformed Hamiltonian $\tilde{\bf H}$.

${\bf ACKNOWLEDGMENTS}$

This work was supported by a RGC grant (HKU 705409P) from the government of
HKSAR.

## References

  * (1) Electronic address<EMAIL_ADDRESS>  * (2) Molecular Electronics, editd by G. Cuniberti, G. Fagas, and K. Ricchter (Springer, Berlin, 2005).
  * (3) S. Kurth, G. Stefanucci, C.-O. Almbladh, A. Rubio, and E. K. U. Gross, Phys. Rev. B 72, 035308 (2005);
  * (4) Y. Zhu, J. Maciejko, T. Ji, H. Guo, and J. Wang, Phys. Rev. B 71, 075317 (2005).
  * (5) X. Qian, J. Li, X. Lin, and S. Yip , Phys. Rev. B 73, 035408 (2006).
  * (6) J. Maciejko, J. Wang, and H. Guo, Phys. Rev. B 74, 085324 (2006).
  * (7) N. Sai, N. Bushong, R. Hatcher, and M. Di Ventra, Phys. Rev. B 75, 115410 (2007).
  * (8) X. Zheng, F. Wang, C.Y. Yam, Y. Mo, and G.H. Chen, Phys. Rev. B 75, 195127 (2007).
  * (9) G. Stefanucci, S. Kurth, A. Rubio, and E. K. U. Gross, Phys. Rev. B 77, 075339 (2008).
  * (10) J. Taylor, H. Guo and J. Wang, Phys. Rev. B 63 245407 (2001); ibid, 63 121104 (2001).
  * (11) A. D. Stone, M. Ya. Azbel, and P. A. Lee, Phys. Rev. B, 31, 1707 (1985);
  * (12) P. Johansson, Phys. Rev. B, 41, 9892 (1990); T. Kwapiński, R. Taranko, and E. Taranko, Phys. Rev. B 66, 035315 (2002); T. Kwapiński, Phys. Rev. B 69, 153303 (2004).
  * (13) Q.-F. Sun and T.-H. Lin, J. Phys.: Condens. Matter 9, 3043 (1997); idib. 9, 4875 (1997); Q.-F. Sun, B.-G. Wang, J. Wang, and T.-H. Lin, Phys. Rev. B 61, 4754 (2000); J. Q. You, C.-H. Lam, and H. Z. Zheng, Phys. Rev. B 62, 1978 (2000); Z.-G. Zhu, G. Su, Q.-R. Zheng, and B. Jin, Phys. Rev. B 68, 224413 (2003).
  * (14) S. Kohler, J. Lehmann, P. Hnggi, Phys. Rep. 406, 379 (2005); L. E. F. Foa Torres, G. Cuniberti, C. R. Physique 10, 297 (2009).
  * (15) M. Plihal D. C. Langreth and P. Nordlander, Phys. Rev. B 61, R13341 (2000);
  * (16) L. Y. Chen and C. S. Ting, Phys. Rev. B 41, 8533 (1990).
  * (17) W. R. Frensley, Phys. Rev. B 36, 1570 (1987).
  * (18) F. B. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 (2005).
  * (19) G. Stefanucci and C.-O. Almbladh, Europhys. Lett. 67, 14 (2004); G. Stefanucci and C.-O. Almbladh, Phys. Rev. B 69, 195318 (2004); N. Bushong, N. Sai, and M. Di Ventra, Nano Lett. 5, 2569 (2005).
  * (20) A. Schiller and S. Hershfield, Phys. Rev. B 62, R16271 (2000); Y.X. Xing, Q.F. Sun, and J. Wang, Phys. Rev. B 75, 125308 (2007).
  * (21) By solving the time dependent density matrix using the equation of motion method, Zheng et alchen1  have calculated the transient current using the wideband limit.
  * (22) B. Wang, Y. Xing, L. Zhang, and J. Wang, Phys. Rev. B 81, 121103 (2010).
  * (23) N. S. Wingreen, A.-P. Jauho and Y. Meir, Phys. Rev. B 48, 8487 (1993).
  * (24) A short account of this method can be found in Ref.WangBin, .
  * (25) A.-P. Jauho, N. S. Wingreen and Y. Meir, Phys. Rev. B 50, 5528 (1994).
  * (26) M. H. Pedersen and M. Büttiker, Phys. Rev. B, 58, 12993-13006 (1998); B. G. Wang, J. Wang, and H. Guo, Phys. Rev. Lett., 82, 398-402 (1999).
  * (27) M. Buttiker, J. Phys. Condens. Matter 5, 9361 (1993).
  * (28) Y. D. Wei and J. Wang, Phys. Rev. B. 79, 195315 (2009).
  * (29) The exact solution of transient current on the orthogonal basis set has been derived in Ref.Maciejko, .
  * (30) J. K. Viljas, J. C. Cuevas, F. Pauly, and M. Häfner, Phys. Rev. B, 72, 245415 (2005).
  * (31) D. H. Lee and J. D. Joannopoulos, Phys. Rev. B 23, 4997 (1981); ibid, 23, 4988 (1981).
  * (32) M. Buttiker, A. Pretre, and H. Thomas, Phys. Rev. Lett. 70, 4114 (1993).
  * (33) B. G. Wang, J. Wang, and H. Guo, Phys. Rev. Lett. 82, 398 (1999).
  * (34) A. Erdélyi, Asymptotic Expansions (Dover, New York, 1956).
  * (35) M. Brandbyge, J. L. Mozos, P. Ordejon, J. Taylor, and K. Stokbro, Phys. Rev. B 65, 165401 (2002).
  * (36) Master thesis by T. Markussen, Quantum transport calculations using wave function propagation and the Kubo formula.
  * (37) B. Wang, Y. Xing, Y. Wei, L. Wan, and J. Wang, Carbon, 43, 2786 (2005); B. Wang, Y. Xing, L. Wan, Y. Wei, and J. Wang, Phys. Rev. B 71, 233406 (2005).
  * (38) S. Sanvito, C. J. Lambert, J. H. Jefferson, and A. M. Bratkovsky, Phys. Rev. B 59, 11936 (1999).