Charge and heat transport through quantum dots with local and correlated-hopping interactions

The transport properties of junctions composed of a central region tunnel-coupled to external electrodes are frequently studied within the single-impurity Anderson model with Hubbard on-site interaction. In the present work, we supplement the model with an important ingredient, namely the charge-bond interaction, also known as correlated or assisted hopping. Correlated hopping enters the second-quantised Hamiltonian, written in the Wannier representation, as an off-diagonal many-body term. Using the equation of motion technique, we study the effect of the correlated hopping on the spectral and transport characteristics of a two-terminal quantum dot. Two different Green functions (GFs) appear: one of them describes the spectral properties of the quantum dot, the other the transport properties of the system. The calculation of the transport GF requires the knowledge of the spectral one. We use decoupling procedures similar to those which properly describe the standard Anderson model within the Kondo regime and outside of it. For an arbitrary ratio $x$ between the amplitudes of correlated and single-particle hopping terms, the transport GF fulfils the $x \leftrightarrow 2-x$ symmetry of the model. The average occupation of the dot also obeys this symmetry, albeit the spectral function of the quantum dot, calculated within an analogous decoupling scheme as for the transport GF, does not. We identify the physical reason for this behavior, and propose a way to cure it. Since the correlated-hopping term breaks the particle-hole symmetry of the model and modifies all transport characteristics of the system, the detailed knowledge of its influence on measurable characteristics is a prerequisite for its experimental detection. Simple, experimentally feasible methods are proposed.


I. INTRODUCTION
The study of quantum dots [1] continues to enjoy a high popularity. Transport properties of nanostructures consisting of quantum dots (QDs) or molecules placed between two or more external electrodes have been intensively investigated in the past few decades [2,3], with the goal, e.g., to achieve efficient heat to electricity conversion at the nanoscale. Nanodevices with ferromagnetic or/and superconducting leads may be relevant as sources of a pure spin current [4] or entangled electrons [5], needed for spintronics and quantum information technology [6]. For example, spin valves have a high potential for controlling spin currents [7].
The properties and functionalities of such structures strongly depend on the state of the leads, their coupling to the central region [8], the interactions of the electrons on the central region, and on external conditions like temperature, magnetic field, etc. The experimental control of the relevant parameters and the theoretical understanding of their effect on measurable characteristics of devices is at the heart of their application potential.
The standard theoretical modelling of such systems is based on the Anderson Hamiltonian: where H LR describes the external leads (L: left lead, R: right lead), H C the central region typically containing the Hubbard [9] repulsion, and H tunn the coupling between the leads and the central QD. The coupling is visualised as the tunneling of electrons between electrodes and QD, and is described by the following term in the Hamiltonian: The amplitude V λkσ is a single-particle transfer proportional to Ψ λkσ |h(r)|Ψ cσ , where Ψ cσ and Ψ λkσ are the wave functions of the central region and the extended states in the leads (λ = L, R), respectively, and h(r) the single-particle part of the first-quantised Hamiltonian of the system. It turns out that besides the single-particle term there may exist another term promoting the transfer of electrons from the electrodes to the central region and vice versa. This term, with the amplitude denoted by K λkσ , has a many-body origin, with In second-quantised representation, the corresponding Hamiltonian is given by whereσ denotes spin opposite to σ, i.e.,σ = −σ. H ass is known as assisted or correlated hopping [10]. Due to its dependence on the charge state of the central site, it is also called charge-bond interaction [11]. Apparently, this term describes the transfer of a spin-σ electron between the dot and the electrode, provided another electron with opposite spin occupies the dot. Both terms have been of considerable interest in studies of strongly correlated bulk materials. The singleparticle hopping V has been intensively investigated in the context of heavy fermions [12], where it describes the coupling between, e.g., d and f orbitals. On the other hand, the presence of the correlated-hopping term K in solids has been found to affect collective properties of materials; it was mainly studied in the context of hightemperature superconducting cuprates as the interaction promoting the appearance of the superconducting instability [13], and/or explaining the asymmetry between the superconducting domes of electron-doped and hole-doped compounds [14,15].
In the context of nanodevices, correlated hopping is expected to be present in most cases. However, this term has not attracted the attention it probably (in our opinion) deserves. In fact, studies of correlated hopping in the context of transport via nanostructures are sparse [16][17][18][19][20][21][22], and mainly by numerical techniques. Correlated hopping has been proposed to explain anomalous features observed inter alia in transport through quantum point contacts [23] and single-electron molecular transistors [24], but-as far as we know-no compelling evidence exists on its experimental relevance. This calls for detailed theoretical studies, in order to find and quantify possible ways for its experimental detection.
In this work, we present a systematic analysis of the role of correlated hopping in transport via QDs, i.e., the modifications it introduces to the standard behavior of the single-impurity Anderson model with Hubbard-only interaction. To solve the problem analytically, and to gain a deeper insight into the physics of correlated hopping, the equation of motion (EOM) technique is employed. We return to this aspect in the concluding section.
While the EOM method does not provide the exact GFs (except for non-interacting systems), as it relies on decouplings and projections of higher-order GFs onto lower order ones, it is easy to implement in quite arbitrary situations and for (almost) arbitrary Hamiltonians, in the linear (small voltage) regime and beyond. Nevertheless, it has to be noted that the near-equilibrium results of [20,21] are based on the numerical renormalization group (NRG) technique, which is known to capture correlation effects in an essentially exact manner, as discussed recently in some detail [25].
The goals of the paper are: (i) to generalise the EOM method, earlier applied to the Anderson Hamiltonian (where it properly describes Kondo correlations [26,27]), to the model with correlated hopping; (ii) to analyze the role of this contribution, which breaks particle-hole symmetry, on the Kondo peak (its width, temperature dependence, etc.) and the transport characteristics of a two-terminal system; and (iii) to determine kinetic and transport coefficients of the two-terminal QD, and iden-tify experimental signatures of the extra term. Last but not least, we demonstrate the power of the EOM approach by presenting selected results beyond the linear regime.
In the next section, Sec. II, we describe the model and its parametrization used throughout the paper, and express the charge and heat currents in terms of the appropriate GF. The full set of equations for the transport GF is introduced in Sec. III. The results are presented and discussed in Sec. IV and Sec. V. The summary and conclusions are given in Sec. VI. Some technical details and lengthy calculations are relegated to the appendices. The Supplementary Material [28] contains additional details.

II. THE MODEL AND BASIC DEFINITIONS
The Hamiltonian to be studied is similar to the standard single-impurity Anderson model, albeit modified to include correlated hopping: where n λkσ = c † λkσ c λkσ and n σ = d † σ d σ denote particle number operators for the leads and the dot, respectively. The operators c † λkσ (d † σ ) create electrons in respective states λkσ (σ) in the lead λ (on the dot). The energies of the leads are measured from their chemical potentials µ λ , so ε λk = ε 0λk − µ λ , with the dependence of ε 0λk on λ allowing for a different spectrum in each of the leads. The spin is σ = ±1 (↑, ↓), and ε σ = ε d + σµ B B, with B the magnetic field, µ B the Bohr magneton, and ε d the dot electron energy level. The Hubbard parameter U describes the repulsion between two electrons on the dot. The operator D σ = d σ (1 − xnσ) takes care of the occupation dependence of the hopping.
The state-dependent hopping has been parameterised by x which is minus the ratio between K λkσ and V λkσ , The parameter x, in principle, may be complex, and even in case it is real it may have both positive and negative values [20,21]. For simplicity, it is assumed not to depend on λkσ, i.e., to have the same constant, spin and wave vector independent value for both leads. This assumption should hold provided the two leads are composed of similar (or even identical) materials. Generally we also expect the k dependence to be of lesser importance, as usual in Fermi liquid theory. However, the spin dependence may become relevant for magnetic leads. Here, following [20,21], we assume x to be real, and focus on the interval 0 ≤ x ≤ 2.
A. Currents in the two-terminal system The charge current and energy current flowing out of the electrode λ are calculated as the time derivative of the average charge, N λ = kσ n λkσ , respectively average energy H λ = kσ ε λkσ n λkσ of lead λ. The derivation is sketched in App. A. Application of those results to the two-terminal QD we are interested in here, provides I = I L = −I R , which expresses current conservation in the system: describing the coupling between the dot and the electrode are assumed to be independent of energy E, which corresponds to the wide-band limit. Similar expressions can be derived for the heat current flowing from the left, and right electrodes: (9) It can also be verified thaṫ Q + (µ L − µ R )I = 0 (10) in agreement with energy conservation: hereQ = J L +J R is the total heat current leaving the leads. We emphasize that the transport GF, differs from the spectral one, the latter inter alia describing the occupation n σ of the quantum dot. The function is the Fermi distribution function describing the electrons in the lead L/R, assumed to be in equilibrium at temperature T L/R and chemical potential µ L/R . The above formulae for the currents are valid for arbitrary voltages V = (µ R − µ L )/e, where e is the electron charge. In particular, Eq. (6) allows the calculation of the conductance beyond the linear regime. In the general (V = 0) case, we define the differential conductance as B. Linear transport coefficients Assuming the temperature difference between the right and left electrode ∆T = T R − T L as well as the voltage V to be small parameters, we can expand the formulae (6) for the current across the system and the similar one for the heat fluxQ = J L + J R , in order to obtain the (symmetric) Onsager matrix of linear kinetic coefficients L ij , as well as the related set of transport parameters: the conductance (G), the Seebeck coefficient (S), and the thermal conductance (κ) [3,29,30]. In the present model, the linear coefficients are given by the moments M n of the imaginary part of the transport GF, ImG r σ (E) (14) where f = ∂f (E, T )/∂E; here we set µ = µ L = µ R , and T = T L = T R . The linear conductance and Seebeck coefficient read It has to be noted that the transport density of states, N tr (E) = (−1/π)Im D σ |D † σ r E , may have features narrow on the scale of k B T . In such a case the Sommerfeld low-temperature expansion [31] is not valid, and one has to use the above expressions to calculate the linear Seebeck coefficient. For parameters such that N tr (E) is a smooth function of energy on the scale k B T around the chemical potential, one finds the approximate formulae where the prime means the derivative with respect to energy.

III. CALCULATION OF THE GREEN FUNCTIONS
To calculate the spin-dependent retarded GFs [32] G r σ (ω) and g r σ (ω), we use the EOM method [33] and the approximation scheme known as Lacroix approximation [34][35][36], with some important extensions proposed recently [26]. In this section and the following ones, we shall (mostly) work in units such thath = k B = 1. We also use frequency as argument of all GFs below, and omit the r subscript with the understanding that we shall first calculate the retarded GFs, and advanced and lesser GFs will be obtained from them by known relations [28].

A. Transport Green function
The transport GF D σ |D † σ ω has been calculated in App. B, and we only quote the final formula here: where and For various definitions and details, see App. B, in particular, Eq. (B43) in conjunction with Eqs. (B42) and (B44)-(B46).

B. Spectral Green function
The following expression for the spectral GF d σ |d † σ ω has been obtained by employing the same decoupling scheme as above [28]: where and It has to be stressed again that both GFs, i.e., the spectral and transport one are coupled together. They both have to be calculated simultaneously as various quantities they depend on require the knowledge of both of them. Needless to say that for x = 0 the transport GF (19) reduces to the spectral one (25) as it has to be, and that the result agrees with the formula found earlier by Lavagna [26].
The symmetry of the Hamiltonian suggests that both GFs are symmetric with respect to x = 1. The approximate transport GF given in Eq. (19) is indeed symmetric, and calculated for x = 0 even analytically is the same as for x = 2. However, this is not correct for the spectral function, Eq. (25): its value at x = 2 is not the same as that for x = 0. This issue, which, however, does not affect the symmetry of the transport coefficients, will be discussed later on.
In a recent paper [27] the GF for the x = 0 model has been studied in the context of a three-terminal QD, which in the strongly non-equilibrium limit works as a heat engine. We have shown that the above formula is quantitatively correct in describing the spectral and transport properties of the QD in the Kondo regime, even in the particle-hole-symmetric case which is notoriously difficult to capture by the EOM technique. This is true in equilibrium as well as far from it. In Sec. IV we shall discuss the transport characteristics of the two-terminal quantum dot with correlated hopping, but with focus on the linear regime.

C. Lifetimes: second order calculations
When writing the expressions for various self-energies, we have introduced the parametersγ σ 1 andγ 2 , which replace the infinitesimal parts γ = 0 + in the self-energies; clearly, they represent the inverse lifetimes of singly and doubly occupied states on the dot, respectively. These decay rates take into account higher-order processes neglected at the present level of approximation. The importance of including such decay rates for the proper description of the Kondo resonance has been observed in [37], and found to result from higher-order processes. Lavagna argued later [26] that they can be calculated perturbatively using Fermi's golden rule. She also noted that fourth order contributions vanish for systems in equilibrium and without magnetic field.
The direct usage of Fermi's golden rule, with V I being the tunneling part of the Hamiltonian (5), leads to the following expressions valid to second order in the dot coupling V λkσ : The last equation shows that the contribution from the doubly occupied states vanishes for x = 1, as expected: for this value of x the doubly occupied state is totally decoupled from the system.

IV. SPECTRAL AND TRANSPORT GFS: NUMERICAL RESULTS AND SYMMETRY DISCUSSION
This section is devoted to the presentation of the results for the system in equilibrium, and for vanishing external magnetic field and spin-independent tunnelings. The transport coefficients can be calculated in the linear regime via Eqs. (14)- (16). In the following, all energies are measured in units of Γ 0 = Γ L ↓ = Γ L ↑ . We start the discussion with the imaginary parts of the transport and spectral GFs, i.e., the transport and dot's densities of states.

A. The Green functions and their symmetry
A brief inspection of the formulae (19) and (25) for the transport and spectral GFs suggests that the former is symmetric with respect to changes of x by (2 − x), while the latter is more difficult to judge due to a more complicated x dependence. Thus we resort to numerical calculations and postpone further discussion of these formulae to the end of this section. In Figs. (1)-(5) we show the imaginary parts of both GFs as a function of energy for a number of x values. The numerically perfect symmetry of the transport GF with respect to x = 1 is visible in Fig. (1a). The imaginary part of the transport GF is plotted as function of energy for a number of x and symmetry-related (2 − x) values. The (hardly visible) dashed black curves for the parameter (2 − x) overlap with the colored curves for x. Panel (b) of the figure is a magnification of the part of panel (a) close to the maxima, where the differences are largest. The observed differences are typically smaller than 1%.
The systematic evolution of the transport density of states with changing x is shown in Fig. (2). We present the results for x ∈ [0, 1] only, as the curves for the other values of x are related by symmetry. For illustration we have assumed a particle-hole symmetric situation with ε d = −4, U = 8, and the same temperature for both leads, T = 0.3. For these parameters and x = 0, both transport and spectral GFs are identical, and the quantum dot is in the Kondo regime. Hence one observes the Abrikosov-Suhl resonance, also known as Kondo resonance, at the chemical potential µ = 0, and two Hubbard bands located symmetrically around zero energy. The occurrence of the Kondo effect in the particle-hole symmetric Hubbard model shows the power of the present version of the EOM technique supplemented with lifetime effects [26].
The increase of x results in distinctive modifications of the transport density of states. First, one notices that the curves (for x = 0) are no longer particle-hole symmetric. Concomitant with this observation is the modification of the Kondo resonance, which develops a strong asymmetry with respect to the chemical potential (E = 0). The lower Hubbard band changes rather weakly with x. Its height slightly increases, and the position slightly moves towards the chemical potential. Most dramatic changes are apparent in the upper Hubbard band, which strongly decreases with increasing x, getting narrower and finally vanishing completely for x = 1. For the considered parameters the center of the upper Hubbard band moves slightly to the right in the figure. The result for x = 1 requires additional comments. Let us note that the operator D σ for the doubly occupied state with nσ = 1 vanishes. Under these conditions, the upper Hubbard band composed of the doubly occupied states does not contribute to transport.
In Figs. (1) and (2) we have shown the evolution of the transport density of states with x for the particlehole symmetric model for which δ = ε d + U/2 = 0. Non-zero values of x break the particle-hole symmetry  of the model. It turns out that for arbitrary values of δ the transport GF is symmetric with respect to x, albeit the differences are slightly larger than those in One can see that for this set of parameters the transport density of states N tr (µ) at the chemical potential (µ = 0) strongly changes with x. Both the absolute value and the slope are affected. Considering the formulae (15) and (16) for the transport parameters, which strongly depend on the transport density of states close to E = 0, one expects for these parameters a noticeable changes of both conductance and thermopower with x. On the other hand, for the set of parameters used in the Fig. (3b), the transport GF for energies close to the chemical potential hardly changes with x; thus both conductance and thermopower are expected to vary only slightly with x. The symmetry of the transport GF with respect to x ↔ 2 − x ensures, as we shall see in the next section, the same symmetry of the transport coefficients.
. The other parameters are ε d = −4, U = 8, and T = 0.3. The inset shows the x dependence of the average occupation n per spin for ε d = −2, −4, −8. For this quantity, the departures from perfect symmetry are very small. In order to increase their visibility, we plot the dependence n vs. x as solid lines, and the dependence n vs. (2 − x) as points of the same color.

B. The x vs. (2-x) symmetry in further detail
As argued in detail in [21], the model at hand is symmetric under the transformation x ↔ 2 − x. However, we find that the (approximate) spectral GF derived above does not obey this symmetry, in contrast to the (approximate) transport GF. In particular, this deficiency is already apparent in the analytical formula for the spectral GF, Eq. (25). The departures from the x ↔ 2 − x symmetry are clearly visible in Fig. (4), where we show the density of states for x = 0.2, 0.5, and the symmetry related values (2 − x) = 1.8, 1.5. The comparison of the curves obtained for the pairs of various x (0.2 and 1.8, and 0.5 and 1.5, respectively) shows that the differences are largest for energies close to the chemical potential µ, i.e., the point where the Kondo effect is expected. The differences at other energies seem to be related to those at µ by the sum rule, +∞ −∞ dEN (E) = 1, which is always fulfilled with an accuracy better than 1%.
To obtain the above results, Fig. (4), we have used the formulae (25) and (26), which have been obtained, see the Suppl. Material [28], using the decoupling scheme I. As discussed there, we have tried several different decouplings. The others, i.e., II and III, overall lead to the same behavior with small quantitative changes only, hence we are not showing the results for them here. For the discussion of decouplings II and III, and also a calculation scheme different from that presented in Sec. III, see App. D and App. E below.
Interestingly, despite the asymmetry of the spectral GF, the average charge density per spin is symmetric un- We see that the symmetry is obeyed with an accuracy of ≈ 0.01, which is only slightly larger than the accuracy of the iterative computation: In the iterative process, n is used as a check of the accuracy of the solution, and we terminate the iteration when the change in n is less than 0.001 in consecutive steps.
Anticipating the discussion in the next subsection, we expect in fact that the calculation of the spectral GF in the interval 0 < x < 1 is, in the present approximation, more reliable than the results obtained for 1 < x < 2. Hence we focus on the former regime, and illustrate in Fig. (5) in more detail the changes of the density of states with increasing x for the particle-hole symmetric case, δ = 0. The general trends in the spectral GF for arbitrary δ are similar to those observed for the transport GF. The modifications of the lower Hubbard band are relatively small, while the Kondo peak and the upper Hubbard band are strongly modified with increasing x. The Kondo resonance disappears, and the upper Hubbard band gets narrower and higher with its center shifting initially towards higher energies.

C. Remarks on the asymmetry of the spectral GF
A careful look at the spectral GF for x = 2 in Fig. (6) shows that a small dip appears at the Fermi energy (µ = 0). Such a dip in the energy dependence of the equilibrium density of states at the chemical potential is a characteristic feature of all previous decouplings [38,39] for the standard Hubbard model, (i.e., for x = 0). The only approach which cures the deficiency is that of Lavagna [26] for the Hubbard model, which is applied here for the correlated-hopping model. Why does the approach fail at x = 2? To elucidate the reason why Lavagna's approach is not effective for the spectral function at x = 2, we have to recall that the existence of the Kondo resonance for the Hubbard model in her approach is intimately related to the lifetime effects, i.e., the use of the parametersγ (1) σ ,γ (2) as discussed earlier [26]. First we note that the explicit dependence on x in the formulae (19) and (20) for the transport GF is through the factors x(2 − x) or (x − 1) 2 . The former vanishes for x = 0 and x = 2, and the latter is symmetric with respect to x = 1. On the other hand, the inspection of the formula (25) for the spectral GF shows that for x = 0 all extra terms vanish, while for x = 2 they do not but rather give a large contribution to the self-energies b 1σ , Σ T 1σ , and b 2σ , Σ T 2σ . These self-energies at low T lead to logarithmically divergent contributions close to the Fermi energy. For x = 0, however, the divergent contributions from those terms which remain in B d and n d eff are cut off by the lifetime effects. On the other hand, this is not the case for x = 2, and a number of diverging self-energies remain.
With this insight, we expect that in order to obtain the correct symmetry of the spectral GF one has to calculate, instead of projecting, those GFs which contain two lead operators. A careful inspection of the decoupling procedures shows, in fact, that the problem lies in the vanishing of certain contributions for x = 0 but not for x = 2. Thus, in order to render the expression for d σ |d † σ ω symmetric, extra terms are needed. These can only result from higher-order contributions to the S sp n,d,c terms. As already noted, the calculations of these "missing" terms can, in principle, be performed in full analogy to previous calculations for the Hubbard model [37], but for the model at hand they are very complicated and will not be pursued here. Thus we conclude that contrary to the Hubbard model, where lifetime effects [26] can mimic the role of the fourth order terms [37], the symmetry of the spectral function of the correlated-hopping model requires calculations up to fourth order in the tunneling amplitude.

V. TRANSPORT CHARACTERISTICS
A. Charge and heat conductance, and thermopower: linear regime We emphasize again that all transport coefficients depend on the transport GF only, and hence fulfil the required x symmetry. We focus the following discussion on the conductance G, the thermopower S, and the Wiedemann-Franz ratio L = κ/(GT ). In agreement with the preceding discussion, the nearly perfect x vs. (2 − x) symmetry is very well visible in G and S. In Fig. (7) we show the linear conductance G and the thermopower S: the comparison of these quantities vs. x (solid curves) and vs. (2 − x) shows that the symmetry is well obeyed. In Fig. (7) the symmetry is illustrated for the model with δ = 0, but it is also valid for arbitrary values of this parameter. For non-vanishing values of δ the functional dependence of G(x) vs. x differs but the symmetry remains intact. Interestingly, the overall dependence of the conductance on x shown in Fig. (8) is similar to that in the previous figure. In both cases the conductance takes on maximal values for x = 0 and x = 2. However, it has to be noted that the behavior of G(x) in the particle-hole symmetric case (shown in Fig. (7)) is due to the destructive effect of x on the Kondo peak at the chemical potential. An increase of x destroys the Kondo resonance, which leads to a smaller conductance. On the contrary, for the parameters of Fig. (8) it is the upper band of the transport density of states which is located close to the chemical potential that gives the largest contribution to G. The modification of this part of the transport spectrum (as visible in Fig. (3a)) determines the x dependence of G. The thermopower dependence on x in the above two models is more complicated. In the first case, S attains zero values for the perfectly symmetric transport densities of states at x = 0 and x = 2, and decreases for x tending towards 1. The complicated sign changes of S for the model with ε d = −8, shown in Fig. (8a), can be approximately read off from the slope of the transport density of states shown in Fig. (3a). This is due to the fact that in the linear regime S is proportional to the derivative of N tr (E) taken at the chemical potential; cf. Eq. (18). We have also calculated the thermal conductance, and its dependence on x traces the charge conductance. Thus we are not showing the plots here. Instead in Fig. (8b) the Wiedemann-Franz ratio L = κ/(T G) normalised to the Lorenz number L 0 = (π 2 /3)(k B /e) 2 is presented (κ denotes the thermal conductance). One observes that for the parameters used the ratio is smaller than unity. This indicates a non-Fermi-liquid behavior of the system with hampered charge transport.
One of the most interesting questions related to the present study is the identification of measurable consequences of correlated hopping. In this context, it should be emphasized that one has experimental control over virtually all parameters of the devices under discussion. In particular, the gate bias independence of the parameters Γ λ demonstrated recently [40] for QDs fabricated using the electromigration technique, supports the hope to achieve this goal.
Naturally one would like to measure some transport characteristics of the system and infer the information about the actual value of x. The symmetry with respect to changing x cannot serve the purpose, as this parameter most likely is beyond experimental control. However, other characteristics of the transport coefficients come to mind: namely the existence of the distinctive peaks in the conductance, and the concomitant saw-tooth shape of the thermopower. The answer depends on whether one considers the linear or the strongly non-linear transport regime. In the linear case (we are interested in here), and for temperatures low on the scale of Γ 0 , the transport parameters probe, as visible from Eqs. (14), the region of energy of the width of a few k B T in the vicinity of the chemical potential (here µ = 0), so the observed changes in the transport density of states are directly measurable (17). We claim that it is the x dependent relative height of the two conductance peaks of the single-level quantum dot which gives direct information on the correlatedhopping term. The observed maxima of G, measured as function of gate bias, are related to the corresponding maxima in the transport density of states: one peak builds up when the on-dot energy band around ε d crosses the Fermi level (µ = 0), and the other when the upper Hubbard band centered around 2ε d + U sweeps through µ. For x = 0 the two peaks are identical as visible in Fig. (9a). Similarly the corresponding "anti-symmetry" is observed for the thermopower, as visible in panel (b) of Fig. (9). This argumentation is applicable for a symmetrically coupled quantum dot, i.e., for Γ L = Γ R .
In Fig. (10a) we show the dependence of the conductance on gate bias, i.e., on ε d , for a few values of x, namely x = 0, 0.1, 0.2, 0.3, and 0.5. One observes a change of the relative height between the lower and the upper conductance "bands" with increasing x. The lower conductance peak decreases with x while the upper one stays constant. The observed decrease is faster than linear as shown in the inset to the figure. The proportionality factor in the linear fit, here equal to 0.04, is not universal, but depends on the details. However, the decrease of the lower peak height with x is a universal effect for a symmetrically coupled quantum dot, and gives immediate information on the very presence of the correlated-hopping contribution. We emphasize that the observation of the different heights of the two consecutive conductance peaks in the two-terminal quantum dot provides a unique proof of the existence of correlated hopping. As we are studying a single-level quantum dot, the main condition related to experimentally studied devices is that the distance between consecutive levels in the dot has to be larger than the Coulomb repulsion U . Otherwise, the consecutive conductance peaks would correspond to singly occupied levels. This probably is the most serious condition to fulfil. Additional information can be drawn from the analysis of the thermopower. However, the gate bias de- pendence of S is slightly more complicated, as is visible from Fig. (10b): An increase of x leads to an increase of the amplitude of S for δ values corresponding to the lower conductance peak, δ ≈ −12, while S remains virtually unchanged for δ ≈ +10, corresponding to the upper conductance peak. The effect of non-symmetric couplings on the conductance and the thermopower is shown in Fig. (11) for Γ R /Γ L = 2, 3, 4, 6. The anisotropy only weakly affects the lower conductance peak. Its effect on the upper peak is appreciable and includes a decrease of the magnitude and an increase of the width. The former effect masks the asymmetry in the peak heights induced by finite x, hence its unique identification becomes more difficult. However, the thermopower (Fig. (11b)) reacts in a more complicated way. Its overall amplitude diminishes in comparison to the symmetrically coupled dot, but both low and high δ parts are modified. The decrease of the magnitude of the thermopower variations can be understood by noting that S is proportional to the slope of the conductivity at the corresponding energy, which is known as Mott relation.  Our EOM results quantitatively agree with those obtained by the NRG technique [20,21]. In particular, an increase of x induces similar modifications of the conductance and thermopower in both methods. This is well seen by comparing, e.g., our Fig. (10) with panels (a) and (b) in Fig. 3 of [21]. It is more difficult to directly compare our results with those presented in [20], as these authors concentrate on such aspects as spin conductance and temperature dependencies. However, some curves shown in their Fig. 3 are close to our results for the corresponding set of parameters.
As a brief intermediate résumé, we note that all measurable characteristics exhibit the required symmetry properties. In particular, we have argued that a detailed experimental analysis of the gate bias dependence of both conductance and thermopower, in devices without orbital degeneracy and such that an adequate control of the symmetry of the couplings is feasible, may elucidate the role played by the correlated hopping, and may even allow for the extraction of x.

B. Non-linear conductance
The non-equilibrium Green function approach is well suited to treat finite voltages, since the EOM captures, albeit in a not well controlled way, higher-order scattering processes including those analysed in [41]. However, the full analysis of the conductance and other transport characteristics in the non-linear regime is beyond the scope of the present paper: these quantities not only depend on x, δ, and V , but also on temperature, Coulomb interaction, and the anisotropy of the couplings.
Here we focus on the differential conductance, G d (x, δ, V ) = ∂I/∂V . We present results for the dependence of G d on x for δ = −4 and a number of voltages V (Fig. (12), panel (a)), and the dependence of G d on δ for x = 0. 3 (panel (b)). In both cases U = 8 and T = 0.3. The linear conductance, formally corresponding to V = 0, is shown by red pluses. For a small voltage, V = 0.1, the differences between small-voltage and zerovoltage results are small, but they strongly increase with V . The departures from the linear regime are more pronounced for small x (and 2−x), and close to the resonant values of the gate bias when the conductance is maximal.
For small values of x and 2 − x, the changes of G d with V are relatively large, but decrease for x approaching x = 1, see panel (a) in Fig. (12). These variations can be understood by recalling the full formula (6) for I, and noting that the spectral Green function entering it also depends on the voltage V in a rather complicated way. For a finite voltage, the energy integration interval depends on temperature, but generally is of order V at low T , thus the actual value of the current and the differential conductance depend on the detailed behavior of the transport spectral density in the considered energy interval. For x = 1, the differences are very small, due to the rather smooth dependence of the transport spectral density on energy in the region between µ L and µ R (cf. Fig. (3)). It is worth noting that the x ↔ 2−x symmetry of the conductance is valid also in the non-linear regime.
The asymmetry between the conductance peaks for positive and negative values of δ, visible in panel (b) of Fig. (12), are related to the particle-hole symmetry breaking by the correlated hopping. The asymmetry depends on the actual value of x and may thus be utilized in precise experiments to obtain information on its value.

VI. SUMMARY AND CONCLUSION
The transport of charge and heat via nanostructures consisting of a quantum dot (QD) coupled to two external normal electrodes has been studied by the nonequilibrium Green function (GF) technique in combination with the equation of motion (EOM) approach. The system is described by the Anderson Hamiltonian containing not only the standard single-particle tunneling term, but also an additional one of many-body origin, known in the literature as assisted or correlated hopping. This term, often neglected in the analysis of transport measurements of QD nanostructures, breaks the particlehole symmetry of the model and thus may be important or even decisive for the interpretation of various experiments. The correlated-hopping term, which we have characterized by the parameter x (0 ≤ x ≤ 2), modifies the tunneling part of the single-impurity model. Employing the non-equilibrium GF method to describe charge and heat transport, it becomes apparent that two different Green functions are needed: one of them, which we denote transport GF, enters the formulae for the charge and heat flux, while the other is related to the dot density of states and hence the dot's average occupation. However, the equations for both GFs are coupled to each other via the occupation and certain self-energies. To obtain each of the GFs within the EOM technique, one has to project higher-order GFs, arising in the chain of EOMs, onto lower order ones. The simplest decoupling leads to a transport GF fulfilling the x symmetry of the model. However, the spectral GF is found not to be symmetric under x ↔ 2 − x. Attempting to cure this deficiency, we have tried three different decoupling schemes (see [28]), but failed to achieve the required symmetry. Even the calculation of the full matrix GF with various decouplings did not restore the x-symmetry. We argue that the symmetry restoration in the spectral GF g r (E) requires the inclusion of higher order GFs, thereby introducing higher powers of x. In contrast to the spectral GF, the dot occupation and all transport coefficients preserve the proper symmetry. Hence we are confident that the transport properties calculated in this work, and their parameter dependencies, are reliable.
The analytical approach employed here clarifies explicitly that for the description of the model with correlated hopping two characteristic Green functions are needed. One of them defines the transport through the system, and the other the thermodynamic properties, like the ondot density of states or the occupancy of the dot. The intimate coupling between both GFs is realized via the on-dot occupancy of the opposite-spin electrons and various self-energies. This, in conjunction with the general formulae for the currents, clearly shows the various ways the voltage V , the magnetic field B, and the temperature difference ∆T enter the transport characteristics.
However, one should note that within the EOM method not only many of the leading, but also some veryhigh-order contributions-in the sense of perturbation theory-are included, hence an interpretation of the results in terms of low-order processes (which can be quite illuminating, if applicable, see, e.g., Ref. [41]) is beyond reach.
In order to elucidate the role of correlated hopping, and find possible ways to infer its existence (and maybe even its value) from transport experiments, we studied in detail the spectral and transport properties of the QD system. The data presented in Fig. (9) show that the model for x = 0 leads to a conductance symmetric and a thermopower anti-symmetric with respect to δ = 0. A sizeable value of x implies a distinct asymmetry, namely different heights of the conductance peaks and clear changes in the thermopower. Thus we conclude that a thorough symmetry analysis of the consecutive peaks in the conductance, and of the related features in the thermopower, can provide information on the very existence of the correlated hopping term as well as its magnitude. One of the conductance peaks of the single-level quantum dot is not affected by a change in x, while the other decreases in height. The increase of [G(0) − G(x)]/G(0), where G(0) is the conductance maximum at the upper peak and G(x) its value at the lower peak, is faster than linear with x, i.e., faster than the function y = a · x (a = 0.04 for the parameters in Fig. (10)). Measuring the peak variation and calculating the value of a for the known parameters of the experimental setup hence allows for a conservative estimation of x and thus the correlated-hopping term. Note that Γ L also can be estimated from experiment, as it is approximately given by the half-width of the upper conductance peak.
The anisotropy of the couplings affects the conductance and thermopower, as shown in Figs. (11a) and (11b), respectively. This figure exhibits the transport co-efficients for a few values of the asymmetry (i.e., Γ R /Γ L ) for x = 0.5. An increase of this parameter mainly affects the half-width and the height of one of the peaks. In our setup, the impaired peak corresponds to large δ. However, the anisotropy may mask the effect of the x parameter, thus making its unique identification uncertain. However, as discussed above for a strongly asymmetric coupling the simultaneous measurement of the gate bias dependence of the thermopower provides additional information from which the very existence of a non-zero x can be inferred.
The non-equilibrium Green function technique in conjunction with the EOM allows the study of transport coefficients beyond the linear approximation, as exemplified above. In this paper, we have limited ourselves to calculations of the differential conductance: the dependence of G d on x and δ for a number of voltages is shown in Sec. V B. The departures from the linear (small voltage) results depend on x and δ in a complicated way. However, for the studied parameter values the conductance decreases with increasing V , except in the limits of a nearly empty or doubly occupied dot, i.e., at the outer wings of the conductance peaks. This is well visible in panel (b) of Fig. (12) for δ < ∼ −6 and δ > ∼ 6.

ACKNOWLEDGMENTS
The work reported here has been supported by the M. Curie-Sk lodowska University, National Science Center grant DEC-2017/27/B/ST3/01911 (Poland), the Deutsche Forschungsgemeinschaft (project number 107745057, TRR 80), and the University of Augsburg. This sum rule for the correlated-hopping model, which is exact in the wide-band limit, is an important formal result of our paper. Its proof is given in the Suppl. Material [28]. The sum rule (A8) extends that found earlier [26] for the standard single-impurity Anderson model.
Using the above result for the lesser GF, one finds the currents flowing out of the λ electrode as follows: These expressions can be used for calculating the currents in an arbitrary system consisting of the central dot and several terminals. Formally the above manipulations are similar to those arising in the calculation of the currents in the standard Anderson model [43]. However, here we deal with completely different GFs. Moreover, as we shall see below, to calculate the transport GF one also needs the spectral one. Note that the kinetic and transport coefficients are expressed through the imaginary part of the transport GF only. Due to this fact, we shall denote the imaginary part of the transport GF as "transport density of states".

Appendix B: Calculation of the transport Green function
Before calculating the relevant GF, let us note the following identities: which will be occasionally used in various formulae below. The above identities show that the point x = 1 is a special one. Indeed, the model at hand is symmetric with respect to x = 1 for 0 ≤ x ≤ 2. For x = 0 the only hopping is that of single-particle type V λkσ , while for n −σ = 1 and x = 2 one gets the effective hopping equal −V λkσ . This together with the fact that V λkσ enters all formulae as |V λkσ | 2 explains the equivalence of the model at these two limiting points. We remark in passing that a similar change of sign of the effective hybridization is also observed in the periodic Anderson model [44]. The symmetry of the present model goes beyond these two points, x = 0, 2, and is valid for arbitrary x ∈ [0, 2] as discussed earlier [21]. It has to be stressed that within the present approach only the transport GF fulfils this symmetry, while the spectral one does not, as discussed in Sec. IV.
To find the transport GF, we apply the EOM technique to two-time GFs and perform the appropriate decoupling. The quality of the solution in this method depends on the decoupling procedure. Before proceeding, let us recall that the decouplings in the EOM technique typically are not well controlled. We shall benchmark the proposed approximation scheme by checking the symmetry of the solution with respect to changing x ↔ 2 − x. We start with the calculation of the transport GF, but as will be evident higher-order GFs are needed to solve the system of equations. The coupling between various GFs is provided by some correlation functions, inter alia including the average occupation of the dot nσ .
In Zubarev notation [33] for fermionic operators A and B, the equation for the two-time GF written in frequency ω space reads Application of the above EOM to operators A = D σ and The EOM for the next GF, allows one to write The factor in front of the GF on the rhs of the preceding equation defines the self-energy: In the wide-band limit one approximates (B9) by its imaginary part: typically assumed to be energy independent. The higher-order GF which multiplies U in Eq. (B6) reads Equation (B11) suggests that the calculation of the transport GF, independently of the forthcoming decouplings, requires the knowledge of nσ and thus of the spectral GF, dσ|d † σ ω . The remaining GFs entering the rhs of Eqs. (B6) and (B11) fulfil With the auxiliary notation We shall not calculate the GFs containing two c λkσ operators but approximate the GFs in question, avoiding the appearance of functions which describe spin-flip pro-cesses. Thus we project higher-order GFs as follows: and (B24) As already alluded to, the above decouplings are analogous to those which in the context of standard Hubbard model are known as Lacroix decouplings [34]. The decoupling of the following GF, introduces a novel GF which has not appeared hitherto, namely d σ |D † σ ω . To obtain this one, we use an exact relation: deduced from the operator identity D σ = d σ − xd σ nσ. If the GF at hand is multiplied by (1 − x), we can express it by the functions appearing on the lhs of Eqs. (B18) and (B19): (B27) This closes the system of equations for D σ |D † σ ω , except that we still need the spectral GF to calculate nσ . It is worth noting in advance that the spectral GF d σ |d † σ ω turns out to be coupled back to the transport one and the function nσd σ |D † σ ω . The solution of (B18) and (B19) is a relatively easy task. First, using the presented decouplings, one calculates the parameters S tr n , S tr d , and S tr c . For S tr n one finds where vanishes in the wide-band limit and for energy independent coupling, the approximation assumed to be valid here. Thus we end up with where we have omitted the frequency dependence of the self-energy Σ 0σ (ω). Occasionally we shall use this convention in the following. Let us note that the proposed decouplings of the GFs containing two operators describing the electrons on the leads provide a simple expressions for S tr n,d,c in terms of D σ |D † σ ω and nσD σ |D † σ ω only. Later on we shall apply analogous decouplings to find the spectral GFs d σ |d † σ ω and nσd σ |d † σ ω . The remaining two auxiliary parameters read where the novel symbols denote the various "summed" correlation functions or self-energies. For example, Using the definition of the operator D σ , the last correlation function can be split as, e.g.,b 1σ (ω) = b 1σ (ω) − xN 1σ (ω), with obvious definitions of b 1σ (ω) and N 1σ (ω). The other self-energies entering (B31) and (B32) are defined as where the energies ε 1(2) are defined as The symbolsγσ 1 andγ 2 refer to the inverse lifetimes of the singly (doubly) occupied states on the dot. For the standard Anderson model they have been found [26] to play an decisive role in assuring the proper Kondo behavior at low temperatures. They are also important here due to the same reasons. The correlation functionb 2σ (ω) is given by the following combination: Some of the self-energies are expressed by the transport and other by the spectral GF as shown in the Suppl. Material [28], but to calculateb 2σ (ω) one needs nσD σ |D † σ ω . On the other hand, the calculation of b 2σ requires the knowledge of the closely related function nσd σ |D † σ ω . Thus the whole set of GFs required to solve the self-consistent set of equations comprises functions of diagonal: D σ |D † σ ω , d σ |d † σ ω , and offdiagonal: nσd σ |d † σ ω , nσD σ |D † σ ω character. The solution of Eqs. (B18) and (B19) for the transport GF is now written in a closed form. Defining the auxiliary function which for x = 0 reduces to I d (cf. Eq. (26) in the main text) [27], and finds and The related GF nσD σ |D † σ ω is given by Using the relations (B2) and (B41), we get the last required GF: The various symbols used above are summarised in App. C, where they are also expressed in terms of the transport GFs and the auxiliary functions like (B47) and (B48). The calculation of the average occupation of the dot nσ requires the knowledge of the spectral GF (25). Also the self-energyb 1σ (ω) requires the knowledge of the spectral GF and the related function nσd σ |d † σ ; cf. Eq. (C10).
Appendix C: Self-energies in terms of the Green functions Here we list for completeness all self-energies entering the solutions expressed self-consistently in terms of the relevant GFs. The self-energyb 1σ (ω) has been obtained in the previous section. It depends on the transport GF only:b The calculation ofb 2σ (ω) requires the knowledge of b 2σ (ω) and N 2σ (ω), which are expressed as and At first glance the calculation of b 2σ (ω) and N 2σ (ω) requires two new GFs. However, it turns out that these quantities enter the formulae for the transport GFs in the combinationb 2σ (ω) = (1 − x) 2 b 2σ (ω) + x(1 − x)N 2σ (ω). Using their definitions and the relation (B26), one arrives atb The remaining self-energies read The following self-energies: and do not depend on the GFs and take on the limiting values Σ 0σ if the inverse lifetimesγσ 1 andγ 2 are positive infinitesimals 0 + . We have already seen the relationb 2σ = (1 − x)b 2σ . Even so it is possible to find the latter from the former, it is much more convenient to calculateb 2σ directly. The result reads and shows that its calculation requires the GF n σ dσ|D † σ r ω . It can be easily obtained from (B47), and is given by (B48).
Other self-energies are expressed in terms of the spectral GF: In a similar manner one finds: The definition D σ = d σ − xnσd σ can be used to express b 1σ in terms of the spectral GF, dσ|d † σ a ω , and the related one, n σ dσ|d † σ a ω .
of the seemingly more involved decouplings, presented in the Suppl. Material [28], leads to an improvement of the results with respect to the x symmetry.

Appendix E: Matrix formulation
As discussed in the Suppl. Material [28], one may have another look at the spectral and transport GFs, stemming from the definition of D σ = d † σ (1 − xnσ). This allows to write and shows that to get both GFs one has to calculate a matrix GF formally consisting of d σ and d σ nσ operators only. The matrix GF reads G σ = φ σ |φ † σ ω where φ = {d σ , nσd σ } T , and T denotes the matrix transpose operation. The knowledge of G σ is enough to get both spectral and transport GF. This formulation seems to be more symmetric compared to that used previously, and thus could, in principle at least, lead to the required symmetry not only of the transport but also the spectral GF.
The details of the calculations and the decouplings are presented in the Suppl. Material [28]. However, we have to stress again that the higher-order GFs appearing in all the entries of G σ (ω) have to be approximated. It turns out that the decoupling I applied to all four components of the matrix GF leads to results which are most symmetric and closest to the direct calculations presented in App. B. This is illustrated in panel (a) of Fig. (6) of the main text for x = 0.5, and ε d = −5, U = 8, T = 0.03. The differences between each set of curves obtained by the direct formulae (marked with I), and those obtained with the help of the matrix formulation (marked with "m"), are small, hardly visible for x = 0.5; but they do depend on x and are largest for x = 2.
Panel (b) of of Fig. (6) shows the spectral and transport densities of states for two x values, namely x = 0 and x = 2, calculated by the matrix method. One notices that in the matrix formulations the transport density of states for x = 2 differs from that for x = 0 also. This is related to the fact that in this method not only the charge density n σ but the full spectral GF enters the formulae for the transport GF introducing some asymmetry. The comparison of the curves for N (E) for x = 0 and x = 2 best illustrate the lack of the required symmetry in the spectral function. One of the main problems here is the complete suppression of the Kondo resonance in g r (ω) for x = 2. While the transport GF for x = 2 is only quantitatively different from that calculated for x = 0, the spectral functions at those two x values are completely different with a minimum (for x = 2) at the chemical potential where the Kondo resonance should appear. The source of the asymmetry is discussed in Sec. IV C. The results shown in Fig. (6b) have been obtained from the matrix formulation, using decouplings of all four GFs analogous to decoupling I. The other two decouplings produce qualitatively similar results. In summary, none of the studied decoupling schemes leads to the appearance of the Kondo resonance in the spectral function for x = 2.