Full-counting statistics of transient energy current in mesoscopic systems

TL;DR

This paper develops a non-equilibrium Green's function approach to analyze the full-counting statistics of transient energy currents in mesoscopic quantum dot systems, revealing oscillatory behaviors and scaling laws.

Contribution

It introduces an exact formalism for calculating the full-counting statistics of energy transport in transient regimes of quantum dots, including effects of finite lead bandwidths.

Findings

Transient energy current shows damped oscillations.

Oscillation frequency in single QD is bandwidth-independent.

Double QD exhibits Rabi oscillations influenced by interdot coupling.

Abstract

We investigate the full-counting statistics (FCS) of energy flow carried by electrons in the transient regime. Based on two measurement scheme we formulate a non-equilibrium Keldysh Green's function theory to compute the generating function for FCS of energy transport. Specifically, we express the generating function using the path integral along Keldysh contour and obtain exact solution of the generating function using the Grassmann algebra. With this formalism, we calculate the transient energy current and higher order cumulants for both single and double quantum dot (QD) systems in the transient regime. To examine finite bandwidth effect of leads to FCS of energy transport, we have used an exact solvable model with a Lorentizian linewidth where all non-equilibrium Green's functions can be solved exactly in the time domain. It is found that the transient energy current exhibits damped oscillatory behavior. For the single quantum dot system the frequency of oscillation is independent of bandwidth of the leads while the decay rate of the oscillation amplitude is determined by the lifetime of resonant state which increases as the bandwidth decreases. At short times, a universal scaling of maximum amplitude of normalized cumulants is identified for the single QD system. For the double QD system, the damped oscillation of energy current is dominated by Rabi oscillation with frequency approximately proportional to the coupling constant between two quantum dots. In general, the transient energy current increases when the coupling between two QDs is stronger. However, when the interdot coupling is larger than half of the external bias the transient energy current is suppressed significantly. All these results can be understood analytically.

Equations

$$Z(\lambda ,t)\equiv ⟨e^{i\lambda \Delta ϵ}⟩=\sum _{\Delta ϵ}P(\Delta ϵ,t)e^{i\lambda \Delta ϵ}.$$

$$\langle\langle(\Delta\epsilon)^{j}\rangle\rangle=\frac{\partial^{j}\ln Z(\lambda)}{\partial(i\lambda)^{j}}\bigg{|}_{\lambda=0}.$$

$$H=\sum_{k\alpha}\epsilon_{k\alpha}c^{\dagger}_{k\alpha}c_{k\alpha}+\sum_{n}\epsilon_{n}d_{n}^{\dagger}d_{n}+\sum_{k\alpha n}\Big{(}t_{k\alpha n}c^{\dagger}_{k\alpha}d_{n}+\mathrm{h.c.}\Big{)},$$

$$H_L=\sum _k{ϵ}_{kL}{c}_{kL}^{†}{c}_{kL}.$$

$$Z(\lambda ,t)=\mathrm{Tr}\left[\rho (0)e^{i\lambda H_L(0)}{e}^{-i\lambda H_L^h(t)}\right].$$

$$H_L^h(t)=U^{†}(t,0)H_L(0)U(t,0),$$

$$Z(\lambda ,t)=\mathrm{Tr}\left\{\rho (0)U_{\lambda /2}^{†}(t,0)U_{-\lambda /2}(t,0)\right\},$$

$$U_{\gamma }(t,0)=\mathcal{T}\exp \left[-\frac{i}{\hslash }{\int }_0^t{H}_{\gamma }(t^{\prime })d{t}^{\prime }\right],$$

$$H_{\gamma }$$

$$e^{i\gamma H_L}{c}_L(0)e^{-i\gamma H_L}=\sum _n\frac{{\hslash }^n{\gamma }^n}{n!}[{\partial }_t^n{c}_{kL}(t){]}_{t=0}=c_L(t_{\gamma }).$$

$$Z(\lambda ,t)=\int D[\bar{\varphi }\varphi ]e^{iS[\bar{\varphi }\varphi ]},$$

$$S[\bar{\varphi }\varphi ]$$

$$S[\bar{\Psi }\Psi ]={\int }_0^{t}d\tau {\int }_0^{t}d{\tau }^{\prime }{\bar{\Psi }}^T(\tau )M(\tau ,{\tau }^{\prime })\Psi ({\tau }^{\prime }),$$

$$M=\left(\begin{array}[]{ccc}g^{-1}_{kk^{\prime}L}(\tau_{\gamma},\tau^{\prime}_{\gamma})&-t_{kLn^{\prime}}\delta&0\\ -t^{*}_{k^{\prime}Ln}\delta&g^{-1}_{nn^{\prime}}(\tau,\tau^{\prime})&-t^{*}_{k^{\prime}Rn}\delta\\ 0&-t_{kRn^{\prime}}\delta&g^{-1}_{kk^{\prime}R}(\tau,\tau^{\prime})\\ \end{array}\right),$$

$$Z(\lambda ,t)=\det (G{\tilde{G}}^{-1}),$$

$$G^{-1}$$

$${\tilde{G}}^{-1}$$

$${\tilde{\Sigma }}_L(t,t^{\prime })$$

$${\tilde{\Sigma }}_L^r$$

$${\tilde{\Sigma }}_L^k$$

$${\tilde{\Sigma }}_L^{\bar{k}}$$

$$\ln Z(\lambda ,t)=\mathrm{Tr}\ln [I-G({\tilde{\Sigma }}_L-{\Sigma }_L)],$$

$$I^{E}_{L}(t)=2\mathrm{Re}\int dt^{\prime}\mathrm{Tr}\big{[}G^{r}(t,t^{\prime})\breve{\Sigma}^{<}_{L}(t^{\prime},t)+G^{<}(t,t^{\prime})\breve{\Sigma}_{L}^{a}(t^{\prime},t)\big{]},$$

$${\breve{\Sigma }}_L^{\chi }(t^{\prime },t)=\sum _k{ϵ}_{kL}{\Sigma }_{kL}^{\chi }(t^{\prime }-t),$$

$$\ln Z(\lambda ,t)=-\mathrm{Tr}[G({\tilde{\Sigma }}_L-{\Sigma }_L)],$$

$$\ln Z(\lambda ,t)=\mathrm{Tr}[G^{<}({\tilde{\Sigma }}_L^{>}-{\Sigma }_L^{>})+G^{>}({\tilde{\Sigma }}_L^{<}-{\Sigma }_L^{<})].$$

$$\ln Z(\lambda ,t)=(n_d-1)M_{L1}+n_d{M}_{L2},$$

$$M_{L1}$$

$$M_{L2}$$

$$A_0(E)=\frac{4{\Gamma }_L(E)}{\pi }\frac{{\sin }^2[(E-{ϵ}_0)t/2]}{(E-{ϵ}_0{)}^2}.$$

$$(-i{)}^n{\partial }^{n}Z(\lambda ,t)/\partial {\lambda }^n{|}_{\lambda =0}=x(t),$$

$$C_n=\frac{{\partial }^n\ln Z}{\partial (i\lambda {)}^n}={\int }_{-\infty }^{{\Delta }_L}dE{A}_0(E)E^n,$$

Full text

††preprint: APS/123-QED

I Introduction

With the rapid development of nanotechnology, the transport of single electron can now be studied in electronic devices Lu ; Bylander ; Fuji , leading to renewal interest on statistical distribution of electrons and the energy it carries. Full-counting statistics is a new methodology to characterize full probability distribution of electron and energy transport by calculating the corresponding generating function Belzig ; Bagrets ; Pilgram ; Gogolin ; Scho ; Saito ; Flindt ; Urban ; T1 . Experimentally the real-time counting of electrons has been carried out in a quantum dot (QD) which direct measures the distribution function of current fluctuations Bylander ; Gust ; Flindt2 . In the experiment, the measurement of higher order cumulants up to 15th has also been reported in quantum point contact systems.Flindt2 Such a stochastic process shows universal fluctuation for noise and higher order moments which is a common feature in mesoscopic systems. Full-counting statistics offers a superior way to study the noise and those all higher order correlations. Furthermore, quantum entanglement and quantum information have been also related to quantum transport in terms of full-counting statisticsKlich ; Song ; Les .

One of the important issues in the non-equilibrium transport process is the energy transport which gives information on how energy is dissipated and its correlations for any working electronic devices. Recently, the energy transport through the one-dimensional system, such as the trapped ion chains, has been measured experimentally Ramm . The energy dissipation and fluctuation can be characterized by the energy current, which can be investigated by Landauer-Büttiker type of formalism theoretically in the dc transport for non-interaction electrons But ; sivan ; kearney . The energy current $I^{E}_{\alpha}$ is also related to the heat current $I^{h}_{\alpha}$ by $I^{h}_{\alpha}=I^{E}_{\alpha}-\mu I_{\alpha}$ where $\mu$ is the chemical potential and $I_{\alpha}$ is the particle current. It is known that the Joule heating $J$ due to the leads is related to the heat current by $\sum_{\alpha}I^{h}_{\alpha}=J$. The heat current driven by external bias and temperature gradient plays a central role in studying the efficiency of nano-scaled heat engineengine1 . Recently ac heat current, its shot noise, and relaxation resistance have been investigated in mesoscopic systemssanchez1 ; sanchez2 ; mosk1 ; mosk2 ; jchen . Using the non-equilibrium Green’s function transient heat current has been studied in mesoscopic systemsMichelini . Moreover, A first-principles calculation for transient heat current through molecular devices has also been carried outYu . It would be interesting to further study the FCS of energy transport in ac regime.

We note that FCS of energy transport of phonon has been studied extensively. An exact formula for cumulant generating function of heat transfer has been derived in harmonic networks to study non-equilibrium fluctuations SD1 ; SD2 . Moreover, energy fluctuations in a driven quantum resonator has also been studied by full-counting statistics C1 ; C2 . Using the phonon non-equilibrium Green’s function, the generating function has been obtained for phonon transport in the transient regime as well as steady statesW1 ; W2 ; W3 . Various cumulants of thermal current and entropy production have been studied numericallyW1 ; W3 ; W4 . So far most of investigations of FCS of energy transfer focus on the phonon transport. However, less attention has been paid on FCS of the electron energy transferfoot5 . It is the purpose of this paper to address this issue.

In this paper, we develop a Keldysh non-equilibrium Green’s function (NEGF) theory to study FCS of transferred energy in transient regime. Based on two measurement scheme we derive the expression of generating functions for FCS in terms of non-equilibrium Green’s functions. This allows us to calculate $n$th cumulant $C_{n}$ of transferred energy of mesoscopic systems in the transient regime as well as in long-time limit. We then apply our formalism to investigate the energy transport in the transient regime for both single and double QD systems. To study the finite bandwidth effect on the cumulants of transferred energy, an exact solvable model is used with Lorentzian linewidth so that the non-equilibrium Green’s function can be obtained exactly. As expected, the cumulants of transferred energy show linear characteristics in time in the long-time limit for both systems. For the single QD, the transient energy current exhibits damped oscillatory behavior. The oscillation frequency is found to be independent of bandwidth of the leads and the decay of oscillation amplitude is proportional to the lifetime of resonant state of QD which decreases as the bandwidth increases. At short times, the maximum amplitude $M_{n}$ of the normalized $n$th cumulant $C_{n}(t)/C_{1}(t)$ show a universal behavior for the single QD system. Specifically, we find $M_{2k}=a_{1}e^{\kappa k}$ and $M_{2k+1}=a_{2}e^{\kappa k}$ for different system parameters where $a_{1}$ and $a_{2}$ are non-universal constants. The universal slope $\kappa$ is found to be close to 3. For the double QD system, we find that the transient energy current shows the damped Rabi oscillations with frequency approximately proportional to the interdot coupling constant $v$ between two QDs. A threshold of interdot coupling $v_{c}$ is found below which transient energy current increases with the increase of $v$ while for $v>v_{c}$ the transient current is suppressed significantly. These interesting results can be understood analytically.

This paper is organized as follows. In Sec. II, the formalism of generating function for studying full-counting statistics of transferred energy in the transient regime is first presented. In Sec. III, we apply the formalism obtained to both single and double QD systems and show numerical results of various cumulants, transient energy current and the corresponding higher order cumulants. Finally, the discussion and conclusion are given in Sec. IV.

II Theoretical Formalism

To study full-counting statistics of energy transport, we need to obtain the probability distribution $P(\Delta\epsilon,t)$ of the transferred energy carried by electrons $\Delta\epsilon=\epsilon_{t}-\epsilon_{0}$ between an initial time $t_{0}$ (for simplicity, we set $t_{0}=0$) and a later time $t$ which can be calculated from two-time quantum measurement. Denoting $\epsilon$ the eigenvalue of the Hamiltonian of the left lead $H_{L}$ where we measure the energy flow. Taking measurement at $t$ gives $\epsilon_{t}$ which is a stochastic variable. The generating function $Z(\lambda,t)$ with the counting field $\lambda$ can be obtained by the Fourier transformation of the probability distribution as Scho ,

[TABLE]

The $j$th cumulant of transferred energy $\langle\langle(\Delta\epsilon)^{j}\rangle\rangle$ is defined by,

[TABLE]

We now derive the generating function using NEGF theory for a general QD system coupled with two semi-infinite leads in the transient regime. For this purpose, we assume that the couplings between the QD and leads are turned on at $t=0$. The Hamiltonian of the whole system can be written as,

[TABLE]

where, $c^{\dagger}(c)$ and $d^{\dagger}(d)$ are the creation (annihilation) operators of leads and QD, respectively. $\epsilon_{n}$ is the energy level for the QD and $\epsilon_{k\alpha}$ is the energy levels of the lead $\alpha(\alpha=L,R)$. $t_{k\alpha n}$ is the coupling constant between two leads and the QD.

To investigate the energy current through the left lead where the measurement is made, we focus on the energy operator of the left lead

[TABLE]

Since we study the behaviour in the transient regime we assume that the bias is applied to the leads at $t=-\infty$ while the leads and QD are disconnected. All the couplings are switched on at $t=0$ giving rise to a transient energy current. The switching of coupling between QD and leads can be done by a quantum point contact that is controlled by a gate voltage. Since the system is disconnect before $t=0$, the initial density matrix of the whole system at time [math] is the direct product of the subsystems expressed by $\rho(0)=\rho_{L}\otimes\rho_{D}\otimes\rho_{R}$. Similar to the cases of phonon and electron charge transportW2 ; T2 , the generating function of transferred energy can be expressed as,

[TABLE]

Here, $H^{h}_{L}(t)$ denotes the energy operator of the lead $L$ in the Heisenberg picture, which is related to the energy operator in the Schrödinger picture $H_{L}(0)$ (Eq. (4)) by

[TABLE]

where $U(t,0)$ is the evolution operator.

In terms of the modified Hamiltonian $H_{\gamma}$ given in Eq. (9), the generating function can be rewritten as,

[TABLE]

where the modified evolution operator is,

[TABLE]

with

[TABLE]

with $t_{\gamma}=\hbar\gamma$ and $\gamma=\lambda/2$. In deriving Eq. (9), we have used the following relation,

[TABLE]

Using Grassmann algebra the generating function becomesT2 ,

[TABLE]

where $D[\bar{\phi}\phi]=\Pi_{x\sigma}d\bar{\phi}^{\sigma}_{x}d\phi^{\sigma}_{x}$ with $x\in k\alpha,n$ and the action $S[\bar{\phi}\phi]$ is given by,

[TABLE]

where $\phi$ and $\bar{\phi}$ are the Grassmann variables which are two independent complex numbers gras and $\sigma=+,-$ denoting the upper and lower branches of the Keldysh contour, respectively.

After the Keldysh rotationKamenev we rewrite the action in Eq. (II) in a matrix form,

[TABLE]

with $\bar{\Psi}^{T}(\tau)=[\bar{\psi}_{kL}^{T}(\tau_{\gamma}),\bar{\psi}_{n}^{T}(\tau),\bar{\psi}_{kR}^{T}(\tau)]$ and

[TABLE]

where $\delta$ is a unit matrix in Keldysh time space.

Using the functional integration of the Gaussian integral for the Grassmann fields, the generating function can be expressed by the Keldysh non-equilibrium Green’s function as T2 ,

[TABLE]

with

[TABLE]

Here, $g=g_{nn^{\prime}}(\tau,\tau^{\prime})$ is the Green’s function of the isolated QD in Keldysh space. $\Sigma_{R}=\sum_{kk^{\prime}}t^{*}_{k^{\prime}Rn}g_{kk^{\prime}R}(\tau,\tau^{\prime})t_{kRn^{\prime}}$ is the self-energy of the right lead in the Keldysh space in the time domain. $\widetilde{\Sigma}_{L}=\sum_{kk^{\prime}}t^{*}_{k^{\prime}Ln}g_{kk^{\prime}L}(\tau_{\gamma},\tau^{\prime}_{\gamma})t_{kLn^{\prime}}$ is the self-energy with the counting field, meaning that the two-time measurement is done in the left lead. It is easy to find the expression of $\widetilde{\Sigma}_{L}$ given by,

[TABLE]

where

[TABLE]

with $\Sigma_{L}^{k}=2\Sigma_{L}^{<}+\Sigma_{L}^{r}-\Sigma_{L}^{a}$.

Note that the expression of Eq. (15) is the same as that of FCS of charge transport in Refs. T2, and ep1, . The difference lies in the expression of $\widetilde{\Sigma}_{L}$. For the charge transport the counting field appears as an extra phase while for the energy transport the presence of counting field is to shift the time from $t$ to $t_{\gamma}$. In Eq. (15) the generating function is expressed in terms of a determinant in both time and space domains. As a result, the calculation of generating function is very time consuming for realistic systems. To calculate $Z(\lambda,t)$ for fixed $\lambda$ and $t$ for a system with $N$ degrees of freedom, we need to discretize the time $t$ into $N_{t}$ uniform mesh and calculate the Green’s function as a function of both position and time. The computational complexity to evaluate the determinant is of order $N^{3}N_{t}^{3}$ for each $t$ foot3 . For this reason, the transient calculation of FCS is limited to a single QD or double QD with $N=1$ or $2$ where the exact solution of time-dependent non-equilibrium Green’s function is available. If one wishes to study a realistic system with $N=100$ (for instance), the amount of calculation increases by six orders of magnitude. However, if we are interested in the first a few cumulants, we can first find their expressions using Eqs.(2) and (15) and then calculate them numerically using these expressions rather than calculating cumulant generating function numerically.

Now we examine the limiting cases of generating function defined in Eq. (15). First of all, we look at the energy current in the transient regime. From Eq. (15) and using the relation $\ln\det A=\mathrm{Tr}\ln A$, the cumulant generating function is given by

[TABLE]

According to Eq. (2), the transient energy current can be expressed as (see derivation in appendix A and we have set $\hbar=e=1$),

[TABLE]

where

[TABLE]

with $\chi=<,a$ and $\Sigma^{\chi}_{kL}(t^{\prime}-t)$ being the self-energy of the left lead in the absence of the counting field. This expression of transient energy current agrees with that obtained directly by the Green’s function method Yu .

We now consider the short time behavior of the generating function. Using the fact that $G(\widetilde{\Sigma}_{L}-\Sigma_{L})$ is of order $t^{2}$, we find for small t

[TABLE]

where the quantities in trace are in Keldysh space. It is straightforward to show

[TABLE]

In the limit of weak coupling case for a single quantum dot system with a single level, we can use Green’s function of isolated quantum dot to replace $G^{<}$ and $G^{>}$ and we find

[TABLE]

where $n_{d}$ is the initial occupation number of the isolated QD and

[TABLE]

with

[TABLE]

where $\Gamma_{L}(E)$ is the linewidth function of the left lead. Here, $\epsilon_{0}$ is the energy level of the QD and $\alpha(E)=1$ or $E$ for charge transport and energy transport, respectively. In the wideband limit (WBL) and $\alpha=1$, Eq. (27) recovers the result of Ref. ep1, from which a universal behavior of $n$th cumulant of charge transport $C_{n}$ has been derived that was first demonstrated experimentally in Ref. Flindt2, . Note that an important relation holds for charge transport when $n_{d}=0$, i.e.,

[TABLE]

which is independent of $n$ with $n>0$. This allows one to obtain an analytic expression for $C_{n}$ for charge transport in the short time limit and very weak coupling regime leading to this universal behavior.

For the energy transport with $\alpha=E$, Eq. (31) does not hold anymore. Although no analytic expression is available for $n$th cumulant of energy transport $C_{n}$, some asymptotic behavior can be derived. Assuming $n_{d}=0$, i.e., initially there is no electron in the QD. From Eq. (27), it is straightforward to find the expression of $n$th order cumulant of energy transport ($\alpha=E$) at zero temperature

[TABLE]

where $\Delta_{L}$ is the bias voltage of the left lead. Now we use WBL such that $\Gamma(E)$ is a nonzero constant only for $|E|<W$Urban . In the regime $W=\Delta_{L}\gg|\epsilon_{0}|>0$, we have

[TABLE]

where $u=E-\epsilon_{0}$. Since $A_{0}(u)$ is an even function, this integral can be evaluated in the large W limit. For even $n=2k$, the major contribution comes from $\int duA_{0}(u)u^{2k}$ which gives

[TABLE]

with $a_{1}=\frac{4\Gamma_{L}}{(2k-1)\pi}$, the next order $W^{2k-2}$ term depends on $t$. For $C_{2k+1}$, it is dominated by $\int duA_{0}(u)(2k+1)\epsilon_{0}u^{2k}$ from which we find

[TABLE]

with $a_{2}=\frac{4\Gamma_{L}(2k+1)}{(2k-1)\pi}\epsilon_{0}$. Denoting $F_{n}=\ln(C_{n}/C_{1})$, we have $F_{2k}=\ln(a_{1}/C_{1})+(2k-1)\ln W$ and $F_{2k+1}=\ln(a_{2}/C_{1})+(2k-1)\ln W$. This suggests that both $F_{2k}$ and $F_{2k+1}$ depend linearly on $k$ with the same slope but different intercepts. Since the slope depends only on the bandwidth $W$, it is universal. Obviously, the above discussion is qualitative under certain limits, the detailed numerical study on the universal behavior at short time will be presented in the next section.

Now we investigate the long-time behavior of the generating function in the transient regime. When $t$ goes to infinity, the Green’s function and self-energy in the time domain become invariants under the time translation T2 . Therefore, the cumulant generating function in Eq. (22) in the long-time limit in the energy space becomes,

[TABLE]

In the next section, we will give numerical result of FCS of transferred energy in the transient regime. We will study FCS for two systems, single QD and double QD systems. In calculating the generating function numerically from Eq. (15), the Green’s functions defined in Eq. (16) for an occupied single QD can be expressed as,

[TABLE]

with the Heaviside step function $\theta(\tau_{1}-\tau_{2})$. For a double QD system, the Green’s function with the counting field in Eq. (17) should be written as,

[TABLE]

Here, $g^{-1}_{1}$ and $g^{-1}_{2}$ are the Green’s function for the first and second isolated QD, respectively. $v$ is the coupling constant between two QDs.

In order to calculate the self-energy $\Sigma_{L(R)}$ and $\widetilde{\Sigma}_{L}$ in Eqs. (16) and (17), we use the Lorentzian linewidth function to describe the self-energy so that the equilibrium energy dependent self-energy can be written using a finite band width $W$,

[TABLE]

with the linewidth amplitude $\Gamma_{\alpha}$. This is a special model that allows us to find the Green’s function exactly while still going beyond the WBL. Note that one can not tune the bandwidth experimentally. Then the retarded self-energy can be given by,T2

[TABLE]

where $\Delta_{\alpha}$ is the external bias applied on the lead $\alpha$. In order to calculate the lesser Green’s function analytically, we focus on zero temperature. We find $\Sigma^{<}_{\alpha}(\tau_{1}-\tau_{2})=\frac{i}{8}\Gamma W$ for $\tau_{1}=\tau_{2}$, and otherwise T2

[TABLE]

with $\tau=\tau_{1}-\tau_{2}$ and $\mathrm{Ei}(x)=-\int_{-x}^{\infty}\frac{e^{-t}}{t}dt$. Note that the diagonal element of lesser self-energy diverges for large $W$. It has been confirmed in Ref.joseph, that the transient charge current at WBL can be obtained as follows: calculating transient current as a function of $W$ and then taking large $W$ limit. We have confirmed that transient energy current at WBL can be obtained similarly.

Before we end this section, we mention that the approach presented in this paper is suitable only for non-interacting problems. Under a special situation where electrons couple with a single phonon mode, this type of approach can be generalized (see Ref. Urban, ).

III Numerical results

In this section, we first apply our formalism to a single QD system which is assumed to be half occupied at $t=0$. The dependence of cumulants on the occupation will be examined later. The linewidth amplitude in Eq. (40) is set to be $\Gamma_{L}=\Gamma_{R}=\Gamma/2$ and the bandwidth $W$ is also set to be the same for both leads. The energy level of the QD is assumed to be $5\Gamma$ and a bias with amplitude $\Delta_{L}=10\Gamma$ is chosen for this system. In the following numerical calculations, we set $e=\hbar=\Gamma=1$ for simplicity.

Figure 1 shows the 1st to 4th cumulants of transferred energy counted from the initial time $t=0$ in the left lead of the system for different bandwidths $W=10\Gamma$, $20\Gamma$, $50\Gamma$, and $80\Gamma$. For the 1st and 3rd cumulants, they decrease immediately once the system turns on and increase after reaching a minimum. In this region, the 1st and 3rd cumulants with smaller bandwidths have larger value until the crossover occurs at the time around $0.45$ and $2.5$, as shown in Fig. 1(a) and 1(c), respectively. After the crossover, situation reverses, i.e., the 1st and 3rd cumulants with the larger bandwidth $W$ becomes smaller than those with smaller $W$ in the self-energy. From Fig. 1(b) and 1(d) we see that the 2nd and 4th cumulants show a sharp rise first when the system turns on and then increases almost linearly after the transient regime. Roughly speaking the larger the bandwidth $W$, the larger the values of the 2nd and 4th cumulants.

We find that all cumulants of transferred energy increase with small oscillations when the time increases and exhibit linear characteristics at the long time which agrees with the long-time limit of the cumulant generating function in Eq. (36). The small oscillation of cumulants with time can be seen clearly from their derivative with respective of time, as shown in Fig. 2. Fig. 2(a) presents the time derivative of the 1st cumulant, namely, the transient energy current of the left lead for different bandwidths $W$. We see that they first drop down exhibiting dips with negative value once the system is connected. After that they increase to maximum values and then decay in an oscillatory fashion. In the long-time limit, they reach the values of dc energy current. We note that the transient energy current behaves similarly to the transient charge current obtained in Ref. joseph, . For the time derivative of the 2nd cumulant which is related to the shot noise of the system, it increases immediately exhibiting peaks when the system turns on, and drops down to the long-time limit very quickly with tiny oscillations. The time derivative of the 4th cumulant shows similar behavior to that of the 2nd cummulant besides that it drops to negative values and finally approaches to the positive long-time limit, while that of the 3rd cumulant shows opposite behavior. The short time behavior of $n$th cumulants of energy transport $C_{n}$ in Fig. 2 can be qualitatively understood from Eq. (33) where $C_{n}$ depends on $W$ through $\Gamma(E)\sim W^{2}/(E^{2}+W^{2})$. Hence at short times, a large $W$ gives a large $|C_{n}|$. The sign of $C_{n}$ can also be understood from Eq. (32) for $t$ approaching zero. For $n=2k$, we have from Eq. (30) and (32)

[TABLE]

which is positive definite. For $n=2k+1$, since $\int_{-\Delta_{L}}^{\Delta_{L}}dEE^{2k+1}=0$, we have

[TABLE]

which is negative in agreement with the results of Fig. 1 and Fig. 2.

From the numerical results we find that the frequency of the oscillations is independent of the bandwidth of leads while for the transient energy current, it is found that the time-dependent energy current with larger bandwidth $W$ decays faster than that with smaller $W$. This oscillatory behavior can be understood analytically. For a QD under an upward pulse of bias within the WBL the transient energy current can be expressed in terms of the spectral function $A(\epsilon,t)$ as,

[TABLE]

with

[TABLE]

Clearly the oscillatory behavior of transient energy current is due to the oscillatory term $\exp[i(\epsilon-\epsilon_{0}+\Delta_{L})t-(\Gamma/2)t]$ in the spectral function $A(\epsilon,t)$. We note that the period of oscillation of the transient energy current is only dependent on the energy level $\epsilon_{0}$ of the QD and the applied bias. In addition, the damping of this oscillation is dominated by life time of the resonant state of the QD which is proportional to $\Gamma$ in the WBL. In our numerical calculation, finite bandwidth $W$ is used which affects the lifetime of the resonant state. To find the influence of $W$ on the lifetime, we investigate transmission coefficient versus energy to look for resonant behavior and find the resonant state (transmission peak) and its lifetime (inverse of peak width). From the inset of Fig. 2(d), we find that the transmission peak is mediated by the resonant state of the single QD system is broadened by increasing the bandwidth $W$ in self-energy. This shows that the life time of the resonant state is proportional to bandwidth and explains why the transient energy current decays faster for the system with larger bandwidth (see Fig. 2(a)). From the transmission coefficient shown in the inset of Fig. 2(d), it is clear that larger $W$ corresponds to a large dc charge or energy current which is consistent with our dc limit of transient energy current.

From the discussion in Section II, we see that there may exist a universal behavior for $C_{n}$ at short times. In this section, we provide numerical evidence to show that this is indeed the case. Denoting $M_{n}$ as the maximum amplitude of the normalized transient energy cumulants $C_{n}(t)/C_{1}(t)$. In Fig. 3, we show the logarithmic plot of $M_{n}$ versus $n$ for different system parameters. We see that both $\ln(M_{2k})$ and $\ln(M_{2k+1})$ depend linearly on $k$ with the same slope $\kappa$ but different intercepts. Varying system parameters such as $W$ and $\epsilon_{0}$ will change the intercepts while the slope remains unchanged. Hence we have $M_{2k}=a_{1}e^{\kappa k}$ and $M_{2k+1}=a_{2}e^{\kappa k}$ where $a_{1}$ and $a_{2}$ are constants. The universal slope $\kappa$ is found to be close to 3.

In order to study the effect of initial occupation number of the single QD on the transient energy current, the time-dependent energy currents calculated by the time-derivative of 1st cumulant with $W=50\Gamma$ for different initial occupation number are plotted in Fig. 4(a). To understand this behavior, we note that in the transient regime the lesser Green’s function can be expressed by zwm ; leizhang ,

[TABLE]

Substituting this expression into Eq. (23) we find the transient energy current consists of two terms $I^{E}_{L}(t)=I^{E}_{L,un}(t)+I^{E}_{L,in}$ where $I^{E}_{L,un}(t)$ is the transient energy current for a system which is initially unoccupied while the transient energy current due to the initial occupation is

[TABLE]

In Fig. 4(a), we plot $I^{E}_{L,in}(t)$ for the single QD system which is initially full occupied (defined as $I^{E}_{L0}$), namely, $g^{<}(0,0)=1i$. Therefore, the transient energy current for a single dot system with initial occupation of $\alpha$ is $I^{E}_{L}(t)=I^{E}_{L,un}(t)+\alpha I^{E}_{L0}$. We have checked that three curves (black, red, and blue lines) in Fig. 4(a) indeed satisfy this relation.

As a second example, we consider a double QD system with the Hamiltonian of \left(\begin{array}[]{cc}\epsilon_{1}&-v\\ -v^{*}&\epsilon_{2}\\ \end{array}\right). We set the energy levels of the first and second QD to be $\epsilon_{1}=4\Gamma$ and $\epsilon_{2}=6\Gamma$ which connect with the left and right lead, respectively. The bandwidth of the lead and the bias are set to be $W=10\Gamma$ and $\Delta_{L}=10\Gamma$, respectively. The bias is applied to the isolated leads at $t=-\infty$ and the couplings between the leads and QDs are switched on at $t=0$. For the double QD system, we assume that the initial state for $\epsilon_{1}$ is occupied initially with the lesser Green’s function given in Eq. (38), while $\epsilon_{2}$ is unoccupied so that its occupation number (proportional to its lesser Green’s function) is 0.

Figure 5 presents the 1st to 4th cumulants of transferred energy of electrons counted in the left lead with different coupling constants of $v=1\Gamma$, $2\Gamma$, $3\Gamma$, and $5\Gamma$ between two QDs. Generally speaking, as the time increases, all calculated cumulants increase with oscillation except for the 3rd cumulant with $v=5\Gamma$ which starts to decrease slowly when $t$ exceeds 4.5. The oscillations of the 3rd and 4th cumulants with a large interdot coupling constant $v=5\Gamma$ decay much more slowly than those of other cases. Similar to the cumulants for the single QD system, the long-time limit of different cumulants for the double QD system is a linear function of time as shown in Eq. (36). It is also found that in the long-time limit, the $n$th cumulant with coupling constant of $v=3\Gamma$ becomes the largest while $v=\Gamma$ is the smallest compared with those with other coupling constants for all $n$. We also note that the behaivor of cumulants of transferred energy for the double QD system resembles that of cumulants of transferred charge reported in Ref. T2, .

The time derivative of the cumulants for the double QD system are then presented in Fig. 6. The time derivative of all cumulants show oscillations with their amplitudes decreasing with the increase of time and their frequencies approximately proportional to the coupling constant, especially for large coupling constant. To understand the behavior of oscillation, we note that in the resonance regime the electron oscillates between two QDs with different energy levels and coupling constants, leading to different oscillation frequencies related by the Rabi frequency defined as,

[TABLE]

which is actually the difference between two eigenvalues of the Hamiltonian for the double QD system. In our case, the period of oscillation for the transient energy current is $T=\frac{2\pi}{\omega}=\frac{\pi}{\sqrt{1+v^{2}}}$. Therefore, higher frequency is obtained for the time derivative of cumulants plotted in Fig. 6 for system with larger coupling constant. Specifically, for the transient energy current, namely, the time derivative of the 1st cumulant, it is found that for the systems with coupling constants of $v=1\Gamma$, $2\Gamma$, and $3\Gamma$, larger energy current is obtained for system with larger coupling constant in the long-time limit since the transmission peaks become higher and wider when the coupling constant is increased , as shown in the inset of Fig. 6(d). However, when the coupling constant is further increased to $5\Gamma$, the transmission peaks start to shift out of the bias windows [0,10] which dramatically reduces the energy current in the long-time limit. For higher order cumulants, we observe similar behaviors (which can be understood using the augments discussed above): (1). the larger the interdot coupling, the larger the oscillation frequency becomes; (2). in the long-time limit, the value of cumulants increases with $v$ as long as $v<4.5\Gamma$. (3). For $v>4.5\Gamma$, the value of cumulants is the smallest in the long-time limit.

Moreover, the transient energy currents calculated by the time-derivative of 1st cumulant of the double QD with $v=2\Gamma$ for different initial occupation condition are plotted in Fig. 4(b). Similar to the case of single QD, we also calculate the contribution of transient current from the initial occupation condition by Eq. (48) for the case that only $\epsilon_{1}$ ($I^{E}_{L1}(t)$) and $\epsilon_{2}$ ($I^{E}_{L2}(t)$) is fully occupied, respectively, as shown in Fig. 4(b). Therefore, for a double QD system in which the occupation number of $\epsilon_{1}$ and $\epsilon_{2}$ are $\alpha$ and $\beta$, respectively, the transient current can be calculated simply by $I^{E}_{L}(t)=I^{E}_{L,un}(t)+\alpha I^{E}_{L1}+\beta I^{E}_{L2}$.

IV Conclusion

We have investigated the FCS of transferred energy in the transient regime. Two time measurement scheme was used to derive the generating function of FCS of transferred energy in the transient regime using the Keldysh non-equilibrium Green’s function. Our formalism was then applied to both single and double QD systems to study the 1st to 4th order cumulants of transferred energy in the transient regime. Oscillations are observed in the transient energy current for both single and double QD systems. At short times, universal scaling was found for maximum amplitude of normalized cumulant of energy current for the single QD system. For the single QD system, we find that the frequency of oscillation in the transient energy current is independent of the bandwidth of the self-energy while and for the double QD system the frequency is proportional to the coupling constant between two QDs for large coupling constant.

Acknowledgements.

This work was financially supported by the Research Grant Council (Grant No. HKU 705212P), the University Grant Council (Contract No. AoE/P-04/08) of the Government of HKSAR, NSF-China under Grant No. 11374246.

Appendix A Transient energy current

In this appendix, we present the details of obtaining the transient energy current from the cumulant generating function. Accoding Eqs. (2) and  (22), the transient energy current, namely, the 1st cumulant of transfered energy, can be written as ($\hbar=e=1$),

[TABLE]

where,

[TABLE]

and

[TABLE]

from the definition of $\widetilde{\Sigma}_{L}$ in the Keldysh space in Eq. (II). For the lesser self-energy $\widetilde{\Sigma}^{<}_{L}$ with the counting field, it can be expressed in terms of the lesser green’s function of left lead,

[TABLE]

where,

[TABLE]

with $\epsilon^{0}_{kL}$ the bare energy levels of the left lead in the absence of external bias.

It is easy to find

[TABLE]

from which we have,

[TABLE]

where

[TABLE]

with $\chi=<,a$ and $\Sigma^{\chi}_{kL}(t^{\prime}-t)$ being the self-energy of the left lead without the counting field. Similarly, we can obtain,

[TABLE]

By substituting Eqs. (56) and (58) into Eq. (52), we have,

[TABLE]

Therefore, the transient energy current can be obtained from Eq. (50),

[TABLE]

Since $G^{k}=2G^{<}+G^{r}-G^{a}$, it can be finally written as,

[TABLE]

For the case of heat current, the cumulant generating function for the counting of the heat is found to be,

[TABLE]

with

[TABLE]

where ${\tilde{\Sigma}}_{L}$ is defined in Eq. (II).

Therefore, the transient heat current can be expressed as,

[TABLE]

which is the known result. We also wish to comment that this formalism is valid for heat transport driven by either temperature gradient and bias voltage.