Transient charge and energy flow in the wide-band limit
F. Covito F. G. Eich R. Tuovinen M. A. Sentef A. Rubio
cond-mat.mes-hall







[]florian.eich@mpsd.mpg.de




Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany


Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany
Center for Free-Electron Laser  Science, Luruper Chaussee 149, 22761 Hamburg, Germany
Center for Computational Quantum Physics (CCQ), The Flatiron Institute, 162 Fifth avenue, New York NY 10010
Nano-Bio Spectroscopy Group, Departamento de Fisica de Materiales, Universidad del País Vasco, 20018 San Sebastián, Spain

 
  The wide-band limit is a commonly used approximation to analyze transport through nanoscale devices. In this work we investigate its applicability to the study of charge and heat transport through molecular break junctions exposed to voltage biases and temperature gradients. We find that while this approximation faithfully describes the long-time charge and heat transport, it fails to characterize the short-time behavior of the junction. In particular, we find that the charge current flowing through the device shows a discontinuity when a temperature gradient is applied, while the energy flow is discontinuous when a voltage bias is switched on and even diverges when the junction is exposed to both a temperature gradient and a voltage bias. We provide an explanation for this pathological behavior and propose two possible solutions to this problem.


Transient charge and energy flow in the wide-band limit
    A. Rubio
    
=======================================================





§ INTRODUCTION
 
Over the last decades great effort has been spent to miniaturize electric circuits. The goal is to realize the fundamental building blocks of electronic circuits, such as transistors, on the scale of single molecules. There has been a great success in shrinking electronic devices down experimentally. In order to understand the properties of molecular break junctions a quantum mechanical description of the device is required. Perhaps the most successful and wide-spread theory to describe how charge flows through a nanoscale junction is the so-called Landauer-Büttiker approach <cit.>, which describes the charge transport as a scattering problem. Essentially, the flow of charge through a molecular junction is determined by the transmission function of the device–describing how impinging electrons are scattered–and the occupation function of the electrons in the (metallic) leads connected to the junction. 

In recent years there has been renewed interest in addressing not only the charge flow, but also the energy (or heat) flow through nanoscale devices. Understanding how charge and energy flow depend on voltage and temperature biases across the device provides crucial insight for the development of thermoelectric circuits, which could be used to convert waste heat into useful electric energy <cit.>. Furthermore, recent experiments demonstrate that local temperatures in nanoscale conductors can be measured with a spatial resolution of tens of nanometers <cit.>. A common path to address the effect of temperature gradients across the nanoscale device is to allow for different temperatures in the occupation functions characterizing the leads in the Landauer-Büttiker formula. Conceptually this can only be justified if the leads are considered to be disconnected from the device initially (partitioned approach). This artificial partitioning of the system, however, is problematic, for it assumes that it is possible to perfectly decouple the leads from the molecular junction–a rather optimistic assumption if one considers atomic-scale devices. For times much larger than the typical time-scale of molecular break junctions, which are on the order of femtoseconds <cit.>, the assumption on a decoupled initial state does not play a crucial role. However, for transient dynamics the initial state matters. As pump-probe experiments are now able to investigate phenomena happening at this timescale <cit.> it is important to properly describe the initial state.

An alternative to the partitioned approach is to couple the device and leads at all times and trigger the charge flow by switching a potential bias <cit.>. This partition-free approach leads to the same steady state as the partitioned approach, but the transient dynamics of the device will, in general, be different <cit.>. The advantage of the partition-free approach is that the transient charge and energy/heat flows are not spoiled by the dynamics induced by connecting leads and device, because the electronic states in the device are allowed to hybridize with the leads before any temperature or voltage bias is applied. Importantly, it is also possible to take into account temperature differences in the leads within the partition–free approach: We consider a thermo-mechanical potential, which couples to the local energy density of the system–much like the usual electric potential couples to the charge density <cit.>. This thermo-mechanical potential acts as mechanical “proxy” <cit.> for local temperature variations. An intuitive way to understand this is to consider the occupation function, which is determined from the ratio of the energy and the temperature. Accordingly, a change in occupations due to a change in temperature can alternatively be viewed as a change in energy keeping the temperature fixed. The thermo-mechanical potential allows to rescale the energy locally, thereby mimicking a locally varying temperature. Applying this idea in the context of transport means that different temperatures in the leads are described by rescaling the bandwidth of the leads <cit.>.

A wide-spread simplification used to describe transport through nano junctions is the so-called wide-band limit (WBL). The WBL assumes that the detailed structure of the density of states in the leads is not important for the description of transport, which substantially simplifies computations. The WBL for charge transport is justified when the bandwidth is large compared to the applied bias <cit.>.

In this work we investigate whether the WBL can be employed in conjunction with the thermo-mechanical potential. An immediate question that comes to mind is: What is the meaning of rescaling an infinite band? In the following we will show that the steady state is well described in the WBL, provided the WBL is taken properly. The transient currents, however, exhibit peculiarities at short times. Specifically, we see that the charge current jumps at the initial time when the device is exposed to a temperature gradient and, similarly, the heat current behaves discontinuously when a voltage bias, but no temperature bias, is switched on. Even more dramatically, the heat current diverges as (t-t_0)^-1, with t_0 being the time at which a temperature and charge bias is applied to the system. By comparing the WBL transient charge and heat currents to results obtained at finite bandwidth, we highlight that this pathological behavior of the WBL can be attributed to the fact that–at short times–the natural cut-off, provided by the finite bandwidth, plays a crucial role for the dynamics.
 


§ MODEL AND METHOD
 

We consider a simple tight-binding model Hamiltonian to describe a molecular break junction. A single molecular level is connected to two metallic leads (cf. sketch in Fig. <ref>). The Hamiltonian reads

    Ĥ    = ϵ_c ϕ̂^†_cϕ̂_c + 
      ∑_α kϵ_α kϕ̂^†_α kϕ̂_α k
       = + ∑_α k( ϕ̂^†_α k V_(α k) cϕ̂_c + ϕ̂^†_c V_c (α k)ϕ̂_α k)  ,

where ϵ_c is the energy of molecular level, ϕ̂^†_α k and ϕ̂_α k are the field operators of the leads, with α=L(eft),R(ight) and k labels the basis functions in the leads, and ϕ̂^†_c, ϕ̂_c represent the field operators associated to the molecular level. The matrix elements V_(α k) c = [V_c (α k)]^⋆ take the coupling between the molecular level and the leads into account. The leads are modeled as non-interacting one-dimensional tight-binding chains, i.e., the dispersion of the electrons in the leads is given by

    ϵ_α k=-2t_αcos(k) + c_α ,

where t_α is the nearest neighbor hopping in lead α, yielding a bandwidth of 4t_α. The energy c_α corresponds to the center of the band of the lead α, i.e., it determines the alignment of the band with respect to the chemical potential, which we take to be at zero energy. Finally the hopping to the central site is V_α k = V_αsin(k). The embedding self-energy due to lead α is then given by

    Σ^R/A_α(z)     =  ∑_k V_(α k) c g^R/A_α k(z)
      V^⋆_(α k) c
        = |V_α|^2/t_α S(z-c_α/2 t_α)  ,

with g^R/A_α k(z) being the retarded/advanced Green's function of the isolated lead α. The function S(z) is given by

    S(z) = z - √(z-1)√(z+1) ,

where the character of the function S(z), i.e., whether it is the advanced or retarded self-energy, is determined by the sign of the imaginary part of z [The function S(z) has a branch cut on the real axis from z=-1 → z=1.]. The Green's function for the molecular level is then simply given by

    G^R/A(z) = [z - ϵ_c
      - ∑_αΣ^R/A_α(z) ]^-1 .

The inverse of the imaginary part of the self-energy yields a finite lifetime for the quasi-particles in the molecular junction, and the real part of the self-energy shifts the energy of the quasi-particles. 

The WBL is defined as the limit t_α→∞ (infinite bandwidth) while keeping the ratio |V_α|^2/t_α, which corresponds to the decay rate into lead α, constant. Expanding the expression of the self-energy for large t_α we obtain

    [Σ^R/A_α(z)]_WBL = ∓ i |V_α|^2/t_α = ∓ i Γ_α/2 ,

where the ∓ sign refers to the retarded/advanced self-energy, respectively. As we can see from this expression, the only effect of the leads is to provide a decay-mechanism for the quasi-particles.

Expressing the field operators ϕ̂ in the Heisenberg picture and using their equations of motion, the charge and heat currents are given by <cit.>:
 
  
    I_α    = - ∂_t ∑_k ⟨ϕ̂^†_α k(t)
        ϕ̂_α k(t)⟩ , 
    
        Q_α    = -∂_t [∑_k ϵ_α k⟨ϕ̂^†_α k(t) ϕ̂_α k(t)⟩.
         
       = .  + 1/2∑_k n( V_(α k) n⟨ϕ̂^†_α k(t) ϕ̂_n(t)⟩ + h.c.) ].


Note that we define the heat current Q_α as the temporal change in the energy within the leads plus half the coupling energy<cit.>.



§ TRANSPORT SETUP
 

In this work, we will investigate the validity of the WBL in the case of dynamical heat and charge transport in the junction described in the previous section. Once a non-equilibrium situation is created by applying a potential bias and/or temperature gradient, transient dynamics will take place and electrons will move, resulting in charge and heat currents flowing across the junction. We focus on the specific cases of quenches, i.e., the electric and thermo-mechanical potentials suddenly change at a certain time t_0. Transient dynamics, induced by changing the potentials, occur on the order of a characteristic time scale τ given by the inverse of the decay rate provided due to the leads, i.e., 

    τ^-1 = ∑_αV_α^2/t_α .

For times t ≫τ the junction will reach a steady state. We choose the hopping V_L = V_R = V as our unit of energy: the molecular energy level is taken to be at ϵ_c = 0.2 V, the nearest neighbor hopping of the leads t_L=t_R=5V, the chemical potential defines the zero of the energy, the center of the bands of the leads are aligned with it (c_α=μ=0), and the (inverse) temperature β = (k_B T_0)^-1 = 100 V^-1. 

For t<t_0 the system is taken to be in thermal equilibrium at temperature T_0. In order to induce a charge current through the junction the left lead is shifted up in energy by U=2V for t≥ t_0, i.e., the energy dispersion of the left lead is given by

    ϵ_L k = 
    -2t_Lcos(k)    for t<t_0, 
     -2t_Lcos(k) + U    for t ≥ t_0.

In order to describe a temperature gradient across the junction–in addition to the potential bias–we apply a thermo-mechanical potential ψ = T_α - T_0/T_0 = 1 in the left lead, which rescales the bandwidth for t≥ t_0. This thermo-mechanical potential effectively doubles the temperature in the left lead:

    ϵ_L k = 
    -2t_Lcos(k),    for t<t_0 
     (1+ψ)(-2t_Lcos(k) + U)    for t ≥ t_0.

Figure <ref> sketches of the molecular junction in the initial equilibrium and in the steady-state limit, showing that in the steady-state the energy dispersion of the left lead is broadened by a factor of two. In Fig. <ref> we depict the time-dependent charge and heat currents through the molecular junction. In this calculation, both, a charge bias and a temperature gradient is applied across the junction and we observe fast transient oscillation of the currents on the timescale τ followed by a saturation to a steady current.




§ RESULTS
 
In order to test the WBL we compute the time-dependent charge and heat currents flowing from the leads into the impurity in the WBL and compare the results to calculations taking the full frequency dependence of the lead self-energy [cf. Eq. (<ref>) and Fig. <ref>] into account. Specifically, we rescale the bandwidth of the leads, making it effectively wider, while keeping the ratio |V_α|^2/t_α constant. Hence, we use

    V^λ_α = √(λ)V_α   ,   
      t^λ_α =λ t_α ,

with a rescaling factor λ, which allows us to approach the WBL as λ→∞. 
We focus on two different scenarios: 1) A situation where only a potential bias is applied to the left lead [cf. Eq. (<ref>)], 2) A situation where, both, a potential bias and a temperature difference are applied across the junction [cf. Eq. (<ref>)]. The numerical algorithm to compute the transient currents–taking the full frequency dependence of the embedding self-energy into account–has been already discussed in Ref. EichVignale:16. In the following we refer to these results as the “full” calculation. Very recently progress has been made in evaluating the time-dependent currents in tight-binding models within the WBL analytically <cit.>. It turns out that this is also possible if a thermo-mechanical potential–describing temperature gradients–is present. Accordingly, all WBL results are obtained analytically. The explicit derivation of the analytical expression will be presented elsewhere.

Steady state currents. For times much longer than the characteristic lifetime τ the system reaches a steady state. In general we find that the steady state currents obtained in the WBL coincides with the results of the full calculation when the scaling factor λ is increased. However, there is a subtle point in the evaluation of the heat current in the steady state: it turns out to be crucial to take the WBL at the end of the calculation and not inside the integral defining the steady state current. The difference between taking the WBL inside the integral and taking the WBL after performing the integral is only present if two leads at different temperatures are connected to the same state in the device. This is trivially the case for a molecular junction modeled by a single site. We present a careful derivation in App. <ref> showing that the order of limits matters.


Transient charge current. The transient charge induced by a potential bias alone is nicely reproduced in the WBL [cf. App. <ref> for the corresponding plots]. If a temperature gradient–in addition to the potential bias–is applied the charge current exhibits a jump at the initial time, but otherwise represents the full calculation for times t ≳τ. In Fig. <ref> we depict the charge current for t ≪τ (a plot of I_α for t ≳τ is provided in App. <ref>). It can be shown analytically that the jump, Δ I_α, at the initial time is proportional to

    Δ I_α∝V̅_α( V_α/t_α
      - V̅_α/t̅_α)  ,

where the hopping amplitude inside the leads for t<t_0 is denoted by t_α and for t≥ t_0 by t̅_α. Similarly, we could write the coupling between lead α and the molecular region as V_α before t_0 and V̅_α afterwards. In our setup the couplings are held constant at all times, i.e., V̅_α = V_α, and the temperature gradient is mimicked by changing the hopping inside the leads as discussed in Sec. <ref>. Specifically, from Eq. (<ref>) we have t̅_L = (1 + ψ) t_L, which means that the charge current has a finite jump when a temperature gradient is applied across the molecular junction. Equation (<ref>), however, suggests that the jump can be avoided if the temperature bias is mimicked by scaling the couplings V_α in the same way as the hopping inside the leads. This would imply that V̅_α / t̅_α = V_α / t_α, which is sufficient to make Δ I_α vanish even in the presence of a temperature gradient. We stress that this cannot be achieved in the partitioned approach, because V_α is zero by definition for t<t_0 if the system is initially decoupled.


Transient heat current. Turning to the transient heat current we find that if only a potential bias is applied the heat current of the full calculation is reproduced in the WBL for times t ≳τ, but exhibits a jump at t_0. The short time behavior, t ≪τ is depicted in Fig. <ref> (cf. App. <ref> for t≳τ). Similar to the case of the charge current induced by a temperature gradient, we can see that the WBL approximates a discontinuity at t_0 in the limit λ→∞ in the full calculation. Again, the heat current in the full calculation always vanishes as t → t_0, but the WBL leads to a finite step in the heat current already in the presence of only a potential bias. In contrast to the charge current it is not possible to extract a simple expression as Eq. (<ref>), but instead the jump depends on the details of the molecular junction, i.e., on the quasi-particle energy levels.


If a temperature gradient is applied across the junction the heat current flowing from the left lead into the molecule depends strongly on the bandwidth for short times, even in the full calculation. In Fig. <ref> we can see that the heat current oscillates strongly for t ≲τ with a frequency proportional to the bandwidth (which, in turn, is proportional to λ). These oscillations correspond to transitions between the band edges of the leads and have been already observed in Ref. EichVignale:16. In the WBL these oscillations are absent since there are no band edges, but instead the heat current diverges as

    Q_α∼V̅_α( V_α/t_α
      - V̅_α/t̅_α) 1/t-t_0 .

In the partition-free approach this divergence can be tamed by rescaling the couplings V_α in the same way as t_α, i.e., by applying the thermo-mechanical potential not only inside the leads but also on the boundary of the junction. It turns out, however, that the sub-leading order for the heat current exhibits a logarithmic divergence as t→ t_0.

This is shown in Fig. <ref>, where we compare the transient charge and heat currents in the WBL scaling only the hopping inside the leads (I) and scaling also the coupling to the impurity (II). We can see that the charge current starts from zero for both leads if also the coupling V_L is rescaled with the temperature, but it exhibits a finite jump in the lead where the temperature is changed if the temperature only rescales the hopping inside the lead. For the heat current we see that currents in all leads exhibit a logarithmic divergence as t → t_0 if, both, t_α and V_α is rescaled due to the change in temperature. If only t_α is rescaled we see the aforementioned (t - t_0)^-1 divergence.



§ DISCUSSION AND CONCLUSION
 
In this work we have carefully examined the WBL for the transient and steady-state charge and heat currents through a molecular break junction. While we find that the long-time dynamics are faithfully captured in the WBL, at short times the WBL deviates considerably from a calculation taking the full frequency dependence of the embedding self-energy into account. This can be understood intuitively by considering that short times implies a wide spread in energy, and, therefore, the dynamics will be sensitive to whether the self-energies have a high frequency cut off (or decay) or whether they are constant for all frequencies. Specifically we have shown that the charge current induced by a temperature gradient, and the heat current induced by a potential bias, exhibit an unphysical jump at the initial time, when the system is suddenly quenched. Even more dramatically, the heat current diverges shortly after the quench if in addition to the potential bias also a temperature gradient is applied to the system. We have shown that these unphysical behaviors of the charge and heat current due to a temperature quench can be mitigated by considering that the temperature change not only affects the metallic leads, but also the boundary between the leads and the molecular junction. Since, in practice, this boundary is not sharply defined, we consider this a legitimate fix for the WBL. We stress that this fix can only be applied in the partition-free approach to the transport problem, i.e., when the coupling between the metallic lead and the molecular junction is already taken into account in the initial state of the system (before the quench). While this fix renders the charge current physical, in the sense that the initial current vanishes, the divergence in the heat current remains, but is only logarithmic. 

In order to address this, we see two possible solutions: 

1) In an actual experiment temperature gradients and potential biases will never be switched on infinitely fast, so a description as a sudden quench is questionable–to say the least–considering short time transient dynamics. It seems plausible that any kind of continuous switching will lead to a physical result (zero initial charge and heat currents). 

2) The second possible “solution” concerns the very definition of the energy or heat current between the leads and the device <cit.>. In this work, the heat current Q_α has been defined as the change in time of the internal energy of the leads plus half of the coupling energy. Alternatively, the energy current J_α from a certain lead α could be defined excluding the energy associated with the coupling to the device. The divergence at small times is due to the internal energy of the leads, which occurs in both the energy (J_α) and the heat current (Q_α). However there is yet another possible definition of the energy flowing between the leads and the device, i.e., we can define an energy current, E_α, via the change of the energy stored in the device itself excluding the coupling. This leads to an expression ∂_t ⟨ H ⟩ = ∑_α E_α, where ⟨ H ⟩ is the expectation value of the energy inside the junction. For a simple one-site model we have trivially E_α = ϵ_c I_α, i.e., the new energy current is proportional to the charge current. However, for a multi-state device region, this is not necessarily the case.

We point out that from a numerical point of view it would be highly desirable to employ the WBL to compute transient charge and heat currents, because it affords an analytical solution in terms of the quasi-particle states and energies in the molecular device. Harnessing this analytical solution would allow for an efficient simulation of mesoscopic devices. An approach somewhat intermediate between taking the full frequency dependence of the self-energy into account and the WBL, could be to approximate the self-energies by Lorentzians, which provide a self-energy with the proper decay at high frequencies while allowing also for a (semi-)analytic solution of the transport problem <cit.>.

In this work we only consider the non-interacting case. Through the Keldysh formalism <cit.> a generalization of the Landauer-Büttiker formula for interacting electrons is also possible <cit.>. There is a very interesting alternative approach for tackling the interacting transport problem using time-dependent density-functional theory (TD-DFT) <cit.>, where the interacting problem is mapped onto a fictitious noninteracting problem. This implies that the Landauer-Büttiker formula applies. The effect of the electron-electron interaction is taken into account via an effective potential, which renormalizes the effective bias driving the charge flow <cit.>. Furthermore, the coupling to a thermo-mechanical potential can be used to generalize TD-DFT to allow for a direct description of charge and energy flow <cit.>. We are confident that the combination of these approaches holds promise for studying the transient charge and energy flow in large molecular junctions for interacting systems.



F. G. E. has received funding from the European Union's Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie Grant Agreement No. 701796. R. T. and M. A. S. acknowledge funding by the DFG through the Emmy Noether programme (SE 2558/2-1). A. R. acknowledges financial support from the European Research Council(ERC-2015-AdG-694097) and Grupos Consolidados (IT578-13).






§ ADDITIONAL PLOTS FOR THE TRANSIENT CURRENTS


In this appendix we provide additional plots comparing the full results to the WBL. Figure <ref> shows the transient charge currents for times comparable to the quasi-particle lifetime τ (upper panel) and for very short times (lower panel). It clearly shows that the WBL represents the λ→∞ limit of the full calculation.

Fig. <ref> shows the charge and heat current for times comparable to the quasi-particle lifetime. The charge current is shown for the second scenario, i.e., when a potential bias and a temperature gradient is applied. The heat current is shown for the first scenario, i.e., when only a potential bias is applied at t_0. We see that in both cases the WBL currents correspond to the λ→∞ currents obtained in the full calculation for times t≥ t_0. However, in the WBL both currents approach a finite value for t → 0 as discussed in Sec. <ref>.







§ DERIVATION OF THE STEADY-STATE ENERGY CURRENT
 
Here we present the analytical evaluation for the steady-state energy current discussed in Sec. <ref>. In our derivation we consider a generic Hamiltonian of the form

    Ĥ    = Φ̂⃗̂^†·H·Φ̂⃗̂ + 
      ∑_α kϵ_α kϕ̂^†_α kϕ̂_α k
       = + ∑_α k( Φ̂^†_α kV⃗^†_α k·ϕ̂⃗̂
      + Φ̂⃗̂^†·V⃗_α kϕ̂_α k)  ,

which is the generalization of Hamiltonian introduced in Sec. <ref>, Eq. (<ref>) to multiple states in the molecular junction. We denote vectors in the single-particle state space of the impurity by bold symbols, e.g., Φ⃗ or V⃗_α k, and matrices by underlined bold symbols, e.g., H. In the WBL the retarded (R) and advanced (A) Green's function are given by
 
  
    G^R(z)     = ∑_n R⃗_n A⃗^†_n/z - ω_n , 
    G^A(z)     = ∑_m A⃗_m R⃗^†_m/z - ω^⋆_m .


The vectors R⃗_n are the right eigenvectors of the retarded eigenproblem and the vectors A⃗_m are the right eigenvectors to the corresponding advanced eigenproblem, i.e.,

  
    ( H - i 1/2Γ) ·R⃗_n     = ω_n R⃗_n  , 
    ( H + i 1/2Γ) ·A⃗_n     = ω^⋆_m A⃗_m  .


They are normalized by requiring A⃗_m^†·R⃗_n = δ_nm. The Landauer-Büttiker formula for the energy current explicitly reads

    J_α    = i ∑_α'∑_n mR⃗^†_m ·Γ_α·R⃗_n
      A⃗^†_n ·Γ_α'·A⃗_m 
       ×∫_-∞^∞dϵ/2 π iϵ[ f(ϵ_α) - f(ϵ_α') ]
      1/ϵ - ω_n1/ϵ - ω^⋆_m ,

where we introduced the short hand ϵ_α = β_α/β (ϵ - U_α), with β_α being the (inverse) temperature, and U_α the potential in lead α (note that β_α / β = [1 + ψ_α]^-1). Using the representation of the Fermi function in terms of a Matsubara summation, i.e.,

    f(z) = 1/2 - 1/β∑_f 1/z - i z_f ,

with the Matsubara poles z_f = (2f+1)π / β, we arrive at

    J_α    = -i ∑_α'∑_n mR⃗^†_m ·Γ_α·R⃗_n
      A⃗^†_n ·Γ_α'·A⃗_m 
       ×[ C^α_nm - C^α'_nm]  ,

with 

    C^α_nm = 1/β∑_f 
      ∫_-∞^∞dϵ/2 π iϵ/ϵ_α - i z_f1/ϵ - ω_n1/ϵ - ω^⋆_m .

Using

    ϵ/(ϵ - ω_n)(ϵ - ω^⋆_m)
      = ω_n/ω_n - ω^⋆_m1/ϵ - ω_n
      + ω^⋆_m/ω^⋆_m - ω_n1/ϵ - ω^⋆_m

we can further decompose

  
    C^α_nm    = 1/ω_n - ω^⋆_m[F^α_n + (F^α_m)^⋆]  , 
    
        F^α_n     = ω_n 1/β∑_f 
        λ^2 ∫_-∞^∞dϵ/2 π i1/ϵ_α - i z_f1/ϵ - ω_n .


The integral can be evaluated by closing the integration contour in the upper half of the complex plane, leading to

    F^α_n     = ω_n 1/β∑_f>01/i z_f - ω_n_α .

The sum (<ref>) does not converge, but it can be combined with the corresponding sum from C^α'_nm, i.e.,

    F^α_n - F^α'_n     = ω_n 1/β∑_f>0( 1/i z_f - ω_n_α
      - 1/i z_f - ω_n_α') 
        = ω_n ∑_f>0ω_n_α - ω_n_α'/(i z_f - ω_n_α ) (i z_f - ω_n_α') .

Expression (<ref>) can be summed explicitly using

    1/β∑_f>01/i z_f - x1/i z_f - y
      = D(-x) - D(-y)/x - y ,

where we defined

    D(z) ≡ -1/2 π iψ_0(1/2 - i z β/2 π)
       ,

in terms of the Digamma function ψ_0(z). This leads to the final result

    C^α_nm - C^α'_nm    = 
      ω^⋆_m/ω^⋆_m - ω_n[ D(ω^⋆_m_α) - D(ω^⋆_m_α') ] 
        + ω_n/ω_n - ω^⋆_m[ D(-ω_n_α) - D(-ω_n_α') ]
       ,

where we used that [D(z)]^⋆ = - D(- z^⋆), which follows from [ψ_0(z)]^⋆ = ψ_0(z^⋆), we can see that the expression for the energy current is a real number.

In the derivation presented above we have replace the frequency dependent self-energies by the frequency independent WBL approximation inside the integrand. In the following we will repeat the calculation keeping a “minimal” frequency dependence, i.e.,

    Σ^R/A_α(ϵ)
      ≈1/2Γ_αλ/ϵ± i λ ,

which reduces to the WBL as λ→∞. Only at the end of the calculation we will take the limit λ→∞. Equation (<ref>) implies that in the integral for F^α_n we have an additional factor of

    λ^4/(ϵ-iλ)^2(ϵ+iλ)^2
      = . ∂_x ∂_y λ^2/(ϵ-iλ x)(ϵ+iλ y)|_x=y=1 .

Accordingly, we have

    F^α_n = ω_n λ^2 ∂_x∂_y
       (1/β∑_f1/i λ x_a - i z_f1/iλ(x+y)(iλ x - ω_n)
        + 1/β∑_f>0( β_α/β)^2/(i z_f - iλ x_α)(i z_f + i λ y_α)
      (i z_f - ω_n_α))_x=y=1 .

Now we use

  
    ∑_f 1/z - i z_f = D(z) - D(-z)  , 
       ∑_f>01/i z_f - a1/i z_f - c1/i z_f - c
        = D(-a)/(a-b)(a-c) + D(-b)/(b-a)(b-c)
        + D(-c)/(c-b)(c-b) ,


to arrive at

    F^α_n     = ω_n λ^2 ∂_x∂_y( D(iλ x_a)/iλ(x+y)(iλ x - ω_n)
        +D(iλ y_a)/iλ(x+y)(iλ y - ω_n)
        + D(-ω_n_α)/(ω_n - i λ x)(ω_n + i λ y))_x=y=1 .

From the asymptotic expansion of the Digamma function, ψ_0(z) ∼log(z), it follows that

    D(iλ z_α) ∼ -log[zβ_α/2 π]/2 π i ,

which, in turn, leads to the asymptotic expansion 

    F^α_n    ∼ω_n ∂_x ∂_y 
      ( D(-ω_n_α)/xy
        + log[λ x β_α/2 π]/2 π i(x+y)x 
      + log[λ y β_α/2 π]/2 π i(x+y)y)_x=y=1 .

Combining the terms due to the different leads yields

    F^α_n - F^α'_n    ∼ω_n (D(-ω_n_α)
      - D(-ω_n_α') 
       ∼ω_n ( - i 1/2πlog[β_α / β_α'] )  .

Equation (<ref>) does not depend on λ anymore and we can safely take the limit λ→∞, because the neglected terms in the asymptotic expansion are of order λ^-1. Plugging Eq. (<ref>) into Eqs. (<ref>) and (<ref>) finally gives

    J_α    = -i ∑_α'∑_n mR⃗^†_m ·Γ_α·R⃗_n
      A⃗^†_n ·Γ_α'·A⃗_m 
       ×( ω_n/ω_n - ω^⋆_m[D(-ω_n_α) - D(-ω_n_α')] 
       ×( + ω^⋆_m/ω^⋆_m - ω_n[D(ω^⋆_m_α) - D(ω^⋆_m_α')] ) 
        + 1/2 π∑_α'Tr[ Γ_α·Γ_α']
      log(T_α/T_α')  .

The first term corresponds to the result when taking the WBL inside the integral. The second term is the correction if, instead, the WBL is taken after the integration. This correction term vanishes if there is no temperature gradient between the leads. Furthermore, it vanishes if the coupling matrices do not overlap in the single-particle state space of the molecular Hamiltonian. In the results presented in this work this correction is crucial in order to reproduce the steady-state heat currents in the WBL.





