Charge redistribution in correlated heterostuctures within nonequilibrium real-space dynamical mean-field theory

We address the steady-state behavior of a system consisting of several correlated monoatomic layers sandwiched between two metallic leads under the influence of a bias voltage. In particular, we investigate the effect of the local Hubbard and of the long-range Coulomb interactions on the charge redistribution at the interface. We provide a detailed study of the importance of the various system parameters, like Hubbard $U$, lead-correlated region coupling strength, and the applied voltage on the charge distribution in the correlated region and in the adjacent parts of the leads. Our results are obtained within non-equilibrium (steady-state) real-space dynamical mean-field theory (R-DMFT), with a self-consistent treatment of the long-range part of the Coulomb interaction by means of the Poisson equation. The latter is solved by the Newton-Raphson method and we find that this significantly reduces the computational cost compared to existing treatment. As impurity solver for R-DMFT we use the auxiliary master equation approach (AMEA), which addresses the impurity problem within a finite auxiliary system coupled to Markovian environments.


I. INTRODUCTION
Correlated systems out of equilibrium and especially electronic transport through heterostructures made from different materials, has attracted increasing interest due to the recent impressive experimental progress to fabricate correlated heterostructures 1-6 with atomic resolution and, in particular, growing atomically abrupt layers with different electronic structures [1][2][3] .
From a theoretical perspective, investigating and understanding the physical processes which govern the behaviour of such systems is a great challenge in the field of theoretical solid state physics. For instance, it was shown that, due to the proximity effect, any finite number of Mott-insulating layers become metallic when sandwiched between semi-infinite metallic leads. [7][8][9][10][11][12][13][14] For such a geometry the effect of impact ionization in periodically driven Mott-insulating layers was studied 15,16 as well as resonance phenomena in a system consisting of several correlated and non-correlated mono-atomic layers 17 . Another challenging aspect of such systems that was investigated is the capacitance of multilayer systems made from correlated materials. [18][19][20] Due to the local Hubbard and long-range Coulomb interaction present in these systems, charge redistribution takes place. [20][21][22] The equilibrium situation was addressed, e.g., in Refs. 21 and 22. In particular, Ref. 21 studied the charge redistribution and the corresponding thermo-electric properties for a metalstrongly-correlated barrier-metal device where the onsite energies of the correlated region are shifted compared to the metals, while Ref. 22 investigated the behavior of the correlated thin film in a transverse electric field. Finally, Ref. 20 considered correlated layers described by the Falicov-Kimball model, where one spin-species is immobile, with emphasis on the non-equilibrium situation arising due to an applied bias voltage.
Here, we investigate a system of correlated layers sand- FIG. 1: (Color online) A Schematic representation of the system consisting of Lc = 4 correlated interfaces (red) sandwiched between two semi-infinite metallic leads (blue). In addition to the local Hubbard interaction, present only within the correlated layers, we also take into account long-range Coulomb forces extending into the leads. We take them into account by solving the Poisson equation in an extended region including part of the lead layers (L lead = 23 for each side).
Here LIL, CIL, and CML stands for lead interface layer, correlated interface layer, and correlated middle layer, respectively.
wiched between two metallic leads, see Fig. 1 for an illustration. Similar to Ref. 20 we take into account longrange Coulomb interactions, but here we use the Hubbard model where both spin-species are mobile. The goal of the current work is to investigate the influence of local Hubbard and long-range Coulomb interactions on the charge redistribution in a non-equilibrium steady state situation produced by an applied bias voltage.
The first two terms of the Hamiltonian (1) describe nearest-neighbor intra-layer and inter-layer hoppings, with hopping amplitudes t z and t zz , respectively. The third term introduces the local Hubbard interactions U z , which are nonzero only for the correlated region. The last term describes the onsite energies, whereby v (0) z is chosen such that we obtain the required bulk filling in the z-th layer with the special case of v (0) z = −U z /2 at half-filling (HF). Furthermore, v z describes the Hartree shift of the onsite energies due to the long-range Coulomb interaction (LRCI) produced by the charge inhomogeneity and has to be determined self-consistently. In contrast to the local Hubbard interaction, the LRCI affects not only the correlated region, but the leads as well. Therefore, we incorporate parts of the leads, namely L lead layers per side, into the region. Here, L lead has to be chosen large enough such that the (self-consistently determined) electron density v z converges to the bulk filling of the leads far away from the correlated region. To summarize, the extended central region contains L = L c + 2L lead layers. The corresponding indices vary from − L−1 2 to L−1 2 . |z| < L c /2 describes the correlated region, L c /2 < |z| < L c /2+L lead corresponds to the left (z < 0) and right (z > 0) leads which we treat explicitly, while |z| > L c /2 + L lead corresponds to the semi-infinite lead layers (z < 0 left lead and z > 0 right lead). Here we note that the chosen labeling convention leads to half-integer indices for even L c considered throughout the manuscript.
We take the Hubbard interaction to be uniform within the correlated region, i.e. U z = U for |z| < L c /2 and U z = 0 on the lead layers (|z| > L c /2). We assume isotropic nearest-neighbor hopping parameters within the correlated region as well as in the leads, respectively. This amounts to the choice t zz = t z ≡ t c for the correlated region (|z| < L c /2) and t zz = t z ≡ t α=l,r for the leads (|z| > L c /2). Finally, the lead correlated region junction (LC-junction) coupling is the same on both We work in units where e = = k b = a = 1, with a denoting the lattice spacing and take t c = 1 as unit of energy.
The non-equilibrium situation is reached by applying a bias voltage V = v l − v r . Here v l and v r are the onsite energies far away from the correlated region (v l/r = v z=±∞ ). Notice that, in general, V is not equal to the difference between the chemical potentials of the leads ∆µ = µ l − µ r due to the contribution from backscattered electrons. 72 .
In order to investigate steady-state properties of our system, we work within the Keldysh Green's function formalism 66,68,69,73,74 and use real-space Dynamical meanfield theory (R-DMFT) combined with the Poisson equation to treat the Hubbard interaction and long range coulomb forces, respectively. The shaded area corresponds to the Poisson loop and we use the vector notation for z-dependent quantities, v, n, Σ(ω) and G loc (ω) introduced in Sec.II C. Moreover, χ Φ and χ ∆ are cost functions for the convergence criteria defined in Eqs. (18) and (19), respectively.

B. Real-space Dynamical Mean-Field theory
Here, we give only a brief overview of the nonequilibrium real-space DMFT approach 12,13,17,[57][58][59][60][61] together with the employed impurity solver, namely the Auxiliary Master Exquation Approach (AMEA) 62,63,71,75 In the non-equilibrium situation, the model remains translationally invariant along the xy plane (parallel to the layers), which allows to introduce the corresponding momenta k = (k x , k y ). Moreover, since the steady state Green's functions depend only on the time difference, it is convenient to transform them to the frequency domain ω.
The Green's function for the extended central region, which consists of L = L c +2L lead layers, can be expressed via Dyson's equation Here, boldface indicate L × L matrices, while γ stands for retarded (R), advanced (A) and Keldysh (K) components. G A and G R are related via G A = (G R ) † , while G K , in general, is independent of G R and needs to be determined separately.
The inverse of the non-interacting Green's function reads Where E z (k) is the dispersion relation for the z-th layer of the the extended central region and describes the hybridization between the semi-infinite leads and the extended central region. g γ l (ω, k) and g γ r (ω, k) denote the Green's functions for the interface layers of the semi-infinite leads disconnected from the extended central region. Their retarded component can be expressed as 26,27,76 where v α=l/r + v (0) α=l/r and E α=l/r (k) denote the onsite energies and the dispersion relation for the left/right lead, respectively. The sign of the square-root for negative argument in (6) must be chosen such that the Green's function has the correct 1/ω behavior for |ω| → ∞. Since the disconnected leads are separately in equilibrium, we can obtain their Keldysh components from the retarded ones via the fluctuation dissipation theorem 73 Here, f α (ω) is the Fermi distribution for the chemical potential µ α and temperature T α .
Finally Σ γ zz (ω) = δ zz Σ γ z (ω) stands for the self-energy matrix, which due to the DMFT approximation is diagonal and k-independent. To determine it, we map each correlated layer z to a (non-equilibrium) single impurity problem (SIAM) with Hubbard interaction U z and onsite energy v z +v (0) z , coupled to a self-consistently determined bath. The latter is specified by its hybridization function obtained as (see e.g. Ref. 17,24) where the local Green's function is defined as To calculate the diagonal elements of the matrices G γ (ω, k) from Eq. (2) we use the recursive Green's function method 16,17,77,78 which we generalize to the present situation of Keldysh Green's functions 17 .
To describe the lattice structure of the isolated layers we use a Bethe-lattice density of state (DOS). Due to this choice, we can replace E z (k) by t z ε and dk (2π) 2 by dερ(ε), where ε is a dimensionless parameter characterizing the energy and ρ(ε) = 1 The corresponding impurity problems are then solved with AMEA which is a state-of-the-art impurity solver particularly suited to address the steady state. AMEA is based upon mapping 62,75 the SIAM to an open quantum system of finite size, which includes one correlated site, N B non-interacting bath sites and two Markovian environments, whose dynamics is governed by a Lindblad master equation. The resulting open quantum system is then solved by numerical many-body techniques such as Krylov-space based 63,71 methods, Matrix Product States (MPS) 79 or so called stochastic wave function algorithm 80,81 .

C. Charge reconstruction
To take into account long range Coulomb forces on a mean-field level, we calculate the onsite energies v z selfconsistently by solving the corresponding Poisson equation It is convenient to adopt von Neumann boundary conditions, which in discretized form amounts to setting the Coulomb potential of the two bulk semi infinite leads equal to the one of the boundary layers of the extended central region : Here c z ≡ 1 ε0εr,z , ε r,z is the relative permittivity of layer z and ε 0 is the permittivity of free space. Moreover is the electron density at layer z obtained from nonequilibrium R-DMFT and n bulk is the bulk electron density, which we set equal to 1 (half-filling) throughout this paper. 82 One way to proceed would be to fix the bias voltage V and in the present particle-hole symmetric case v l = −v r = V /2. In this case, one should adjust the asymptotic chemical potentials µ l and µ r of the leads in order to obtain the correct asymptotic charge neutrality n z→±∞ → n bulk = 1. This is numerically demanding. Another alternative is to carry out the calculations for given µ l = −µ r = ∆µ/2 and update the values of the onsite energies in the semi-infinite leads after each iteration according Eq. 12. The bias voltage is then determined by V = v l − v r a posteriori. Here we follow the second strategy as it is numerically more convenient. In fact, we find that the difference between ∆µ and V is quite small in most of the calculations presented in this paper (1% or smaller), except for weak to intermediate U at large LC-junction coupling, as we will discuss below.
For better readability we introduce a vector notation for the z-dependent quantities, namely Obviously, the elements of Σ are zero outside of the correlated region. The electron densities depend, through G K loc,z (ω) in Eq. (13), on the onsite energies as well as on the selfenergy. The self-energy in turn is, through the selfconsistency in R-DMFT, a functional of the onsite energies and of itself, i.e. Σ = Σ( v, Σ). Thus, we have to solve Eqs. (11)-(13) together with the R-DMFT equations in a self consistent manner.
For a fixed self-energy Σ(ω), we solve Eq.(11)-(13) by formulating it as a root searching problem which we treat by the Newton-Raphson method. To this end, we define the function of which we seek the zero. Following the Newton-Raphson scheme, we expand For the technical details about the discretization of the Poisson equation and the expression for the matrix elements M ij we refer to appendix A.

D. Self-consistency loop
Here, we describe the self-consistency loop used to determine the self-energies Σ(ω) together with the onsite energies v as self-consistent solution to the R-DMFT equations coupled, through the electronic number densities n( v, Σ(ω)), with the Poisson equation, Eq. 11. An illustration of the algorithm is presented in Fig. 2. In short, the iterative solution of the Poisson equation constitutes an inner loop to the R-DMFT self-consistency and is done for fixed self-energies Σ(ω) before the determination and solution of the impurity problems, which is more time demanding.
In more detail, we start with an initial guess of the selfenergies Σ(ω) and onsite energies v. Next, the Poisson loop is performed by calculating the electronic densities n, Eq.(13), and updating the onsite energies according to Eq. (16). These two steps are then iterated until convergence 83 is reached, for which we require where Φ is the required accuracy. For each converged Poisson loop we proceed, with the corresponding onsite energies v, to the R-DMFT iteration which consists of computing the bath hybridization functions, Eq.(8)- (9), and solving the corresponding impurity problems thereby obtaining a new set of selfenergies Σ(ω). The alternate solution of the Poisson equation and the impurity problems is then iterated until convergence of the R-DMFT loop. We quantify the accuracy of the latter by the weighted difference between the hybridization functions of two consecutive loops 84

III. RESULTS
As mentioned in the introduction, the emphasis of the present work lies on the influence of electronic correlations on the charge redistribution in a non-equilibrium situation. To this end, we consider the heterostructure sketched in Fig. 1 which is driven out of equilibrium by an applied bias voltage.
In order to understand the behaviour of the charge distribution, for finite LC-junction coupling (t lc > 0) it is instructive to begin with a qualitative discussion of the expected behaviour in the limit in which the correlated region is isolated from the leads (t lc = 0), but still capacitively coupled to them via the long range Coulomb interaction. In that case, when the correlated region is metallic, i.e. for weak to intermediate Hubbard interactions, the system consists of two capacitors (one at each LC-junction) connected in series. On the other hand, when the correlated region is insulating, i.e. for large values of the Hubbard interaction, it can be viewed as one capacitor with a dielectric material placed between two conducting materials. Applying a bias voltage will cause in both cases opposite charging of the facing surface layers of the lead and the correlated region, which can be viewed as dipole-like layers. For definiteness, we will refer to them as lead interface layer (LIL) and the correlated interface layer (CIL), respectively (see Fig. 1).
We perform calculations for L c = 4 and L c = 40 correlated layers, with a homogeneous local Hubbard interaction U z = U . For L c = 4 (L c = 40), we explicitly consider L lead = 23 (L lead = 30) non-interacting, U z = 0, layers for each lead, to allow for proper charge redistribution in the leads as well. Therefore, in total, the extended central region, where the long range Coulomb interaction is accounted for, contains L = 50 (L = 100) layers. The infinite region outside of this range is treated exactly, whereby we take the charge and the Coulomb potential to be equal to its asymptotic bulk values. This   All calculations are performed at ambient temperature T l = T r = 0.025 and we consider an isotropic Coulomb parameter with the moderate value c z = c = 1.5. Due to particle-hole symmetry, properties of the z-th and (−z)-th layer are connected by a particle-hole transformation. For the self-energies, the relation reads Consequently, we need to calculate the self-energies only for half the system, i.e. z < 0. Finally, all results for L c = 4 are obtained with N b = 6 bath sites in the AMEA, while for L c = 40 we considered N B = 4 due to the increased numerical effort. 85

A. Effect of the bias voltage
First, we investigate the effect of an applied bias voltage for intermediate, U = 4, and strong, U = 8, Hubbard interaction, as well as small (t lc = 0.2) and large (t lc = 1) coupling strengths between the leads and the correlated region.
Our calculations show that at the LC-junction the system still hosts dipole-like layers for small but non-zero LC-junction coupling strengths. Fig. 3(a) indeed shows for t lc = 0.2, that the charge density deviations from halffilling, ∆n z = n z − 1, for the CIL and the LIL have opposite signs and their absolute values increase with bias voltage V (see Fig. 3(b)) for both considered Hubbard interactions. So, similar to t lc = 0, also for t lc = 0.2 LIL and CIL can be viewed as dipole-like layers.
On the other hand the behaviour is qualitatively different for large values of the LC-junction coupling (t lc = 1) and in particular sensitive to the value of the Hubbard interaction. For strong interaction (U=8), we obtain that ∆n z of the LIL and CIL have the same sign and their absolute values increase with the bias voltage. When considering a weaker interaction (U=4) this stays true for ∆n LIL (charge density deviation from half-filling for the LIL), while ∆n CIL (charge density deviation from half-filling for the CIL) shows non-monotonic behaviour and a sign change as a function of the bias voltage. So, in contrast to small values of the LC-junction coupling strength, for large ones dipole-like layers are only present at the LC-junction for weak to intermediate U and low bias voltages.
Let us remind that the bias voltage V = v l −v r and the difference between the chemical potentials ∆µ = µ l − µ r differ from each other. As we have already discussed in Sec. II D, it is numerically more convenient to perform calculations for fixed ∆µ and evaluate V a posteriori.
For weak values of the LC-junction coupling or for large value of U , the difference between V and ∆µ is negligible (1% or smaller). However, there is a significant deviation for the case of t lc = 1 and U = 4, see Fig. 3(c).

B. Effect of the LC-junction coupling strength
We further investigate the effect of the LC-junction coupling strength t lc between the leads and the correlated region. We perform calculations for several values of t lc , fixing ∆µ = 2 and again considering U = 4, 8.
When the LC-junction coupling strength is increased, the current through the heterostructure rises. Thus, we expect that more charge is transferred from the left lead to the correlated region. Indeed our results, Fig. 4   interestingly, for the CIL, ∆n CIL changes sign at some U -dependent value t * lc . Furthermore, we find that t * lc decreases with increasing U and for non-interacting correlated region (U = 0) ∆n CIL is negative for all values of t lc we have considered. From here, it follows then that correlations lead to an earlier disappearance of dipole-like layers with respect to the LC-junction coupling strength. This can be understood by the following: For U = 0, the behavior of the system can be intuitively understood by the Hydraulic Analogy, where a fluid takes over the role of the electric charge and pipes represent wires. In this picture larger t lc translates into a bigger diameter of the "LC-junction-pipe". For the behavior of the LIL, this means that less fluid gets jammed at the interface. When thinking about the behaviour of the left-CIL in the hydraulic picture it is easiest to consider the jam created at the right-CIL, since the two are connected by particle-hole symmetry, which will also get decreased with increasing t lc . This means, that the trends observed in Fig. 4(b) are consistent with the Hydraulic analogy.
Coming back to the reason why for stronger Hubbard interaction t * lc is lowered, we can thus interpret the slope of ∆n CIL (t lc ), for low t lc , to originate from the U = 0 behaviour and thus the value of t * lc is mainly influenced by the starting value ∆n CIL (t lc = 0) which is suppressed by the Hubbard interaction leading to the decrease of t * lc as a function of U .

C. Effect of the local interaction
Finally we investigate the effect of the interaction U for small (t lc = 0.2) and larger (t lc = 1) values of the LC-junction coupling strength. We consider differences between the chemical potentials, ∆µ = 2 and ∆µ = 0.5, respectively. These values are chosen such that for small interactions the opposite charging of the LIL and CIL is most pronounced, see Fig. 3(b). Furthermore, in order to better resolve the charge distribution, we also present results for a system with a larger correlated region (L c = 40), in addition to the case with L c = 4. When studying the charging dependence as a function of U , we should expect that in the limit of large U , ∆n z vanishes for the correlated region, since in this limit any double occupation is extinguished.

Small correlated region (Lc = 4)
First, we discus the effect of the interaction for weak LC-junction coupling (t lc = 0.2) and ∆µ = 2. Fig. 5(a) and 5(e) show that the opposite charging of the interface layers is suppressed by the Hubbard interaction. Further, ∆n z for LIL converges monotonically to some finite value for U → ∞, while for the CIL it converges to 0 as expected. In order to investigate the behavior of the boundary charge, we fit them (for U ≥ 2) with exponential functions (see Fig. 5(e)), namely ∆n CIL = A 1 exp(−A 0 U ) and ∆n LIL = B 2 + B 1 exp(−B 0 U ). The resulting fit parameters are given in the figure caption. Notice that both fits give approximately the same exponent, that is A 0 ≈ B 0 .
For small LC-junction coupling (t lc = 0.2) and increasing interaction strength U , as we already mentioned above, the charging of the LIL and CIL is exponentially suppressed, but these layers still have opposite sign and for any finite U , while being reduced, the dipole-like layers are still there.
On the other hand, this is no longer the case for stronger LC-junction coupling strength (t lc = 1) and ∆µ = 0.5, as can be anticipated based on the results presented in previous subsections. Indeed, from Fig. 5(b) and Fig. 5(f) we can see that ∆n z is non-monotonic for both surface layers and in addition the CIL displays a sign change at U ≈ 5 which approaches zero only for higher values of the interaction. In order to understand this behavior, it is important to recall that the results presented in Figs. 5(b) and 5(f) are performed for fixed ∆µ = 0.5, which corresponds to different bias voltages V (see inset of Fig. 5(f)). When examining Fig. 5(f) more closely, one can see that the shape of ∆n LIL (U ) for U > 4 resembles that of V (U ) from the inset. Moreover, from Fig. 3(b) we know that ∆n LIL (V ) is just proportional to V and almost insensitive to U . Based on that, to exclude the dependence on the bias voltage we plot n(U ) = 0.04V (U ), where the coefficient of proportionality is extracted from Fig. 3(b), see green line in Fig. 5(f). One indeed finds that the behaviour of ∆n LIL for U > 4 is controlled by the V (U ) dependency. We thus expect that the curve of ∆n LIL vs. U for fixed V would continue its downward trend also for U > 4 and converge to some value as in the case of the smaller LC-junction coupling strength t lc = 0.2. In contrast to the behaviour of ∆n LIL , fixing V would not affect qualitatively the behavior of ∆n CIL versus U . As a matter of fact, taking the dependence on V into account, one would expect an even more pronounced maximum in the behavior of ∆n CIL (see red curve in Fig. 5(f)).
To investigate the reason for the different behavior of the system for weak and strong values of the LC-junction coupling t lc and why ∆n CIL changes sign at U ≈ 5, we analyze the local electric field E z = − dv z dz present in the system due to the charge redistribution (see Figs. 5(c) and 5(d)). A comparison between the two figures clearly shows that the electric field always has a minimum at z = 0 for t lc = 0.2, where the charging of the surface layers is always opposite, while for t lc = 1 the minimum turns into a maximum exactly at the values of the interaction where ∆n CIL changes sign and again flattens when approaching its asymptotic value. A different behavior of the system between the regimes of weak and strong LC-junction coupling strengths can be also seen by considering the steady state spectral functions A z (ω) = − 1 π mG R z (ω) (see Fig. 6). For t lc = 0.2 and ∆µ = 2 the spectral function does not show a Kondolike peak at ω = µ l = 1. We attribute this fact to a combined effect of the width of the Kondo-like peak being so small that we are not able to resolve it as well as the substantial bias voltage present in the system leading to decoherence which suppresses the resonance. In contrast, for large values of the LC-junction coupling strength there is a clear Kondo-like peak for the CIL (at ω = µ l = 0.25) up to interactions as strong as U = 10. This is not surprising, because the width of the Kondo-like peak is proportional to t 2 lc and correspondingly the difference between these two cases is O(100) and in addition the considered ∆µ is a factor of four smaller. Fig. 6 also shows the spectral function for the central middle layer (CML) featuring, as expected 86 due to the increased distance to the leads, a less pronounced Kondo-like peak compared to the CIL which is already destroyed for U = 10.
It appears that the Kondo-like peak in the spectral density occurs whenever ∆n LIL and ∆n CIL have the same sign, which indicates that the mobility within the correlated region is small as compared to t lc .

Large correlated region (Lc = 40)
We now want to investigate how far the charging of the interface region extends into a bulk system. To this end, we enlarge the correlated region to L c = 40. Results are obtained with N b = 4 auxiliary bath sites in the AMEA impurity solver. 85 To check convergence of our results in respect to N b we also perform calculations with N b = 6 for U = 4. We obtain that results obtained with N b = 6 cannot be distinguished from one obtained with N b = 4 on the scale of Fig. 7.
At this point it is worth noting that for a metallic material, one would expect that only the surface is charged with an exponential tail into the bulk since the induced charge on the surface will compensate the electric field in the bulk. Indeed, our results for ∆n z and v z , presented in Figs. 7(a) and 7(b) respectively, show that the charging and onsite energies behave as expected and fall off exponentially into the bulk. Further, we find that the corresponding penetration depth for charging, although increasing with U , depends only weakly on U and that this dependence is more pronounced for the onsite energies. Note that the system is still metallic for all values of the interaction U ≤ 10 and the exponential suppression can therefore be attributed to screening. The trend that the penetration depth increases with U can thus be interpreted as less effective screening due to the lower density of states around ω ≈ 0.
As in the previous results for L c = 4 the main effect of the interaction is to reduce the absolute value of the charging at the interface between the correlated and uncorrelated region. As can be seen from Fig. 7(c) the behaviour agrees qualitatively with the ones observed for L c = 4, see also Fig. 5(e). The fact that the exponential dependence on U is not so obvious in Fig. 7(c) can be attributed to the lower accuracy due to the increased numerical challenge to converge the self-consistent equations.

IV. CONCLUSIONS
We addressed the steady-state properties of a system consisting of a multilayer correlated region attached to two metallic leads. The model was solved by nonequilibrium R-DMFT whereby AMEA 62,63,71 was used as impurity solver. We studied the charge redistribution in the system induced by the local Hubbard and the long-range Coulomb interactions in the presence of a bias voltage. We find that its behavior is very different for weak and strong LC-junction coupling strengths, especially for strong local interactions. The influence of U on the lead layers is due to the proximity effect and therefore less pronounced in the lead compared to the correlated region.
Our results indicate that the charges (considered with respect to the bulk value) on opposite sides of the LCjunction, can have equal or opposite signs depending on the system parameters. The case of opposite signs can be interpreted as the formation of a dipole-like layer. In particular, these dipole-like layers are present for small but finite LC-junction coupling strengths. In contrast, for stronger values of the LC-junction coupling strength this is only true for intermediate to weak interactions at low bias voltages. For strong interactions, as well as for intermediate to weak interactions at moderate to high values of the bias voltage, the dipole-like layers are destroyed and the charging of the LIL and CIL have the same sign. The dependence of ∆n CIL on the local Hubbard interaction U is quite peculiar, being exponentially decreasing for small t lc while for large t lc it displays a non-monotonic behaviour and even changes sign as a function of U .
This behavior can be understood from the fact that the dipole-like layers are formed if the charges flow faster out of the transition region than they flow in, i.e. of our results is reported in the three-dimensional plot Fig. III C 2.
Finally, we want to emphasize that the results presented in this work obtained for the Hubbard interaction differ from the ones for the Falicov-Kimball model for large values of the LC-junction coupling (t lc = 1) in Ref. 20. In the latter, the LIL and CIL are always oppositely charged. This indicates, that the sign change of ∆n CIL is not a generic feature of strong local correlations paired with long-range coulomb forces. Rather it is a combined effect of strong local Hubbard interactions together with long-range coulomb forces.