Non-equilibrium BBGKY Hierarchy from the Redfield Equation

A BBGKY-like hierarchy is derived from the non-equilibrium Redfield equation. Two further approximations are introduced and each can be used to truncate and solve the hierarchy. In the first approximation such a truncation is performed by replacing two-particle Green's functions (GFs) in the hierarchy by their values at equilibrium. The second method is developed based on the cluster expansion, which constructs two-particle GFs from one-particle GFs and neglects the correlation part. A non-equilibrium Wick's Theorem is proved to provide a basis for this non-equilibrium cluster expansion. Using those two approximations, our method of solving the Redfield equation, for instance, of an N-site chain of interacting spinless fermions, involves an eigenvalue problem with dimension $2^{N}$ and a linear system with dimension $N^2$ in the first case, and a nonlinear equation with dimension $N^2$ in the second case, which can be solved iteratively via a sequence of $N^2$ linear systems. Other currently available direct methods correspond to a linear system or an eigenvalue system with dimension $4^N$ plus an eigenvalue system with dimension $2^N$. As a test of the methods, for small systems with size N=4, results are found to be consistent with results made available by other direct methods. Although not discussed here, extending both methods to their next levels is straightforward. This indicates a promising potential for this BBGKY-like approach of non-equilibrium kinetic equations.


I. INTRODUCTION
A dynamic equation, such as Newton's equation or the Schrödinger's equation, describes only dynamical processes and strictly speaking it does not describe evolution towards thermal equilibrium, although sometimes it is used to do so together with a presumed thermal equilibrium distribution, such as in the derivation of Kubo formula of linear response theory 1 . Starting from the dynamic equation and its corresponding BBGKY hierarchy 2 the uniqueness of a stationary solution at thermal equilibrium, and particularly being the Boltzmann distribution, remains as an assumption but not a proved theorem 2 . The Redfield equation 3 , although it may show an unphysical transient process due to the use of the Markov approximation, on the other hand puts the description of dynamical evolution, thermal evolution towards the equilibrium state and evolution towards nonequilibrium stationary states (NESSs) under a common framework.
The Redfield equation has long been a standard tool in the study of the relaxation process in the theory of nuclear magnetic resonance 4 , optical spectroscopy and chemical dynamical systems 5 . Due to the fact that usually in those studies the central system of interest is modeled by a Hamiltonian with a very low dimension, the lack of an efficient algorithm to solve the Redfield equation, other than direct diagonalization and direct integration, is not a serious problem. However, when the Redfield equation is applied to transport calculations, the size of the system is usually much larger. It is then necessary to have a more efficient way to find the NESSs, without which many physical questions remain unanswered.
A famous example is the validity of the phenomenology law of transport. Advancement in science and technology has made it possible to build nanoscale electronic devices 6 , where classical phenomenological laws of transport such as Fourier's Law and Ohm's Law may not be valid any more 6 . In order to construct a theory of transport for mesoscopic or microscopic systems, and also to check under what circumstances those phenomenological laws hold, one may start from the first principles, i.e. the dynamic equation of classical or quantum systems, and add in as few extra assumptions as possible. The Landauer formula 7 makes use of scattering waves from the Schrödinger's equation but a biased distribution of those waves is assumed to 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 non-equilibrium stationary states (NESSs), and treats coupling between the two systems perturbatively 8 . For interacting systems, usually NEGF is used together with the density functional method 9 , 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 10 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 system and reservoirs within second-order perturbation. This generalization of the Redfield equation to transport studies was first implemented by Saito 11 and more or less followed by others [12][13][14] . 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 6 . 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 increasing 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 12,13 . 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.
A general Redfield equation can be cast into the following form 12 , where L HS ρ = −i [H S , ρ] is the Hamiltonian of the central system. Sometimes we also use where V in is a possible induced potential, for example 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 nonzero. 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 systems 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 12 or an eigenvalue system 14 with dimension 4 N plus an eigenvalue system with dimension 2 N . Currently, for interacting central systems, one usually can only solve the Redfield equation numerically up to N = 10 12-14 while for non-interacting systems, the equation with N around a hundred can be solved in terms of single-particle Green's functions (GFs) 11 . For example, for non-interacting systems, in Ref.11, a closed equation of single-particle GFs, was derived from the stationary Redfield equation Eq(2) and solved.
Our idea is basically to extend this GF based solution of the Redfield equation for non-interacting systems onto 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 single-particle GFs G 1 m † , m ′ , generally denoted as G 1 , will be coupled to the two-particle GFs also generally denoted as G 2 , which are then coupled to three-particle GFs, 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 directly the Redfield equation. The hierarchy has to be truncated first and then solved. In the rest part 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 value of G 2 at equilibrium for the unknown G 2 appearing in the first equation while the second expresses G 2 as combinations of G 1 via cluster expansion. After either one, 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 one is capable of dealing with relatively small systems but with large interaction strength while the second one can deal with much larger systems but with relatively small interaction strength.

II. DERIVATION OF BBGKY-LIKE HIERARCHY
For concreteness in presenting our general formulation, let us start from a Redfield equation describing an Nsite 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 set = 1, k B = 1, the lattice constant a = 1 and 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 Vin term in Eq(1), is neglected.
The corresponding Redfield equation reads 10,12 , 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 −iHS 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 operatorms. Detail is included in Appendix A. There we will also see that change of variable between summation over k and integration over energy in Eq (7) involves the density of states of the baths D α (ω). We combine this density of states together with coupling constant V α k , and set D α (Ω mn ) |V α kmn | 2 (see Appendix A for detail) as an overall constant, which is included in λ 2 .
When only a long-time steady-state solution ρ ∞ is of interest, we may derive a stationary form, Eq(2), from the kinetic Redfield equation. Furthermore, for a physical quantity of the central system with operator A, from Eq(2), we have generally wherem † α (m † α ) is the hermitian conjugate ofm α (m α ). All equations of GFs in the rest of this paper will be derived from this equation. For example the first and the second equation of the hierarchy can be derived from Note that since the set of all polynomials of c l , c † l forms a complete basis of the operator space, operatorsm are certain functions of polynomials of c l , c † l . Therefore, as expected G 1 is coupled to G 2 from Eq(9b), and possibly also G 3 or higher GFs from Eq(9c), and G 2 is coupled to G 3 from Eq(10b), and possibly also G 4 or higher GFs from Eq(10c). Solving such an equation hierarchy is no 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 certain order with some further approximations, such as the molecular-chaos assumption in the classical Boltzmann equation 10 , or replacing high order GFs by cluster expansion of lower order ones 15,16 . In this work, we suggest the following two approximate methods: (1) substitution of certain high-order GFs by their values at equilibrium; (2) expressing high-order GFs as combinations of lower-order ones plus a correlation part via cluster expansion and then ignoring the correlation part at 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 preformed 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 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.

III. SOLVING THE HIERARCHY
In order to solve Eq(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 Eq(9) to be a closed equation by dealing with G 2 terms in the equation differently.
We will first discuss a more accurate even for a large V 0 but computationally costly method, 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 later works only for relatively small V 0 but it can be applied on 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 Eq(2). A measure of relative distance between two matrices A and B, is used to describe the accuracy of our approximations.
A. Method 1: Starting from Equilibrium States As explicitly worked out in Eq(A2) in Appendix A, operatorms can be written in eigenmodes of H S , which can be solved from an exact diagonalization of H S , a 2 Ndimension eigenvalue problem. Then in the language of super-operator space 14 , where operators are treated like vectors -so called super-vectors, super-vectorsm can be expanded under the basis -polynomials of c l , c † l , Here we keeps only the linear polynomial in the present work although further expansion is possible. Using the definition of inner product between super-vectors A|B = tr A † B , we have and operators V 0 D and V 0D are just the rest part of operatorsm andm respectively. Here 2 N −1 is a normalization constant to make d α,l = 1 whenm α = c l . With above expressions of operatorsm, Eq(9) becomes, Notice 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 Eq(14c) by their values at equilibrium, denoted here as G E,(0) 2 where the superscript (0) means to use the thermal equilibrium (denoted by superscript E) as the zeroth order approximation 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 Eq(14d), using the first term as an example, Next let us calculate G from Eq (14), where the superscript (1) means that the approximate calculation takes care of the first equation of the hierarchy, Eq (14). We organize all G E,(1) 1 m † , n as a vector, then Eq (14) for given value of m, n is the equation occupying the (mN + 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 Eq (14) for all m, n then becomes a linear system on g where vector ν comes from ordering Eq(14b) in the same way as g correspondingly from ordering Eq(14c) and Eq(14d), and from Eq(14a) one gets matrix Γ (1) . For example, assuming m, n are not at boundaries, one may read from Eq (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 . Here G Ex m † , n = tr c † m c n ρ ∞ , where ρ ∞ is the exact solution from Eq(2).

results
First, we set V 0 = 0.2 as a constant, and check the accuracy of g  with different values of V 0 . The worst case is d (1) = 0.3% as shown in Fig.1(b). Overall, d E,(1) is always much smaller than d E,(0) . J E,(0) , J E,(1) , J Ex are calculated respectively from g E,(0) 1 , g E,(1) 1 and g Ex 1 . From Fig.1(c) and (d) we see that in both cases, J E,(1) is very close to the exact one J Ex while J E,(0) 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 non-equilibrium stationary states. It is also worth mentioning that this method generates reasonable results for very large V 0 . Furthermore, it is likely the approximation could be improved: from further expansions in terms of higher order polynomials of c l , c † l and substituting 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 the 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 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 λ 2 ∆T term in λ 2 ν in Eq (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 In order to get some information on how such approximation at the next order improves the accuracy, we also want to compare ∆ which is estimated in the same way from the difference between the equations respectively for for g Ex 1 and g E,(0) 1 . See Appendix C for detail of those equations and the estimation. Here we summarize the results that and ,T g Ex Here ∆ We refer 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. Judged from this term, as long as λ 2 V 0 is a small number compared with t then the method under consideration is very reasonable. As of 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 V0 t where this method is still applicable is by 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 -the smaller the larger n. This explains why as we see from Fig.1 that 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 on accuracy has been found. More detail and more systematic analysis will be presented elsewhere.

B. Method 2: Non-equilibrium Cluster Expansion
Another way to make Eq(9) to be 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 on higher-order GFs, for example by a similar expansion on G 3 and setting G 3 = 0. The fact that 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 non-equilibrium GFs, we are going to propose the same expansion and this requires a non-equilibrium Wick's Theorem to be true that G 2 = 0 when V 0 = 0 holds for non-equilibrium GFs. Fortunately, this can be proved (see Appendix B for detail). Setting G 2 = 0 with V 0 = 0 is like 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 could not be applied, we have to expand operators operatorsm in higher-order polynomials of c l , c † l . This can be done as following. Other than the exact direct diagonalization, operatorsm can also be found perturbatively analytically(see Appendix A for the full detail). The basic idea is to start from assuming c l (t) = c and then derive and solve the equations of motion of c (0) l , c (1) l from the 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 . Kept only the first order, operatorsm becomê where definitions of D α;m and D α;m1m2m3 are given in Appendix A.

results
Similarly we define ∆ C,(0) 1 (d C,(0) ) as the absolute (relative) distance between g In both cases, d C,(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 V0, as long as V0 is not too large, J C,(0) already provides a major part. From (d) we find that for relatively larger V0, 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 V0. with different values of ∆T . From Fig.2(a) we can see that the worst case is about d (1) = 1%. Secondly, we set ∆T = 0.4T as a constant, and check the accuracy of g C,(1) 1 with different values of V 0 . The worst case is about d (1) = 2% as shown in Fig.2(b). Overall, d C, (1) , is always much smaller than d C,(0) . From Fig.2(c) we see that for a small V 0 , J C,(0) already provides a major part. However, in Fig.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 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 λ 2 and V 0 are small. Define similarly ∆ and g Ex 1 , and then compare the three equations while ignoring certain higher-order terms such as terms which are proportional to λ 2 V 0 . See Appendix C for detail of those equations and the estimation. We arrived at, and This agrees with the numerical tests that ∆ We refer readers to Appendix C for definitions of all Γ matrices. Most importantly, we see again that ∆ . 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 small number compared with t then 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 diag-onalization of a 2 N -dimension matrix. Although in above examples in order to be comparable with exact solution, only N = 4 is used, on a normal PC a system N = 100 has been successfully tested.

IV. 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, non-equilibrium stationary states of interacting systems are calculated. It is found that they are consistent with results made available by other direct methods. We also estimate accuracy of the two approximations. The difference among the two is also discussed that the first method is applicable for large V 0 while the second method is more efficient so it can be applied to larger systems. We have not tested performance of further orders of both methods, although it seems quite straightforward. Our method can also be applied to the local-operator Lindblad equation 14,17 . It may also be worth pursuing a comparison between results from the Lindblad equation and from the Redfield equation. Besides its application on non-equilibrium stationary states, 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 relations between these methods and the NEGF method may also be interesting.

Appendix A: Perturbation decomposition ofms
From its definition in Eq (7), in the representation of eigenvalues {E m } and eigenstates {|m } of H S , operator m can be written as where we have used ∞ 0 dτ e iωτ = πδ (ω) + iP 1 ω and neglected the principal value part. We have also assumed that it is possible to perform a change of variable on V α k such that it becomes V α (k nm ), where k mn is defined by ω kmn,α = Ω mn , i.e. a bath mode resonant with this transition. This limits the possible forms of V α k and ω k,α . For example, for a given energy Ω mn , there should be a unique value of V α knm . In this work, we take V α k as a constant so this condition is satisfied. D α (ω) is the bath's density of states. We arrived at where Ω mn = E m − E n = −Ω nm . We furthermore set V α knm 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 operatorŝ m perturbatively.