{'fraction_non_alphanumeric': 0.04211179073766102, 'fraction_numerical': 0.009245426614426249, 'mean_word_length': 4.293643688451208, 'pattern_counts': {'":': 1, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transport in Conductors and Rectifiers: Mean-Field Redfield Equations and Non-Equilibrium Green's Functions

We derive a closed equation of motion for the one particle density matrix of a quantum system coupled to multiple baths using the Redfield master equation combined with a mean-field approximation. The steady-state solution may be found analytically with perturbation theory. Application of the method to a one-dimensional non-interacting quantum wire yields an expression for the current that reproduces the celebrated Landauer's formula. Nonlinear rectification is found for the case of a mesoscopic three-dimensional semiconductor p-n junction. The results are in good agreement with numerical simulations obtained using non-equilibrium Green's functions, supporting the validity of the Redfield equations for the description of transport.


I. INTRODUCTION
Open quantum systems present an interesting challenge for theory. One popular approach relies on a quantum master equation (QME) and the Born-Markov approximation [1][2][3][4][5][6] . QMEs can be used to investigate the transport properties in condensed matter systems, for example in spin chains or fermionic or bosonic tightbinding models. The Lindblad equations, a specific kind of Markovian QME, have been frequently used [7][8][9][10][11][12][13][14][15] . The Lindblad equations maintain positivity of the density matrix (all eigenvalues are non-negative), but there are difficulties with the definition of the currents and with thermalization 4,7,13-22 . The Redfield equation is an alternative QME that, despite not respecting the positivity of density matrix for some initial conditions [23][24][25] , may describe the dynamics more accurately 19 since it does not require further approximation beyond the Born-Markov approximation such as the secular (or rotating-wave) approximation. The secular approximation is frequently invoked in derivations of the Lindblad equations but can be problematic under some circumstances 26 and can break down in the limit of large system size, especially for gapless systems. A separate problem for interacting systems is the exponential growth of the many-body Hilbert space and the dimension of the density matrix with system size [27][28][29] . In this paper we treat interactions in a meanfield approximation to avoid this problem.
A parallel approach to QMEs employs non-equilibrium Green's functions (NEGF) [30][31][32] that extend the standard equilibrium formalism [33][34][35] and provides, in principle, a systematic approach to deal with out-ofequilibrium quantum many-body systems. The NEGF approach has been successfully applied to a broad variety of phenomena ranging from electronic transport through semiconductors 36 and nanostructures 37 to ultrafast pump-probe spectroscopies of strongly correlated materials. The combination of NEGF with dynamical mean-field theory (DMFT) provides 38 a powerful tool to describe, for instance, electronic currents through strongly correlated materials or the time-evolution of cold atoms in optical lattices under sudden quenches. The NEGF approach is formally exact. In principle it avoids the Born-Markov approximations. The Kadanoff-Baym equations of motion 31 describing the time-evolution of an open quantum many-body system are causal and make no assumptions on the strength of the interactions in the system nor on the strength of the coupling of the system to the environment. The equivalent Keldysh formulation of NEGF provides a precise description of the steady state of an open quantum system in which the initial electronic correlations are lost. An exact expression for the stationary current through a system coupled to large fermionic leads can be obtained 37 making it useful for benchmarking less computationally intensive methods such as the Redfield approach discussed below.
Here we examine the validity of the Redfield QMEs for the study of non-equilibrium properties, especially transport, of weakly-interacting systems. To do this we formulate a modified Redfield equation (MRE) by dropping the imaginary Cauchy principal value parts (CPVP) of the dissipator. Within a mean-field approximation such MREs yield a closed equation of motion for the oneparticle density-matrix (OPDM) which can be solved. To test the MRE approach we apply it to a one-dimensional quantum wire and to a three-dimensional semiconductor p-n junction, investigate transport analytically and numerically, and show that the results agree well with NEGF methods.
The rest of the paper is organized as follows. A brief introduction to the full Redfield equations is presented in Section II A with details left to Appendix A. Section II B discusses the MRE and the mean-field approximation is made in Section II C. Section II D briefly discusses the numerical method used to find the self-consistent solution of MRE and Section II E discusses the connection between Redfield equations and Lindblad equations. Section III introduces the NEGF method we use. Section IV applies the MRE method to different systems, and includes a comparison to results obtained by NEGF. Conclusions are discussed in Section V.

II. THE REDFIELD EQUATIONS
A. Full Redfield equation In this section we derive the full Redfield equations. For concreteness we consider the specific problem of a one-dimensional tight-binding chain, connected to separate baths (leads) at each end as shown in Fig. 1. The total Hamiltonian including the chain and the baths is: is the non-interacting system Hamiltonian, is the interaction, is the Hamiltonian for the left and right baths, and models the coupling between the system and the baths. Operator c † i,σ (c i,σ ) creates (annihilates) a fermion at chain site i with spin σ, operator f α † λ,σ (f α λ,σ ) creates (annihilates) a spin σ fermion in a bath α eigenstate λ with eigenenergy ǫ α λ,σ , µ i is the on-site potential and V i,σ;j,σ ′ is the interaction matrix element. We have set i α = 1 if α = L and i α = N if α = R. In the following for simplicity we suppress the spin index σ unless specified otherwise.
Treating the combination of the two baths and the chain as a whole, the equation of motion (EOM) for the density matrix is given by: where we have set = 1 and in the following.
To obtain the EOM for the reduced density matrix of the system ρ S = Tr B ρ total , the usual approach is to make the Born-Markov approximation to trace out the degrees of freedom of the baths 3,5 . The approximation assumes that the coupling between system and baths is weak, that the baths are too large to be affected by the system, and that the baths retain no memory of past history. A detailed derivation and discussion may be found in Appendix A and the final result is where L(ρ S ) is called the "dissipator": is the density of states for bath α, and f α (Ω) = 1/(e (Ω−µ α )/kB T α + 1) is the Fermi-Dirac distribution function. T α and µ α are the temperature and the chemical potential of the bath α respectively. The eigenoperators A(Ω) are defined by: where Π(ǫ) = i |ǫ i ǫ i | is the projection operator onto the subspace spanned by all the eigenstates with energy ǫ. When the system is non-interacting, [c i (Ω)] † = λ λ|i c † λ δ ǫ λ ,Ω and c i (Ω) = λ i|λ c λ δ ǫ λ ,Ω , where |λ are the system eigenstates.
Though Eq. (7) and Eq. (8) are derived without further approximation beyond Born-Markov, it is well known that the Redfield equation may violate the positivity of the density matrix in some circumstances and does not necessarily give the Gibbs distribution at equilibrium. Nevertheless, dropping the imaginary CPVP in Eq. (8) yields a modified Redfield equation that does give the desired Gibbs distribution at equilibrium. To see this first note that Application of the Baker-Campbell-Hausdorff identity yields and e βHS c i (Ω)e −βHS = e −βΩ c i (Ω) .
Substituting the equilibrium thermal distribution ρ th = e −βHS /Tr(e −βHS ) into Eq. (8), combined with Eq. (11) and Eq. (12), and dropping all CPVPs, shows that L(ρ th ) = 0 if the bath chemical potential is set to zero, which means the left hand side of the Eq. (7) also vanishes because H S commutes with ρ th . The recovery of the equilibrium thermal distribution thus motivates discarding the CPVP that appears in Eq. (8) as discussed in more detail in Section II B.

B. Modified Redfield equation for the OPDM
Single-particle operators that are bilinear in the fermion creation and annihilation operators are often of greatest interest, and for such observables the full density matrix using Eq. (7) and Eq. (8) is not needed. Therefore, in the following we focus on the OPDM defined as For a non-interacting system coupled to the bath by a single fermion operator as in Eq. (5), the EOM for the OPDM can be closed. Multiplying both sides of Eq. (7) by c † j c i and taking the trace yields: where h = λ |λ ǫ λ λ| is the system Hamiltonian in firstquantized form, ρ λ| is the OPDM of the system with the temperature and chemical potential of bath α, P iα = |i α i α | is the projection operator that projects onto site i α , andĴ α = π λ |λ N α (ǫ λ ) λ| is the matrix related to the spectral density of bath α.F αI = λ |λ F αI (ǫ λ ) λ| andK α = λ |λ P ∞ ∞ dω N α (ω) ω−ǫ λ λ| are two matrices arising from the imaginary CPVP in Eq. (8). In the wide-band limit πN α (ǫ λ ) = J α , and K α is nearly proportional to the identity matrix. The additional Lamb shift Hamiltonian, which is often neglected, is given by H LS = − α P iαK α . The Lamb shift term has its origin in the interaction between system and baths, analogous to the energy shift of electrons due to their interaction with the quantized photon "bath" of quantum electrodynamics. Further neglecting the [P iα , F αI ] term yields a much simpler equation for the OPDM: It is then obvious that the OPDM of the steady state given by Eq. (15) when the system is only connected to one bath α is just ρ (1) α , which implies that Eq. (15) gives correct Gibbs state at equilibrium, consistent with the analysis of Section II A. Equation (15) is the central result that will be used in this paper. The generalization to the three-dimensional case with translational invariance in two of the three dimensions is straightforward and only slight modifications are required. For more details and further discussion see Appendix B.
We now discuss the decision to drop the CPVPs in Eq. (8) and Eq. (14). Ref. 39 observes that the CPVPs come from the renormalization of the system Hamiltonian. By ignoring the CPVPs, the steady state given by Redfield equation due to single bath is the desired Gibbs distribution. In the context of the Lindblad equations, the CPVPs nearly cancel out with the remainder being the Lamb shift (see Section II E) that is often ignored.
We now consider the existence and uniqueness of the steady state solution of Eq. (15). For the steady state, from Eq. (15) we obtain with A = −ih − α J α P iα . This equation is of the Lyapunov, or more generally, Sylvester form. For a unique solution to exist, A and −A † , or h − i α J α P iα and h + i α J α P iα , are forbidden to have common eigenvalues 40 . This is equivalent to requiring The condition is violated when the system has disconnected parts because then the Hamiltonian h can be block diagonalized in the subsystem basis. Numerical solutions of Eq. (16) also become challenging if u i ≈ u * j for some i, j despite the formal existence of a unique solution.
In some instances it is possible to treat J α = J in Eq. (16) as a small parameter and obtain a perturbative solution. The solution should be valid for small systems with large spacing between energy levels such as semiconducting quantum dots; see Appendix D. Insight into the structure of MRE comes from the lowest-order solution to the OPDM in eigenenergy basis: This equation says that the density matrix of the system is simply a superposition of the density matrices set by the different baths.
C. Mean-field approximation For a general system with interaction, it is not possible to obtain closed EOM for the OPDM, the quantum analogue of the classical Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy problem. Following the approach taken in Ref. 41,42 , our strategy here is to treat the system Hamiltonian in a mean-field approximation, as then the system is a collection of independent particles with some self-consistently determined parameters. The full Hamiltonian h is replaced by the mean-field Hamiltonian h MF , and ρ  α (t). The mean-field approximation can also be made from the outset prior to the derivation of the Redfield equation. As the Markov approximation assumes that ρ(t) does not change appreciably over the memory time scale of the baths, the same final result is obtained.

D. Numerical technique
We briefly discuss the numerical technique used to find the non-equilibrium steady-state of the MRE within the Hartree mean-field approximation. The calculation is carried out in the real-space position basis. First, an initial choice for the mean-field Hamiltonian h MF is made. Next h MF is used to obtain ρ (1) α as described in Section II C and the corresponding initial OPDM ρ (1) 1 is found with the use of Eq. (16). The calculation is carried out with a Lyapunov equation solver. The OPDM and the mean-field Hamiltonian are then repeatedly updated using ρ n where 0 < ǫ ≪ 1 is a small parameter that controls the rate of updating, andρ (1) n is the solution to the Lyapunov equation with operator A determined by the OPDM from the previous step, ρ (1) n−1 . The procedure continues until the change in the OPDM and the mean-field Hamiltonian is smaller than a predetermined threshold; at that point a good approximation to the self-consistent solution ρ (1) st is at hand. One possible obstruction can occur when, as already noted in Section II B, u i ≈ u * j at an intermediate step. Such a situation may arise if the on-site potential in space varies too much, and can be avoided by careful choice of the initial Hamiltonian and ǫ.

E. Connection with Lindblad equations
In this section we briefly discuss the connection between the Redfield equations and the more widely used Lindblad equations. The Lindblad dissipator is given by: where L k is called the Lindblad operator. Depending on whether or not the Lindblad operators are local in space, the Lindblad equations are sometimes classified as local Lindblad equations 8-12 or global Lindblad equations 13,14 , and their validity in different circumstances is still under investigation 43,44 .
To derive the global Lindblad equation with Lindblad operator L ∼ c † λ or c λ typically requires the use of the so-called secular (or rotating-wave) approximation and omits all the terms with Ω = Ω ′ in Eq. (A12) because e i(Ω−Ω ′ )t oscillates rapidly in the interaction picture and average to zero at the coarse-grained time scale is the smallest nonzero value. However, the treatment inevitably becomes problematic in the thermodynamic limit of large system size, for example in a metal, as the energy level spacing is infinitesimally small and timescale separation breaks down. Such global Lindblad equations also raise some difficulty in defining the current and it is necessary to include a fictitious current between partitions that are not directly coupled to ensure that the current obeys continuity. To derive the global Lindblad equation from the Redfield equations, Eqs. (7) and (8), operators c iα are replaced by c iα (−Ω) and c † iα by c † iα (Ω) in Eq. (8). Consequently, the projection operators P iα in Eq. (14) are also replaced byP iα = λ |λ | i α |λ | 2 λ| and thus all the off-diagonal elements of P iα in the energy eigenbasis are ignored. The substitution leads to the decoupling of the diagonal and off-diagonal elements of the OPDM in energy eigenbasis, and since [P iα ,F αI ] = 0, the standard global Lindblad equation in terms of OPDM is obtained: where h LS = − αP iα K α is again the Lamb shift Hamiltonian, due to the interaction between system and reservoirs.
The local Lindblad equation with Lindblad operator L ∼ c † iα or c iα can be arrived at by an alternative approximation. For the chain, if the hopping between sites is weak compared to the memory time of the bath τ B , namely t S τ B ≪ 1, to a good approximation the local operators evolve solely due to the local Hamiltonian, namely Inserting this equation into Eqs. (A6) and (A7) after transforming back to the Schrödinger picture, or alternatively replacing c † iα (Ω) by c † iα , c iα (−Ω) by c iα , leads to the following local Lindblad equation for the OPDM: Another way to arrive at this result is to notice that in the simple model bath, the bath relaxation time is of order 1/k B T if the reservoir is in the wide-band limit 19 , and the temperature scale k B T is much larger than the bandwidth of the chain 2t S . In this limit all of the system energy levels may be considered to be the same and Eq. (21) follows directly from Eq. (14).

III. NON-EQUILIBRIUM GREEN'S FUNCTION APPROACH
We review the Keldysh NEGF formalism used later. For concreteness again consider the one-dimensional chain described in Section II A. The main expressions needed to describe transport are described. Similar expressions for the three-dimensional case can be found in Appendix C.
In the stationary state that we are interested in, timedependent Greens functions G(t, t ′ ) are a function only of the time difference t − t ′ , so the Fourier transform can be expressed in terms of a single frequency: G = G(ω). The Green's functions of the chain coupled to the baths needed to evaluate observables can be obtained from Dyson's equation: where G r,a,< (ω) and Σ r,a,< (ω) denote the retarded, advanced, and lesser Green's functions and self-energies re-spectively that describe the coupling of the chain to the baths. The self-energies read: (23) where we again take the wide band approximation for the leads so that J α = πN α (ǫ λ ) is energy independent. The particle current in term of the Green's function is well known and given by: 37 and the occupation along the chain reads: The occupation along the chain can be obtained from Eqs. (22)- (23) and Eq. (25): Within the Hartree approximation we could in principle obtain all of the Green's functions by solving Eqs. (22) and (26) iteratively, but a more accurate and efficient numerical evaluation of Eq. (26) can be achieved 45 by splitting the integration into the equilibrium and nonequilibrium contributions: where the equilibrium part reads: while the non-equilibrium contributions are: Here, we evaluate the equilibrium contribution using Matsubara sums. From the definition of the Matsubara Greens function we find: where G α i,i,σ (iω n ) is the Matsubara Green's function evaluated with the chemical potential of either lead, µ = µ α . The ω n = 0 contribution in the last Matsubara sum is weighted by a factor 1/2. The integrals appearing in the non-equilibrium part of the occupation as well as in the current Eq. (24) are performed using Gaussian quadrature.
It is enlightening to illustrate the connection between MRE and NEGF here. QME can also be formally derived using Keldysh diagrammatic technique, by carefully picking certain relevant diagrams 46 . Here we make the connection from the perspective of Kadanoff-Baym equation 31 . The equation of motion version of the Keldysh formalism, as introduced by Kadanoff and Baym, leads to an equation for G < (t, t ′ ) which reads 36 : in matrix form. Note that in order to solve for G < (t, t ′ ), the time dependence of G r,a (t, t ′ ) is also needed. This dependence can be generally obtained by solving the Dyson equation of motion in integro-differential form for G r,a (t, t ′ ). We now apply the above expression (31) to the same non-interacting one-dimensional tight-binding model used before. Equation (31) above simplifies in the stationary regime in which not only the baths (which are in equilibrium) but also all the Keldysh Green's functions depend only on the time-difference G(t, t ′ ) = G(t − t ′ ). After Fourier transforming the Greens functions and setting t = t ′ we obtain: where we have used the definition: ji . Note that for non-interacting electrons in the absence of baths, the retarded and advanced Greens function can be expressed in terms of the eigenvalues and eigenvectors {ǫ λ , |λ } of the system Hamiltonian H 0 S as: (33) Substituting this expression into Eq. (32) leads to an equation equivalent to Eq. (14) that was obtained within the MRE approach without neglecting the CPVPs (note that we have implicitly assumedK α = 0 in Eq. (23)). Therefore, our MRE approach effectively ignores the direct contribution from the system-bath coupling to G r,a (ω) in Eq. (32). A similar argument shows that the implementation of the mean-field approximation within MRE is equivalent to replacing all the explicit selfenergies in Eq. (31) with the mean-field self-energies plus the system-bath self-energies given by Eq. (23), and only keeping the self-consistent mean-field self-energies in Eq. (22). Note that, due to the wide band approximation, the instantaneous self-energy Σ r (t, t ′ ) ∝ δ(t − t ′ ), which when substituted into Eq. (31) leads to local-in-time dynamics. Non-Markovian dynamics may, however, arise as a result of the edge singularities present in generic fermionic baths. 47

IV. APPLICATION TO SPECIFIC SYSTEMS
A. One-dimensional ballistic transport in non-interacting metallic quantum wire We first study the simplest case of a spinless noninteracting chain of fermions that is connected to reservoirs at each end as shown in Fig. 1. The system Hamiltonian is given by Eq. (2) and for simplicity we set µ i = V i,σ;j,σ ′ ≡ 0, and J α ≡ J. Eq. (15) shows that, as the system approaches the steady state, in real space the OPDM obeys with particle current j equal to where 1 ≤ n ≤ (N − 1)/2. For the uniform chain, we expect that the interior density ρ Similarly we have where we use the fact that a system at equilibrium has no current. Combining Eqs. (35), (36), (37), and (38) we obtain Alternatively we may write Eq. (40) as where Γ(E) = | 1|E | 2 = 2(1−E 2 /4t 2 S )/(N +1), N (E) = N/(2πt S 1 − E 2 /4t 2 S ) is the density of states of the system and we have assumed thermodynamic limit. As the Redfield equations are based on the weak-coupling assumption 20 , care must be taken in using Eq. (41) when J is large. Eq. (39) can be interpreted with Landauer ballistic transport as Eq. (41) is similar to Landauer's formula j L = ∞ −∞ dE f L (E) − f R (E) /2π for only one conducting channel 48 ; however, different energy levels contribute differently to the current and the total current is suppressed by a factor of ∼ J/(1 + J 2 /t 2 S ). This difference originates from the fact that the system-bath interface in our model can reflect electrons that travel from the system to the bath. Fig. 2(a) shows excellent agreement between numerical results obtained by RME, NEGF, and the analytic formula of Eq. (41). As shown in Fig. 2(b), the occupancy determined by MRE and NEGF differs at the boundary of the system but in the bulk of the system the difference between two methods is negligible. These results suggest that RME is reliable for the investigation of transport properties in the thermodynamic limit.

B. Semiconductor p-n junctions under DC bias
Next we consider the three-dimensional face-centercubic (FCC) tight-binding model illustrated in Fig. 3. Hopping is turned on only between nearest-neighbor sites and thus connects sublattice A to B. To open a gap, alternating on-site energies of ±∆ are assigned to the two sublattices. Because the system has translational invariance along the two transverse directions, for each different wavevector k = (k x , k y ) the three-dimensional problem is described by an effective 1-D tight-binding ladder model, as sketched in Fig. 4. The one-dimensional Hamiltonian is (42) witht xy = 2t x cos(k x ) + 2t y cos(k y ) and we have indicated index k implicitly. To model p-n junctions the positive charge background must be included as well as the long-range Coulomb interaction as these ingredients are necessary to model the depletion layer of a standard semiconductor p-n junction. Here we make the Hartree approximation with a mean-field Hamiltonian that ignores correlation and exchange: where n R denotes the charge of the positive charge background at site R, and V RR ′ = V0 |R−R ′ | is the Coulomb interaction. We have restored the index of spin σ for clarity, and take advantage of spin rotational symmetry. For an intrinsic semiconductor without doping, the system is at half-filling, so we set n R ′ = 1 for all R ′ . By tuning n R ′ to be slightly above or below unity we can model either a n-doped or a p-doped semiconductor. Due to the translational symmetry along x and y directions, it is possible to further simplify the Hamiltonian to where i, j are the indices labelling the coordinate of site along z direction, R ⊥ labels the coordinate along x and y direction, n j,σ = n j,A,σ + n j,B,σ /2 and V ij ≈ −4πV 0 |i − j|, which is obtained by approximately replacing R ′ V RR ′ with an integral and ignoring the overall energy shift term. Transforming Eq. (44) to kspace we obtain the effective 1-D interaction Hamiltonian for each different k: where N x and N y are the number of unit cells along x and y direction respectively, s is the index labelling sublattice A and B, and n k,j,σ = n k,j,A,σ + n k,j,B,σ /2. By setting n j differently on the left and right halves of the ladder, a semiconductor p-n junction can be modeled. The corresponding MRE of the OPDM for the 3D system is: Here |k, λ is the eigenstate of h k + h int k , and P k,iα,s projects onto state |k, i α , s .
Note that here we have assumed the generalized spectral density defined in Appendix B to be N α (k; i α , A; i α , A) = N α (k; i α , B; i α , B) = J and vanishes otherwise. Solving Eq. (46) self-consistently yields the OPDM of the non-equilibrium steady state.
The I-V curves of the p-n junction obtained from the essentially exact NEGF and MRE are shown in Fig. 5(a). There is only a small difference between the two approaches, providing support for the MRE approximation. The I-V curves are similar to those of a textbook semiconductor p-n junction and show clear rectification. From Fig. 5(b) and Fig. 5(c), we see that the size of the depletion region changes when the bias changes as expected. However, since there is no disorder or dissipation mechanism in the bulk of the system, the transport is ballistic everywhere except at the p-n interface and the carriers are in a highly out-of-equilibrium state everywhere. This contrasts with textbook descriptions of macroscopic p-n junctions where the carriers are in a quasi-equilibrium state away from the p-n interface and a recombination current is important in the vicinity of the junction. We also note from Fig. 5(c) that the shift in the electrostatic potential at the baths does not exactly equal the external bias. This difference is due to the absence of disorder and dissipation in our model.

V. CONCLUSION
In this paper, we derive the Redfield equation of motion for the one-particle density matrix. By ignoring the Cauchy principle value parts (CPVPs) and making a mean-field approximation, we obtain a closed equation (15) that can be efficiently solved numerically. Application to the one-dimensional metallic quantum wire system shows that the approach captures the essential ballistic transport physics and leads to the Landauer current formula. One difference with the non-equilbrium Green's function (NEGF) approach is apparent in the occupancy near the boundary with the baths. The difference has two sources: (1) correlations between the system and bath, which are neglected in general quantum master equations, may be important near the system-bath interface; and (2) the CPVPs that we ignore are localized near the boundaries. Nonetheless, for mesoscopic and macroscopic systems, transport appears to be only minimally affected (see Fig. 2(a)). Application of the modified Redfield equation (MRE) to a three-dimensional semiconductor p-n junctions yields a good match between the MRE and NEGF I-V curves with strong rectification as expected (see Fig. 5). There are slight differences between our results and those of textbook p-n junctions, as discussed in Section IV B. MREs serve as an alternative method for the investigation of transport in the inhomogeneous heterostructures 49,50 .
MREs may be better suited to describe transport in mesoscopic systems than the Lindblad equations, as the microscopic derivation of the Lindblad equations requires either a secular approximation, which breaks down, or a high-temperature limit. The range of validity of MREs can be explored by direct comparison to NEGF (see Section III), showing that the two methods are closely related and in fact almost equivalent in the weak systembath coupling limit. For mean-field approximations MRE is more numerically more efficient than NEGF as it avoids summations over frequency and the need to solve complicated integro-differential equations. Positive definiteness of the density matrices is not guaranteed for MREs, however, and positivity needs to be checked on a case-by-case basis.
The MRE approach introduced here seems to be limited, in practice, to weakly interacting systems that are adequately described by a mean-field approach. At the mean-field level, the set of MREs decouple forming a closed system of equations that can be solved. Blind application of MREs to strongly interacting systems would, in principle, require the solution of an infinite hierarchy of equations that couple correlation functions at different orders. This closure problem may be avoided by truncating and solving the hierarchy of MREs through a coupled cluster expansions up to a desired order in the interaction 51 . Truncated Redfield equations have already been used for analyzing charge transport through 1D chains with interacting spinless fermions and thermal conduction through XXZ spin chains 42 yielding promising results. It may be possible to study transport through 3D strongly correlated system systems via the MREs in the limit of infinite dimensions by analogy to the NEGF + DMFT approach which neglects non-local spatial correlations. These theoretical routes are beyond the scope of our present work and deserve further exploration. where is the evolution operator, H S SB is the system-bath coupling Hamiltonian in the Schrödinger picture, T is the timeordering operator, ρ S total is the total density matrix in the Schrödinger picture and initially ρ total (0) = ρ S (0) ⊗ ρ E (0). With Eq. (A1) one can obtain Assuming that the coupling between bath and system is weak and the bath evolves only slowly (the Born approximation) the density matrix can be written as ρ I total (t) ≈ ρ I S (t) ⊗ ρ I E (0), therefore after tracing out the degree of freedom of bath and taking derivative with respect to t, Eq. (A4) becomes where we have used the fact that f α λ B = f α † λ B = 0, and ... B denotes the average value taken with respect to the bath. Here we will assume that the system time scale τ S is much larger than the baths' relaxation time scale τ B so that at time scale τ B ≪ τ ≪ τ S the dynamics of the coarse-grained density matrix of the system becomes local in time, which is known as the Markov approximation. The evolution of the system does not depend its history, and the bath is memoryless. Therefore, we can replace all the ρ I S (t ′ ) by ρ I S (t) in Eq. (A5) and extend the lower limit in the integrand of Eq. (A5) to −∞. Introducing the concept of eigenoperator defined in Eq. (9), one can easily show that Using the fact that and the definition of spectral density N α (Ω) = λ |T α λ | 2 δ(Ω − ǫ α λ ), one can obtain Transforming back to Schrödinger picture one obtains Eq. (7) and Eq. (8).

Appendix B: Generalization of the Redfield equation to three dimensions
For a three-dimensional system, we generalize the system-bath Hamiltonian to: where R α in general denotes the coordinate, and other degrees of freedom (such as orbital) of the state connected to bath α. Similar calculations yield a formula analogous to Eq. (A5): We could define the generalized spectral density and another calculation that includes a transformation back to Schrödinger picture yields the three-dimensional Redfield equation in real space: where for example and a similar definition holds for H α (R α , R ′ α , Ω) and H αI (R α , R ′ α , Ω). To simplify the above equation, we assume that the whole model, including both baths and system, has translational invariance along the transverse direction and there is no interaction. We separate R α into two parts: R ⊥ α which represents the coordinate of the unit cell in the transverse direction, and R α along the junction. Due to the translational invariance in transverse direction, all functions depend only on the difference R ⊥ α − R ⊥′ α and we can define the Fourier transform for a general function g( and its inverse where N x and N y are the number of unit cells in x and y direction respectively. Translational symmetry also means that the energy eigenstates can be labelled by the transverse momentum k ⊥ and the longitudinal index λ k ⊥ . Therefore, we have and Substituting these equations into Eq. (B6) and summing over R ⊥ α and R ⊥′ α yields where Ω = ǫ(k ⊥ , λ k ⊥ ) is implied. We may also obtain the EOM for the OPDM. Ignoring the dependence of N α (k ⊥ ; R α , R ′ α , Ω) on Ω and neglecting all CPVPs, following Section II B we obtain where ρ (1) and P k ⊥ ,α = R α ,R ′ α πN α (k ⊥ , R α , R ′ α )|k ⊥ , R ′ α k ⊥ , R α |.

(B16)
Appendix C: Generalization of NEGF to three dimensions With the same assumption and notations as in Appendix B the Green's functions describing the system can now be expressed as G R ,R ′ (k ⊥ , ω) due to translational invariance in x − y direction. The Dyson equations for the Keldysh Green's functions of the chain coupled to the baths now read: where the self-energies of the chain due to the coupling to the baths are given by where we have again made the wide-band approximation for the baths.
In the three-dimensional model, the occupations along the junction read: The occupations are again split into equilibrium and nonequilibrium contributions. The equilibrium part reads: while the non-equilibrium contribution is: Finally, the expression for the particle current through a unit cell becomes: Appendix D: Perturbative solution to the Redfield equation for the OPDM Consider a system that may have some degeneracies and obeys Eq. (15). If the system-bath coupling J is small, or in other words J ≪ ǫ α − ǫ β for any α = β that ǫ α = ǫ β , we can expand ρ = ρ 0 + Jρ 1 + J 2 ρ 2 + ... and solve the equation order by order. To zeroth order, it is easy to show that the OPDM ρ 0 is block diagonalized in the energy eigenbasis and each block corresponds to a degenerate space. To first order, where we have used the shorthand P = α P iα . We may also write the above equation in block matrix form: and where indices u and w label the blocks corresponding to different degenerate eigenspaces. Working in the eigenbasis that diagonalizes P in each block, as long as P does not have opposite eigenvalues, one can solve for ρ 0 in each degenerate eigenspace leading to: In particular, if there is no degeneracy, Eq. (D4) becomes Eq. (17).