Non-equilibrium transport through a model quantum dot: Hartree-Fock approximation and beyond

The finite-temperature transport properties of the spinless interacting fermion model coupled to non-interacting leads are investigated. Employing the unrestricted time-dependent Hartree-Fock (HF) approximation, the transmission probability and the non-linear $I$-$V$ characteristics are calculated, and compared with available analytical results and with numerical data obtained from a Hubbard-Stratonovich decoupling of the interaction. In the weak interaction regime, the HF approximation reproduces the gross features of the exact $I$-$V$ characteristics but fails to account for subtle properties like the particular power law for the reflected current in the interacting resonant level model.


INTRODUCTION
Out-of-equilibrium quantum systems have received much attention in the past few decades, both experimentally and theoretically [1]. A major goal is to understand the transport of charge and energy through molecular or nanoscale systems coupled to reservoirs such that voltages and temperature gradients can be applied. While there exist powerful numerical and analytical methods to calculate ground state and finite temperature properties of isolated interacting quantum systems, the situation becomes much more involved when these systems are coupled to reservoirs and driven out of equilibrium, even in the case when a stationary state is reached [2]. Due to these difficulties most of the previous studies have been restricted to single site models like the spinless interacting resonant level model (IRLM) or the single impurity Anderson model, and the main focus was on the zero temperature I-V characteristics, in particular in the linear regime. A variety of methods, both numerical and analytical, have been applied like, e.g., the time-dependent density matrix renormalization group (tdDMRG) [3][4][5], non-equilibrium Green's functions [6,7], the time evolving block decimation method [8,9], iterated summation of the path integral [10], and time-dependent density functional theory (tdDFT) [11][12][13][14][15]. While in the Meir-Wingreen approach [6] the problem is formally solved, for interacting models it is generally not possible to evaluate the various Green's functions needed as input without further approximations. In the purely numerical methods there exist severe limitations with respect to size and dimensionality of the systems that can be studied, and even for single site models the approaches are computationally very expensive.
In the present study we utilize the Hartree-Fock approximation in order to calculate the I-V characteristics of a spinless fermion model that can be seen as a generalization of the IRLM to several sites. The obvious advantages of this method are the relatively low computational costs compared to numerically exact approaches, and the great flexibility with respect to dimen-sionality, system size and geometry, finite temperature, and type and range of interactions. On the other hand, it is well known that HF calculations for isolated systems in the ground state or in thermal equilibrium tend to overestimate the appearance of spurious broken symmetry phases. Therefore it is most important to benchmark the results against exact solutions. The purpose of this study is to assess the reliability of the HF approach in the non-equilibrium setting in comparison with available exact results.

Spinless fermion model
We consider a one-dimensional model of spinless fermions, where N C central sites (the molecule) with nearest-neighbor interaction U are coupled to a left and a right lead of N L and N R non-interacting sites. The hopping parameter t 0 is, for simplicity, chosen to be the same in the leads and within the molecule, while the coupling between the leads and the molecule is described by a hopping parameter t and a nearest neighbor interaction U . The Hamiltonian readŝ witĥ  wheren l =ĉ † lĉl . The sites a + 1 and b are the first and the last sites of the molecule, respectively (see Fig. 1), and the system size is N = N L + N C + N R . The charge currents between the leads and the molecule can be obtained from the equation of motion in the Heisenberg picture. Assuming that each fermion carries the charge e, the charge current operator from the left lead into the molecule is given bŷ whereN L is the number of particles in the left lead, and is set to one. Correspondingly, the operator of the current from the molecule into the right lead readŝ Due to particle number conservation, in a stationary state the expectation values of left and right currents are identical, Î L = Î R . In order to probe the transport properties of a system we drive it out of equilibrium by applying a voltage and/or a temperature gradient. There are two different schemes that have been proposed in the literature [16,17] to simulate a non-equilibrium condition. In the first setting the molecule is initially considered non-interacting and decoupled from the leads, and each of the three isolated subsystems is assumed to be in grand canonical thermal equilibrium at its own temperature and chemical potential. At time t = 0 the interaction is switched on, and the molecule and the leads are connected by addinĝ H LC andĤ CR . The time evolution of the system is then governed by the full Hamiltonian. For sufficiently long leads one expects that after a transient phase a stationary state arises where the current is time-independent, and the goal is to calculate this steady state current as function of the applied temperature and voltage bias.
In the second scheme one considers a situation where initially the leads and the molecule are coupled and in thermal equilibrium at the same temperature and chemical potential. At time t = 0 a voltage V is applied by adding different potentials ±eV /2 to the leads. Subsequently, the time evolution is governed by the full Hamiltonian including the voltage bias. The stationary currents thus obtained coincide with the ones of the first quench scheme as long as the applied potential is much smaller than the band width of the leads [18]. However, for larger voltages there are significant differences in the I-V characteristics. While in the partitioned scheme the current approaches a finite value in the limit V → ∞, in the second scheme it goes to zero when the potential difference exceeds the band width.
Since it is not clear how to realize a well-defined temperature bias in the second setting, and since we want to retain the possibility to treat both voltage and temperature gradients we use the first scheme for our calculations. The merits and shortcomings of both quench schemes have been extensively discussed in the literature [5,18,19].

Non-interacting system
Generally, the state of a quantum mechanical system is described by the density matrixρ. In our case each subsystem α = L,C,R is initially in thermal equilibrium with inverse temperature β α and chemical potential µ α , and the density matrix readŝ where Z is determined from the condition Tr{ρ(0)} = 1. The time-evolution of the density matrix is determined throughρ In order to calculate time-dependent expectation values of observables, like the current, one has to compute the equal-time Green's function defined as For non-interacting systems, i.e., when the Hamiltonian H contains only bilinear combinations of the Fermi operators, G can be obtained from the equation of motion where the matrix H is defined bŷ Equation (8) is valid also for time-dependent but noninteracting Hamiltonians. Expectation values of arbitrary products of Fermi operators can be expressed in terms of the Green's function G lm by applying Wick's theorem. Therefore, once the matrix H is numerically diagonalized, it is straightforward to calculate timedependent charge currents for arbitrary times.

Hartree-Fock approach
When interactions are taken into account the task of determining the time-evolution of the density matrix becomes much more involved. Exact solutions for interacting systems out of equilibrium are extremely rare and limited to very special points in the parameter space of the considered model, e.g., the self dual point U = 2t 0 of the IRLM [20]. A rather simple and physically intuitive approximation arises from the Hartree-Fock decoupling of the interaction terms in the Hamiltonian, As a result, the time-evolution of the density matrix is governed by a non-interacting but time-dependent HamiltonianĤ HF , and the equation of motion for the one-particle density matrix reads This equation of motion can be solved numerically with arbitrary precision using the short-time propagation and choosing sufficiently small time steps ∆t. So far, we have discussed the time-evolution of the system starting from a given initial state in order to extract the transport properties from the stationary state that emerges in the course of time. As an alternative, one can use an approach that allows one to calculate the stationary currents directly without considering explicitly the time-evolution. Formally this is achieved by the Meir-Wingreen approach [6,7], based on the nonequilibrium Green's function formalism. For an effectively non-interacting system likeĤ HF their result for the I-V relation is equivalent to the Landauer formula [21,22]. It has been shown [23] that the stationary state Green's function G HF can be expressed in terms of scattering states as where |kα is the eigenstate of H HF that corresponds to an incoming wave from lead α = L,R with wave number k, and f α ( k ) = (e βα( k −µα) + 1) −1 with k = −2t 0 cos k is the Fermi function that accounts for the thermal occupation of the states within each lead. Explicitly, up to a normalization constant the plane wave scattering state originating in the left lead is given by and k = (2π/N L )ν with 0 < ν < N L /2. Here we have assumed for simplicity that the wave functions in both leads are continued periodically up to infinity, and that N L = N R . Numerically it is straightforward to connect the plane wave ansatz (14) in the left and right lead by solving the Schrödinger equation for the HF Hamiltonian within the molecule. Since H HF itself depends on the equal-time Green's function, Eq. (13) has to be solved self-consistently. Numerically this is achieved through iteration. Starting with some reasonable guess for G HF one calculates H HF using Eq. (10) and from there again G HF via Eq. (13). The procedure is iterated until convergence is reached. From the self-consistent solution the steady state current can then be calculated using Eq. (3). The expression for the stationary current can also be cast into the form of the usual Landauer formula where T ( k ) = |t k | 2 is the transmission probability. In contrast to the non-interacting case, here T ( ) is not only a function of the energy, but depends also on the voltage and the temperature of the leads due to the selfconsistency condition.

Discrete Hubbard-Stratonovich transformation
If the number of interacting sites is small, ideally if there is only interaction across a single link of neighboring sites as it is the case for N C = 2 and U = 0, it is possible to calculate the time-dependent density matrix ρ(t) exactly with moderate numerical effort. The starting point is to write the time evolution operatorÛ as a product of M short-time evolution operators, with t = M ∆t. For small ∆t one may further split each exponential by dividingĤ =T +V into non-interacting contributionsT and interaction termsV . Using the symmetric Trotter breakup, one obtains e −iĤ∆t = e −i(T /2+V +T /2)∆t e −iT ∆t/2 e −iV ∆t e −iT ∆t/2 (17) with an error of O(∆t 3 ). Finally, each exponential con-tainingV can be replaced by a sum over an Ising variable s = ±1 using Hirsch's discrete Hubbard-Stratonovich decoupling [24] e −iU (n l −1/2)(n l+1 −1/2)∆t = Ising variables, each of them representing the time evolution of a non-interacting system in a time-dependent potential. For not too long leads and a number of time slices M ≤ 10 one can compute these sums numerically without introducing any further source of error. The discrete Hubbard-Stratonovich decoupling has previously been used in combination with an iterative summation scheme in order to calculate transport properties of the single impurity Anderson model [10]. We will use the Hubbard-Stratonovich method later to benchmark the results of the Hartree-Fock approximation. The inset shows the broadening Γ/Γ0 of the peak as function of U compared to the broadening of the exact spectral function of the zero temperature IRLM [28] for the same hopping parameter. The solid curves are polynomial interpolations as guide to eye.

RESULTS
In the following we present numerical results for the model defined in Eq. (1). We restrict ourselves to the cases N C = 1 which corresponds to the IRLM, and N C = 2 which we will refer to as two-site model, for brevity. In the IRLM U is the only interaction parameter, whereas in the two-site model we set U = 0 and vary the interaction U between the two atoms of the molecule.
The interaction dependence of the linear conductance for a two-site model similar to ours was studied in [25] using DMRG. In contrast to our work, in that paper the interaction between molecule and leads was varied, while the interaction on the molecule was kept fixed.
Unless stated otherwise, in all calculations the length of the leads is N L = N R = 100, the inverse temperature is β L = β R = 20/t 0 , and the overall chemical potential of the unbiased system is µ = 0 which corresponds to half-filled band in the leads. The time-dependent current I is calculated as the mean value of left and right currents if they are different. The currents obtained from the plateau value reached in the time evolution are identical to the steady state currents calculated from the selfconsistent HF approach within a relative error typically of the order of 10 −5 , with the exception of the cases where hysteresis occurs in the self-consistent solution.

NC = 1: Interacting resonant level model
In the case where the molecule consists of a single site the Hamiltonian (1) is identical to the IRLM with interaction parameter U . The I-V characteristics of this model for the special value U = 2t 0 is known analytically [20], and is given in closed form by [26] where 3 F 2 is a hypergeometric function, and eV c /t 0 = r(t /t 0 ) 4/3 . The prefactor r ≈ 3.2 is the lattice regularization of the corresponding field theory. In Fig. 2 we compare the I-V characteristics obtained within the HF approximation with the exact analytical result of Eq. (19) for hopping parameters t = 0.3t 0 and t = 0.5t 0 . For small voltages there is a linear relation between current and voltage due to the fact that at half filling the Fermi level is exactly at the resonance, and therefore the conductance is identical to the conductance quantum G 0 = e 2 /h of spinless fermions. Expanding Eq. (19) to leading order in the voltage yields For comparison, the relative deviation ∆ of the HF currents from perfect conduction is shown in the inset of Fig. 2. A linear fit of the data on a double logarithmic scale yields ∆ ∝ V 2 , just like for the non-interacting system, in contrast to the nontrivial analytical result ∆ ∝ V 6 . Perturbative calculations [27] for the IRLM away from the self-dual point U = 2t 0 indicate that there exists a nonvanishing contribution ∆ ∝ V 2 to the backscattered nonlinear conductance although with a reduced weight compared to the noninteracting system. For larger voltages, the exact currents reach a maximum and then decrease slowly but steadily with negative differential conductance. The HF currents, on the other hand, increase up to somewhat larger voltages. Then a sudden drop occurs followed by a range of negative differential conductance quite close to the analytical result. For voltages beyond the band width, eV > 4t 0 , the HF current remains constant while the exact one continues to decrease. The location of the jump is independent of the method how the HF currents are calculated. In particular there is no sign of hysteresis, i.e., solving the self-consistency equations for voltages increasing by small steps by taking the converged result of the previous voltage as starting point for the next iteration leads to the same curves as for incrementally decreasing voltages. The discontinuous behavior, which is obviously an artifact of the HF approximation, can also be observed in the time-evolution of the current displayed in Fig. 3 for t = 0.5t 0 and several voltages close to the jump. While for voltages below the transition (red curves) the stationary state is reached quite fast and the plateau values increase with voltage, there is a quite strong reduction of the stationary current within a very small range of voltages, and it takes much longer for the system to reach the non-equilibrium steady state.
The most striking feature is a pronounced broadening of the central peak with increasing interaction. The zero temperature spectral function of the IRLM has been calculated numerically using the Chebyshev expansion [28], and the broadening of the central peak has been obtained by fitting the data to a Lorentzian, In the inset of Fig. 4 we compare our data with the broadening calculated from the exact spectral function. Up to U ≈ 0.5t 0 the HF results agree reasonably well with the exact ones but for larger interactions there is a quite strong discrepancy indicating the limitations of the HF approximation.

NC = 2: Two-site model
The transport properties of the model (1) in the case N C = 2 are different from those of the IRLM in many respects due to the fact that now the Fermi level lies exactly between the two transmission resonance peaks of the non-interacting system. Therefore, the inclusion of interaction may not only broaden or shift these peaks but also modify the transmission at small energies and thus strongly influence the conductance. Figure 5 shows the HF transmission coefficient of the unbiased two-site model for several values of the interaction U . While for attractive interaction, U = −t 0 , the  Fig. 7. The size of the symbols is larger than the estimated error. The dotted curve is the HF current for U = 1.5t0 and β = 50/t0. resonance peak is shifted to smaller energy values and slightly broadened with respect to the non-interacting one, the opposite is the case for repulsive interaction. Correspondingly, the transmission at zero energy is reduced with increasing U which is expected to result in a decreased conductance. This expectation is confirmed in the I-V diagram (Fig. 6) where the linear conductance is largest for attractive interaction U = −t 0 , and decreases with increasing values of U . Remarkably, at higher voltages there is a strong and for U ≥ 1.5t 0 even discontinuous increase of the current such that all curves nearly collapse in the large voltage regime.
In order to check if the strong increase of the currents in the I-V characteristics is an artifact of the HF approximation, we have calculated the exact time-dependent currents using the Hubbard-Stratonovich approach described in section 2.4 for shorter leads (N L = N R = 20). The number of Trotter time slices M = 9 is chosen such that discretization errors in the data shown in the figure are much smaller than the size of the symbols. To avoid finite size effects due to the discrete spectra of the leads, we have used a somewhat smaller inverse temperature, β = 8/t 0 , for the data shown in Fig. 6 and in Fig. 7. To get an idea about the influence of finite temperatures we also show the HF current for U = 1.5t 0 and β = 50/t 0 as dotted curve in Fig. 6. The deviation from the β = 8/t 0 result is negligible with the exception of the voltage region close to the jump. For small voltages, V ≤ 2t 0 /e, the HF currents agree reasonably well with the exact ones displayed as open symbols of the same color, whereas for larger voltages they are completely off. To elucidate this behavior, the time-dependent HF currents are compared with the exact ones in Fig. 7 for two  voltages, V = t 0 /e and V = 4t 0 /e. For the smaller voltage (upper panel), the HF currents for repulsive interaction nearly coincides with the exact currents and deviate only slightly for repulsive interaction, U = −t 0 . Note that the stationary currents for much longer leads indicated by arrows on the right axis are nearly identical to what one would obtain by averaging over the small oscillations observed in the HF data for N L = 20.
In the case of large voltage (lower panel), the HF currents resemble the exact ones only for short times during the transient phase of the time evolution but converge all to the same current plateau of the non-interacting system. The exact currents on the other hand approach different stationary values depending on the interaction strength U .
It is therefore clear that the collapse of the different curves at large voltages in the I-V diagram of Fig. 6 is an artifact of the HF approach and indicates a fundamen- tal failure of the approximation in this parameter range. This observation is in contrast to what we found for the IRLM where the HF data were in reasonable agreement with the analytical results both for small and for large voltages, see Fig. 2.
Interestingly, the discontinuity of the currents appearing in the I-V diagram for sufficiently strong interaction and low temperature can be associated with the existence of multiple solutions of the self-consistent HF equation. In Fig. 8 the currents obtained from the iterative solution of the HF equation and the stationary current taken from the time-evolution are shown in a small range of voltages close to the jump. Performing the iteration for voltages that increase by small steps one obtains an I-V -curve that jumps from a low to a high current branch at a voltage of V c1 ≈ 3.03t 0 /e while in the opposite direction the jump from high to low currents occurs at V c2 ≈ 2.77t 0 /e. The true transition voltage obtained from the time-evolution lies in between at V c ≈ 3t 0 /e. The existence of multiple solutions in the I-V characteristics within the HF approximation and within the adiabatic local density approximation of density functional theory has recently been discussed for a model with Hubbard-type interaction [29].

CONCLUSION
The time-dependent HF approximation is a computationally cheap and versatile approach to calculate the I-V characteristics of weakly correlated systems at finite temperatures. The time-evolution of the currents until a plateau value is reached as well as an iterative solution of the self-consistent HF equations for the stationary state yield identical results with comparable numerical effort. However, the self-consistent approach sometimes allows for multiple solutions which leads to hysteretic behavior when the voltage is varied adiabatically. This ambiguity can be avoided using the stationary current obtained from the time-evolution approach. For a model of interacting spinless fermions, the HF data agree well with available exact results, with the exception of the large voltage regime of the two-site model where a spurious discontinuous transition is observed within the HF approximation. It is straightforward to generalize the HF approach in many respects. In addition to the charge currents also energy or heat currents can be calculated which is of particular interest when besides the voltage there is also a temperature gradient. Furthermore, dynamical properties, e.g., the response to a time-dependent gate voltage, can be studied without significant additional effort.