{'fraction_non_alphanumeric': 0.0380354307530373, 'fraction_numerical': 0.008456282446969832, 'mean_word_length': 4.231846182564868, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Monte Carlo method for solving the NEGF equations for electron transport

We derive, by introducing restrictions to the lesser self energy, a Monte Carlo scheme that solves the NEGF equations for electron transport. In doing so we formally prove that the Monte Carlo estimator has an expectation value equal to the lead current of the NEGF solution, and we provide a simple test of the Monte Carlo scheme by calculating conductivity in nanowires within Buttiker's approximation of scattering. Good agreement between the Monte Carlo simulations and the alternative approaches are obtained, and we also demonstrate the existence of a regime where the Monte Carlo method is the fastest method. In our tests this regime is to extreme to be of practical use. We discuss various ways in which to speed up our prototype and how it can be extended to include more physics. Although the full applicability range of our assumptions about the lesser self energy remains to be better understood, we argue that they should apply as long as the transport process can be considered stationary.


I. INTRODUCTION
In the study of electron transport, the most commonly applied formalism is the Boltzmann transport equation, which works well in most applications. However, effects that originate in the wave like nature of electrons are difficult to address with the semiclassical Boltzmann equation. In order to include such effects, it is instead preferable to use a transport formalism derived from the Schrödinger equation. An example of such a formalism, is the Non-equilibrium Green's functions (NEGF) 1 , a generalization of many particle perturbation theory to non-equilibrium. Due to reduction in size and increased complexity of combinations of different materials, recent years has brought significant interest in developing NEGF solvers 2,3 .
However, a remaining challenge with NEGF is the computational burden. NEGF is formulated in terms of matrices of size N × N, where N is proportional to the spatial system size, and to the number of included states. In addition, these matrices are functions of energy E, and in many cases also of a Bloch vector k. To obtain acceptable accuracy, the continuous parameters E and k must be discretized in the form of a dense integration grid.
Thus, without the introduction of approximations, NEGF calculations involve the solution of large equation sets.
In fact, similar challenges also apply when solving the Boltzmann equation. The Boltzmann distribution function, without approximations depends both on position and on Bloch vector. Thus, the number of parameters is still proportional to the system size and large integration grid. By introducing the relaxation time approximation, the procedure can be simplified, particularly in the linear response regime. However, in cases where this approximation is not warranted one must instead solve the general Boltzmann equation. A common approach is to utilize Monte Carlo techniques, which involves the explicit simulation of electron motion, with random scattering between k-points 1,4 .  [5][6][7] and is also similar to a Neumann expansion. The DMC thus shares some features with the NEGF approach to be introduced in this work. However, an important difference is that DMC performs a Monte Carlo evaluation also on the self energy. In the case of electrical transport, typical NEGF solvers 2,3 assume the self energies to be given explicitly.
The evaluation of Monte Carlo estimators in quantum transport will generally result in opposite signs or phases. In a large sample one can thus expect partial cancellation 1 , leading to a variance which is greater than the expectation value. The Monte Carlo integral thus converges slowly. In the DMC method this is partially alleviated by performing some of the diagrammatic summations analytically 7 .
In this work, we suggest a Monte Carlo scheme based on a Neumann expansion of the lesser Green's function G < in the NEGF formalism. In the method suggested here, we will assume that the system is stationary. As will be shown, the Neumann expansion can be expressed in terms of positive quantities, and the "sign problem" is avoided. However, this is at the cost of increased computational complexity per estimate.
Briefly, the suggested Monte Carlo method can be described as a process where electrons are randomly scattered from state to state. This is similar to the DMC method and to BMC simulations. However, while states in the BMC simulations are described by a combined position in physical space and k-space, we instead describe states using wave functions |ψ . These wave functions are found by solving a modified Schrödinger equation with a source term |i . The source term is picked randomly among the eigenvectors of the operator Λ(|ψ ′ ψ ′ |), where |ψ ′ is the previous state, and Λ is a function with the property that Λ(G < ) = Σ < . G < being the lesser Green's function and Σ < the lesser self energy. This process is described in more detail in section II D.
As long as sparse matrices can be used to model the Hamiltonian and the scattering operators, the solution of the modified Schröedinger equation will be a linearly scaling oper-ation. Thus, since the calculation of N eigenvectors will always scale at least quadratically with the system size N, the bottleneck of the method will generally be the diagonalization of Λ(|ψ ψ|). In order for the method to be practical, it is likely that approximation must be applied in order to make the diagonalization tractable. Among such approximations, is Büttiker's approximation, where Λ(|ψ ψ|) is diagonal in some basis {|i }, independent of |ψ . In this work, all calculations have been performed within this approximation in order to test the validity of the method.
This work is organized as follows: Section II presents the method. Here we prove that the described Monte Carlo scheme solves the NEGF equations. A discussion of the linear limit of the procedure and a description of Büttiker's approximation is also given. In section III we describe the prototype implementation in more detail.
The results of the initial tests are presented in section IV, where we calculate the conductance of nanowires. In particular, we study how the Monte Carlo computation time required for a fixed accuracy of one percent scales with the length of the wires, and how this compares to iterative and direct schemes of solution.
Our prototype Monte Carlo solver is implemented in Python. Still, we demonstrate the existence of a regime where the Monte Carlo calculations are faster than the alternative approaches. This regime is however not to be expected to be relevant in practical applications.
Thus, in order for the presented Monte Carlo scheme to be useful, steps must be taken to further increase performance. We do expect a significant improvement upon transition to a compiled language. Further suggestions of how it is possible to reduce the computational burden are discussed in section V, where we also discuss how more general physics can be included. Finally, section VI summarizes and presents a few conclusions.

II. THEORY
This section is divided into the the following parts: Section II A gives a quick summary of the NEGF formalism, while in section II B we derive a Neumann expansion of the NEGF equation, which is needed to define the general Monte Carlo integration scheme. This scheme is derived in section II C, where the most general version of the method is described. This requires that we assume some probability distribution g, from which sequences of electron states |ψ can be drawn. In section II D we introduce a concrete procedure which specifies the probability distribution g. Finally, in section II E we discuss the linear limit for elastic scattering, and in section II F Büttiker's approximation.

A. Summary of the NEGF formalism
The Non-equilibrium Green's function formalism for electrons takes its simplest form when expressed in terms of two correlation functions G r and G < . These are respectively referred to as the retarded and lesser Green's functions, and can be expressed as matrices Here a i and a † i are fermionic annihilation and creation operators, operating on electrons in the single particle state i. θ(t) is the unit step function, curly brackets represents the anti-commutator, and the brackets represent the average with respect to the quantum state.
If the Green's functions describe a stationary process, then they will only depend on the difference between the time arguments, so that we may define the Fourier transforms 8 In addition, the formalism makes use of the advanced Green's function G a , which can be found through the relation G a (E) = G r (E) † , as well as the spectral density The literature operates with multiple conventions regarding the definition of A(E), but in this work we follow the line of Datta 8 , and define this without additional prefactors.
In equilibrium there is a simple relationship between the lesser Green's function and the spectral density, namely 1,8 where f (E) is the Fermi distribution. Away from equilibrium however, there is no such simple relation, and accordingly G < and G r must be calculated independently. Within diagrammatic perturbation theory, it can be shown 9,10 that for stationary processes this can be done by solving the equations where H is the single particle Hamiltonian of the electrons, while Σ r and Σ < are respectively the retarded and the lesser self energies. These self energies will in general be functions of G r and G < , and their form will depend on the particular perturbation used. Examples of such perturbative expressions can be found in the literature 1,8,11 , while a more heuristic model is discussed in section II F.
In a typical electron transport problem, we are considering a device or material which is connected to a set of leads {p}. The leads are assumed to be internally close to equilibrium, at least compared to the device, so that their electronic occupations are described by Fermi distributions f p (E). It can then be shown 8 that the effect these leads have on the device can be accounted for in terms of the self energies Here Σ r s and Σ < s are self energies due to scattering, while Σ r p and iΓ p (E)f p (E) are respectively the retarded and lesser self energies due to lead p. Expressions for Σ r p can be found in the literature 8 , while Γ p is defined as where Σ a p (E) = Σ r p (E) † . In addition to Eq. (11), we also make the definitions where Σ a (E) = Σ r (E) † and Σ a s (E) = Σ r s (E) † . It can be shown 8 that Finally, it can also be shown 8 that the electric current I q and heat current Q q at a lead q will be given respectively as I q = −e i q (E)dE and Q q = (E − µ q )i q (E)dE. Here e is the elementary charge, µ q the electrochemical potential of lead q, while Boltzmann distribution function, L is a linear integral operator, and f 0 is an independent source term which arises from the initial conditions 1 . We can put Eq. (8) into a similar form by introducing the following assumption about Σ < s : Assumption (i) The lesser scattering self energy can be written as where the function Λ(E, E ′ , G) is linear in G.
This assumption holds at least in the self consistent Born approximation, except for electron-electron scattering. Further justification for the assumption is discussed in appendix A 1.
Taking into account Eq. (10) and Eq. (16), Eq. (8) becomes where we have defined The integral operator Ξ is now defined as and then we write Eq. (17) as The expression for G < 0 can be simplified by making the definitions so that Eq. (21) can be written Eq. (20) has a form similar to the integral Boltzmann equation. A minor difference is that the source term G < 0 arises from coupling to the leads rather than from initial conditions. A Neumann expansion of G < , similar to what is done with the Boltzmann equation 1 can now be performed.
However, before we proceed we note that Eq. (20) depends on the retarded Green's function G r solving Eq. (7) through Eq. (18). In general, both of the scattering self energies Σ r s and Σ < s will depend on both Greens's functions G r and G < , while H may also depend on G < through a Poisson potential. Thus, Eq. (7) and Eq. (8) are in fact coupled equations. Since a Neumann expansion based on Eq. (20) must involve G r , it will generally not be applicable unless both equations have already been solved. This problem will be discussed further in section II E, and to some extent in section V B 3. Until then we will simply regard G r as known.
To find the Neumann expansion of Eq. (20), we begin by iteratively defining the functions where the last step makes use of assumption (i). Assuming that the series converges, the iteration of Eq. (25) will converge to This function solves Eq. (20), as can be seen by the evaluation Thus, as indicated by the notation, G < (E) is in fact the lesser Green's function, and Eq. (26) its Neumann expansion.
By now applying Eq. (26), we can then find the energy resolved current of Eq. (15) as Since we assume G r and thus the operator Ξ to be known, Eq. (28) could in principle be evaluated directly. However, this is computationally demanding for the general case.

C. General Monte Carlo method
The expansion Eq. (28) is the starting point of the Monte Carlo method we here propose.
There are several different ways in which this expression could be evaluated through a Monte Carlo process. Probably, the simplest such method would be to express the sum in terms of matrix elements as where Γ ij q (E) = i|Γ q (E)|j , G ij p (E) = i|G p (E)|j , and Ξ ij kl (E, E ′ ) = k|Ξ(E, E ′ , |i j|)|l for some basis {|i }. If a set of energies E 0 · · · E n−1 and indices i 0 j 0 · · · i n j n are picked at random from some probability distribution f , then the stochastic variable will have Eq. (29) as its expectation value. Thus, the average of a large number of random values of x could be used as an estimate of the sum in Eq. (28). An advantage of this particular approach is that the calculation of each individual Monte Carlo estimate x is numerically fast. It is simply the product of a series of scalar quantities. However, a major challenge is that x is in general a complex quantity. Consequently one can expect instances of x to cancel each other.
Such "sign problems" as they are often referred to 6 , cause major issues in most quantum Monte Carlo methods. In our case, it makes the use of the estimator of Eq. (30) completely unfeasible, as all cancellations would make the estimate converge very slowly. If possible we should rather try to find an estimator which is always positive in order to remedy the sign problem. From an assumption which will be introduced in the next section, it will follow that the operators in Eq. (28) are positively definite. Thus, one possibility would be instead of x, to use the estimator where the number n is drawn randomly from the probability distribution p n . The estimator y also has an expectation value equal to the sum in Eq. (28), but is always positive. The estimator y however, has another problem. The calculation of each estimate is numerically expensive. This is partially because it requires integration over energy, but more importantly because each step of the calculation handles huge and potentially dense matrices Ξ k G < 0 , for k = 0, . . . , n.
Thus, since Monte Carlo estimators expressed in terms of scalars and operators both have issues, we shall instead try to follow a middle ground. We seek an estimator expressed in terms of vectors. The easiest way to express Eq. (28) in terms of vectors, is to diagonalize the operators. For this purpose we introduce an additional assumption about the scattering operators: Assumption (ii) Λ preserves Hermiticity in G. That is, for all E and E ′ , and any Hermitian operator G, Λ(E, E ′ , G) is also Hermitian.
The justification for this assumption is discussed in appendix A 2.
Since it can be shown that the operators G p (E) = G r (E)Γ p (E)G a (E)f p (E) are Hermitian, it follows from assumption (ii) that all of the operators Ξ(E, E n−1 , Ξ(· · · Ξ(E 1 , E 0 , G p (E 0 )) · · · )) are also Hermitian and diagonalizable. The sum in Eq. (28) can thus be expressed as where λ (n) i and |i (n) are respectively the eigenvalues and eigenvectors of the operator Ξ(E, E n−1 , Ξ(· · · Ξ(E 1 , E 0 , G p (E 0 )) · · · )). It follows from the assumed positive definiteness of these operators that an estimator λ While this estimator is expressed in terms of vectors, handling the operators ) · · · )) for k = 0, . . . , n is still required during intermediate steps of the calculation. However, this issue is easy to eliminate, since we can start the diagonalization at the innermost operator. Using Eq. (18) and the linearity of Λ, we get where λ i k and |i k are now respectively the eigenvalues and eigenvectors of the operators Λ (k) , defined recursively as for k ≥ 1, and as Λ (0) = Γ p (E 0 ).
Inserting Eq. (33) into Eq. (28), we get Now consider a probability distribution g(p, n, φ, E, r) from which we draw an initial lead p, a final lead r, a number n, a vector of indices φ = [i 0 , i 1 , . . . , i n ], and a vector of energies We define a stochastic variable the expectation value of which will be Eq. (37) is seen to be identical to the sum of Eq. (35), except that the integral is also over the final energy E n . Thus, we can estimate the current at lead q as whereX q denotes the average of X q over the Monte Carlo ensemble. The heat current Q q can be evaluated by an analogous expression.
The assumption of positive definiteness leads to a positive estimator X q . We thus avoid the sign problem. In addition, the calculation of X q has the desired property of mainly involving vector operations. The only exception to this is that one must diagonalize the operators Λ (k) . However, given reasonable approximations one can assume these operators to be sparse. This approach is then still significantly more efficient compared to dealing with the operators Ξ k G < 0 , which are likely to be dense for sufficiently high k. In fact, by utilizing an additional approximation in section II F, we can avoid diagonalization all together. However, let us first define an algorithm which determines the probability distribution g.

D. Scattering Monte Carlo method
In the literature 1,4 one will find described a particularly intuitive BMC method, where electrons are scattered from state to state with probabilities proportional to the physical scattering rates. While improvements in terms of convergence can be made using so called weighted Monte Carlo methods 1 , the benefit of the unweighted intuitive approach is that the estimators have particularly simple expressions. Generally, one should expect complicated estimators to have smaller variance than simpler ones only in carefully constructed cases.
So as a starting point, we in this work seek to mimic the intuitive scattering Monte Carlo method, in the hope that this will result in a moderate variance.
We start by interpreting the function Λ as somewhat analogous to the scattering rate Γ in the Boltzmann equation, and seek to express the probability of scattering to new states in terms of Λ. The involved quantities need to be positive and we need to introduce an assumption about positive definiteness. While we have briefly mentioned such an assumption earlier, in relation to discussions of the sign problem, we now introduce it formally.
Assumption (iii) Λ preserves positive definiteness in G. That is, for all E and E ′ , and any Hermitian positively definite operator G, Λ(E, E ′ , G) is also positively definite.
The justification for this assumption is discussed in appendix A 3. We make use of assumption (iii), and define the following Monte Carlo procedure: Step 1. Draw a starting lead p, all leads having equal probability 1/N.
Step 2. Draw a starting energy E 0 from some probability distribution h(E).
Step 3. Diagonalize Γ p (E 0 ) to find its eigenvectors |i and eigenvalues λ i .
Step 4. We construct a probability distribution From the expressions for Γ p found in the literature 8 , we see that it is positive definite.
Furthermore, the spectral density A must also be positively definite due to its relation to the density of states D i (E) = i|A(E)|i . Thus, p i is positive and sums to one. Draw an index i 0 from this probability distribution, and set n = 0.
Heuristic justification: The probability that an electron coming from lead p will enter the system in state |i , is proportional to the projected density of states in |i , and to the strength of interaction between the state |i and lead p. The density of states in |i is Step 5. Find the state vector |ψ n = G r (E n )|i n . It can be found either by explicitly calculating the retarded Green's function G r , or alternatively by solving the modified Schrodinger equation Heuristic justification: This step represents the free propagation of electrons between scattering events.
Step 6. Calculate the quantities Tr Γ r |ψ n ψ n | = ψ n |Γ r |ψ n for all leads r, and also Step 7. From Eq. (34), we recognize From assumption (iii), the eigenvalues λ i are positive, so that Accordingly, S n is also positive, so we can construct a probability distribution Attempt to draw a lead r from this distribution, which is positive, but does not sum to one. Since the distribution does not sum to one, a successful draw is not guaranteed. If the draw is successful and a lead r is drawn, stop the process and choose lead r as the end lead.
Otherwise continue the process.
Heuristic justification: After free propagation, the electron will either scatter to a new state, or exit the system through one of the leads. Following the reasoning behind step 4, the probability of exiting through a lead r should be proportional to Tr Γ r |ψ n ψ n | = ψ n |Γ r |ψ n .
By interpreting A(E) as the density of states, and the function Λ as describing the scattering interaction strength, one could argue that the scattering rate to energy E should be given by . Thus, S n given by Eq. (41) represents the total scattering rate.
Since the probability of scattering is proportional to S n , while the probabilities of exiting through a lead r is proportional to ψ n |Γ r |ψ n , the requirement of total summation to one gives Eq. (44).
Step 8. Draw an energy E n+1 from the probability distribution which integrates to 1, and is positive from Eq. (43).
Heuristic justification: Tr A(E)Λ (E, E n , |ψ n ψ n |) represents the scattering rate to energy E, while S n represents the total scattering rate.
Step 10. Draw an index i n+1 from the distribution .
Heuristic justification: Following again the reasoning of step 4, the probability for an electron to scatter from state |ψ n to |i , is proportional to the projected density of states in |i , and to the strength of interaction between |i and |ψ n . The density of states in |i is given by D i (E) = i|A(E)|i , while the interaction strength is represented by the eigenvalues λ i from step 9. Requiring summation to one we get Eq. (46).
Step 11. Increment n by 1, and return to step 5. Continue the procedure from there.
The steps 1 to 11 define the Monte Carlo procedure. The procedure is truncated only when an end lead r is drawn in step 7. At that point we will have selected start and end leads p and r, and a set of n + 1 energies and n + 1 indices. Thus, the procedure defines a probability distribution g(p, n, φ, E, r).
We can evaluate this probability by taking the product of the probabilities associated with each step of the procedure. Clearly, steps 1 to 4 always give rise to a combined factor The contributions from the remaining steps will depend on whether the electron scatters or exits at a lead in step 7. The probability of scattering is the same as the probability of not exiting at a lead, so that every scattering event will contribute a factor to the total probability, p r being given by Eq. (44). In addition, each scattering event will also contribute a factor p(E) from Eq. (45), and a factor p i from Eq. (46). Thus, assuming that the electron is scattered n times before exiting, the total contribution to g from scattering becomes n−1 The probability of finally exiting at lead r after the n scattering events, is given by Eq. (44). Thus, taking the product of Eq. (47), Eq. (49) and Eq. (44), we obtain Having found an explicit expression for the probability distribution g, we can find an expression for the estimator X q using Eq. (36). Inserting Eq. (50) in Eq. (36), we find that the eigenvalues and the quantity what remains is While these cancellations simplifies the expression, what remains is still fairly involved.
However, as we will see, Eq. (51) can be simplified even further by making an additional assumption about Λ(E ′ , E, G): Assumption (iv) For any Hermitian and positively definite operator G, Assumption (iv) can be interpreted as an assumption of the scattering mechanism being charge conserving. This is more thoroughly discussed in appendix A 4.
Employing assumption (iv), we first note that in step 6 above, we can calculate the quantity S n of Eq. (41) simply as so that we no longer need to perform the integral explicitly. Next, making use of Eq. (14) and Eq. (9), we see that Inserting this into Eq. (51), we see that the remaining terms in the product also cancel, so that X q becomes simply which has the simple form we were looking for.

E. Linear elastic limit
In the preceding discussions we assumed that either the retarded Green's function G r , or at least its inverse is known. In the previous section for instance, an expression for G r or its inverse was required in step 5. But in reality G r is rarely known beforehand, since the scattering self energy Σ r s usually has a functional dependence on G < . In fact, even the Hamiltonian H may depend on G < , for instance through a Poisson potential. Thus, Eq. (7) and Eq. (8) must generally be solved self consistently with each other. The method presented above will thus not be applicable. In section V B 3 we briefly outline a possible approach for how to deal with this in the general case. In this section we restrict the discussion to linear transport.
In the linear limit, the existence of Kubo expressions tells us that the currents should be expressible in terms of equilibrium quantities. If we could replace G r above with its equilibrium value, then Eq. (7) and Eq. (8) would no longer be coupled, and the described method would thus be applicable. While it is difficult to show that this substitution is possible in the general case, it is rather straightforward if we also assume elastic scattering.
In that case, we have so that the operator Ξ is given simply by In addition, we need one final assumption about the scattering mechanism, namely Assumption (v) When the scattering mechanism is elastic, the function Λ(E, G) satis- where F and G are any two Hermitian and positively definite operators.
Intuitively, this assumption can be interpreted as a statement of detailed balance. Its justification is discussed in appendix A 5.
It follows from assumptions (v) and (iv) that the quantity Tr A(E)Γ q (E) can be expanded as The proof of this is given in appendix B. Inserting Eq. (59) into Eq. (28) yields . Assume now that the the system's offset from equilibrium is described by some parameter x, which could for instance be a voltage or a temperature difference between two of the leads. Then all functions in Eq. (60) are functions also of x. Expanding to first order in x yields where all functions on the right are evaluated at equilibrium. But in equilibrium f p (E) = f q (E), so that in fact Thus, all transport coefficients can be evaluated from the transmission functions T Eq qp (E) at equilibrium. The conductance G qp for instance, are given by the expression where Th(E) = −∂f /∂E. This means that G qp can be estimated as where the stochastic variables X pq are found using the Monte Carlo procedure described above with a few modifications; the initial lead p is kept fixed, f (E) is replaced by Th(E), and G r and the function Λ are replaced with their respective equilibrium values. If we also choose the probability distribution as h(E) = Th(E), then Eq. (55) becomes simply

F. Büttiker's Approximation
In the general case, the bottleneck of the procedure described in section II D will be step 9, which requires diagonalization of the operators Λ (n+1) . In order for the procedure to be applicable in practice, one must probably always find some approximation which improves the calculation efficiency of this step. The simplest of such is to assume that Λ (n+1) is already diagonal in the adopted basis. In fact, the procedure can be simplified even further by making some additional assumptions. First, we assume the existence of a single basis set If we in addition to Büttiker's approximation also assume elastic scattering, theñ In that case, we also see that so that assumption (v) is satisfied automatically.
We are now at a point where the procedure of section II D can be greatly simplified. First of all, since Λ (E, E n , |ψ n ψ n |) is already diagonal in the basis {|i }, step 9 can be omitted.
Secondly, the total probability associated with steps 8 and 10 will be Accordingly, we can reverse the order of these two steps by first picking an index i with and then an energy E with probability The probability of picking a particular index i, given that a lead is not drawn, is then given precisely by Eq. (73).
Thus, steps 6 to 10 of section II D, can be replaced by the steps Step 6' Draw an internal index i n+1 with probabilities given by Eq. (75), or a lead r with probabilities p r = ψ n |Γ r |ψ n r ψ n |Γ r |ψ n + j γ j | ψ n |j | 2 .
If a lead is drawn, stop the process. If an internal index is drawn, continue.
Step 7' Draw an energy E n+1 from the probability distribution of Eq. (74).
After step 7', the process is continued from step 11.

A. Model
Our tests of the described method have been performed within Büttiker's approximation, using the specialized method described in section II F. Application to more general scattering mechanisms, using the method of section II D, is postponed for future studies. In addition, the performed tests assume elastic scattering.
In general the lifetime τ could be a function of both E and i, but in our tests it is set to a constant.
All tests are performed with only two leads p = 1 and p = 2. To simplify the procedure we write the lead self energies as where P p is a projection operator projecting onto orbitals on the boundary with lead p, and Γ l is a constant describing the strength of coupling between the system and leads. This is an unrealistic model of Σ r p (E), but the physics relevant to the tests remains correct.

B. Monte Carlo implementation
Our prototype implementation follows the approach described in section II F. It is written in Python, and makes use of the RGF algorithm 11 to calculate G r by solving Eq. (7). Since the basis set {|i } is the tight binding basis set, we only need a single column of G r for each execution of step 5. Thus, the diagonal of G r is calculated in the beginning of each simulation, while the required columns of G r are calculated as needed, upon execution of step 5.
The instantiation of stochastic variables is done by transforming to a homogeneously distributed variables 16 , which are drawn using the Numpy rand function. In particular, the execution of step 6' is done as follows: An array a j is defined, storing the values γ i | ψ n |i | 2 , and ψ n |Γ r |ψ n . The cumulative sum b j of the array a j is calculated. b j is then divided by its last element, which normalizes it to a cumulative probability distribution. A random number x is picked homogeneously in the range [0, 1], and the Numpy argmax function is used to find the first value in b j which is larger than x. If the index j corresponds to an element of a j given by γ i | ψ n |i | 2 , then the internal index i n+1 = i is picked. If the index j corresponds to an element of a j given by ψ n |Γ r |ψ n , then the lead r is picked as an exit point.

C. Direct and iterative methods
For comparison reasons we also implemented a scheme where Eq. (8) is solved directly.
Within the elastic Büttiker approximation, Eq. (8) can be written Eq. (79) can be solved directly or iteratively, and we test both of these against the Monte Carlo method. The direct method employs the Numpy solve function, while the iterative method employs the Scipy implementation of the gmres method.
To evaluate the transport coefficients G pq , we must integrate the quantity i p (E) from Eq. (15). This is done by a midpoint integration scheme, using a regular integration grid of N E points in the range µ ± 10k B T .

IV. RESULTS
As an initial test of the Monte Carlo method we estimate the conductivity of nanowires.
In this case of CdTe, with a square cross section of two by two unit cells. The conductivity is found by fitting against the conductances of wires, with lengths 10 × n unit cells, n being in the range 1 to 10. In all calculations the temperature is set to 300 K, and the chemical potential is set to µ = 4.5 eV, which is about 1 eV above the conduction band edge.
The scattering models use Büttiker's approximation. Different values of γ =h/τ , namely 0.02, 0.04 and 0.06 eV were sampled. This corresponds to life times of τ = 33, 16 and 11 fs, respectively. The coupling to the leads is set to Γ l = 0.2 eV.
The scaling of the Monte Carlo method is compared to that of the direct and iterative schemes described in section III C. We required an accuracy of at least one percent for each method. In the Monte Carlo method this is assured by ending the sampling when a relative deviation of 10 −2 is achieved. In the other methods, the accuracy is controlled by the resolution N E of the energy integration grid. We use a resolution of N E = 13 in all cases. This grid was chosen since it results in an error of one percent when compared to the converged N E → ∞ result in the γ = 0.06 eV case.

A. Transport results
The wire resistance as a function of wire length is shown in figure 1. The Monte Carlo calculations are consistent with the other results, except for γ = 0.02 eV. We located this to be the result of N E = 13 being to coarse for this particular case. To show that integration is in fact the source of this discrepancy, we also include a calculation with N E = 26, where the results at γ = 0.02 eV is indeed in better agreement.
The discrepancy between the γ = 0.02 eV Monte Carlo and N E = 13 result is seen to be largest for short wire lengths. This indicates an increased sensitivity to the energy resolution due to edge effects penetrating further into the wires, which would be caused by the weaker scattering of γ = 0.02 eV. For the longer wires we thus expect N E = 13 to be sufficient, and This means that the Monte Carlo results stay within three standard deviations from the direct/iterative solutions, as expected.
The coefficients of Eq. (80) and Eq. (81), which differ slightly from those of the fits in figure   3, represent the fits obtained with fixed exponents. The number of samples required to obtain an accuracy of 0.01 will be n = ∆X 2 / X 2 /0.01 2 , where X is the estimator. X is given by Eq. (66), which with some approximation can be

Monte Carlo method
where p(x) is the probability that the electron reaches a distance x into the wire, and s(x) is the average number of steps taken by an electron reaching this distance. The probability of reaching at least a distance x, is given by the transmission function T (x), which scales as λ/x. The probability p(x) of an electron reaching a distance of exactly x is then given by the derivative of T , and is thus proportional to λ/x 2 . Finally, the electron is undergoing a random diffusion type motion through the nanowire, which means the number of steps required to move a distance x scales as (x/λ) 2 . By inserting these estimates into Eq. (82) we get s ∼ L/λ.
In fact, not only does the number of samples and the number of steps per sample show linear behavior, both are also proportional to 1/λ. Since the mean free path λ should be inversely proportional to the scattering rate γ, the slopes of the fitting functions in figure 4 should be proportional to γ. Table I also shows the ratios a/γ, which is indeed seen to be roughly independent of γ. The average of all three ratios is respectively 42389 and 3.1 eV −1 for the fits of figures 4a and 4b. The average constant term b n is found to be 70803. From this, we obtain which we assume is valid for all values of γ. We have here ignored the constant term b s , since this is too small to be of significance at the length scales of these calculations. We adopt the following models for the computation times t d and t s : where the constant terms have been ignored due to their small size. The computation time for the total computation time.

Comparison between the methods
The total computation times of all methods are shown in figure 6, where they are compared to their respective fitted models. These models are given by Eq.    Figure 7 also makes reference to a third method, which will be discussed in section V A 5.
By the arguments above, the relevant parameter determining the required computation time of the Monte Carlo method is actually the mean free path λ, rather than the scattering strength γ. We can obtain some very rough estimates of the mean free paths in our calculations by utilizing the fact that edge effects will penetrate roughly a few mean free paths into the wires. Observing figure 1, we see that in the γ = 0.04 eV case, edge effects are causing a deviance from the linear behavior for wire lengths up to 20 nm, while the same edge effects Assuming that edge effects are visible for wire lengths shorter than ∼ 2λ, we estimate the mean free paths in those two cases to be 10 and 20 nm. Since λ should be inversely proportional to the scattering rates, we obtain the model λ = 0.4 eVnm/γ. This would mean the charge carriers in our example nanowire have an approximate average group velocity of 4 eVÅ, which is a reasonable number. Adopting this model, the condition γ < 70 µeV can be translated to the requirement that λ must be larger than approximately 6 µm. This is unlikely to be the case in most practical calculations.

A. Performance
Since the Monte Carlo method requires large mean free paths > 6 µm in order to be faster than the iterative method, our current implementation is not useful for the chosen example. This situation is unlikely to change for other finite size geometries. Thus, in order for the method to be useful, we are in need of an improved implementation. There are multiple ways in which both the method itself, and its implementation can be improved, performance wise. Some of these will be discussed in the following. Their implementation and verification will be postponed to future studies.

Compiled solvers
Our prototype tests make use of a straightforward Python implementation of the Monte Carlo solver. It is well known that there is some overhead, possibly orders of magnitude due to the non-compiled nature of Python. It is important to emphasize that also the other methods in our tests were implemented in Python. However, most of these made use of precompiled solvers, such as gmres in the case of the iterative method. Thus, a direct comparison between the methods is not entirely fair, and we can expect improved performance by writing the Monte Carlo solver in a compiled language.
However, a significant improvement is not guarantied to follow from such conversion alone.
While the RGF algorithm was also implemented in Python, it utilizes a large number of matrix inversions and products, functionality which is already precompiled through NumPy.
Any expected improvements by switching the entire method to a compiled implementation will be limited by the ratio in cost between these matrix operations and the Python overhead.
This will largely depend on the size of the matrices, which in the case tested here was only 32 by 32.

Systems with k-dependence
The relative difference in performance between the Monte Carlo method and the other approaches can be expected to improve in systems with a k-dependence. If the nanowires of section IV were replaced with thin-films of varying thickness, the Hamiltonian and other operators would then have a two dimensional k-dependence. If a 10 × 10 k-point grid was used for the Brillouin zone, then the number of unknowns in Eq. (79) would increase by a factor of 100. Accordingly, the direct and iterative solution schemes would respectively be slower by six and four orders of magnitude. On the other hand, the change in computation time of the Monte Carlo method would be determined by the change in the variance of the estimator. Although the variance can be expected to increase, a change by four orders of magnitude is unlikely.

Parallelization
One large benefit to Monte Carlo calculations, is that they are very easy to parallelize, and that they scale near perfectly with processor distribution. Thus, although what we refer to as the total computation time would not be reduced, the real time requirement of the calculations could be significantly improved by utilizing high throughput parallelization.
In fact, we have already made use of parallelization in our test calculations, which were

Wave function based solver
Switching to a wave function based solver is another possibility. That is, instead of calculating G r explicitly, the wave functions |ψ n = G r |i n is calculated by solving Eq. (40). In the calculations of section IV a wave function based solver would not cause a very significant improvement. This is due to the utilization of the RGF algorithm, which scales linearly with system size along one spatial dimension. However, for the other spatial dimensions the RGF algorithm has cubic scaling. It is thus best suited for long structures, such as thin nanowires.
In more bulk like systems, wave function based methods utilizing efficient sparse linear solvers, such as those employed by the ballistic transport software Kwant 17 are most likely better suited. Alternatively, iterative sparse solvers like the gmres method would be less memory demanding, and thus better suited for parallelization.

Interpolation
The In figure 7 we also show the parameter range in which this interpolating Monte Carlo method is faster than the iterative method. The scaling estimate of the interpolating method was obtained using an analysis similar to that of section IV B 2. Figure 7 shows that the minimal mean free path requirement is reduced to about 200 nm when the interpolating method is used. This is a significant improvement, but still not quite in the range of practical applicability. However, because the Monte Carlo method consists only of very simple operations when the wave function calculations are eliminated, we believe that the computational time of the interpolating method is mainly determined by the Python overhead. We expect further reduction in efficiency when traversing this solution to a compiled language.

Improved estimator
Finally, the efficiency of the Monte Carlo method could also be improved by reducing the variance of the estimators. One way of approaching this is the use of weighted Monte Carlo methods, like those applied in Boltzmann calculations 1 . This would involve modifying the probability distributions of section II D to new probabilities p ′ , and then weighting the estimator with the probability ratios p/p ′ . If this is done correctly, an estimator X ′ will have smaller variance. Less samples will thus be required to achieve a particular accuracy.
There is also another way in which the variance can be reduced for bulk quantities like the electrical conductivity. From Eq. (83) we see that the number of samples N mc is increased significantly by the presence of a large constant term. This term stems from the constant term A in 1/T = A+L/λ, which reduces the binomial probability and increases the variance.
The term A is related to the probability of transmission between the system and the leads.
Thus, we could reduce the variance of the estimator by eliminating the leads from the expression. Intuitively, it should be possible to estimate conductivity from the distribution of electron step lengths |x n+1 − x n |, rather than from expressions like X pq , which involve the leads p and q. Here x n = ψ n |x|ψ n .

B. Generalizations
In addition to the suggested strategies to improve the computational efficiency of the calculations, it is also possible to incorporate additional physics in the Monte Carlo model.

Generalizations of the scattering model
The test calculations of this work utilize two approximations to the scattering model: 1) Büttiker's approximation, and 2) elastic scattering. In calculations intended to model material systems for technology applications, these approximations are likely to be too limiting. In particular, self energies derived from various perturbative schemes, e.g. modeling electron-phonon interaction and disorder scattering will rarely conform to Büttiker's approximation.
If such models are employed, it is necessary to use the more general method of section II D, rather than the simplified algorithm of section II F. This requires diagonalization of the operators Λ (n+1) , which increases the memory requirements and the computation time of the method. However, since scattering is usually a local phenomenon, Λ (n+1) can most likely be approximated as a sparse matrix, possibly ordered, so that the increased computational requirements are manageable.
In addition, both electron-electron and electron-phonon interactions give rise to inelastic scattering. While acoustic phonon scattering can often be approximated as elastic, this is not the case with scattering of optical phonons. Also, in applications where energy relaxation, for instance thermalization of states is important, even acoustic scattering may need to be modeled as inelastic. In high field applications this is merely a question of finding a way to draw energies from the probability distributions of Eq. (45) or Eq. (74), followed by application of the algorithms in sections II D or II F.
However, when one is interested in linear transport coefficients such as G pq , this approach is not well suited. At low fields the currents I q will be very small, so that the signal to noise ratio of the estimators X q will also be low. Thus, in the linear limit we believe a better approach is to find an expression which generalizes Eq. (66) to inelastic scattering. This generalized expression is likely to be considerably more involved.

Generalizations of the Hamiltonian
The Hamiltonian employed in our calculations could also be generalized in various ways.
Within a tight binding framework, one could improve the model either by including more orbitals per atom, or by including interactions of longer range than the nearest neighbor. In addition, one could also include effects such as a Poisson potential or strain dependent model parameters, and interface the Hamiltonian directly to some ab-initio software, employing for instance Density Functional Theory (DFT) or post DFT methods.

Nonlinearities
The Monte Carlo methods described in sections II C, II D and II F are all based on the assumption that the function Ξ(E, E ′ , G) defined in Eq. (18) is linear. There are two ways in which this linearity assumption could fail: 1) The scattering function Λ(E, E ′ , G) could itself be nonlinear, or 2) the Green's functions G r and G a could be functions of G < .
As discussed in appendix A 1, Λ will typically be nonlinear if one makes use of perturbation theory to orders higher than one, or if one includes electron-electron scattering to any order. Also, as discussed in section II E, G r and G a can be functions of G < either through the scattering mechanism, or for instance a Poisson potential. Thus, the function Ξ(E, E ′ , G) will in most realistic cases be nonlinear in G.
We will handle the nonlinearity by introducing an object x = F (G < ), and write the function Ξ as Ξ(E, E ′ , G, x), where Ξ is now linear in the third argument. Note that this is always possible in principle, since we could choose x = G < . The reason for using the more general expression is that it will sometimes be more efficient to choose x differently. For instance, if the nonlinearity comes from a Poisson potential, one could choose x = ρ, where ρ is the charge distribution.
Given any fixed value of x, Ξ is linear, and we can solve Eq. (17) using the Monte Carlo method. The solution is a lesser Green's functionG < (x) which depends on the value used for x.G < (x) will only equal the true solution G < if x is set to its true value x ⋆ . That is, . Since x ⋆ = F (G < ), we can find x ⋆ by solving the equation x = F (G < (x)). In particular, if F (G < (x)) can be estimated in the Monte Carlo method using some estimator X, then x ⋆ is the solution of X|x = x ⋆ = x ⋆ . Here the notation X|x = x ′ means that the estimator is calculated in a simulation where x = x ′ .
One way of solving this equation is to construct a sequence of random estimates X n = X|x = x n , where x n is always the average of all previous estimates. That is, for k > 0, while x 0 equals some initial guess. Intuitively, x n should approach the correct value of x as n → ∞. Such an approach is sometimes employed when Monte Carlo simulations are used to solve the Boltzmann equation in the presence of a Poisson potential. It is referred to as the incident flux approach 4 .

Generalization to phonons
In principle, the Monte Carlo method described here could also be applied to phonon transport. Since equations similar to Eq. (7) and Eq. (8) can be found also in phonon NEGF formalism, most of the expressions in section II should carry directly over to the phonon method. An important exception however, is assumption (iv), which is a consequence of charge conservation. Unless a similar conservation condition can be shown to hold also for phonons, Monte Carlo estimates must be based on Eq. (51) rather than Eq. (55).

VI. CONCLUSIONS
We have shown that a Monte Carlo scheme similar to that used for the solution of the Boltzmann equation can also be used to solve the NEGF equations, as long as certain assumptions about the scattering self energies are satisfied. While the most general Monte Carlo method, described in section II C, requires only assumptions (i) and (ii) to be satisfied, the more specialized method described in section II D requires also assumption (iii). In addition, the simplified estimator of Eq. (55) is based upon assumption (iv). We tested the method with a scattering model conforming to Büttiker's approximation.
The tests show that at the level of accuracy required, the Monte Carlo method is faster than the alternative approaches in a specific regime. However, this regime involves mean free paths larger than what can be found in most applications. Our prototype implementation was performed to demonstrate the concept and present introductory performance numbers to guide future studies. As we have discussed, there are several ways in which the current implementation can be be improved to reduce the computation time. More work is thus needed to demonstrate how successful such improvements will be, and whether the discussed physical generalizations are applicable in practice.
Appendix A: Justification for the scattering assumptions In this appendix, we try to provide some justification for the assumptions (i) -(v) introduced in the main text. A formal proof that these assumptions are always valid, is not the scope of this appendix, and this may not even be the case. Instead, we propose heuristic arguments, in an attempt to justify why these assumptions are reasonable, or at least why they are reasonable as approximations.

Assumption (i), Linearity
In the case of electron-phonon scattering in the self consistent Born approximation, the self energy Σ < s can be expressed as 1,8,11 This expression can be written on the form of Eq. (16), with which is seen to be linear in G as required. Similar expressions apply for scattering on impurities and other forms of disorder, and indeed for self consistent Born approximations of scattering on any other field than the electron field itself.
In the case of electron-electron scattering, the lowest order perturbation would instead take the form and Λ would thus be second order in G. Clearly, other scattering mechanisms will also result in nonlinear Λ if higher order expansions than Eq. (A1) are used. However, in principle we can rewrite Λ as Λ(E, E ′ , G, G < ), where Λ is still linear in the third argument. For instance, Eq. (A3) could be written on the form of Eq. (16) with which is now linear in G. However, this is only useful if the fourth argument of Λ can be approximated in some other way. One way of doing this is to expand it in terms of the unperturbed Green's function g < . Alternatively, one could try the approach suggested in section V B 3.
Since the lead contributions to Eq. (10) are Hermitian, so must Σ < s be. But this just means that the integral in Eq. (16) is anti-Hermitian, and does not necessarily mean that where the potentials U r ij describe perturbations to the electronic Hamiltonian introduced by lattice perturbations x r , and D rs (E) = dt x r (t)x s (t ′ ) e iE(t−t ′ )/h , where the brackets now denotes an average over the lattice state. Eq. (A10) makes use of the assumption of a stationary process, so that x r (t)x s (t ′ ) is a function of t − t ′ alone.
The first step is to show that the matrix D(E) with elements given by Eq. (A10) is positively definite for all E. Consider a set of functions f r (t). We have We can also express the integral in the transformed variables T = (t + t ′ )/2 and ∆ = t ′ − t.
This means that under the assumption of a stationary process, assumption (iii) is at least satisfied for two very common scattering models, namely self consistent Born approximations of disorder and phonon scattering. Generalizing to higher orders is more difficult, but we can make a heuristic argument. By an argument very similar to that of Eq. (A11) to Eq. (A15), one can also show that −iG < (E) is positively definite for stationary processes. Then, since ψ|Σ < |ψ = ψ|(E − H − Σ a )G < (E − H − Σ r )|ψ > 0, −iΣ < must also be positively definite.
Due to the fact that Σ < and G < are related by Eq. (10) and Eq. (16), it seems likely that Λ should preserve positive definiteness. First of all, although it is possible that −iΣ < is positively definite while −iΣ < s is not, this seems rather controversial, since the lead contribution to Eq. (10) can be chosen arbitrarily. Secondly, this same freedom suggests that any Greens' function G < could in principle be a solution of Eq. (8), if the lead self energies were chosen appropriately. Thus, in order for −iΣ < s to always be positively definite, the integral of Λ(E, E ′ , G(E)) must be positively definite for any positively definite function G(E). This must also be the case for G(E) ∼ δ(E − E ′ ), leading to a positive definite Λ(E, E ′ , G) for any G. and accordingly, Λ(E, G) becomes so that Thus, assumption (v) would follow if we could say D ij kl (E) = D lk ji (E). In the case of electron phonon interactions, we have from Eq. (A9) and Eq. (A10) as required. Similar expressions apply for interactions with other fields, while for disorder scattering we have from Eq. (A7) where V is now the disorder potential. Thus, assumption (v) applies at least to these commonly used scattering models. = · · · = TrΞ n−k−1Â q (E)Λ(E, Ξ k A p (E)) = · · · = TrÂ q (E)Λ(E, Ξ n−1 A p (E)) = Tr G a (E)Γ q (E)G r (E)Λ(E, Ξ n−1 A p (E)) = Tr Γ q (E)Ξ n A p (E), and we thus finally obtain Tr A(E)Γ q (E) = ∞ n=0 p Tr Γ q (E)Ξ n A p (E).