Calculation of the current response in a nanojunction for an arbitrary time-dependent bias: application to the molecular wire

Recently [Phys. Rev. B 91, 125433 (2015)] we derived a general formula for the time-dependent quantum electron current through a molecular junction subject to an arbitrary time-dependent bias within the Wide Band Limit Approximation (WBLA) and assuming a single particle Hamiltonian. Here we present an efficient numerical scheme for calculating the current and particle number. Using the Pad\'e expansion of the Fermi function, it is shown that all frequency integrals occurring in the general formula for the current can be removed analytically. Furthermore, when the bias in the reservoirs is assumed to be sinusoidal it is possible to manipulate the general formula into a form containing only summations over special functions. To illustrate the method, we consider electron transport through a one-dimensional molecular wire coupled to two leads subject to out-of-phase biases. We also investigate finite size effects in the current response and particle number that results from the switch-on of such a bias.


Introduction
The Landauer-Büttiker (LB) formalism has for some time provided a reliable method for modelling coherent transport through molecular junctions [1,2]. In the original formulation, this method enabled computation of the current in the leads coupled to a molecular region C, when a bias in at least one of the leads causes a constant shift in the Fermi energies of the leads. Work done with the LB formalism has mainly focused on the calculation of steady-state properties, corresponding to a time regime that consigns the bias switch-on time to the distant past [3]. The LB expressions for the current and quantum noise may be derived from the Nonequilibrium Green's Function (NEGF) formalism, which is a more general alternative to the S-Matrix approach [4,5]. This requires the solution of the Kadanoff-Baym equations for the lesser Green's function of the molecular region [4], from which the electric current in the leads and electron number in the molecule may be extracted. When the bias is time-independent, all GFs are functions of the time difference, and hence it is possible to work purely in the frequency domain when solving for the components that enter the expression for the current.
However, recent studies have opened up the possibility of calculating the time-dependent (TD) current in a multilevel nanojunction using a closed integral formula; expressions that include the transient effects of the switch-on (and hence the system preparation) have been derived for the cases of a constant [4,6,7] and timedependent [8] bias. Common to both approaches is the assumption of the Wide-Band Limit Approximation (WBLA), which neglects the detailed electronic structure of the leads and renders the Kadanoff-Baym equations for the lesser GF analytically tractable. Computationally, the main cost of the method we presented in [8] lies with the evaluation of several integrals in the (ω, t) plane. The purpose of this paper is to present an efficient computational scheme for evaluating these integrals with no loss of accuracy making consideration of complex multi-level and multi-terminal systems within the WBLA numerically possible.
The structure of this paper is as follows: in Section 2, we introduce the model of the generic switchon problem and the basic assumptions involved in the WBLA. In addition we prove a continuity equation for the current formula derived in [8]. In Section 3, the formula for the current is re-expressed using the Padé expansion of the Fermi function to remove all frequency integrals, forming the basis for our numerical implementation. Section 4 presents numerical results for the current through a molecular wire, modelled with a simple tight-binding Hamiltonian. In this implementation, expressions are worked out for a sinusoidal bias whereby Bessel function expansions are used to remove all t−integrals in the TD current. This enables simulations of the transport resulting from any combination of sinusoidal biases in the leads, including those which contain a symmetry-breaking phase. Furthermore it enables us to identify some novel finite size effects on the dynamics and to study the conservation of charge as the transient current dies out.

Hamiltonian and GF Components
The NEGF formalism is a method for calculating the evolution in time of ensemble averages, following an event at a particular time t 0 which breaks the physical symmetry of the system between times t = ±∞. In quantum transport one is primarily concerned with calculations of the time-dependent current following the switch-on of a bias in one of the leads. Here we model this switch-on by adding a spatially constant time-dependent shift to the lead energy levels.
We will in general be concerned with the following Hamiltonian: Here,d kα ,d m andd † kα ,d † m are annihilation and creation operators of leads and central system electronic states, where for simplicity spin degrees of freedom are neglected. We collect elements of this Hamiltonian into a matrix consisting of 'blocks' corresponding to each of the physical subsystems it describes. The first term is a Hamiltonian of the lead states k belonging to each lead α. The second term is the Hamiltonian of the molecule coupled to the lead; it contains inter-molecular hopping matrix elements V mn between the m and n molecular orbitals, defining a molecular Hamiltonian matrix 'block' h CC with elements V mn . The third term describes the coupling of the molecule to the leads, and defines the α − C 'block' of the Hamiltonian matrix, h αC with elements V kα,m . The switch-on problem is specified by assigning values to the matrix elements of H (t) at times before and after t 0 .
We assume that, prior to t 0 , the system has equilibrated with lead-molecule couplings present in the Hamil-tonianĤ (t < t 0 ) ≡Ĥ 0 , so thatĤ 0 is as in Eq. (1) with time-independent lead state energies ε kα (t < t 0 ) = ε kα . The system is therefore described by an equilibrium density operator ρ 0 = Z −1 e −β( H0−µ N ) , where Z is the partition function,N the number operator, β = 1/k B T the inverse temperature and µ the chemical potential. At times following the switch-on, the lead energies acquire a time-dependent shift, ε kα (t > t 0 ) = ε kα + V α (t), where the function V α (t) is arbitrary. This is an example of a partition-free approach [9,6,8] to the switch-on problem, as the system Hamiltonian is not divided into decoupled subregions prior to t 0 . Partitioned approaches [5] add the lead-molecule coupling at the switch-on time and therefore cannot be expected to give a physically meaningful transient.
The problem of calculating the current in response to this kind of switch-on is mapped to the solution of a set of coupled integro-differential equations for Green's Functions (GFs) defined on a complex time contour γ, which we refer to as the Konstantinov-Perel' contour [10]. One then solves for various components of the two-time GF of the central region CC 'block', whose rigorous definition may be found elsewhere [4,8]. Here we simply note that the retarded GF satisfies the following equation of motion: On the right-hand side of this expression, one encounters the retarded component of the embedding selfenergy, We have in this expression introduced the level width matrix which forms a Hilbert transform pair with Λ α,mn (ω). In addition we have defined the phase factor: In the WBLA, one assumes that Γ α,mn (ω) is replaced with a frequency-independent value Γ α = Γ α,mn ε F α , where ε F α is the equilibrium Fermi energy of lead α [5]. This has the effect of mapping Eq. (2) onto an expression which is analytically tractable, resulting in the following exact expression: In (6), Γ α is an effective Hamiltonian of the central region, whose eigenvalues correspond to unstable eigenmodes of the molecular structure, which have acquired a finite lifetime due to the presence of the leads. All other components of the GF can also be explicitly calculated in the time domain [8]. In particular, the greater and lesser GFs can be expressed as where we make the definition and note that all information concerning the time-dependent bias V α (t) of lead α is to be found in the phase factor ψ α (t, t 0 ).

Current and Continuity Equation
In [8], results were presented for the transport through a quantum dot using direct numerical integration to evaluate the current at each time t > t 0 . The current in lead α is defined as the time derivative of the average particle number in that lead, , and may be expressed in the following rather compact form (taking electron charge q = −1 and a spin degeneracy of 2): The molecule which has been sandwiched between macroscopic leads is fundamentally an open quantum system, so the total charge in the molecular region is not necessarily conserved even though the total amount of charge in the molecular region plus that of all the leads is fixed for all times. This is a condition automatically satisfied by the NEGF formalism [4]. The condition for local charge conservation is that the sum of all currents measured in the leads of a multiterminal device is equal to zero: This condition states that the rate of charge entering the molecular region is balanced by the rate of charge leaving it. One should now recall that there is no direct coupling between leads, so that all charge transfer between leads is mediated by the central region, and also that the definition of the current in lead α was given as I α (t) ≡ q dNα(t) dt , i.e. it was given as the rate of increase of charge in this lead. Eq. (10) is simply a statement of the condition that if charge increases in one lead connected to the molecule, then this increase is compensated exactly by a net decrease of charge from the other leads making up the junction. The total amount of charge in the molecular region at any given time is directly related to the lesser GF: The time derivative of the charge number dN C (t) dt in the central region has been referred to as the displacement current in the literature [4]. To make precise the relation between this quantity and the currents measured in the leads, we take the time-derivatives of the matrix S β ≡ S β (t, t 0 ; ω) and its complex conjugate, e.g.
This expression, along with its complex conjugate, is substituted into the derivative of N C in Eq. (11) to give If we then sum over lead indices in the expression (9), the following theorem is established: This is a statement of global charge conservation in the nanojunction. When the bias in every lead is constant, (11) reduces to the formula for the particle number given in [6]. In the long time limit the number of charges on the molecule is constant in time, i.e. an ideal stationary state is reached; hence, trivially the condition for local charge conservation (10) is satisfied in this case. Although the condition of global charge conservation (13) is never violated, physical situations exist in which the condition (10) for local charge conservation in the molecular (central) region of a nanojunction is violated, i.e. dN C (t)/dt = 0.
3 Pade Expansion of the Time-Dependent Current

Exact Formula for the Current: Hurwitz-Lerch Functions
In previous work [8], numerical integration was used to compute the current formula (9) at different times for a single level system. We now seek to extend this work to a numerical method that facilitates time-dependent transport calculations in systems of real experimental interest, for instance carbon nanotube transistors [11].
To proceed with a numerical implementation of Eqs. (9) and (11), one can introduce the right and left eigenproblems for the renormalized Hamiltonian matrix h ef f CC [12]: The eigenenergiesε j contain an imaginary part that is strictly negative (as Γ is positive-definite), and the same value ofε j corresponds to each of the left and right eigenvectors. Using the idempotency property the formula (9) can be put into the form: is the Fermi function and c.c j↔k denotes the complex conjugate of the preceding term with indices j and k exchanged. We therefore reduce I α (t) to a sum of integrals over scalar functions in the (x, t) plane. One integral whose structure is repeated throughout (16) can be performed analytically, using a contour integral over a semi-circular contour in the upper-half of the complex plane: Here we have introduced the so-called Hurwitz-Lerch Transcendent Φ [13]: The formula (17) enables us to replace several frequency integrals with a fast-converging series expansion. We note in passing that, if one follows the steps in the formal integration of the lesser GF G < CC (t 1 , t 2 ) which enabled the derivation of Eq. (9) in Ref. [8], one uses a well-known [4] transformation of Matsubara summations into frequency integrals when evaluating a set of convolution integrals taken along the vertical branch of the Konstantinov-Perel' contour. However, these steps can be omitted altogether, and a direct route can be taken to an expression for I α (t) in terms of Φ which retains the Matsubara summation at all stages.

Removal of Frequency Integrals: Padé Expansion
The Padé approximation to the Fermi function is a series expansion whose terms possess a simple pole structure [14,15,16]: This is identical in structure to the Matsubara expansion, which has parameter values η l = 1 and ζ l = π (2l − 1), but it can be shown to converge much more quickly than the Matsubara expansion as N increases [15,16]. Unlike the Matsubara expansion, in which the poles are spaced at constant intervals of 2π β along the imaginary axis, the poles iζ l /β in the Padé expansion are spaced unevenly along the imaginary axis, and the prefactors η l are all real and positive-valued. In practice, the expansion f N (x) truncated at some finite value of N is used, with N chosen such that for values of L ε, where ε denotes the typical energy scale of the problem, the deviation δf N (L) ≡ |f (x) − f N (x)| βx=L < 10 −p , where p can be chosen to give arbitrary accuracy [15,16,17]. In practice we find good convergence for N ∼ 20 when the Padé parameters are used.
The expansion (19) is then substituted into (16), resulting in x-integrals that can be evaluated analytically with the residue theorem, for example: Using this method and the identity (17) one finally obtains an expression for the current in lead α that is asymptotically exact (as N → ∞), and which has all frequency integrals eliminated: Here we have introduced the digamma function Ψ, defined as the logarithmic derivative of the complex gamma function [5]. Importantly, the difference of two digamma functions appearing in the expression for I α (t) can be converted into a well converging numerical series thanks to the property: The general result (21) forms the basis for all subsequent numerical work involving the calculation of the current response to an arbitrary time-dependent bias. In practice, the time integrals may be performed numerically or removed analytically given the assumption of a particular functional form for V α (t). Several of the decaying modes in the current can be directly identified by inspection of (21), although the exact nature of its long-time behaviour depends on the bias chosen. In particular, if we separate out the real and imaginary parts of the complex eigenvalues,ε j = λ j − iγ j , we see that modes appearing in the single summation decay as e −γj (t−t0) , whereas the prefactor governing the decay of modes in the double sum is e −(γj +γ k )(t−t0) . This defines a timescale of τ ≡ M ax{1/γ j } as the time taken for transient behaviour to vanish following the switch-on.

Results
We shall now apply the method which was outlined in Section 3 to a molecular wire coupled to two leads, illustrated schematically in Fig. 1(a). We describe this system using a tight-binding model for the wire with nearest-neighbor hopping, as described in [18,19]. Specifically, we assume that the leads are connected by a wire of 5 sites, with one state per site, and that the interaction of the chain orbitals with those in the leads is only via the end sites. Within the WBLA, this translates into the condition that only the Γ 11 and Γ 55 elements in the level-width matrix are non-zero and are given by Γ 11 = 2π i∈L |V i1 | 2 δ ε F L − ε iL and Γ 55 = 2π j∈R |V N j | 2 δ ε F R − ε jR , respectively. For simplicity, we assume that for each site in the chain there is a single energy level, so that effectively we are describing a chain of quantum dots coupled to each other with a nearest-neighbor hopping, similarly to the system studied in [20], resulting in a tridiagonal molecular Hamiltonian matrix h CC , with elements [h CC ] k,k ≡ E and [h CC ] k,k+1 = [h CC ] k+1,k ≡ τ . Additionally, we choose the chemical potential µ to be zero.
We now focus on the case of a sinusoidal AC bias turned on in each lead, with a lead-dependent constant part, amplitude, frequency and phase shift: This enables us to remove analytically all the remaining integrals in (21), by use of the identity: where J r denotes a Bessel function of the first kind of order r. When Eq. (24) is substituted into (21), the single and double time integrals in (21) can be removed analytically and the problem of calculating the time-dependent current response is reduced to a simple embedded summation over left/right eigenvectors, the Padé parameters and Bessel functions. For the special case of the sinusoidal bias (23), this summation can be expressed entirely in terms of the special functions (18) and (22). We confirm that in the case of a single quantum dot our method reproduces exactly the I − t characteristic curves published in [8], which were obtained by direct numerical integration of (9). In this work, we apply the bias (23) to the leads L and R in the system shown in Fig. 1(a). The amplitudes A α and constant bias shifts V α are identical for each lead, as is the driving frequency of the bias, Ω α . The physical symmetry of this transport problem is then broken by a relative phase-shift of φ R = −π/2 of V R (t) with respect to V L (t), whose phase is fixed at φ L = 0. The tridiagonal central region Hamiltonian has all diagonal terms equal to 1, and all off-diagonal terms τ = 0.1. Other parameters in the model are V L = 5.0 = V R , A L = 4.0 = A R , Γ 11 = 0.5, Γ 55 = 0.5, and the inverse temperature β = 10. The resulting current in each lead is plotted in Fig. 1(b) (thick lines), along with the suitably rescaled bias applied to the leads (dotted lines). In addition, we plot the steady-state value of the current at each time, using the LB formula. This serves as a comparison between the best predictions of the steady-state formalism and our fully time-dependent method, from which it is instantly apparent that the TD current amplitude is an order of magnitude greater than that of the LB current. There is a much sharper peak in the transient current in lead L than in R, because the bias applied in L has the value V L + A L at the switch-on time t = 0, whereas the corresponding bias in lead R has the value V R . In the long-time regime, the current signal in each lead has the same form but the signal in R is phase-shifted by −π/2 with respect to the signal in L. Additionally, 'ringing' oscillations in the current occurring at higher frequencies than the driving frequency Ω α can be seen, as also observed in [5] and [8]. We identify the physical source of these frequencies in the energy gap |V α + µ − E| between onsite energies and the constant bias shift in the leads.
One can also see from Fig. 1(b) that, whereas the LB formalism preserves the condition of charge conservation, I L (t) + I R (t) = 0, at all times, this condition is not satisfied by the two TD currents. We understand this as follows: electrons propagate in the lead at the finite Fermi velocity v F α , so that the effects of a bias applied locally will take a finite time to propagate to other parts of the junction [4]. If a constant bias is switched on in each lead, then following a time delay equal to the transmission time of this signal across the junction, the rate of flow of electrons from C to L (I L ) must equal the rate of electron transport from R to C (-I R ), as in this case the number of electrons in the molecular region N C (t) = −iTr [G < CC (t, t)] is constant in time. If the bias is time-dependent, then the Hamiltonian is no longer symmetric under time translations and, in general, N C (t) also varies with time. Therefore local charge conservation does not apply in this case, a fact that is reflected in the current characteristics of Fig. 1(b).
It is a point of concern with Fig. 1(b) that it appears to indicate an unbroken net flow of charge into the molecule, because after the usual transient 'ringing' regime is over, the average value of I L (t) + I R (t) per cycle appears to be positive. This is a property which, if sustained indefinitely, means that there will always be a charge increase in the central region of the junction -clearly an unphysical situation. To investigate what is occurring here, we compute in Fig. 2 the transport properties over a much longer time range of the same 5-site model as was considered in Fig. 1(b). The long time plots of both the sum of the currents (blue curve) and the particle number N C (black) are displayed in Fig. 2, for this system, where to compute N C Figure 2: The electron number (black line) and sum of currents (blue line) for the sparse tight-binding hamiltonian is compared with the electron number (red line) and sum of currents (orange line) for a hamiltonian in which each lead is coupled with the same strength to every site in the chain.
we have expanded (11) into the left/right eigenbasis in a similar fashion to (21). Certain qualitative features of the plots stand out, in particular a periodic 'sloshing' of the charge in the molecular region as the number of particles increases there, before levelling off at a value below the maximum charge of 10. Simultaneously with this dynamic filling of the molecular levels, the average value of I L (t) + I R (t) drops until the integral under this curve evaluates to zero over a full cycle. Thus in addition to the fast decay of the initial transient, there is a secondary underlying decay which occurs over a much longer time scale. Note that this decay over two time scales is not possible with a single-level model, so that it is fundamentally a finite-size effect.
In Fig. 2 we compare the results obtained for the sparse tight-binding level-width matrix with the same plots for level-widths that are relatively non-sparse. Specifically, we take [Γ L ] ij = δ ij 0.5 = [Γ R ] ij , giving a level-width with a constant term along the main diagonal. The timescale for decay of all transient modes in formula (21) is τ ≡ M ax{1/γ j }, where the γ j are defined as imaginary parts of the complex eigenvalues, ε j = λ j − iγ j . For the sparse tight-binding level-width τ = 61.9, whereas in the non-sparse case we have a decay time of τ = 2.0. This order-of-magnitude difference is well illustrated by the extremely fast relaxation of the non-sparse number density (red curve) to a steady oscillation below N C = 10. This implies that the integral of I L (t) + I R (t) over one cycle is zero, implying no long-term change in the sum of currents, as shown by the orange curve in Fig. (2). Additionally, we note that despite the qualitatively similar 'sloshing' between two values of the electron number in the molecule after t ∼ τ has elapsed, there is a difference in the amplitude of the oscillating signal for the sparse and non-sparse cases. This is simply due to the fact that the overlap integral prefactors (for example ϕ L j Γ α ϕ R j ) are smaller on average when Γ α is sparse.

Conclusions
In this work we have shown that the general formula for the current in a multi-terminal junction derived in [8] can be manipulated into another form which is very convenient for numerical calculations. Using the Padé approximation for the Fermi function it is possible to convert infinite frequency integrals into fastconverging summations. In the case of a sinusoidal bias applied to the leads all the time integrals can be removed analytically and expressed in terms of special functions. To illustrate this formalism, we considered a molecular chain in a two-terminal junction with a −π/2 relative phase shift in the sinusoidal biases applied to both leads. The calculated current demonstrates similar transient and ringing effects to those observed previously [5,8] and reveals novel effects due to this symmetry-breaking. Specifically, we demonstrate the violation of local charge conservation in the molecular region. We anticipate that in future we will be able to apply this formalism to the study of quantum 'pumping', where symmetry-breaking parameters in the periodic driving signal lead to a net transfer of charge through the system [21]. We also find that for extended systems there can be two distinct decay processes in the current characteristics: an initial 'ringing' transient following the bias switch-on, and also an underlying decay of the entire signal implied by global charge conservation. The first type of decay is a well-known effect, and has been studied before for static [7] and TD [8] biases. The latter type of decay is characterized by the fact that the signal retains its shape throughout, and also by a decay time which is related to sparseness of the level-width matrix. We anticipate that this double decay effect will be observable in more realistic extended systems than the one considered here, and that the formalism developed here provides a fast way to access transient properties in those systems.