{'fraction_non_alphanumeric': 0.034103317933641325, 'fraction_numerical': 0.007335853282934342, 'mean_word_length': 4.818833496252851, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time-dependent transport in open systems based on quantum master equations

Electrons in the active region of a nanostructure constitute an open many-body quantum system, interacting with contacts, phonons, and photons. We review the basic premises of the open system theory, focusing on the common approximations that lead to Markovian and non-Markovian master equations for the reduced statistical operator. We highlight recent progress on the use of master equations in quantum transport, and discuss the limitations and potential new directions of this approach.


I. INTRODUCTION
The term master equation traditionally refers to differential equations that describe the time evolution of the probability that a given physical system will occupy a set of allowed states. These equations are typically first order in time and often, but not necessarily, linear in the probabilities. Nowadays, the term master equation is used more broadly: in the theory of few-level open systems, it refers to equations that describe the evolution of the open system's statistical operator (usually called the reduced statistical operator or the reduced density matrix) in the presence of coupling with an environment that has a large, perhaps infinite, number of degrees of freedom and is responsible for the irreversibility in the open system's evolution [1,2]. Master equations are sometimes obtained as an "educated guess" (i.e. phenomenologically), but can often be derived from the general framework of the open system evolution and a few reasonable assumptions. Master equations can be divided into Markovian master equations, in which temporal evolution of the reduced statistical operator does not depend on its past but only its current state, and non-Markovian equations, in which the so-called memory effects play an important role and involve information about the evolution of the environment. How to best quantify non-Markovian effects in open quantum systems [3][4][5][6] and how to experimentally control the information flow between the system and the environment, potentially driving the system between the Markovian to the non-Markovian regime [7], are currently very active areas of inquiry.
Electronic systems in semiconductor nanostructures are open quantum systems, exchanging particles and information with the rapidly dephasing reservoirs of charge, often referred to as contacts, and possibly interacting with phonons or photons as well [8]. In the open system theory, environments are commonly considered to be bosonic, which is fine for the interaction of electrons with light or lattice vibrations. However, electronic transport in the presence of contacts is a case of a fermionic open system coupled to fermionic reservoirs, which is a largely unexplored problem [9,10]. When referring to electronic transport calculations, the use of the term * knezevic@engr.wisc.edu master equation falls in two camps: on the one hand, we have few-level models (e.g. the resonant-level models used for quantum dots [11]) for which master equations continue to refer to the dynamics of the reduced statistical operator. On the other hand, if we strive to account for the generally continuous single-particle energy spectrum of an electron in a nanostructure (e.g. when capturing current in structures without resonances) and the fact that many electrons are available to populate them, then calculating the full many-body reduced statistical operator becomes both intractable and unnecessary, as a great deal of information can be obtained from the singleparticle quantities. In this case, master equation can refer to the equations for the time evolution of the single-particle density matrix (e.g. Redfield-type equations [12]) or just its diagonal terms (e.g. the Pauli master equation [13]).
In this paper, we review the basic premises and recent progress on the use of Markovian and non-Markovian master equations in the description of quantum electronic transport. In Sec. II, we present the basics of the open system formalism, including the concept of complete positivity of a dynamical map. We discuss microscopic derivations of Markovian master equations in Sec. III, focusing on the weak-coupling limit, and follow with examples from quantum transport in Sec. IV. In Sec. V, we overview general features of non-Markovian master equations and present the Nakajima-Zwanzig projection operator technique. Examples of non-Markovian master equations for time-dependent quantum transport are given in Sec. VI. A summary and outlook conclude this paper in Sec. VII.

II. GENERAL OPEN SYSTEM FORMALISM
Consider a quantum-mechanical system S interacting with an environment E. This composite S + E system is generally described by the full statistical operator that, like other operators, lives in the Liouville space, which is (in the case of finite systems) isomorphic to the square of the composite Hilbert space, H 2 = (H S ⊗ H E ) 2 . Here, H S and H E are the Hilbert spaces while H 2 S and H 2 E are the Liouville spaces of the system and environment, respectively. Operators acting on the Liouville space are often called superoperators. If S + E is closed, the dynamics of its statistical operator ρ is given by the Liouville equation (in the units of = 1) L is the Liouville superoperator and H the total S + E Hamiltonian, generally of the form H = H S ⊗I E +I S ⊗H E +H int . The integral form of the Liouville equation is with T denoting time ordering. In the case of a time-independent H, U (t, t 0 ) = exp (−iH(t − t 0 )). Assume that we are interested only in the evolution of S. Its statistical operator can be obtained by tracing out the E degrees of freedom, i.e., S is often referred to as the reduced system and ρ S as the reduced statistical operator or the reduced density matrix. We will use the term reduced statistical operator, because the term density matrix is usually reserved for the single-particle quantity in quantum transport studies. While the dynamics of S + E, a closed system, is unitary, the dynamics of S is not. If the environment has a large number of degrees of freedom, on the timescales accessible in experiment the evolution of the reduced system S effectively becomes irreversible. Quite generally, the dynamics of the reduced statistical operator is given by The central goal of the open system theory is to obtain the evolution of the (relatively small) reduced system S while minimizing the information that has to be gathered about the (relatively large) environment E. This quest is understandably very difficult, and approximations must be employed to yield tractable equations. A first major simplification is assuming that, at some point in the S + E evolution, the S and E were decoupled. It is assumed that, up until a certain point in time, usually designated as t = 0, S and E were mutually isolated, non-interacting, and therefore the initial ρ is of the uncoupled, tensor-product form, Thereafter, the interaction is turned on, presumably adiabatically. (A reader interested in the field of nanoelectronics can immediately ask if this assumption is ever satisfied in electronic systems, and the answer is "sometimes." For instance, it can be considered true when we have high tunnel barriers between the active region (S) and contacts (E) in a nanostructure, and have let the active region and environment each relax on its own, with minimal tunneling between them.) The assumption of an uncorrelated initial state is a very appealing one to adopt, because it guarantees [14] the existence of a subdynamics (also known as reduced dynamics), i.e. it guarantees that the evolution of the reduced statistical operator can in principle be fully described within H 2 S . In other words, the existence of a subdynamics means there exists a generally nonunitary evolution superoperator W(t, 0), such that W(t, 0) is often referred to as a dynamical map. In general, there exists a non-Hermitian generator of the dynamical map, K(t), which satisfies The differential equation of motion for a subdynamics can be written as Clearly, it is practically impossible to obtain K and W from first principles, and approximations are commonly made to the structure of K. Generator K should be such that the resulting evolution does not violate the unit trace or the positivity of ρ S . While the former is quite easily satisfied (any generator that acts as a commutator or a sum of commutators will preserve the trace of ρ S ), the latter is generally a tall order and is, in fact, not fulfilled by many common approximations. Complete positivity of a dynamical map [15] is a stronger criterion than positivity (i.e. requiring that the map preserve the positivity of the statistical operator). Namely, if you have two systems whose evolution is such that the density matrix of each remains positive (i.e. the evolution operators W for each subsystem are positive maps), it is not guaranteed that the composite map (their tensor product) will be a positive map, i.e. it is not guaranteed that the composite statistical operator will remain positive throughout evolution. If, however, each one of the evolution maps is completely positive, then the tensor product is completely positive. In essence, complete positivity of the evolution map is a stricter criterion than positivity and necessary for a successful description of composite systems [1].
A time-independent generator K corresponds to Markovian approximations. Evolution operators W(t, t ′ ) generated by a time-independent K form a semigroup, with K then referred to as the semigroup generator. It has been shown by Lindblad [16] that the most general case of a generator of a completely positive Markovian evolution must be of the form (given in the Schrödinger picture): where γ k are nonnegative coefficients. The last term on the right-hand side is often referred to as the dissipator. The dynamical map W(t, t ′ ) generated by K from Eq. (9) is a completely positive Markovian map. Lindblad's form of K is very useful because it enables development of physically reasonable approximate forms of semigroup generators.

III. MICROSCOPIC DERIVATIONS OF MARKOVIAN MASTER EQUATIONS
In the most general terms, completely positive Markovian equations for the subsystem dynamics can be obtained in the weak-coupling limit [17][18][19][20], singular coupling limit [2], and by coarse graining over time [21] (discussed in more detail in Sec. VI A). In electronic systems, there is also the large bias limit with contacts that have a constant density of states [11] (also referred to as the wide-band limit [22,23]), which we will discuss separately. The weak-coupling limit is of particular importance, being applicable to electron-phonon interaction and electron-contact coupling in the case of tunnel barriers, so we discuss it in more detail.

A. The weak-coupling limit
The total dynamics in the interaction picture can be written as (in differential and integral forms) where H I (t) is the interaction Hamiltonian in the interaction picture. Putting the integral form in the right-hand side of the differential form results in dρ(t) . Tracing out this equation over the environment degrees of freedom, we obtain It is commonly assumed that Tr E [H I (t), ρ(0)] = 0. This assumption is often satisfied: for instance, if the initial statistical operators of the system and environment are grandcanonical or canonical equilibrium ensembles, they will contain pairs of the creation and annihilation operators associated with the single-particle spectra for S and E. The interaction Hamiltonian is usually linear in these operators, i.e. it is commonly assumed to be of the form (in the Schrödinger picture) where A α are the system and B α the environment operators. Tracing out the product of this interaction Hamiltonian with the environment statistical operator over environmental states gives zero.
In electronic systems, the above approximation is satisfied for the electron-phonon interaction (the interaction Hamiltonian is linear in phonon creation and annihilation operators) as well as for typical model Hamiltonians that describe the interaction of the device with the contacts (Hamiltonian linear in the contact and device creation/annihilation operators) if the evolution starts from equilibrium. However, if contactactive region electron-electron interaction is deemed important and is part of the interaction Hamiltonian, then the term So the equation we are focusing on, which is still exact provided approximation Tr E [H I (t), ρ(0)] = 0 holds, is: The Born approximation assumes that the interaction is weak, so that the environment is negligibly affected by it and no considerable S-E correlations arise due to it over time on the timescales relevant to S. As a result, we can write Equation (13) has memory. The Markov approximation assumes that the interaction magnitude is such that the evolution will depend only on the present state of the system, not its prior evolution, so ρ S (s) is replaced by ρ S (t). This assumption is valid on timescales coarser than the decay time of environmental correlations [2]. We can switch from s to t − s, with s now denoting the temporal distance from t, and integrate over all values of s, because we expect the integrand to be negligible for large values of s (i.e. environmental correlations decay rapidly with increasing s), finally arriving at Equation (14) is the Redfield equation [24] and it still has memory.
For an interaction Hamiltonian of the form α A α ⊗ B α , we can define the Fourier transforms of A α and B α based on the system and environment spectra, where P(ε) projects onto the eigenspace of H S corresponding to eigenvalue ε. As a result, In that case, the interaction Hamiltonian in the interaction picture becomes H Finally, the evolution of the reduced statistical operator becomes where . When the typical timescales for the system evolution, proportional to |ω − ω ′ | −1 , are much shorter than the expected relaxation timescales for the system, the so-called secular approximation (also known as the rotating-wave approximation or RWA) can be applied: all terms with ω−ω ′ = 0 are considered as varying too fast, so that their average contribution on the timescales relevant to S can be neglected. As a result, we obtain the weak-limit Markovian equation of motion Let us define χ corresponds to the so-called Lamb shift, an effective correction to the system Hamiltonian of the form H LS commutes with H S , so it shares the eigenvectors with H S and simply corrects the H S energy levels, and is therefore not a true dissipative term. γ defines the coefficient of the true dissipator, This dissipator is of Lindblad form, which can be shown after proving that γ αβ is positive definite and diagonalizing it.
If the system Hamiltonian is diagonalized in a basis |n as H S = n ε n |n n|, then we can derive an equation of motion for the populations of the eigenstates ρ S (n, t) = n|ρ S (t)|n as S(n, n ′ ) = αβ γ α (ε ′ n − ε n ) n ′ |A α |n n|A β |n ′ being the transition rates obtained from Fermi's golden rule [2]. Equation (21) is known as the Pauli master equation.

A. Pauli master equation for electron-phonon interaction
An example of the Pauli master equation in the treatment of electron-phonon interaction in devices is the work of Fischetti [13,25]. He has shown that, in the Born-Markov approximation and the van Hove limit (time tends to infinity while the coupling strength tends to zero, so that interaction squared times time remains constant and nonzero during the limiting procedure) [17], the master equation for the fermionic active region will include the exclusion principle, thus generally becoming non-linear for high population of states. Scattering states |µ that diagonalize the single-electron Hamiltonian in the active region can be obtained from the solution of the coupled Schrödinger and Poisson equations with open boundary conditions. In order to accurately compute spatially resolved quantities, such as charge density and potential, in the numerical implementation, an appropriately dense set of scattering states is obtained through a mapping onto standing-wave-type solutions (details can be found in [25]).
In the work, the active region -contact interaction is treated through a boundary injection/collection term that acts as a source to the equation. The Pauli master equation in the basis of scattering states reads , where the source term for contact j is given by Here, f (j) ( k µ,j ) is the distribution function in contact j and υ ⊥ ( k µ,j ) is the perpendicular component of velocity associated with state µ and normal to the active region/contact j boundary. Figure 1 shows a comparison between the Pauli master equation and ensemble Monte Carlo simulation of a silicon nin diode.
There is a concern that the Pauli master equation does not conserve current outside of the steady state. It has been shown that current is conserved as long as coupling to the contacts is local [26].

B. Markovian equations for system-environment coupling
One of the early contributions aimed specifically at the treatment of transport in electronic systems via master equations was the paper by Gurvitz and Prager [11]. In their work, the approximation of high bias has enabled the Markov approximation. They discuss resonant transfer in mesoscopic devices, focusing on resonant states as the only relevant eigenstates of the electronic Hamiltonian in the systems of interest. The resonant level model is commonly adopted [22,27,28]. The open system has two terminals and is coupled only to the left and right reservoirs, such that the resonant levels are comfortably inside the transport window (the range of energies between the Fermi levels of the two contacts) and the density of contact states is constant throughout. Markovian evolution can be obtained in the form of the density matrix in the basis of the resonant states, with off-diagonal terms making it different from the phenomenological rate equations [29]. Around the same time, Stoof and Nazarov [30] investigated time-dependent resonant tunneling via two discrete states in the presence of resonant-frequency irradiation based on a phenomenological Markovian master equation for the full statistical operator of this two-level system. An enhancement to the work of Gurvitz and Prager [11] was put forth by Li et al. [31] for a system with multiple resonances, such as a quantum dot, that is connected to the reservoirs via barriers through which tunneling is relatively weak. Starting from the Born approximation and working with conditional density matrices that correspond to a fixed number of electrons getting onto the dot at a given time, the authors are able to derive a Markovian equation of the Lindblad form that does not require the wide-band limit [11].
The work of Harbola, Esposito, and Mukamel [32] uses projection operators (see more in Sec. V A below) to derive a hierarchy of quantum master equations for the many-body statistical operators representing the system with a given number of electrons. They show that Fock-space coherences between states with different populations do not contribute to transport to second order in system-environment coupling, but coherences between different many-body states with the same n are appreciable.
Espostio and Galperin [33] derived a time-local Markovian master equation for molecular transport based on the Redfield equation, which is nonlocal in time, and supplanting it with a kind of time-reversed Redfield evolution that enables a selfconsistent procedure for deriving the generator.
Pedersen and Wacker [34] worked in the basis of the  Fig. 2). The evolution they describe is numerically tractable and contains considerably more information than the rate equations. The approach is referred to as the second-order von Neumann approach (2vN for short), indicating that the correlations between two tunneling events are included [35]. Based on a diagrammatic expansion, Karlström et al. recently showed the equivalence between the 2vN approach and the resonant tunneling approximation, and discussed the limitations of the technique in the calculation of higher order cumulants [36].

V. MICOSCOPIC DERIVATIONS OF NON-MARKOVIAN EQUATIONS
It is known that the general, completely positive, non-Markovian evolution of an open system that started in an uncorrelated state (5) can be written as This form is usually referred to as the operator-sum representation or Kraus representation [38], where R i (t) are the Kraus operators. Approximate Kraus maps based on physically reasonable assumptions have been constructed [39]. However, in contrast to Markovian evolution, where the Lindblad equation (9) specifies the required form for a generator of a completely positive dynamical map, there are no similarly compact criteria to determine if an approximate non-Markovian map is completely positive or not.

A. Nakajima-Zwanzig and time-convolutionless (TCL) projection operator techniques
A general and widely applied technique for the derivation of non-Markovian master equations up to a given order in the S − E interaction is the Nakajima-Zwanzig projection operator technique [40,41]. Commonly, terms up to the second or fourth order in the interaction are retained, but complete positivity of the resulting master equations is generally not guaranteed.
In the S + E Liouville space H 2 , any environment density ρ E matrix generates a projection operator P whose action is given by Px = Tr E (x) ⊗ ρ E , x ∈ H 2 . P is a projector, meaning that P 2 = P. The range (space of images) of P is isomorphic to the system Liouville space H 2 S . The complementary projector is Q = 1 − P.
By projecting the Liouville equation (1) onto the ranges of P and Q, we obtain two equations of motion If the interaction Hamiltonian, in the interaction picture, is of the form ǫH I (t), where ǫ is a unitless number characterizing the smallness of the interaction, we can formally solve the equation for Qρ as where G(t, s) = T exp ǫ t s ds ′ QL(s ′ ) . Substituting this equation into (25a) above, we obtain the Nakajima-Zwanzig equation: ds PL(t)G(t, s)QL(s)Pρ S (s). (27) Commonly, in the case of an uncorrelated initial state (5), the initial environment density matrix is chosen to generate the projection operator P, which means that Pρ(0) = ρ(0) and Qρ(0) = 0. Alternatively, the projector may be chosen so as to annul the odd moments of the interaction Hamiltonian. The choice of P generally depends on the application in mind, and P is often assumed to be associated with the equilibrium canonical or grand canonical statistical operator for the environment.
The time-convolutionless (TCL) projection operator technique, originally due to Shibata et al. [42], writes the Nakajima-Zwanzig equation in a form that depends only on the instantaneous ρ S (t), and all the memory effects are relegated to certain evolution operators, which opens doors to systematic approximations, even if the operators are still quite unwieldy and a partial trace does technically need to be taken over the equation after everything. (Partial-trace-free timeconvolutionless equations of motion and the related concept of memory dressing have been proposed in [43,44].) Here, we quote the TCL equation in the form without the inhomogeneity, i.e. for Qρ(0) = 0. where Obviously, there is an assumption that 1 − Σ(t) is invertible [43]. Upon performing a Taylor expansion of 1 − Σ(t) in terms of ǫ, we can get a series K(t) = n ǫ n K n (t), where the TCL equation yields the Redfield equation (14). A number of examples of TCL equations with second and fourth order coupling can be found in [2]. Timm [45] discusses a diagrammatic expansion of time-convolutionless equations.
It is also worth noting that a projection operator technique can be used to derive the well-known semiconductor Bloch equations [46].

VI. EXAMPLES OF NON-MARKOVIAN TRANSPORT EQUATIONS IN QUANTUM TRANSPORT
One of the early non-Markovian approaches to electron transport in nanostructures was put forth by Bruder and Schoeller [47]. Time-dependence was introduced either by periodic modulation of the Fermi energy or by time-dependent perturbations to the quantum states in the dot. The authors focused on the effects of the Coulomb interaction in the limit of low tunneling rates but finite level spacing.
Vaz and Kyriakidis [48][49][50] calculated the full Redfield tensor in Fock space for a two-level system (Fig. 3). The authors find that Fock-space coherences between states with different particle numbers are robust and may be preserved even in the presence of tunneling into and out of the dot. The authors also note that, while Redfield dynamics could potentially violate positivity of the statistical operator, they have not observed it in practice [48].
Recently, Gudmundsson and co-authors [51] used a non-Markovian transport equation to analyze time-dependent transport in a few-mode nanowire containing a localized region and focused on the effect of nontrivial geometry. The authors pay attention to the fact that arbitrary decisions where the active region ends and contacts begin lead to inconsistencies, and that an effective overlap between the S and E wave functions will yield effective interaction Hamiltonian matrix Zedler et al. [53] present an interesting analysis of non-Markovian versus Markovian equations in the weak coupling limit on the example of a quantum dot coupled to contacts with a Lorentzian density of states (i.e. contacts with a finite electron lifetime), thereby going beyond the high bias limit, and conclude that one must be careful with non-Markovian master equations as they do not necessarily perform better than their Markovian counterparts when non-Markovian effects are strong, and are not in general guaranteed to conserve positivity. The authors compare the exact solution for a single level system with dynamical coarse graining [54], non-Markovian master equation, and the Markovian master equation limit [53].  [55,56] emphasizes the continuous spectrum in the open active region, with forward-and backward-propagating scattering states, whose asymptotic forms are plane waves (a combination of injected and reflected waves in the incoming contact, transmitted wave in the outgoing one). The model interaction Hamiltonian couples each scattering state only with the plane wave with the same wave number k from the injecting contact, i.e.
c † k,L (c k,L ) and c † −k ′ ,R (c −k ′ ,R ) create (destroy) an electron with a wavevector k in the left and −k ′ in the right contact, respectively, d k and d † k do the same for active-region states, and k ′ 2 − k 2 = 2m * eV / 2 (k and k ′ are the wave numbers corresponding to the same energy in the two contacts separated by bias V ). The hopping coefficients ∆ k and ∆ −k are proportional to the current I k carried by each mode, ∆ k = I k eT k , where T k is the transmission coefficient of mode k [56].
To obtain a tractable theoretical approach, the full dynamics is coarse grained over the momentum-relaxation time of the contacts. Contact relaxation occurs on timescales of order 10 1 − 10 2 femtoseconds [57,58], owing to fast electron-electron scattering that results in a drifted Fermi-Dirac distribution [59].
To coarse grain, we partition the time axis into intervals of length τ , t n = nτ , so the environment interacts with the system in approximately the same way during each interval [t n , t n+1 ] [21], where is the averaged value of the map's generator over any interval [t n , t n+1 ] (K is reset at each t n ). If the coarse-graining time τ is short enough, then the short-time expansion of K can be used to perform the coarse-graining [55].
Each term in the short-time expansion of K turns out to be a sum of independent contributions over single-particle states, so in reality we have a multitude of two-level problems, one for each |k , where the two levels are a particle being in |k ("+") and a particle being absent from |k ("-"). Each such 2-level problem is cast on its own 4-dimensional Liouville space, with ρ k = ρ ++ k , ρ +− k , ρ −+ k , ρ −− k T being the reduced statistical operator that describes the occupation of |k and evolves according to a master equation The equations for The above equations describe non-Markovian evolution, because drifted Fermi-Dirac distribution functions in the contacts depend on time through the drift wave vector (related to current). As the transient progresses, the current and the charge density in the structure change, which in turn changes the potential profile, the scattering states available to electrons, the transmission coefficients, and, to a small degree, the interaction matrix elements ∆ ±k , as well as the aforementioned contact distribution functions. Moreover, there may be well-like confined states that cannot be populated by tunneling but only by scattering in the active region. These considerations have been addressed in detail in [56,60]. Figure 5 depicts the potential, charge density, and current density for a single ellipsoidal valley in an nin silicon diode at room temperature. The left and right contacts are doped to 10 17 cm −3 , whereas the middle region is intrinsic (undoped). In the three main panels, the momentum relaxation time in the contacts is taken to be τ =120 fs, based on the textbook mobility values for the above doping density. The characteristic response time of the current is of order hundreds of picoseconds, three orders of magnitude greater than τ . The transient duration can be thought of as the inverse of a typical ∆ 2 k τ among the k's participating in the current flow; shorter τ means weaker coupling and a slower transient (inset).

VII. CONCLUSION
Electrons in the active region of nanostructures constitute an open many-body quantum system, coupled with reservoirs of charge, as well as interacting with phonons and photons. We overviewed the basics of the open system theory, with special focus on the approximations that lead to Markovian and non-Markovian master equations for the reduced statistical operator, and highlighted some recent applications of both types of master equations in quantum transport theory and simulation.
It should be noted that this review did not discuss other widely applied techniques for time-dependent quantum transport, such as the Wigner function simulation [61,62], nonequilibrium Green's functions [28], time-dependent density functional theory [37], Bohmian trajectories [63], or full quantum statistics [64], which will receive due attention in other reviews in this special issue. Also, we did not discuss semiconductor Bloch equations [46,65], which are often employed to address ultrafast optics in semiconductors, and which deserve much more space than available here.
We conclude with some thoughts on the limitations of the master equation framework, as well as potential avenues for further developments.
Active region/contact partitioning. An obvious question is where the active region ends and the contacts begin; there is no a good answer to this question, especially for structures that have no resonances. In large and complex physical systems it is impossible to treat all degrees of freedom quantummechanically, so a boundary between the quantum and the classical (rapidly dephasing) parts has to be adopted, but a boundary should be moved until convergence is reached and the physics no longer varies with its position [66,67]. Rossi [52] has argued that, in the Wigner function simulations, this seemingly arbitrary introduction of the contact/active region boundary results in artifacts that have conceptual, rather than computational origin.
A related issue is that the reduced statistical operator formalism requires that we be able to write the total many-body Fock space as a tensor product of the Fock spaces of the system and environment, and that we write down an interaction Hamiltonian between the two. With S and E containing electrons, we can try to split the total S + E single-particle Hilbert space into S and E subspaces spanned by specific eigenvectors of the position operator, then construct Fock spaces based on these spatially separated single-particle spaces, and finally form a tensor product of said Fock spaces. Unfortunately, this framework artificially makes the interaction local and is not a good choice for capturing current flow that the full S + E Fock space can describe. Rossi [52] shows that consistency requires that the effective interaction depend on the overlap between contact states and active region states, where both contact and active region states in principle extend throughout the whole coordinate space.
Validity of the RWA approximation. The usual secular or RWA approximation -assuming that the system energy levels are so large that the spacing between them is much greater than the system relaxation rate -works well for optical systems and is amply applied in the derivations of master equations for electronic transport, but may not necessarily hold. In fact, in nanostructures with a continuum of states, the spacing between relevantly coupled levels is small and easily smaller than the expected system relaxation rate, especially in the case of strong coupling with the contacts. Therefore, the opposite limit, that of quantum Brownian motion [2,68] may be more applicable in electronic systems with densely spaced system states strongly coupled to the environment. This is a direction in which quantum master equations may have a lot to offer to quantum transport studies [69].
Uncorrelated initial state. Considering that, in reality, the contacts and active region share a Fock space, once we partition it into spatially-determined subspaces and if there are no tunnel barriers, it is not easy to justify the approximation of an uncorrelated initial state. Taking a close look into correlated initial states [70] can be a very fruitful direction of research, one where a tight coupling between approaches that do not adopt contact/active region partitioning, such as TDDFT, with master equations would likely be necessary.
High-frequency transport. Another direction in which the master equation approaches can grow is to look into systems with continua of states and realistic fermionic reservoirs, with a more complete account of intra-reservoir dephasing. This work has opportunities to interface with modern experimental work on GHz-frequency response of nanostructures [71,72].
Deriving single-particle techniques from statistical operator nonunitary dynamics. Capturing the entire statistical operator is feasible only in very small systems. With the statistical operator being the "parent" concept from which single particle quantities such as the density matrix and Green's functions can be derived, it is reasonable to expect that a good non-Markovian approximation for the many-body statistical operator of the electronic system would come first, and from its non-unitary evolution one can further derive single-particle techniques [73]. An open direction of research is to look at single-particle kinetic approaches that originate from non-Markovian approximations for the evolution of the reduced statistical operator. Time-convolutionless non-Markovian equations, thus far underutilized in quantum transport theory, could enable systematic development of singleparticle non-Markovian formalisms that are of a fixed order in the interaction.

VIII. ACKNOWLEDGEMENT
This work has been supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, under award DE-SC0008712.