Dynamic Coulomb blockade in single-lead quantum dots

We investigate transient dynamic response of an Anderson impurity quantum dot to a family of ramp-up driving voltage applied to the single coupling lead. Transient current is calculated based on a hierarchical equations of motion formalism for open dissipative systems [J. Chem. Phys. 128, 234703 (2008)]. In the nonlinear response and nonadiabatic charging regime, characteristic resonance features of transient response current reveal distinctly and faithfully the energetic configuration of the quantum dot. We also discuss and comment on both the physical and numerical aspects of the theoretical formalism used in this work.


I. INTRODUCTION
Dynamics of open electronic systems far from equilibrium has attracted considerable interest in the past few years, since the dynamic processes are closely related to and thus can be used to explore the system properties of interest. It is well known that electron-electron interaction plays a significant role under realistic experimental conditions [1]. Therefore, understanding the electronic dynamics in presence of electron-electron interaction is of fundamental importance to the emerging field of nanosciences.
So far theoretical work has focused mostly on steady or quasi-steady state dynamics, which however may contain only partial information on the system of study. Consider, for example, a quantum dot (QD) as depicted in Fig. 1(a), with two Zeeman levels of ǫ ↑ < ǫ ↓ , in contact with a single electrode. Under the quasi-steady state charging condition, the applied voltage raises the electrode chemical potential (µ) adiabatically. The first charging electron occupies exclusively the energetically favorable spin-up level, while the second one that occupies the spin-closed double occupation state requires an excess energy to overcome the static on-dot Coulomb blockade; see Fig. 1(b). Apparently, the steady-state study offers rather incomplete information, as the single occupation of spin-down state is dark here.
In this work, we resort to transient dynamics of the system, by applying a time-dependent external voltage to the system, i.e., µ varies explicitly in real time, and analyzing the resulting nonadiabatic charging process. It is expected that all transitions within the range of applied voltage will be manifested through the system dynamic properties. Specifically, we consider a single-lead Anderson impurity QD driven by a family of ramp-up external voltages, by exploiting the hierarchical equations of motion (HEOM) theory for open electronic dynamics and transient current [2]. In Sec. II we briefly describe the * Electronic address: chxzheng@ust.hk † Electronic address: yyan@ust.hk theoretical and computational aspects. We demonstrate in Sec. III that the transient dynamics reveals sensitively and faithfully the QD energetic configuration. Section IV concludes this work.

II. METHODOLOGY
We implement a HEOM formalism, which is formally exact for dynamics of an arbitrary non-Markovian dissipative system interacting with surrounding baths [2,3,4,5]. From the perspective of quantum dissipation theory, the electronic dynamics of an open system is mainly characterized by the reduced system density matrix ρ(t) ≡ tr B ρ T (t) and the transient current through the coupling leads. Here ρ T (t) is the total density matrix of the entire system, and tr B denotes a trace over all lead degrees of freedom. The final HEOM is cast into a compact form of [2] Here, L · ≡ [H(t), · ] is system Liouvillian (setting ≡ 1). The basic variables of Eq. (1) are ρ(t) and associated auxiliary density operators (ADOs) ρ n (t), where n is an index set covering all accessible derivatives of the Feynman-Vernon influence functional [6]. For a single-level QD coupled to a single lead, the lead correlation function can be expanded by an exponential series via a spectral decomposition scheme [7,8,9]. In this case n involves anñ-fold combination of (σ, s, m), which characterizes the exponential series expansion with σ = ±, s the spin, and m the index of exponents, respectively. Therefore, in Eq. (1) ρ n |ñ =0 = ρ(t) with γ n |ñ =0 = ρ {−} n |ñ =0 = 0; ρ n |ñ >0 is an ADO at theñ th -tier; γ n (t) collects all the related exponents along with external voltages, to ρ n (t); and ρ {−} n and ρ {+} n are the nearest lowerand higher-tier counterparts of ρ n , respectively. In particular, the 1 st -tier ADOs, ρ n (t)|ñ =1 = ρ σ sm (t), determine exclusively the transient current of spin-s through the lead as [2] N s is the lead total-occupation operator of spin s; a s is the system annihilation operator of spin s; and tr T and tr denote trace operations over the total system and reduced system degrees of freedom, respectively. Since there is only one coupling lead, the displacement current at every time t is simply − s I s (t); thus the Kirchhoff's current law retains [10,11]. The present HEOMbased quantum transport theory admits classical geometric capacitors or capacitive coupling to a gate electrode. The associated displacement current can then be evaluated readily, leading to the conservation of total current [10,11,12]. Gauge invariance [10,11,13] is also guaranteed at all time, as it can be seen from the form of the 0 th -tier of the HEOM [2]. It has been proved [2] that for a noninteracting QD, the HEOM (1) is exact with a finite terminal tier ofñ max = 2; while for an interacting QD, in principle an infinite hierarchy is required for the exact transient dynamics. In practice a truncation scheme is inevitable. In this work a straightforward truncation scheme is adopted; that is to set all the higher-tier ADOs zero: ρ n |ñ >Ntrun = 0, with N trun the preset truncation tier. This scheme leads to a systematic improvement of results as N trun increases, i.e., by inclusion of higher-tier ADOs, as verified by extensive numerical tests. For a moderate system-lead coupling strength Γ < 5 T , quantitatively accurate (if not exact) dynamics is obtained with N trun = 2, where T is the temperature in the unit of Boltzmann constant (i.e., k B = 1). A higher N trun means a drastic increase in computational cost. Therefore, considering the tradeoff between numerical accuracy and computational efforts, all the following calculations are carried out with N trun = 2 by confining the ratio Γ/T less than 5; see Sec. IV for further discussions.
The composite system of interest is described by the Anderson impurity model [14], with its Hamiltonian H T expressed as H represents the interacting QD of our primary interest: Here, ǫ s is the energy of the spin-s (↑ or ↓) state, n s ≡ a † s a s the system occupation-number operator of spin s, and U the on-dot Coulomb interaction strength. In Eq. (3), h B = k s ǫ ksnks describes the noninteracting coupling lead, with ǫ ks being the energy of its single-electron state k of spin s, andn ks ≡ d † ks d ks the lead occupation-number operator. The QD-lead coupling is given by H SB = k s t ks d † ks a s + H.c., where t ks is the coupling matrix element between the lead state k and the QD-level, both with spin s. The widely used Drude model is adopted to capture the coupling lead, i.e., a Lorentzian spectral density function of J(ǫ) = ΓW 2 /(ǫ 2 + W 2 ), with W characterizing the lead bandwidth.
We consider explicitly the scenario depicted in Fig. 1(a). A ramp-up voltage V (t) is applied to the coupling lead at t = 0, which excites the QD out of equilibrium. The lead energy levels are shifted due to the voltage, ∆(t) = −eV (t), with e being the elementary charge, and so is the lead chemical potential µ(t) = µ eq + ∆(t). µ eq is the equilibrium lead Fermi energy, which is set to zero hereafter, i.e., µ eq = 0. Under a ramp-up voltage, ∆(t) varies linearly with time until t = τ , and afterwards is kept at a constant amplitude ∆: By tuning the duration parameter τ , the family of rampup voltage covers three distinct regimes: (A) The adiabatic limit at τ → ∞; (B) the intermediate range with a finite τ ; and (C) the instantaneous switch-on limit corresponding to τ → 0 + . These three regimes are individually explored and elaborated in the subsections of Sec. III. The work flow of our numerical procedures are briefly stated here, while details to be published elsewhere. (a) The lead correlation function is expanded by exponential functions, and hence establishes Eq. (1). (b) The equilibrium reduced density matrix and ADOs in absence of external voltages, {ρ eq n }, are obtained by setting both sides of Eq. (1) to zero. The linear sparse problem is solved by the biconjugate gradient method [15]. (c) The real-time evolution driven by V (t) is solved via a direct integration of Eq. (1) by employing either the Chebyshev propagator [16,17] or the 4 th -order Runge-Kutta algorithm [15]. The outcomes of step (b) are adopted as the initial conditions for the time evolution in step (c), i.e., ρ n (t = 0) = ρ eq n . For every time t, the transient current I s (t) is evaluated via Eq. (2). It normally takes about 1 days for a typical time evolution to 100 ps with a time step of 0.02 ps with a single Core 2 Duo processor.

A. The adiabatic limit
As τ approaches ∞, the infinitely slow application of external voltage results in quasi-static electronic dynamics for the open QD system. Consequently, at every time t the QD can be deemed as in the quasi-equilibrium determined by µ. This adiabatic charging scenario is introduced in Sec. I, which provides only partial information on the QD energetic configuration.
The QD physical subspace is completely spanned by the following four Fock states, i.e., |0 (vacancy), |↑ (single spin-up occupation), |↓ (single spin-down occupation), and |d (spin-closed double occupation). The QD gains an energy of ǫ s as it transits from |0 to |s , or ǫs + U from |s to |d , respectively. Heres labels the spin direction opposite to s. Each value of µ defines a steady state of zero current for the QD, together with the number of electrons, N e , residing on it. Shown in Fig. 1(b) is N e as a function of µ for a specific case of µ eq < ǫ ↑ < ǫ ↓ . The two plateaus corresponding to quantized N e are separated by the magnitude of U , and the vertical step prior to each plateau represents a Fockstate transition which adds one more electron to the QD. The two steps at ǫ ↑ and ǫ ↓ + U correspond to the transitions |0 ↔ |↑ and |↑ ↔ |d , respectively; while the other two transitions, |0 ↔ |↓ and |↓ ↔ |d associated with the excitation energies ǫ ↓ and ǫ ↑ + U , are missing from Fig. 1(b). This is because µ is considered to increase adiabatically; therefore, the first electron occupying the QD exclusively enters spin-up level that is energetically more favorable, and consequently, the second electron occupation requires an excess energy of U to overcome the static Coulomb blockade. We note that steady states offer rather incomplete information on the system of interest. To acquire further knowledge, it is inevitable to go beyond the adiabatic charging limit and explore the transient regime.

B. The finite switch-on duration regime
We now turn to cases of a finite switch-on duration 0 < τ < ∞, in which the lead level shift ∆(t) is exemplified in the inset of Fig. 2(b). The response current emerges at t = 0 immediately after the QD is excited out of equilibrium, and then fluctuates as the lead chemical potential µ sweeps over the various QD Fock states. After t = τ , the current vanishes gradually as the reduced system is drawn towards the steady state with µ(t > τ ) = ∆. In Fig. 2(a)   labeled by numbers from 1 to 4, which are responsible for all the four resonances consecutively activated by the ramp-up voltage with an increasing excitation energy of ǫ ↑ , ǫ ↓ , ǫ ↑ +U , and ǫ ↓ +U , respectively. Figure 2(b) depicts I s (t) for both spins with τ = 50 ps. Within the period of t < τ , the major peaks in I ↑ (t) (peak-1) and I ↓ (t) (peak-4) start to form as the time-dependent µ(t) matches ǫ ↑ and ǫ ↓ +U , respectively; and they remain discernible even for a rather large τ ; see the line for τ = 90 ps in Fig. 2(a). However, the amplitudes of the two minor peaks, peak-3 in I ↑ (t) and peak-2 in I ↓ (t), diminish rapidly as the ramp-up duration τ is lengthened. In the adiabatic limit (τ → ∞), the minor peaks become completely invisible, since the spin-up and spin-down QD-levels are charged sequentially due to the Coulomb blockade; see Fig. 1(b). Influence of Coulomb interaction on the transient current is illustrated by Fig. 3. Cases of different U ranging from a noninteracting scenario to a strongly blockaded QD are studied. For U = 0, there is only one peak in I s (t) for either spin, since the excitation energies ǫ s and ǫ s + U are indistinguishable. The peak is split into two with an increasing displacement as U is augmented. For either spin direction, the first peak always sticks to its original position as in U = 0, while the second one emerges later in time with a larger U . For U as large as 5 meV, the two peaks in I ↓ (t) are almost completely sep- Other parameters are same as in Fig. 2(b).
arated. This confirms our attribution of the peaks to the various resonances: the first peaks in I ↑ (t) (peak-1) and I ↓ (t) (peak-2) correspond to the resonances µ = ǫ ↑ and µ = ǫ ↓ , respectively; and the second ones in I ↑ (t) (peak-3) and I ↓ (t) (peak-4) are associated with µ = ǫ ↑ + U and µ = ǫ ↓ + U , respectively. Thermal effect is also explored. I s (t) under different lead temperatures are plotted in Fig. 4. A series of satellite oscillations appears following each resonance peak in I s (t) at a low T , which is resulted from a time-varying phase factor, exp[ i t τ dt ∆(t)], of the nonequilibrium lead correlation function. The oscillation frequency is thus determined by the magnitude of ∆(t). For a ramp-up voltage, ∆(t) increases monotonically during 0 < t < τ . It is thus reasonable to find the frequency of satellite oscillations grows with time, and the decaying amplitude of the oscillations is due to the dissipative interactions between the QD and the semi-infinite lead.
To summarize, with a finite switch-on duration τ , the QD energetic configuration is resolved distinctly (sometimes completely) through the resonance peaks of response current in real time, in obvious contrast to the adiabatic charging limit.

C. The instantaneous switch-on limit
As τ approaches zero, a ramp-up voltage becomes asymptotically a step pulse, i.e., ∆(t) = ∆Θ(t), with Θ(t) being the Heaviside step function. Upon the instantaneous switch-on of external voltage, both the spin-up and down QD-levels are activated simultaneously, provided the voltage amplitude ∆ is sufficiently large. It is intriguing to see in this limit how the response current would be related to the QD energetic configuration.
In Fig. 5 (a) and (b) we plot the calculated I ↑ (t) and I ↓ (t) in response to a step voltage of ∆ = 10 meV, respectively. Three QDs of different ǫ ↑ are studied. The energetics scenario for the three QDs are µ eq = ǫ ↑ < ǫ ↓ , µ eq < ǫ ↑ < ǫ ↓ , and ǫ ↑ < µ eq < ǫ ↓ , respectively. Note that for all cases are in a large applied voltage regime since ∆ > ǫ s and ∆ > ǫ s + U , which ensure all the four QD Fock states accessible energetically. For both spin directions, the transient current shows an initial overshooting, and then oscillates on top of a decay as the reduced system evolves into a steady state with the lead chemical potential µ = ∆. These rapid oscillations are not artifacts. Their presence reflects the electronic dynamics of QD driven by an external voltage of large amplitude. Upon a small perturbation on ǫ ↑ of 0.1 meV, I ↑ (t) scarcely changes as those with ǫ ↑ = ± 0.1 meV almost overlap with the ǫ ↑ = 0 case shown in Fig. 5(a); while I ↓ (t) displays drastic deviations among the three cases, in terms of the overshooting amplitude, the oscillation frequency, as well as the decaying rate.
To understand the physical origin of the above observations, frequency-dependent current spectrum is calculated, i.e., I s (ω) ≡ F[I s (t)] with F denoting the conventional Fourier transform. The results corresponding to the three QDs are depicted in Fig. 6. For all three cases, two characteristic frequencies are revealed in Re[I ↑ (ω))], i.e., a major kink at ∆ − ǫ ↑ and a minor wiggle around ∆ − ǫ ↑ − U . The line shape of Re [I ↑ (ω)] around the frequency ω ↑ ≡ ∆ − ǫ ↑ much resembles that of a singlelevel QD [18]. For ǫ ↑ < µ eq , Re [I ↑ (ω)] exhibits a peak at ω ↑ ; while for ǫ ↑ > µ eq , a dip shows up. In the case of ǫ ↑ = µ eq , Re [I ↑ (ω)] around ω ↑ is described by the function form of ln(x 2 +1)/x, with x = 2(ω−ω ↑ )/Γ. However, the profiles of Re[I ↓ (ω)] are significantly different among the three cases. To be specific, for ǫ ↑ = µ eq , two prominent dips show up at frequencies ∆ − ǫ ↓ and ∆ − ǫ ↓ − U , respectively; for ǫ ↑ > µ eq , the dip at ∆ − ǫ ↓ − U is considerably suppressed, which is associated with the transition |↑ ↔ |d ; and for ǫ ↑ < µ eq , it is the dip at ∆ − ǫ that is almost completely quenched, which indicates an inactive transition |0 ↔ |↓ . Hence it is demonstrated the response current spectrum varies sensitively to the QD energetic configuration.
The presence (or absence) of a resonance signature in Re[I s (ω)], as well as the issue of peak versus dip at ω ↑ , are closely related to the equilibrium occupancy of the spin-up state prior to the switch-on of the voltage. Specifically, in the case of ǫ ↑ = 0.1 meV (> µ eq ), the initial population of spin-up electrons on the QD is as small as 0.08, and the subsequent incoming spin-down electrons driven by the step voltage would feel scarcely any Coulomb repulsion due to the almost vacant spin-up level. This thus leads to the prominent dip at ∆ − ǫ ↓ and the barely recognizable one at ∆−ǫ ↓ −U in Re[I ↓ (ω)]; see Fig. 6(b). Whereas in the case of ǫ ↑ = −0.1 meV (< µ eq ), the spin-up level is nearly fully occupied at t = 0, i.e., N ↑ (0) = 0.92. Therefore, the majority of spin-down electrons injected afterwards need to acquire an excess energy of U to overcome the Coulomb blockade, and this is why the dip at ∆ − ǫ ↓ becomes trivial; see Fig. 6(f). For ǫ ↑ = µ eq , the spin-up level is half occupied (N ↑ = 0.5) prior to the switch-on of voltage. In such a circumstance, both the resonant transitions |0 ↔ |↓ and |↑ ↔ |d contribute equally to the response current, and hence give rise to the two dips of similar depth at frequencies ∆ − ǫ ↓ and ∆−ǫ ↓ −U ; see This thus confirms that in the instantaneous switchon limit, the ω-dependent response current spectrum resolves sensitively and faithfully the QD energetic configuration through the characteristic resonance signatures. To our surprise, Re [I ↓ (ω)] resolves a third peak at ω ↑ ; see the tiny bump in Fig. 6(b) and its inset. It possibly arises due to the higher-order response, and the mechanism of its presence requires further understanding on the tran-sient dynamic processes. At a lower temperature, all the resonance signatures in the current spectrum are accentuated (not shown), i.e., the peaks (dips) become higher (deeper) and sharper. In the real time, this amounts to an enhanced oscillation amplitude of transient current.

IV. CONCLUDING REMARKS
To conclude, we have investigated the transient electronic dynamics of a single-lead interacting QD driven by a family of ramp-up voltages. We have found that (a) in the adiabatic charging limit, the quasi-steady dynamics provide only incomplete information on the QD energetic configuration; (b) With an appropriately chosen switch-on duration, all transitions among QD Fock states can be resolved distinctly through the resonance peaks of transient current in real time, and hence provide comprehensive (sometimes complete) information on the QD energetic configuration; (c) In the instantaneous switch-on limit, the frequency-dependent responses current spectrum reveals sensitively and faithfully the QD energetic configuration through the characteristic resonance signatures. This work thus highlights the significance and versatility of transient dynamics calculations.
In the sequential tunneling regime where Γ ≪ T (e.g., Figs. 1−3), the present HEOM with truncation tier of N trun = 1, which amounts to the conventional quantum master equation methods with memory [9], will be sufficient. When Γ ∼ T , calculations with N trun = 1 would produce qualitatively correct static and dynamic Coulomb blockade for a single-lead QD, because the basic physics does not involve the cotunneling mechanism. However, to obtain quantitatively accurate results, N trun = 2 is needed. Extensive numerical tests indicate that N trun = 2 can support up to Γ ∼ 5T (e.g., Figs. 5 and 6). It has been shown [2] that the present HEOM theory with N trun = 2, which properly treats the cotunneling dynamics, is equivalent to a real-time diagrammatic formalism [19]; thus the Kondo problems can be addressed. However, a higher truncation tier is numerically needed for convergency when T ≪ Γ. The quantitatively accurate results presented in this work (for Γ < 5T ) are useful for further development of efficient but maybe approximate computational schemes.
Investigating and/or manipulating transient electronic dynamics provide(s) sometimes a unique way to achieve physical properties and novel functions of mesoscopic many-body systems. However, exact theoretical results for transient quantum transport are rare. Weiss et al [20] have recently proposed an iterative real-time path integral approach, which is numerically exact but expensive especially when the applied bias voltage is beyond some simple forms such as a step function. The present HEOM theory is formally exact and provides an alternative approach. It is applicable to arbitrary time-dependent external fields applied to leads and/or to the system, without extra computational cost. However, computational effort increases dramatically if the converged results require a high truncation tier of the hierarchy [2].