Non-equilibrium Green's functions and their relation to the negative differential conductance in the interacting resonant level model

TL;DR

This paper investigates the non-equilibrium spectral properties and current-voltage behavior of the interacting resonant level model, revealing how spectral changes relate to negative differential conductance and impurity decoupling at high bias.

Contribution

It provides nonperturbative, accurate Green's function calculations for the IRLM under bias, linking spectral evolution to NDC and impurity decoupling mechanisms.

Findings

Spectral properties change drastically between low and high bias regimes.

Negative differential conductance correlates with impurity decoupling from leads.

Distribution functions show anomalous evolution near the impurity.

Abstract

We evaluate the non-equilibrium single particle Green's functions in the steady state of the interacting resonant level model (IRLM) under the effect of an applied bias voltage. Employing the so-called auxiliary master equation approach, we present accurate nonperturbative results for the non-equilibrium spectral and effective distribution functions, as well as for the current-voltage characteristics. We find a drastic change of these spectral properties between the regimes of low and high bias voltages and discuss the relation of these changes to the negative differential conductance (NDC), a prominent feature in the non-equilibrium IRLM. The anomalous evolution of the distribution function next to the impurity shown by our calculations suggests a mechanism whereby the impurity gets effectively decoupled from the leads at voltages where the NDC sets in, in agreement with previous renormalization group approaches. This scenario is qualitatively confirmed by a Hartree-Fock treatment of the model.

Figures

Equations

$$H_{\text{IRLM}}=$$

$$H_{\text{L}}=$$

$$H_{\text{R}}=$$

$$H_{\text{dot}}=$$

$$+U\sum _{r=\pm 1}\left(c_r^{†}{c}_r-\frac{1}{2}\right)\left(c_0^{†}{c}_0-\frac{1}{2}\right)$$

$$I(V)=\frac{V}{2\pi }{}_2{F}_3\left[\left\{\frac{1}{4},\frac{3}{4},1\right\},\left\{\frac{5}{6},\frac{7}{6}\right\};-{\left(\frac{V}{V_c}\right)}^6\right]$$

$$I(V\ll V_c)\approx \frac{V}{2\pi }\left[1-\frac{24}{170}{\left(\frac{V}{V_c}\right)}^6+O{\left(\frac{V}{V_c}\right)}^{12}\right]$$

$$I_{r,r+1}$$

$$j(\omega )$$

$$\text{Im}{G}_{r{r}^{\prime }}^K(\omega )$$

$$\text{Im}{G}_{rr}^K(\omega )$$

$$I_L(V)\equiv I_{-2,-1}=2\pi (J{)}^2\int d\omega$$

$$\times$$

$$J^{\prime }⟶{\mathcal{J}}_{\pm }=J^{\prime }+U⟨c_{\pm 1}^{†}{c}_0{⟩}_{HF}.$$

$${\tilde{f}}_{-1}(\omega ;V)$$

$$\alpha (\omega ,|\mathcal{J}{|}^2(V))=|\mathcal{J}{|}^4{\left|\frac{1}{{\omega }^{+}-|\mathcal{J}{|}^2\frac{1}{{\omega }^{+}-J^2{g}_{TB}^R}}\frac{1}{{\omega }^{+}-J^2{g}_{TB}^R}\right|}^2,$$

Full text

I Introduction

Transport through nanodevices such as molecular junctions or quantum dots has become of great interest in the past due to the potential application of these systems as new type of electronic componentsra.13 ; ko.sm.14 . Generically, the working principle of such components is entailed in their current/voltage (I/V) characteristic. In some situations this can display non-monotonic behavior, usually referred to as negative differential conductance (NDC), a peculiar effect that is intriguing by itself but also most useful in potential applicationsce.re.1999 ; ce.wa.2000 ; kr.ir.02 ; xu.su.1999 ; hu.wu.14 . Therefore, a thorough understanding of the NDC is highly desirable.

A prototypical model exhibiting a NDC is the so-called interacting resonant level model (IRLM), a simplistic model featuring a two level quantum dot connected to leads used to study the interplay of quantum fluctuations and electronic correlations in the setting of quantum impurity problems. Introduced by Vigman and Finkelsteinwi.fi.78 in the (equilibrium-) context of the Kondo-problem, the IRLM in non-equilibrium has received increasing attention over the last decade after the discovery of an analytic expression for the I/V characteristicbo.sa.08 in the so-called scaling regime and for a special value of the interaction, referred to as the self-dual point of the IRLM.

Previous works on the IRLM in non-equilibrium extended the analytic treatment of the self-dual pointgo.sc.18u , considering also higher order statistics of charge transportbr.bo.10 ; ca.ba.11 , and provided further validation by numerical treatment of increasing accuracybi.mi.17 . Away from the self-dual point, Perfetto et al.pe.st.12 studied the transport properties of the IRLM employing non-equilibrium Green’s functions (NEGF) focusing on the effect of long-range interactions. In addition, a perturbative treatment within NEGFvi.sc.14 as well as renormalization group (RG) approachesan.me.10 ; ke.me.13 ,valid for weak interactions, provide further insight for small interactions. In particular, it is found that the NDC within RG arises due to a renormalization of the hopping rate into the leads which gets suppressed for higher voltagesbo.vl.07 ; ca.an.10 ; ka.pl.10 . In contrast, less is known about the physical mechanism of the NDC at the self-dual point, i.e., for intermediate values of the interaction. Related to the NDC, but also very interesting in itself, is the spectral function of the IRLM which in equilibrium was numerically studied by Braun and Schmitteckertbr.sc.14 but, to our knowledge, has not been considered so far in a non-equilibrium situation within a non-perturbative treatment.

In this work, we evaluate NEGF of the IRLM in order to investigate their connection with the NDC and how the spectral and effective distribution functions evolve in terms of the bias voltage. Our results are obtained within the so-called auxiliary master equation approach (AMEA) - a numerical method to treat non-equilibrium quantum impurity problems and evaluate their NEGF with considerable accuracy. For simplicity, our calculations refer to the self-dual point, but can be readily carried out for other values of the interaction. Finally, we complement our discussion of the AMEA results with an Hartree-Fock (HF) treatment in order to help with the interpretation.

We find that in the regime of the NDC, the spectral function evolves from a peak at finite frequencies into a dominant central peak and that the NDC can be traced back to the behavior of the effective distribution functions on the first lead sites. We interpret this behavior as an effective decoupling of the impurity from the leads, which is confirmed be the HF calculations.

II Model and Method

II.1 Model

The IRLM is a well-known impurity model of spinless fermions. It features an impurity site connected to two semi-infinite tight-binding chains together with a density-density interaction term coupling the impurity site to the neighboring chain sites, see Fig. 1. The Hamiltonian is defined as,

[TABLE]

where $c_{r}^{\dagger}/c_{r}$ denote the fermionic creation/annihilation operators at site $r$. Here, $H_{\text{L/R}}$ describe the semi-infinite tight-binding chains and $H_{\text{dot}}$ introduces the hopping to the impurity as well as the interaction term. A non-equilibrium steady state situation is induced in the system via an applied bias voltage $V$ simulated by shifting the chemical potentials of the leads symmetrically, that is $\mu_{l}=-\mu_{r}=\frac{V}{2}$. We use $J$ as unit of energy and work in units where $\hbar=e=k_{B}=1$.

II.1.1 The continuum limit and the scaling regime of the IRLM

Here, we want to summarize some well known facts about the IRLM in the so-called scaling regime, which are important for the present work. A nice overview in the non-equilibrium context can be found in the recent works bi.mi.17 ; ca.ba.11 ; wi.fi.78 ; bo.sa.08 and references therein.

When the bandwidth of the leads, $W=4J$, is the dominant energy scale in the system the lattice model, Eq. 1, becomes equivalent to its continuum limit111also referred to as wide-band limit, allowing for a field theoretic description. In this scaling regime of the IRLM, the physics becomes universal with the emergence of a Kondo energy scale $T_{B}\sim(J^{\prime})^{\frac{4}{3}}$. The constant of proportionality is the lattice regularization of the corresponding field theory relating results from the continuum limit to the lattice model.

The continuum model can be solved analytically for the special value of the interaction $U_{\text{c}}^{*}=\pi$, which corresponds to $U_{\text{lat}}^{*}\cong 2$ in the lattice model, where the IRLM exhibits a certain self-duality. Most notably, there is a closed form expression for the steady state current at $T=0$

[TABLE]

with $V_{c}=c(J^{\prime})^{\frac{4}{3}}$, where $c\approx 3.2$ 222In more detail, $V_{c}=\frac{\sqrt{3}}{4^{2/3}}\frac{4\sqrt{\pi}\Gamma(2/3)}{\Gamma(1/6)}T_{B}$ and $T_{B}=\tilde{c}(J^{\prime})^{\frac{4}{3}}$ with $\tilde{c}\approx 2.7$ from bo.sa.08 . and ${}_{2}F_{3}(a,b;z)$ is the generalized hypergeometric function. From Eq. 2, it immediately follows that $I/V=f(V/V_{c})$ depends only on the rescaled voltage and thus has a universal form set by the energy scale $T_{B}$. As is best seen by expanding the hypergeometric function up to leading order

[TABLE]

the current is linear for small voltages $V<V_{c}$. The most prominent feature of the current arises for $V>V_{c}$ where the model exhibits a negative differential conductance, see Fig. 2.

II.2 Method

In this work, we use the Auxiliary Master Equation Approach (AMEA) ar.kn.13 ; do.nu.14 ; do.so.17 to investigate the IRLM under the influence of an applied bias voltage. AMEA is a method to treat non-equilibrium correlated impurity problems which is particularly efficient to target the steady state. It is based upon mapping the noninteracting bath onto an auxiliary open quantum system whose dynamics is described by the Lindblad equation. This mapping becomes exponentially accurate by increasing the number of sites in this auxiliary system. This open quantum systems effectively mimics a system with infinite volume, so that one can reliably reach the steady state. Correlation functions are then obtained by time evolution of the many-body density matrix starting from the steady state.

The dynamics of the auxiliary open quantum system can be solved numerically exact by available approaches. Here, we employ the so-called Stochastic-Wavefunctionsda.ca.92 ; br.ka.97 ; br.ka.98 , whose application to AMEA is presented in Ref. so.fu.18u . Within the mapping, the central interacting region $|r|\leq 1$ described by $H_{\text{dot}}$ (cf. Fig. 1) remains unchanged 333We note that the mapping needs to be performed for each bias voltage separately but is independent of the parameters in $H_{\text{dot}}$, i.e. $J^{\prime}$ and $U$. In total the auxiliary open quantum system, thus, consists of $L=3+2N_{B}$ sites, where $N_{B}$ denotes the number of auxiliary dissipative bath levels used to replace the left(right) semi-infinite leads. For details, we refer to previous publications so.fu.18u ; ar.kn.13 ; do.nu.14 ; do.so.17 .

II.2.1 Steady-state Current

The current $I_{r,r+1}$ across a bond connecting site $r$ and $r+1$, which is clearly independent of $r$ in the steady state, can be expressed within the Keldysh Green’s function (GF) formalism as ha.ja

[TABLE]

provided the interaction self-energy is zero across the bond. Here, a capital $G_{r,r}(\omega)$ denotes the local GF of the full system, while the lower case $g_{r,r}(\omega)$ is the one when the system is disconnected at the bond connecting the sites $r$ and $r+1$. A convenient choice is the bond from one noninteracting bath to the interacting region, i.e. $r=-2$ to $r=-1$.

In equilibrium, $V=0$, the Keldysh and retarded GF are not independent and connected by the fluctuation-dissipation theorem, which for the GF’s appearing in Eq.(4) reads

[TABLE]

where $f(\omega)$ denotes the fermi-dirac distribution function. In analogy, one can define an effective local non-equilibrium energy distribution function ${\tilde{f}}_{r}(\omega)$ via the relation (cf. also e.g. ts.ok.09 ; mu.ts.17 ; ha.ji.13 ; ca.ko.12 ; ca.we.13 )

[TABLE]

With Eq.(5), we can express the current from the left lead into the central region as444Which reduces to the well known Landauer formula for $U=0$ involving only the difference of the left/right Fermi-functions. To see this consider Eq.(8).

[TABLE]

where $A_{r}\equiv-\frac{1}{\pi}\text{Im}G_{rr}^{R}$ is the local density of states, $A_{TB}(\omega)$ denotes the DOS of the disconnected left lead, that is the DOS of a semi-infinite tight binding chain, and $f_{L}$ is the fermi-function of the left lead. Here, for convenience, we have indicated any possible dependence on the bias-voltage. In Eq.(6), the frequency integrand contains the difference between the effective distribution function at the first correlated site $r=-1$ and the one deep into the left lead weighted with the corresponding DOS.

III Results

In this section, we present results for the non-equilibrium spectral properties of the IRLM. We are not aware of previous numerically accurate results for the non-equilibrium Green’s function of this model from the literature. We consider the self dual point $U=2$ and compute results for two different values of the hybridization strength $J^{\prime}=0.2$ and $J^{\prime}=0.5$ at a finite temperature $T=0.025$. The size of the auxiliary system, which controls the accuracy of the bath hybridisation function (see  do.nu.14 ; do.so.17 ) is fixed to $L=13$. Both the steady state as well as the Green’s functions are obtained by time evolution by stochastic wave functions, see  so.fu.18u for technical details. In order to illustrate the accuracy of the approach, we first plot the steady state current as function of the bias voltage in Fig.2. Specifically we compare data from the present $L=13$ auxiliary-system calculation with the ones of the more accurate approach of Ref. so.fu.18u , where the current is obtained via an extrapolation for values of $L$ up to $L=19$. The analytic solution of the continuum model at $T=0$ is also shown for comparison. In this paper we use smaller values of $L$ because a full Green’s function calculation for $L=19$ would be computationally too expensive. These results show that also $L=13$ provides quite accurate results555see also refso.fu.18u for a benchmark of the $L=13$ GF in equilibrium. and in particular reproduces the NDC.

III.1 Spectral properties at the central impurity site

Fig. 3 shows the density of states at the impurity site, $r=0$, for different bias voltages 666We only show the result for positive frequencies since even at finite bias the IRLM is still particle-hole symmetric at the impurity site..The equilibrium ($V=0$) system is characterized by a pronounced peak at $\omega=2$. Upon increasing the bias voltage, the spectral weight is removed from the $\omega\approx 2$ in favor of a second peak at zero frequency, which quickly becomes dominant for large bias voltages. At the same time, the equilibrium resonance develops sidebands at $\omega=2\pm V/2$. This effect is more pronounced for the case of low $J^{\prime}=0.2$ since a stronger $J^{\prime}$ broadens all peak features. At large voltages, $V\gtrsim 3.2$ for $J^{\prime}=0.2$, the left satellite merges with the central peak.

Out of equilibrium the fermionic effective distribution function obviously deviates from the Fermi Dirac distribution and acquires an anomalous, position-dependent shape. In Fig. 4, we plot the effective distribution function, ${\tilde{f}}_{r}(\omega)$ defined in Eq.(5), at the impurity site, $r=0$, for different bias voltages. We find, that the latter is dominated by a double Fermi-step, $2{\tilde{f}}_{r=0}(\omega)=f_{L}(\omega)+f_{R}(\omega)$, for small bias voltages and changes drastically its shape for bias voltages where the NDC sets in.

III.2 Sites next to the impurity ($r=\pm 1$) and relation with the current integrands

To make contact with the current integrands, Eq. 6, we now consider the spectral properties on the sites next to the impurity, see also Sec. II.2.1. 777We only show results for the site $r=-1$ as the properties at $r=1$ are connected by particle-hole symmetry. Fig. 5 displays the local density of states for different bias voltages. It shows two main peaks around $\omega=\pm 2$,888We note that for other values of the interaction $U$, the peaks are not fixed at $\omega=\pm U$ nor do they stay at $\omega\approx 2$ (not shown). and a featureless spectrum in between. For both hybridization strengths, $J^{\prime}=0.2$ and $J^{\prime}=0.5$, the peaks become sharper and higher with increasing voltage. In addition, for $J^{\prime}=0.5$ spectral weight accumulates for negative frequencies up to the lower band edge at $\omega=-2$.

A more interesting behavior can be seen in the corresponding effective distribution function for $r=-1$ presented in Fig. 6. Similarly to the central impurity site, a double Fermi-step persists in the linear regime, while for higher bias voltages the effective distribution function becomes more similar to the effective distribution function of the isolated left lead for which all states for frequencies smaller than its chemical potential $\mu_{l}$ are occupied. More specifically, the plateau in the positive frequency region $0<\omega<\mu_{l}$ rises in the regime of the NDC.

To elucidate the effect of the bias dependent spectral and effective distribution functions on the current (cf. (6)), we display in Fig. 7 for $J^{\prime}=0.5$ the difference in the effective distribution functions entering Eq.6 as well as the current integrand999We do not show the corresponding plots for $J^{\prime}=0.2$ since they are qualitatively the same and can be readily extracted from Fig. 5 and Fig. 6. which can be seen to be dominated by the behavior of the effective distribution function. The difference of the effective distribution functions has a Fermi-window form of amplitude $1/2$ for small voltages which is considerably distorted in the NDC regime. For negative frequencies $\omega\lesssim-1$ the amplitude quickly vanishes due to the corresponding states being filled at larger voltages, see Fig. 6 whereas at positive frequencies, the amplitude gets suppressed with increasing bias voltage which technically leads to the NDC. Outside of the Fermi window the difference of the effective distribution functions becomes slightly negative. One should not overemphasize this negative region, since the negative differential conductance does not depend on this101010That is, the NDC is present even if these negative contributions to the current are ignored..

IV Discussion and Interpretation of the results

In order to understand the behavior of the spectral and effective distribution functions presented above, we discuss the probabilities of certain characteristic many-body configuration states on the correlated sites. These are displayed in Fig.8 and ranked according to their energy for zero voltage. Notice that the configurations in each pair are related to each other by a particle-hole+inversion (PHI) transformation 111111We note that while at zero bias the model is independently particle-hole and inversion symmetric, at finite voltage the model is symmetric under a simultaneous PHI transformation. and, thus, have the same probability. In addition, the states $(II_{a})$ and $(II_{b})$ have the same probability at zero bias voltage. The corresponding probability is given by the diagonal terms of the reduced (many-body) steady state density matrix, which is plotted in Fig. 9 as a function of the bias voltage.

One can see that the lowest-energy state, type $(I)$, initially slightly gains weight as the bias voltage is increased. This occurs approximately until the point where the NDC sets in. In the NDC regime, $V>V_{\text{max}}$, the $(I)$ state looses weight and eventually crosses with the state $(II_{a})$ which becomes the dominant state at high voltages. Further, the states of type $(II_{a(b)})$, which are degenerate in equilibrium, get their degeneracy lifted by the bias voltage favoring the $(II_{a})$ state since it is the one showing more occupation on the left in accordance with the chemical potentials, $\mu_{L}>\mu_{R}$. On the other hand, the weight of the highly suppressed high energy states $(III)$ stays roughly constant for all bias voltages.

IV.1 Impurity Spectral function

As discussed above, in equilibrium ($V=0$), the configuration $(I)$ has a large overlap with the ground state. Adding a particle at the impurity site to $(I)$ leads to the state $(III)$. Since the energy difference, in the atomic limit $J^{\prime}=0$, between these two states is $\Delta E=2$, this process can be associated to the $\omega\approx 2$ resonance. The suppression of the $\omega=2$ resonance for higher voltages immediately follows from the loss of the weight of the $(I)$ state, cf. Fig. 9. It remains to explain the development of the dominant central peak for high voltages. In general, a resonance at zero frequency occurs when two low-lying states differing by one particle, at the corresponding cite, are almost degenerate. This is the case for the states of type $(II)$. The development of the central peak is then readily explained by the increased weight of the state $(II_{a})$ at high bias voltages.

IV.2 Negative differential conductance

In Sec.III.2, we discussed that on the level of NEGF’s the NDC at large voltage in the IRLM arises due to the effective distribution function on the site next to impurity resembling the Fermi-function of the corresponding lead. This can be seen as an effective decoupling of the impurity from the leads at large bias voltage.

Refs. sc.dz.15 ; vi.sc.14 showed that the NDC in the IRLM is already obtained at the Hartree-Fock (HF) level. Therefore, it is interesting to investigate if the mechanism leading to the NDC obtained from our results is qualitatively similar to the one in the HF approximation.

IV.3 Comparison with Hartree-Fock

We will not present the details of the HF calculations, but we will only underline the connection to the AMEA results. For an alternative discussion of the NDC arising already within HF we refer to the work of Vinkler-Aviv et al.vi.sc.14 . Within HF for the particle-hole symmetric case, which we are discussing in this paper, the Hamiltonian is the same as the non-interacting one with the only exception that we have a renormalized complex hopping between the central impurity and the $r=\pm 1$ sites:

[TABLE]

The computation of the GF’s can be taken from the $U=0$ case, keeping in mind that the hopping ${\cal J_{\pm}}$ is complex and has to be determined self-consistently. It occurs that the local NEGF’s within HF depend only on $|\mathcal{J}|^{2}$ and the expression for the distribution function on the site $r=-1$ has the form

[TABLE]

where $\alpha(\omega)$ depends on the bias voltage only through $|\mathcal{J}|^{2}$ and is proportional to121212The precise form is

$\alpha(\omega,|\mathcal{J}|^{2}(V))=|\mathcal{J}|^{4}\left|\cfrac{1}{\omega^{+}-|\mathcal{J}|^{2}\cfrac{1}{\omega^{+}-J^{2}g_{TB}^{R}}}\frac{1}{\omega^{+}-J^{2}g_{TB}^{R}}\right|^{2},$

with $\omega^{+}\equiv\omega+i0^{+}$ and $g_{TB}^{R}$ denoting the GF of a semi-infinite tight binding chain. $|\mathcal{J}|^{4}$.

In Fig.10, we display the (scaled) current and the squared effective hopping amplitude as a function of the scaled bias voltage within HF. The HF current is qualitatively the same as in the exact solution but instead of a smooth transition from the linear regime to the NDC it shows a cusp and a sudden drop 131313Reminiscent of a first order phase transition as is typical within a mean field treatment. at $V/V_{c}\approx 2$. The drop in the current is accompanied by a drop in the squared effective hopping, which becomes small for voltages outside the linear regime. This behavior of $|\mathcal{J}|^{2}$ can be interpreted as an effective decoupling of the impurity from the $r=\pm 1$ sites in the NDC regime, consistent with the interpretation of the AMEA results.

In the regime in which $|\mathcal{J}|^{2}$ is small, i.e. large $V$, the impurity is weakly coupled to the reservoirs. Its spectral function, thus, consists of a single central peak. It follows that the spectral function at site $r=-1$ will be given by the DOS of a semi-infinite tight binding chain. In addition, from Eq.8 it is clear that the effective distribution function ${\tilde{f}}_{-1}(\omega)$ will resemble the one of the left lead since $\alpha$ is strongly suppressed. In the opposite case, when $|\mathcal{J}|^{2}$ is not small, ${\tilde{f}}_{-1}(\omega)$ will be close to a double Fermi-step and the spectral functions, independent of $r$, will resemble the DOS of an infinite tight binding chain. This means that $A_{-1}^{(\text{HF})}(\omega)$ changes between two different shapes in the large and small $V$ regions in contrast to the AMEA results. Similar to the AMEA results, the NDC within HF is also caused by the change in the effective distribution function since the spectral density $A_{-1}^{(\text{HF})}(\omega)$ in the NDC regime has more spectral weight inside the transport window compared to the solution just before the cusp in the current.

V Summary and Conclusion

We calculated the non-equilibrium single particle Green’s functions (GF), as well as the (many-body) steady state density matrix, of the Interacting Resonant Level Model (IRLM) in the presence of an applied bias voltage employing the auxiliary master equation approach (AMEA). We find developments of sidebands in the impurity spectral function which transforms into a single peak at zero energy for high bias voltages in the regime of the negative differential conductance (NDC). Further, on the level of the non-equilibrium spectral and effective distribution functions, the negative differential conductance in the IRLM arises due to the behavior of the effective distribution functions at the sites next to the impurity. In more detail, they feature a double fermi-step which persists in the linear regime of the current and resemble their equilibrium form of one separated lead for high bias voltages which we interpret as an effective decoupling of the system for voltages in the NDC regime. Supplementing our results with a Hartree-Fock (HF) treatment makes the decoupling explicit and shows that the spectral features resulting in the NDC are shared by both approaches.

In conclusion, our results suggest, in accordance with previous results for small interactions, an effective decoupling of the impurity from the leads as origin of the NDC in the IRLM also at the self dual point.

Acknowledgements.

We would like to thank Delia Fugger, Irakli Titvinidze, Gerhard Dorn, and Volker Meden for fruitful discussions. A special mention goes to Antonius Dorda who originally developed most of the code used for the numeric computations. This work was partially supported by the Austrian Science Fund (FWF) within Projects P26508 and F41 (SFB ViCoM), as well as NaWi Graz. The calculations were partly performed on the D-Cluster and L-cluster Graz and on the VSC-3 cluster Vienna