Green functions and density matrix

#tl;dr#
The single particle density matrix is directly related to NEGF as shown here, I wish to find a way to relate NEGF also to density matrices which describe probability distribution of many body states, such as those encountered by solving Liouville Von Neumann based equations for many body problems.

#General question#

I’m wondering regarding the relation between Green functions and the density matrix. I wish to know if one can form a link between density matrix approaches based on the Von Neumann equation, and its open system counter-parts, i.e. Lindblad or Redfield equations.

###Motivation###

• Master equation approaches (such as Lindblad or Pauli rate equations) require some approximations: Born, Markov, and secular. These aren’t necessarily order consistent, and they have some well known subtle effects.
• I think I found a system in which some of these approximations yield undesirable side effects, and I’d like to try and preform the same calculation via NEFG for comparison.
• Even though I could compare expectation values for observables such as single site occupations, which are directly obtained in the GF formalism. A more general comparison would be to compare to the raw results of Lindblad’s equation. To that end I wish to find the density matrix which describes the probability distribution of Fock states via GF as this is what I get from Lindblad.
• As can be seen in several sources master equations can be obtained via diagrammatic expansion, which suggests to me that the diagrammatic technique and perhaps the NEGF formalism contains in some way the information I desire. For completeness such sources are: 1, 2, 3, 4, though I’m sure that there are more. Also the source I referenced below discusses the relation and shows that the relation between the two approaches is trivial when discussing a single particle problem, but how can one generalize it?

###Main Question###
The intro above hinted that I’m interseted in the reduced density matrix of a sub-system coupled to a thermal bath, and is out of equilibrium. Even though I’m actually interested in such a system, for the sake of simplicity let his now discuss a system with a finite number of sites at thermal equilibrium.

The occupation probability of a specific site (and hence the single particle density matrix) can written in terms of equal time lesser Green functions as shown here, by:

$$
P_i=langle d_i^dagger (t) ;d_i (t)rangle= -iG_{ii}^<(t,t)
$$

I wish to know:

1. Can I define a coherence (in the density matrix sense) in terms of lesser Green functions by $rho_{ij}=langle d_i^dagger (t) ;d_j (t)rangle= -iG_{ij}^<(t,t)$? Is this equivalent to the off-diagonals of the density matrix?

2. Can these definitions be generalized to many particle states? The standard definition of density matrix defines the probability distribution over the Hilbert space regardless of whether the states under discussion are single or many particle states. However I’m not sure how to generalize the Green function definition of probability to learn something about the occupation many particle states, it only teaches me about the occupation of a specific site? Can it be done, and if so how?

#Specific Example#
Let us discuss a specific example of a $N=3$ site system, of fermions. The fact that I choose to discuss a fermionic system helps me as I cannot have more than 3 particles in my system and my Hilbert space is finite dimensional: $dim(V)=2^3=8$.

Any linear operator which maps this space to itself can be expressed as on $8times8$ matrix, and that includes the Green function for the system. As discussed in the comments this includes also the creation and annihilation, adn projection operators onto states.

First, in equilibrium I may just write $G^r(omega)=(I(omega+ieta)-H)^{-1}$.

If on the other hand I wish to relate to NEGF for future use, I may write the equations of motion for the NEGF and solve them.

At this point I’m already confused because the matrix form will yield an $8times8$ matrix with 36 independent quantities (even though $G^r$ isn’t hermitian the entires above and below the diagonal aren’t really independent). However if I think of in terms of two point correlation functions: $G_{ij}^r(t)=-itheta(t)langle{c_i(t),c_j^dagger(0)}rangle$, I don’t have that many options.
What went wrong?

As a side note, in response comments (now in chat) the the various Fock states can be written as multiplication of creation operators acting on the vacuum $|Omegarangle$. For instance two of the 8 states will be:

$$
|1,1,0rangle=c_2^dagger c_1^dagger|Omegarangle
|0,0,1rangle= c_3^dagger|Omegarangle
$$

and the projection operators onto these states would be:

$$
P_{|1,1,0rangle}=|1,1,0ranglelangle1,1,0|
P_{|0,0,1rangle}= |0,0,1ranglelangle0,0,1|
$$
If one wishes to write these projectors as matrices and chooses the basis outlined above as the standard basis ,the trivial result is that the matrices are filled with $0$‘s except for a single $1$ somewhere on the diagonal according to the specific projector.