Non-equilibrium stationary states from the equation of motion of open systems

Non-equilibrium stationary states are solved from the equation of motion of open systems via a Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY)-like hierarchy. As an example, we demonstrate this approach for the Redfield equation, which is derived from the whole dynamic equation of motion of a central system coupled to two baths through the projector technique. Generalization to other equations of motion is straightforward. For non-interacting central systems, the first equation out of the hierarchy is closed, whereas for interacting systems it is coupled to higher-order equations in the hierarchy. In the case of interacting systems, two systematic approximations, in the form of perturbation theories, are proposed to truncate and solve the hierarchy. A non-equilibrium Wick's theorem is proved to provide a basis for the perturbation theories. As a test of reliability of the proposed methods, we apply them to small systems, where it is also possible to apply other exact direct methods. Consistent results were found.

3 calculate physical quantities. The non-equilibrium Green's function (NEGF) method starts from two decoupled systems each at their own thermal equilibrium states, which could be very far from the expected NESSs, and treats coupling between the two systems perturbatively [9]. For interacting systems, usually NEGF is used together with the density functional method [10], from which the whole spectrum and corresponding effective wavefunctions are calculated to construct the non-equilibrium density matrix. Both the perturbation and the effective wavefunctions introduce further approximations to the calculation.
The Redfield equation approach to studying transport phenomena is to explicitly couple the system of interest, H S , to reservoirs and then use the projector technique [11] to derive an effective equation of motion for the system. One then solves this Redfield equation to get NESSs. The Redfield equation requires validity of the Markovian approximation and treats the coupling between the system and reservoirs within second-order perturbation. This generalization of the Redfield equation to transport studies was first implemented by Saito et al [12] and more or less followed by others [13]- [15]. In such cases, one quite often expects to deal with systems with size N ≈ 100 in order to be comparable to other methods, such as the Landauer formula and NEGF [7]. This then leads to a Hilbert space with dimension 2 100 , when we, for example, consider an N -site chain of spinless fermions with N = 100, meaning a matrix with dimension 4 100 in the corresponding Redfield equation. Due to this exponential increase of the problem size and the lack of a more efficient alternative method, currently one can only discuss physical behaviors of very small systems with N ≈ 10 [13,14]. Other propagator-based direct methods [6] make it possible to analyze only slightly larger systems. Conclusions drawn from numerical results for such small systems are not regarded as reliable enough. In this work, we will present a new approach, using the idea of a BBGKY hierarchy to study this kinetic equation. This allows us to develop very efficient and systematic approximate methods to find NESSs. In fact, besides these direct methods, the stochastic wave function method [16] has been generalized and used to solve the Redfield equation [17]. Such stochastic methods are capable of dealing with slightly larger systems, say N of around 20. It will be an interesting topic of future work to compare the performance of our method against the stochastic methods.
A general Redfield equation can be cast into the following form [13], where L H S ρ = −i[H S , ρ] is the Hamiltonian of the central system. Sometimes we also use H S = H 0 + V to separate H S into a non-interacting part H 0 and an interaction V . Finally, where V in is a possible induced potential; for example, the electric potential due to charge distribution in the case of charge transport. L B (T, µ) comes from coupling to baths with coupling strength λ, and L P ( T, δµ) exists when baths have different temperatures and/or chemical potentials. This equation describes a dynamical process when λ = 0, thermal relaxation towards equilibrium when λ = 0, T = 0, µ = 0 and evolution towards NESSs when they are all non-zero. If we are interested in the long-time steady-state solution ρ ∞ , then which is sometimes called a stationary Redfield equation. Here, L is a matrix with dimension d 2 , where d is the dimension of the system's Hilbert space; for example, d = 2 N for the N -site chain of spinless fermions mentioned above. Solving this equation is very costly computationally.
It involves solving a linear system [13] or an eigenvalue system [15] with dimension 4 N plus an eigenvalue system with dimension 2 N . Currently, for interacting central systems, one can usually only solve the Redfield equation numerically up to N = 10 [13]- [15], whereas for non-interacting systems, the equation with N equal to about 100 can be solved in terms of single-particle Green's functions (GFs) [12]. For example, for non-interacting systems, in [12], a closed equation of single-particle GFs, , was derived from the stationary Redfield equation (2) and solved.
Our idea is basically to extend this GF-based solution of the Redfield equation for non-interacting systems to interacting systems. We consider G 1 (m † , m ), GFs in the lattice basis, and derive an equation for these GFs. In the presence of interaction, the singleparticle GFs G 1 (m † , m ), generally denoted as G 1 , will be coupled to the two-particle GFs n c m c n , also generally denoted as G 2 , which are then coupled to three-particle GFs, n c l c m c n , also generally denoted as G 3 and so on. We arrive at a BBGKY-like equation hierarchy [2].
In principle, solving the whole hierarchy is as hard as solving the Redfield equation directly. The hierarchy has to be truncated first and then solved. In the rest of this paper, after deriving the hierarchy from the Redfield equation, we will then present two such methods. In the example calculations, we will only truncate the hierarchy at the first equation of the hierarchy. The first method substitutes the value of G 2 at equilibrium for the unknown G 2 appearing in the first equation, whereas the second expresses G 2 as combinations of G 1 via cluster expansion. After either method, the equation is closed and then solved. We will see in the following that both methods are significantly more efficient than the direct methods. The first is capable of dealing with relatively small systems but with large interaction strength, whereas the second can deal with much larger systems but with relatively small interaction strength.

Derivation of a Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY)-like hierarchy
For concreteness in presenting our general formulation, let us start with a Redfield equation describing an N -site chain of spinless fermions coupled with two fermionic baths. Our system of interest is defined by H S , The two heat baths are collections of fermionic modes, where α = L, R indexes the left-and right-side baths and we seth = 1, k B = 1, the lattice constant a = 1 and the hopping constant t = 1. The system-baths coupling is chosen as where the left (right) bath is coupled to the first (last) site: c L = c 1 and c R = c N and so on. Bath parameters, including temperature and chemical potential, are chosen to be (T L , µ) and (T R , µ) with T L/R = T ± T 2 . In the present work, the induced L V in term in equation (1) is neglected.
The corresponding Redfield equation reads [11,13] ∂ρ(t) ∂t Here, n(ω k,α ) = (e β α (ω k,α −µ α ) + 1) −1 is the Fermi-Dirac distribution with the bath temperature T α = 1/β α and chemical potential µ α . If U (t) = e −iH S t is known, then so is c α (t) = U † (t)c α U (t) and therefore the operatorsm. This requires a full diagonalization of H S . Using eigenvectors of H S , one can perform the above integrals to get the operatorms. More details are included in appendix A. There we will also see that a change of variable between summation over k and integration over energy in equation (7) involves the density of states of the baths D α (ω). We combine this density of states together with the coupling constant V α k , and set D α ( mn )|V α k mn | 2 as an overall constant, which is included in λ 2 (see appendix A for more details).
When only a long-time steady-state solution ρ ∞ is of interest, we may derive a stationary form, equation (2), from the kinetic Redfield equation. Furthermore, from equation (2), for a physical quantity of the central system with operator A, we generally have wherem † α (m space, operatorsm are certain functions of polynomials of {c l , c † l }. Therefore, as expected, G 1 is coupled to G 2 from equation (9b), and possibly also G 3 or higher GFs from equation (9c), and G 2 is coupled to G 3 from equation (10b), and possibly also G 4 or higher GFs from equation (10c). Solving such an equation hierarchy is not easier than directly solving the Redfield equation, unless V 0 = 0, so that the above equation of G 1 is closed and is not coupled to G 2 .
We may, however, solve these equations by truncating the hierarchy at a certain order with some further approximations, such as the molecular-chaos assumption in the classical Boltzmann equation [11], or replacing high-order GFs by cluster expansion of lower-order ones [18,19]. In this work, we suggest the following two approximate methods: (i) substitution of certain high-order GFs by their values at equilibrium; and (ii) expression of high-order GFs as combinations of lower-order ones plus a correlation part via cluster expansion and then ignoring the correlation part at a certain order. Specifically in the following example calculation, the first-order form of both approximations (i.e. only the first equation of the hierarchy) is used and substitution or cluster expansion is performed on G 2 . One can do such a substitution or cluster expansion of GFs at further-order GFs and make use of further equations in the hierarchy. A general discussion of the accuracy of such substitutions at different orders will be presented elsewhere. In this work, we focus on the potential of this BBGKY-like formulation and discuss briefly the topic of performance of the two approximations in their first-order forms.

Solving the hierarchy
In order to solve equation (9) explicitly, we will first have to find explicit forms of operatorsm in terms of operator {c l , c † l }. In the following, we will present one exact numerical calculation and one perturbative calculation of those operators. Correspondingly, based on these two methods of finding operatorsm, we will discuss in this section two ways of making equation (9) a closed equation by dealing with G 2 terms in the equation differently.
We will first discuss a more accurate method, even for a large V 0 , but computationally costly: perturbation based on two-particle GFs at equilibrium. Next we will discuss a relatively less accurate but computationally much cheaper method, the non-equilibrium cluster expansion. The latter works only for relatively small V 0 , but it can be applied to much larger systems. The unknown G 1 (m † , n) solved from both methods will be compared against G Ex 1 (m † , n), the exact solution of equation (2). A measure of the relative distance between two matrices A and B, is used to describe the accuracy of our approximations.

Method 1: starting from equilibrium states
As explicitly worked out in equation (A.2) in appendix A, operatorsm can be written in eigenmodes of H S , which can be solved from an exact diagonalization of H S , a 2 N -dimension eigenvalue problem. Then in the language of super-operator space [15], where operators are treated like vectors (so-called super-vectors), super-vectorsm can be expanded on the basis of Here we keep only the linear polynomial in the present work, although further expansion is possible. Using the definition of the inner product between super-vectors A|B = tr(A † B), we have and operators V 0 D and V 0D are just the remaining part of operatorsm andm, respectively. Here, 2 N −1 is a normalization constant to make d α,l = 1 whenm α = c l . With the above expressions for operatorsm, equation (9) becomes Note that every c 0 , c † 0 , c N +1 and c † N +1 that appears in the equation should be recognized as 0. First, let us replace all G 2 s in equation (14c) by their values at equilibrium, denoted here as G E,(0) 2 where the superscripts indicate the usage of thermal equilibrium (denoted by superscript E) as the zeroth-order approximation (denoted by (0)) of the non-equilibrium G 2 . Using the first term as an example, where ρ eq (H S ) = 1 Z e −H S /T . This requires eigenstates of H S . Similarly, one can define G E,(0) D from equation (14d), using the first term as an example, Next, let us calculate G E,(1) 1 from equation (14), where the superscript (1) means that the approximate calculation takes care of the first equation of the hierarchy, equation (14). We organize all G E,(1) 1 (m † , n) as a vector, then equation (14) for given values of m, n is the equation occupying the (m N + n)th row and totally there are N 2 such equations. After substituting G E,(0) 2 and G E,(0) D for the exact but unknown G 2 and G D , the whole set of equation (14) for all m, n then becomes a linear system on g E,(1) where the vector ν comes from ordering equation (14b) in the same way as g E,(1) 1 . The same holds for g E,(0) 2 and g E,(0) D correspondingly from ordering equations (14c) and (14d), and from equation (14a) one gets the matrix (1) . For example, assuming m and n are not at boundaries, one may write from equation (14) ν Next, we calculate single-particle equilibrium GFs, G E,(0) 1 , and organize it in the same way into a vector denoted as g E,(0) 1 . In order to set a reference of the accuracy, we compare d E,(0) , the difference between the exact solution g Ex 1 and the zeroth order g E,(0) 1 , and d E, (1) , the difference between the exact solution g Ex 1 and the first order solution above, g E,(1) where ρ ∞ is the exact solution from equation (2). with different values of V 0 . The worst case is d (1) = 0.3%, as shown in figure 1(b). Overall, d E, (1) is always much smaller than d E,(0) . J E,(0) , J E, (1) and J Ex are calculated, respectively, from g E,(0) 1 , g E,(1) 1 and g Ex 1 . From figure 1(c) and (d) we see that in both cases, J E, (1) is very close to the exact, J Ex , while J E,(0) , the current in the equilibrium state, is always zero. Very high accuracy is found, especially for small T . This indicates that the approximation captures the essential part of the NESSs. It is also worth mentioning that this method generates reasonable results for very large V 0 . Furthermore, it is likely that the approximation could be improved: from further expansions in terms of higher-order polynomials of c l , c † l and substitution of their values at equilibrium for the higher-order unknown GFs in higher-order equations of the hierarchy. Stopping the expansion of operatorsm at linear order of V 0 is compatible with solving only the first equation of the hierarchy. If further equations of the hierarchy are used, then one should also expand operatorsm in further orders of V 0 .
In order to estimate the accuracy of the first-order form of this approximation, and also to get an overview of the accuracy of possibly the next order, let us study the leading order of residues in terms of λ 2 and T T , which are assumed to be small in the following. Hence λ 2 V 0 V 0 ; therefore we know that g D is relatively smaller than the other g 2 term, so we drop it. This in fact requires λ 2 V 0 T , which we assume to be true. Similarly, for the same reason, since λ 2 T T , we drop the λ 2 T term in λ 2 ν in equation (18), and keep only the major term, λ 2 ν 0 (T ), which is independent of T . Here, ν ,T denotes formally a derivative of T on ν. The general idea is then to write down, respectively, equations for g Ex 1 and g E,(1) 1 , and then compare the two equations to estimate E,(1) In order to get some information about how such an approximation at the next order improves the accuracy, we also want to compare E,(1) which is estimated in the same way from the difference between the equations, respectively, for g Ex 1 and g E,(0) 1 . See appendix C for further details of those equations and the estimation. Here we summarize the results that and Here, E,(0) We refer the readers to appendix C for definitions of all matrices. Most importantly, here we see that E,(0) 1 is multiplied by a small number λ 2 V 0 and then becomes a part of E,(1) 1 . Furthermore, this relation holds generally for higher-order forms of this approximation. Judging from this term, as long as λ 2 V 0 is a small number compared with t, the method under consideration is very reasonable. As for the other two additional terms, they can be regarded as (V 2 0 g Ex 3 + V 0 g Ex 2 ) T . Therefore, the limit of V 0 t where this method is still applicable goes as |g Ex 2 | −1 or |g Ex 3 | −1/2 , which could be much larger than 1 since roughly |g Ex n | = |g Ex | n -smaller for larger n. This explains why, as we see from figure 1, this method is applicable even for V 0 larger than t. We have also tested several systems with larger N (up to N = 8) and no qualitative difference in accuracy has been found. More details and a more systematic analysis will be presented elsewhere.

Method 2: non-equilibrium cluster expansion
Another way to make equation (9) a closed equation is to use cluster expansion. In the case of equilibrium GFs, it proposes, for example, at the level of two-particle GFs, and then sets It can be applied to higher-order GFs, for example, by a similar expansion on G 3 and setting G 3 = 0. The fact that the equilibrium Wick's theorem shows that indeed G 2 = 0 when V 0 = 0 makes this expansion plausible for equilibrium GFs. Here, in the case of NEGFs, we are going to propose the same expansion, and this requires a non-equilibrium Wick's theorem to be true so that G 2 = 0 when V 0 = 0 holds for NEGFs. Fortunately, this can be proved (see appendix B for complete details). Setting G 2 = 0 with V 0 = 0 is similar to using the Hartree-Fock approximation, so that, depending on the system and the physical problem under investigation, one may quite often need to go beyond that to the next level of approximation, i.e. keeping G 2 but ignoring G 3 , and truncating the equation hierarchy at the second equation instead of the first equation. In this work, we will use the first-level approximation, i.e. ignoring G 2 .
However, with operators D defined in the previous section, cluster expansion cannot be applied; instead we have to expand operatorsm in higher-order polynomials of {c l , c † l }. This can be done as follows. Other than the exact direct diagonalization, operatorsm can also be found perturbatively analytically (see appendix A for full details). The basic idea is to start by assuming and then to derive and solve the equations of motion of c (0) l , c (1) l from Heisenberg's equation. In this way, one avoids the direct diagonalization of H S so that it simplifies the calculation, but its accuracy depends on the order of V 0 at which the expansion stops. Stopping at the linear order of V 0 is compatible with the cluster expansion at G 2 . If cluster expansion at higher-order GFs is applied, then operatorsm should also be expanded in higher orders of V 0 . Keeping only the first order, operatorsm becomê where definitions of D α;m and D α;m 1 m 2 m 3 are given in appendix A.
With the first-order cluster expansion and the above expansion of operatorsm plugged into equation (9), we have where the five terms are, respectively, the five subequations in equation (27); for example, This equation can be solved iteratively, where we start from g (0) 1 = g C,(0) 1 , which is the exact solution of equation (28) when V 0 = 0. Through the iteration defined above, we get the solution g C,(1) 1 , which in practice stops at large enough n such that g (n) 1 − g (n−1) 1 is small enough.

Results.
In a similar manner, we define C,(0) 1 (d C,(0) ) as the absolute (relative) distance between g C,(0) 1 and g Ex 1 , and C,(1) (1) ) as the absolute (relative) distance between g C,(1) 1 and g Ex 1 . Firstly, we set V 0 = 0.2 as a constant, and check the accuracy of g C,(1) 1 with different values of T . From figure 2(a) we can see that the worst case is around d (1) = 1%. Secondly, we set T = 0.4T as a constant, and check the accuracy of g C,(1) The worst case is around d (1) = 2%, as shown in figure 2(b). Overall, d C, (1) is always much smaller than d C,(0) . From figure 2(c), we see that for a small V 0 , J C,(0) already provides a major contribution. However, in figure 2(d) when V 0 becomes larger, the difference between J C, (1) and J C,(0) becomes more important. We should also note that for larger V 0 , J C, (1) starts to deviate from J Ex . This indicates that the approximation captures the essential part of the interaction, but it is more accurate for small V 0 . Furthermore, it is likely that the approximation could be improved: by keeping G 2 but ignoring correlations in higher-order GFs and calculating perturbatively further terms in V 0 in operatorsm.
In order to estimate the accuracy of this approximation, let us assume that λ 2 and V 0 are small. Define similarly C,(0) n = g C,(0) n − g Ex n and C,(1) n = g C,(1) n − g Ex n . Again, we start from the equations of the three, g C,(0) 1 , g C,(1) 1 and g Ex 1 , and then compare the three equations while ignoring certain higher-order terms, such as terms that are proportional to λ 2 V 0 . See appendix C for more details of those equations and the estimation. We arrive at and This agrees with the numerical tests that C,(0) 1 is proportional to V 2 0 . We refer the readers to appendix C for definitions of all matrices. Most importantly, we see again that C,(0) 1 is multiplied by a small number λ 2 V 0 and then becomes a part of C,(1) 1 . The other term, V 2 0 g Ex 3 ∼ V 2 0 |g Ex 1 | 3 , since roughly |g Ex n | = |g Ex | n , is also much smaller than C,(0) 1 ∼ V 0 |g Ex 1 | 2 . However, for large enough V 0 , the other approximation used in this method, the perturbation expansion of operatorsm, will be invalid. Therefore, as long as V 0 is a  (1) is always much smaller than d C,(0) . J C,(0) and J C, (1) are compared against J Ex in (c) and (d). From (c) we see that for a given value of V 0 , as long as V 0 is not too large, J C,(0) already provides a major contribution. From (d) we find that for a relatively larger V 0 , the difference between J C, (1) and J C,(0) becomes more important. At the same time, the difference between J C, (1) and J Ex also becomes larger for larger V 0 .
small number compared with t, the method under consideration is very reasonable. It should be noted that this method is capable of dealing with large systems since it does not involve a direct diagonalization of a 2 N -dimension matrix.

Conclusion and discussion
To conclude, a BBGKY-like equation hierarchy is derived from the Redfield equation and two systematic approximations are suggested to solve the hierarchy. Using the first-order form of the two methods, NESSs of interacting systems are calculated. It is found that they are consistent with results made available by other direct methods. We also estimate the accuracy of the two approximations. The difference between the two is also discussed, showing that the first method is applicable for large V 0 whereas the second method is more efficient, so it can be applied to larger systems. We have not tested the performance of further orders of both methods, although it seems quite straightforward. Our method can also be applied to the local-operator Lindblad equation [15,20]. Also, it may be worth pursuing a comparison between results from the Lindblad equation and from the Redfield equation. Besides their applicability to NESSs, these methods may also be valuable for perturbation theory on equilibrium states. There, using the second method, the equilibrium interacting g C,(1) 1 can be calculated starting from the equilibrium non-interacting g C,(0) 1 by setting T L = T R and µ L = µ R . The computational cost, being a sequence of linear systems with dimension N 2 , is obviously cheaper than direct diagonalization. It will be interesting to see a further comparison of accuracy and efficiency between these methods and other perturbative methods on equilibrium states. The non-equilibrium equal-time GFs calculated by the proposed methods are also objects of the NEGF method. Investigating the relation between these methods and the NEGF method may also be interesting. where mn = E m − E n = − nm . We furthermore set V α k nm D α (ω) as a constant and absorb it into λ 2 . This procedure involves a direct diagonalization of the isolated system H S . One can avoid this by finding such operatorsm perturbatively.
Next, assuming V 0 is small, we want to express operatorm α in terms of {c m } and V 0 . When V 0 = 0, the system is a tight-binding open chain, the following basis transformation Therefore, c α (t) is a linear function of all c m , Hence,m α is also a linear combination of all c m s. One can imagine that for small V 0 ,m L should not be too far from a linear combination. Denote c l (t) when V 0 = 0 as c (0) l (t). Starting from treating this as the zeroth order to the full dynamical c l (t), and expanding we may derive a perturbative equation of c (n) where the shorthand notation, for non-negative integers n, n 1 , n 2 and n 3 , Then the solution of the above equation can be written as Here the initial condition that c (n) (0) = 0 (∀n 1) is used. Plugging this general solution into equation (7), and after straightforward but tedious algebra, we arrive at the decomposition of m α andm α in equation (26) with expansion coefficients defined as × sin k 2 π(m + 1) N + 1 sin k 3 π(m + 1) N + 1 + sin k 2 π(m − 1) N + 1 sin k 3 π(m − 1) N + 1 . (A.11c)