Full-counting statistics of energy transport of molecular junctions in the polaronic regime

We investigate the full-counting statistics (FCS) of energy transport carried by electrons in molecular junctions for the Anderson-Holstein model in the polaronic regime. Using two-time quantum measurement scheme, generating function (GF) for the energy transport is derived and expressed as a Fredholm determinant in terms of Keldysh nonequilibrium Green's function in the time domain. Dressed tunneling approximation is used in decoupling the phonon cloud operator in the polaronic regime. This formalism enables us to analyze the time evolution of energy transport dynamics after a sudden switch-on of the coupling between the dot and the leads towards the stationary state. The steady state energy current cumulant GF in the long time limit is obtained in the energy domain as well. Universal relations for steady state energy current FCS are derived under finite temperature gradient with zero bias and this enables us to express the equilibrium energy current cumulant by a linear combination of lower order cumulants. Behaviors of energy current cumulants in steady state under temperature gradient and external bias are numerically studied and explained. Transient dynamics of energy current cumulants is numerically calculated and analyzed. The universal scaling of normalized transient energy cumulants is found under both temperature gradient and external bias.

We investigate the full-counting statistics (FCS) of energy transport carried by electrons in molecular junctions for the Anderson-Holstein model in the polaronic regime. Using two-time quantum measurement scheme, generating function (GF) for the energy transport is derived and expressed as a Fredholm determinant in terms of Keldysh nonequilibrium Green's function in the time domain. Dressed tunneling approximation is used in decoupling the phonon cloud operator in the polaronic regime. This formalism enables us to analyze the time evolution of energy transport dynamics after a sudden switch-on of the coupling between the dot and the leads towards the stationary state. The steady state energy current cumulant GF in the long time limit is obtained in the energy domain as well. Universal relations for steady state energy current FCS are derived under finite temperature gradient with zero bias and this enables us to express the equilibrium energy current cumulant by a linear combination of lower order cumulants. Behaviors of energy current cumulants in steady state under temperature gradient and external bias are numerically studied and explained. Transient dynamics of energy current cumulants is numerically calculated and analyzed. The universal scaling of normalized transient energy cumulants is found under both temperature gradient and external bias.

I. INTRODUCTION
Rapid experimental development in the field of nanotechnology makes fabrication of the single-molecule junctions possible [1,2], which could push the limit of Moore's law further. In the electronic quantum transport though nanodevices, the electron-phonon coupling plays an important role. One of the mechanisms that induces electron-phonon coupling is due to the charging of the molecule leading to elastic mechanical deformations. This in turn causes an interaction between electronic and the quantized mechanical degrees of freedom giving rise to the electron-phonon coupling. A variety of intriguing transport properties, such as phonon-assisted current steps and Franck-Condon blockade [5] have been found in the polaronic regime [3,4] when this kind of electron-phonon coupling in molecular junctions is strong. Theoretically, these phenomena could be understood using a quantum dot described by the Anderson-Holstein model [6,7] coupled to two electrodes.
To understand quantum transport in the polaronic regime, many methods have been used, such as the master equation method [8][9][10][11], diagrammatic quantum Monte Carlo method [12], numerical renormalization group method [13], as well as the nonequilibrium Green's function (NEGF) technique [20] that is particular useful in describing time dependent non-equilibrium processes. Perturbation method is applicable when the electron-phonon coupling strength is weak [14][15][16] and it fails in the strong electron-phonon coupling system. Other approximation has to be made in order to deal with the strong and intrinsically nonlinear electron-phonon interaction in the Anderson-Holstein model. In order to decouple the phonon cloud operator in the polaronic regime, dressed tunneling approximation (DTA), in which the leads' self-energies are dressed with the polaronic cloud, has been proposed to eliminate the noticeable pathological features of the single particle approximation (SPA) at low frequencies and polaron tunneling approximation (PTA) at high frequencies [17][18][19][20].
The key in FCS is to obtain the generating function (GF) which is actually the Fourier transform of probability distribution of the related physical quantity. Using the NEGF technique [37][38][39][40] and the path integral method under the two-time quantum measurement scheme [27,[41][42][43], GF was formulated as a Fredholm determinant in the time domain for both phonon [29][30][31] and electron [27,[32][33][34][35] transport. This formalism enables one to study the transport properties in the transient regime providing more information on the short time dynamics [32]. Recently transient dynamics of particle current transport in the molecular junctions has been studied by Schmidt et al. [44,45] in the case of weak and strong electron-phonon couplings and has been reported by Maier et al. using PTA [46] and by Souto et al. using DTA [20] in the polaronic regime.
The transport study of energy flow in the nonequilibrium system could reveal information on how energy is dissipated and its correlation for electronic devices and can be investigated theoretically by Landauer-Büttiker type of formalism for noninteracting systems [47][48][49]. Energy transport in trapped ion chains has been measured experimentally by Ramm et al. [50]. The heat current I h α in the α lead is related to the energy current I E α by the expression I h α = I E α −µ α I α with the particle current I α and the chemical potential µ α in the α lead, and this quantity is quite important in characterizing the efficiency of thermoelectric devices [51]. So far, FCS of energy transfer mostly focuses on phonon transport both in the transient regime and steady states [29][30][31] and less attention has been paid to the FCS of energy transfer carried by electrons in the electronic transport problems. In our previous work, we investigated the transient FCS of energy transfer in the non-interacting system [34]. It would be important and interesting to study of FCS of energy transport carried by electrons of molecular junctions with electron-phonon coupling in the polaronic regime for both transient dynamics and steady states, and this is the purpose of this work.
In this paper, FCS of energy transport carried by electrons in molecular junctions for the Anderson-Holstein model in the polaronic regime is investigated both in the steady states and transient regime. Within the DTA, GF for the energy current is derived from the equation of motion and could be expressed as a Fredholm determinant in the time domain using NEGF. Numerical calculation is performed which allows us to analyze the time evolution of the energy flow towards the steady state for a sudden switch-on of the coupling between the quantum dot and the leads. The cumulant GF of energy current in the steady state is obtained analytically in the energy domain. Universal relations for cumulants of energy current under finite temperature gradient with zero voltage bias are established. In addition, we also calculate and analyze steady state solution for various order of cumulants (from the first to the fourth order) under temperature gradient or external bias.
The rest of the paper is organized as follows. In Sec. II, the model Hamiltonian of a molecular junction is introduced and GF of energy flow in the transient regime is determined in terms of NEGF in the time domain. Sec. III is devoted to the steady state investigation of FCS of energy current, both theoretically and numerically. In Sec. IV, transient dynamics of energy current is investigated under a sudden switching-on of external bias. Finally, a brief conclusion is drawn in Sec. V.

II. MODEL AND BASIC THEORETICAL FORMALISM
Considering only the lowest electronic orbital, the single-molecule is simplified as a single electronic level of a quantum dot (QD) being coupled to localized vibrational mode, which is the simplest spinless Anderson-Holstein model [52]. The QD then is coupled to the left and right electrode so that the system is driven to a nonequilibrium state when the external bias or temperature gradient is applied (Fig. 1). The corresponding Hamiltonian reads as with the Hamiltonian of the central dot (in natural units, = k B = e = m e = 1) where ǫ 0 is the bare electronic energy level, and ω 0 is the frequency of the localized vibron. d † (a † ) denotes the electron (phonon) creation operator in the QD. The localized vibron modulates the QD with the electron-phonon coupling constant t ep . The Hamiltonians of the leads is given in a compact form where the indices kα = kL, kR are used to label the different states in the left and right leads. H T is the Hamiltonian describing the coupling between the dot and the leads with the tunneling amplitudes t kα , The tunneling rate (linewidth function) of lead α is assumed to bear the Lorentzian form and could be expressed as with the linewidth amplitude Γ α and bandwidth W , and one can denote Γ = Γ L + Γ R . The electron-vibron coupling term can be eliminated by applying the Lang-Firsov unitary transformation [53] given byH which leads toH where the bare QD electron energy is changed toǭ = ǫ 0 − g 2 ω 0 . The tunneling Hamiltonian is transformed as with the phonon cloud operator X = exp[g(a − a † )], while Hamiltonians of isolated leads remain unchanged.
In the present work we study the transient dynamics in which the interaction between the leads and the QD is suddenly turned on at t = 0 and afterwards the system evolves to the steady states. The turning on process could be facilitated by a quantum point contact which is controlled by a gate voltage. The initial density matrix of the whole system at t = 0 is the direct product of each subsystem and expressed by ρ(0) = ρ L ⊗ρ S ⊗ρ R . The statistical behaviors of the energy current in a specific lead are all encoded in the probability distribution P (∆ǫ, t) of the transferred energy carried by electrons ∆ǫ = ǫ t − ǫ 0 between an initial time t = 0 and a later time t. The GF Z(λ, t) with the counting field λ is defined as, The kth cumulant of transferred energy (∆ǫ) k could be calculated by taking the kth derivative of cumulant generating function (CGF) which is ln Z(λ) with respect to iλ, One can further define the energy current cumulants which tend to the steady state energy current cumulants in the long time limit t → ∞. The second energy cumulant could be expressed as C 2 (t) = t 0 dt 1 t 0 dt 2 δI E (t 1 )δI E (t 2 ) , so that the second energy current cumulant is ( . One should note that the second energy current cumulant (I E ) 2 is not an average of a squared quantity. To investigate statistical behaviors of the energy current through the left lead, we could focus on the energy operator which is actually the free Hamiltonian of the left lead H L . Under the two-time measurement scheme, GF of transferred energy in the left lead can be expressed over the Keldysh contour as [27,31,33], with the modified evolution operator (γ = ±λ/2 depending on the branch of the contour, see Fig. 2), Here the modified evolution operator is expressed by the modified Hamiltonian, with t γ = γ, and c kL (t γ ) = e iγHL c kL (0)e −iγHL . GF for the transferred charges in transient regime has been expressed by NEGF in the time domain for the noninteracting case [33] and in the polaronic regime using the DTA [18,20]. GF for the energy current has expressed by NEGF and higher-order cumulants has been investigated by Yu et al. for the non-interacting case [34]. We now generalize the GF for the transferred energy to the interacting case in the polaronic regime following the derivation of the GF for transferred charges [20]. Following the procedure outlined in Ref. 54, one can get GF from the derivative of the logarithm of Eq. (12) with respect to the counting field, where we take ′ − ′ in the first part and ′ + ′ in the second for the forward time contour, while inversely for the backward contour (see Fig. 2). The average T C · · · denotes Tr ρ(0) The equation of motion of the three point Green function on the contour which could be written in an integral form [38] T Under DTA, one has the following decoupling [20] T with Λ(t, t 1 ) = T C X † (t 1 )X(t) being the phonon cloud propagator which will be discussed later. Then we have The self-energies due to the coupling to the leads under the DTA could be expressed as, where a, b = +, − denote different Keldysh components and Note that the counting field enters the self-energy in absence of the phonon cloud operator and the modified self-energy can be expressed by [34] where Tr K indicates the trace is over the Keldysh space. Using the fact that Z(λ = 0, t) = 1, the GF could be expressed in the Fredholm determinant by the Keldysh NEGF in the time domain as [19,20], with where G 0 denotes the Green's function of the uncoupled QD, and the tilde indicates the inclusion of the counting field in the self-energy Σ α,D . Note that the Green's functions and self-energies without counting field possess the Keldysh structure, The phonon cloud operator Λ ab δ (t 1 , t 2 ) that is coupled to lead δ = L, R is given by [7], with and I m being the modified Bessel function of the first kind, and Bose factor n Bδ = 1/(e β δ ω0 − 1), β δ = 1/k B T δ . We should mention that the temperature of the phonon cloud operator is dependent on which self-energy it multiplies with, and in the next section we will see that this will ensure the important fluctuation symmetry relation. In the work by Y. Utsumi et. al., a third thermal probe electrode due to the thermal bath was added to determine the temperature of the vibrations [55]. In our work, we only consider the energy flow carried by electrons, and the fluctuation symmetry relation is already satisfied for the two-terminal system in Eq. (43). At zero-temperature α m = α mL = α mR could be simplified as, The remaining components of Λ δ could be calculated by the relations, The Dyson equation bearing a Keldysh structure under DTA is where Utilizing the Dyson equation, Eq. (23) could be written as, so that CGF has the form, by using the relation det B = exp[Tr ln B]. Taking the first derivative of GF and noting that Σ +− t kL , energy current in the transient regime is found to be, whereΣ and we have similar definition forΣ −+ (t ′ , t). The transient current expression formally agrees with the one which was obtained directly by NEGF method [56].

III. STEADY STATE ENERGY TRANSPORT FCS
In the long-time limit, the system goes to steady state, and the Dyson equation Eq. (30) bearing the Keldysh structure in the energy domain could be expressed by so that [18] with The dressed retarded self-energy in frequency domain could be obtained by the Fourier transformation of the time domain counterpart with the form Σ r The real and imaginary part could be obtained using Plemelj formula 1/(E ± i0 + ) = P (1/E) ∓ iπδ(E) which will be used in the numerical calculation. One can verify that the real part and imaginary part satisfies respectively [18]. In the long-time limit, the Green's function and self-energy in Eq. (32) become time translation invariant so that scaled cumulant generating function (SCGF) F (λ) = lim t→∞ ln Z(λ)/t could be expressed in the energy domain as In this expression T mn (ω) is the transmission coefficient involving m and n vibrational quanta in the left and right lead, respectively, with the form, Taking the first order derivative of SCGF with respect to λ, we can get the expression of energy current, Now we consider the universal relations for energy current cumulants under finite temperature gradient with zero bias which is in analogy with the universal relation for particle current cumulants [57,58]. Using the relation α −m = e −βLmω0 α m , α −n = e −βRnω0 α n and f R (1 − f L ) = exp(∆βω)f L (1 − f R ) with ∆β = β L − β R for ∆µ = 0 in Eq. (40), we have the fluctuation symmetry relation with iλ being replaced by ξ for convenience. One can verify that the fluctuation symmetry can only be satisfied by considering the dependency of phonon temperature with respect to the specific lead. In the linear response regime ∆β → 0, we can expand both sides as Taylor series around ∆β = 0 and ξ = 0, which leads to, where we have written the dependence of ∆β of SCGF explicitly out in both sides. Since F (ξ = 0, ∆β) = 0, Eq. (43) gives F (∆β, ∆β) = 0, from which we find that the LHS of Eq.(44) vanishes. The last term in the summation of Eq. (44) is the kth derivative of the SCGF with respect to the counting field ξ, which is actually (I E ) k at equilibrium. Then we have the relation in which the energy current cumulant at equilibrium is expressed by a linear combination of lower order energy current cumulants. This is similar to the case that the particle current cumulant could could be expressed by a linear combination of lower order particle current cumulants in the presence of small voltage bias [57,58]. We now show numerical calculations regarding steady state energy current cumulants under temperature gradient and external bias of molecular junctions in the polaronic regime. The energies are measured in the unit of vibron frequency ω 0 , and the linewidth amplitude is chosen to be Γ = 0.05ω 0 which indicates weak coupling. In addition, WBL is taken in our steady state calculation.
The first to fourth energy current cumulants for increasing g versus temperature gradient ∆T = T L − T R with the left lead warmer and temperature of right lead fixed at k B T R = 0.2ω 0 are shown in Fig. (3). The chemical potentials in both leads are set to be zero and the renormalized energy level of the QD isǭ = 0. The energy current cumulants become smaller with the increasing of g because of the suppression of transport due to electron-phonon interaction. The second energy current cumulant with zero temperature gradient is finite due to the thermal noise in the leads, and it is reduced with increasing g.
In Fig. (4), energy current cumulants with different renormalized energy levels of the QD with g = 1 are plotted. We can see that the first to fourth cumulants and SCGF as well are even functions ofǭ. This can be understood as follows. Since the chemical potentials of both leads are zero, one can set and verify that, using the relation f L+m (ω) = 1 − f L−m (−ω). In the WBL, from Eq. (39), the real and imaginary part of the dressed retarded self-energy are the odd and even function of ω, respectively, so that we have the following symmetry with respect to the transmission coefficient in the polaronic regime where the dependency ofǭ has been written explicitly. Then, we have the following symmetry of SCGF with respect toǭ, with µ L = µ R = 0 in the WBL. One can also see from Fig. (4), I E (ǭ = 2ω 0 ) is smaller than I E (ǭ = ω 0 ) under small temperature gradient, and this is also for the second energy current cumulant. Since the linewidth amplitude Γ = 0.05ω 0 is small, so that the transmission coefficient which is centered aroundǭ is narrow. As a result, the main contribution to the transport process is coming from energy nearǭ. When the temperature gradient across the junction is small, the difference of Fermi distribution functions between left and right lead f L (ω) − f R (ω) is smaller nearǭ = 2ω 0 than nearǭ = ω 0 . When T L increases, the difference of Fermi distribution functions between left and right lead f L (ω) − f R (ω) near ω = 2ω 0 could exceed the difference near ω = ω 0 , so that the first and second cumulant with largerǭ is larger than the ones with smallerǭ. The first to fourth energy current cumulants for increasing g versus external bias ∆µ with µ L = ∆µ/2 and µ R = −∆µ/2 are shown in Fig. (5). Temperatures of both leads are chosen to be very small with k B T L = k B T R = 0.04ω 0 which is almost in the regime of zero temperature. The renormalized energy level of the QD isǭ = 2ω 0 . For the non-interacting case, the energy current and second cumulant are almost zero when bias is below ∆µ = 4ω 0 = 2ǭ and display plateau structures when the external bias exceeds 2ǭ. The width of transmission coefficient is small due to the small linewidth amplitude Γ = 0.05ω 0 . When ∆µ = 2ǭ, chemical potential of the left lead is equal to the renormalized energy of QD, µ L =ǭ, in which energy the transmission coefficient experiences a sharp increase and reaches its largest value as indicated in Fig. 1(b). From the Fig. 5, we observe electron-phonon coupling enables the plateau height to become smaller, however creates smaller steps at ∆µ = 2ǭ + 2nω 0 with n = 1, 2, 3 · · · . This is due to the presence of sidebands in the leads and could be understood as follows. In the presence of the polaronic regime, from Eq. (42), we can approximately write the energy current in presence of bias voltage at zero temperature as, (ignore the terms with product of Fermi distribution function) The energy current is written as a sum of a series, with each term coming from a different sideband in the leads. The first plateau of the energy current in the polaronic regime is mainly due to the first term in Eq. (50), and the second plateau due to the contribution from the second term in Eq. (50) with one polaron involved in the transport process, and etc.. We find that T m /T m−1 = g 2 /m is responsible for the ratios between plateau heights. One can see that when g = 0.5, T 1 /T 0 = 0.25, so that the height of the second plateau is a quarter of that of the first plateau at zero temperature, which explained what we see in Fig. 5. This is also applicable to the case g = 1.0 with T 1 /T 0 = 1.0 and the case g = 1.5 with T 1 /T 0 = 2.25. One should note that the temperature of the system in Fig. 5 is very small.
The plateau structures disappear in the third and fourth energy current cumulants. Instead a dip occurs at ∆µ = 2ǭ for both the third and fourth energy current cumulants with fourth cumulant larger for both non-interacting and interacting cases. Polaronic regime creates smaller dips at ∆µ = 2ǭ + 2nω 0 with n = 1, 2, 3 · · · which could also be identified in Fig. 7. Increasing g reduces the amplitude of the dip at ∆µ = 2ǭ but increases the amplitude at ∆µ = 2ǭ + 2nω 0 . The explanation is as follows. For the non-interacting case under zero-temperature, we have with the ranges of integration from −∆µ/2 to ∆µ/2. We can further take derivative of (I E ) k with respect to external bias ∆µ, ∂ I E /∂∆µ and ∂ (I E ) 2 /∂∆µ is always positive definite since the transmission coefficient for non-interacting case has the form T (ω) = (ω−ǭ) 2 +Γ 2 /4 . However the derivative of the third and fourth cumulant with respect to external bias change sign around ∆µ = 2ǭ and also the transmission coefficient experiences an abrupt change because of small linewidth amplitude. This leads to the the dips of third and forth cumulant of energy current as shown in Fig. 5.
The influence of temperature on cumulants under external bias is depicted in Fig. 6, one can see that both the plateaus and dips get smoothed or even disappeared when temperature increases. In Fig. 7, energy current cumulants with differentǭ with g = 1 are plotted. We can see that the first and third cumulants are odd functions ofǭ, while the second and fourth cumulants are even functions ofǭ. The reason is as follows. Under zero temperature, the transport is unidirectional and Fermi-Dirac distribution function f L(R) has a step-wise form, since the transmission coefficient is peaked around the resonant levelǭ with a very small linewidth amplitude (say δǫ), the energies of electron which mainly contribute to the energy transport are very close toǭ. So if we change the sign ofǭ from positive to negative, then most of electron energies will reverse their signs if δǫ <ǭ. Since the energy current is proportional to energies of electron, this will lead to the energy current reversal.

IV. TRANSIENT DYNAMICS OF ENERGY TRANSPORT
We first investigate the behaviors of energy current at very short time. To do that, we expand the GF to the lowest order in time, The expressions of Green's function for isolated QD and self-energy are given in the Appendix. Under the wide-band limit W → ∞, we can obtain the GF in the short time limit in a compact form as, where with From now on we use f α±n to denote f α (ω ± nω 0 ). We can see from the expression of short time limit of the GF that the transport process is unidirectional in the short time limit. We can get the current expressions in the short time limit as We apply the formalism to perform numerical calculation with respect to the transient dynamics of energy current under temperature gradient and external bias, respectively. The energies are measured in the unit of Γ and 1/Γ is the unit of time. We only consider the case where the QD is initially unoccupied n d = 0 and the linewidth amplitude in Eq. (5) is set to be Γ L = Γ R = Γ/2 and the bandwidth is also set to be the same for both leads with W = 10Γ.
First to fourth transient energy current cumulants, (I E ) k for k = 1, 2, 3, 4, in the left lead for increasing g under temperature gradient and external bias are shown, respectively, in Fig. 8 and in Fig. 9. Increasing g corresponds to the increasing of the electron-phonon coupling strength. The frequency of the localized vibron is ω 0 = 6Γ. The renormalized energy level of the QD isǭ = 2Γ for case under temperature gradient andǭ = 1.5Γ for the case with external bias. The left lead is assumed to be warmer with the temperatures of the two leads to be k B T L = 1.5Γ As a general feature for both the non-interacting (g = 0) and interacting cases, the transient amplitudes of (I E ) k increase with cumulants order. This behavior is universal and will be investigated in detail in Fig. 10. The second and fourth energy current cumulants may even oscillate to negative values at short times. The negativity of the second energy current cumulants can be explained as follows. The energy cumulant C 2 (t) must be positive at all times from a statistical view, however it can oscillate at short times so that the second energy current cumulants which is the derivative of C 2 (t) may not be positive at short times. (I E ) 2 at steady state (long time limit) is positive and can be identified from the figures. The amplitudes of oscillation in the evolution and the asymptotic values of the cumulants are suppressed with the increasing of g. The first and third energy current cumulants in the stationary limit are positive under temperature and external bias, since we put the normalized energy level of QD above the Fermi energy of the both leads so that the electrons with positive energy contribute to the transport process. However, in short times, the energy current and third cumulant oscillate to negative values with a minimum. This could be understood as follows. Since the QD is prepared initially empty, once the system is connected, the contribution to the transport process mainly comes from electron of the left lead which could be seen from Eq. (53). The contribution of energy current cumulants from the energy window [0, µ L ] cancels with the contribution from [−µ L , 0], so that energy below −µ L in the left lead will contribute to the energy transport process which leads to the negativity of the first and third energy current cumulants in the short times. The cumulants of transient energy current approach to their steady state values in the long time limit.
We also plot the logarithm of maximum amplitude of the normalized transient energy cumulants M k = max|C k /C 1 | under temperature gradient [ Fig. 10(a)] and external bias [ Fig. 10(b)]. Different lines with respect to different bandwidths W are plotted, while the other parameters are same as in Fig. 8 and Fig. 9. Maximum amplitudes M k for different interaction parameter g = 0, 1.0 and 1.5 coincide. We can see from the figure that both ln(M 2k ) and ln(M 2k+1 ) are linear with cumulants order k with the slope close to 3 but they have different intercepts. This universal scaling of normalized transient energy cumulants is found under both the temperature gradient and external bias, and it is the result of the universality of the GF in the short time which was also reported in the charge cumulants [20,59]. Theoretical understanding of this behavior for the noninteracting case was reported in our previous work [34]. Interestingly, turning on the electron-phonon interaction does not affect this behavior.

V. CONCLUSION
Both steady state and transient behaviors of energy transport carried by electrons in molecular junctions for the Anderson-Holstein model in the polaronic regime have been investigated using FCS. Using two-time measurement scheme and equation of motion technique, GF for the energy current could be expressed as a Fredholm determinant in the time domain using NEGF. The DTA decoupling scheme [17] which could provide a good description in dealing with the phonon cloud operator has been adapted in obtaining GF. This formalism allows us to analyze the time evolution of energy transport dynamics after a sudden switch of the coupling between the dot and the leads towards the stationary state. The amplitudes of oscillation in the evolution and the asymptotic values of the cumulants are suppressed with the increasing of g. The universal scaling of normalized transient energy cumulants is found under external bias.