Non-interacting causal Green function in localized-orbitals representation Say I have a one-partilce hamiltonian 
$$\hat{h}=\sum_{\alpha} \epsilon_\alpha \hat{n}_\alpha+\sum_{\alpha \neq \beta, } t_{\alpha,\beta} \hat{c}^\dagger_{\alpha }\hat{c}_{\beta }$$
(I will ignore spin for simplicity, since it does not play a relevant part in my doubt), where the $\hat{c}_\mu,\hat{c}_\mu^\dagger$ are anihilation/destruction operators of electrons in states $\vert \chi_\mu \rangle .$
Quick version of my question: 
What is the causal Green function i frequency space for this problem?
Detailed question and problems:
I want to write down explicitely what are the Green functions $g(\mu,t)= \frac{-i}{\hbar} \langle T [\hat{c}_\mu(t) \hat{c}_\mu^\dagger(t^\prime)]\rangle$ in frequency space, $g(\mu,\omega),$ and I want to do it by two methods: 1) Applying the direct definition and explicitely computing the Fourier transform and 2) By using that, if we call $h$ to the matrix of elements $\langle \chi_\nu \vert \hat{h} \vert \chi_\mu \rangle,$ it holds that $g(\mu,\omega)=(\omega \ \text{Id}-h)^{-1}(\mu,\mu)$ (this holds for one-electron hamiltonians and shows that, in that case, the notion of Green function in many-body theory coincides with the concept of Green function in pure mathematics).
1) First method: For simplicity, we assume $t^\prime =0$ and I will drop the hats from the operators, since there's no possible confussion We have:
$\langle T [\hat{c}_\mu(t) \hat{c}_\mu^\dagger(0)]\rangle= \Theta(t)\langle c_\mu(t)c_\mu^\dagger\rangle-\Theta(-t) \langle c_\mu^\dagger c_\mu(t) \rangle.$
We insert in both terms the resolution of the identity, $\text{Id}=\sum_{m} \vert m \rangle\langle m \vert$, where $m$ is an index which runs over all eigenstates of the system. 
On the other hand, recall that $c_\mu(t)=e^{i h t}c_\mu e^{-i h t}$ (setting $\hbar=1$) and that $e^{-i h t} \vert m \rangle = e^{-i E_m t} \vert m \rangle.$ Keeping this in mind is easy to get to:
$$\Theta(t)\sum_{m}e^{i(\omega_0-\omega_m)t}\langle 0 \vert c_\mu \vert m \rangle \langle m \vert c_\mu^\dagger \vert 0 \rangle-\Theta(-t)\sum_{m}e^{i(\omega_m-\omega_0)t}\langle m \vert c_\mu \vert 0 \rangle \langle 0 \vert c_\mu^\dagger \vert m \rangle.$$
Next, to Fourier transfrom this expression we notice that for the first term we need a convergent factor $i\eta$ (that is, we'll have pole sin the upper semi-plane), and for the second term we need to plug in a convergent factor $-i\eta$ (so the poles lie in the lower semi-plane). 
After computing those Fourier transforms and replugging the prefactor $-i$, we get:
$$\sum_ {m} \dfrac{\langle 0 \vert c_\mu \vert m \rangle \langle m \vert c_\mu^\dagger \vert 0 \rangle}{\omega-(\omega_m-\omega_0)+i\eta}+\sum_{m}\dfrac{\langle m \vert c_\mu \vert 0 \rangle \langle 0 \vert c_\mu^\dagger \vert m \rangle}{\omega-(\omega_0-\omega_m)-i\eta}.$$
If the Ground state $\vert 0 \rangle$ is a state of $N$ electrons, in the first term the sum runs over eigenstates $\vert m \rangle $ of $N+1$ particles and in the second term the sum runs over states with N-1$ particles. 
But I don't know hot to proceed any further. Intuitively speaking, it seems that in the first term, the term in the denominator $\omega_m-\omega_0$ should be $\epsilon_\mu$ in case all the $t_{\alpha,\beta}=0,$ and perhaps $\epsilon_mu+\sum_\beta t_{\mu,\beta}$ in the general case (but I'm not sure), since this should be the energy difference between $\vert 0 \rangle $ and $\vert m \vert.$ For the denominator in the second term, the reaosoning is similar and leads (or it would lead, would it be clear) to the same result, that's for sure. Also, the numerator in the first term is about "when an electron in state $\vert \chi_\mu \rangle $ is not in $\vert 0 \rangle,$ and in the second term the numerator is about "when an electron in state $\vert \chi_\mu \rangle$ is in $\vert 0 \rangle $", so it seems that the final result should be (probably) something like (here we call $n_\mu=\langle 0 \vert \hat{n}_\mu \vert 0 \rangle$ to the occupation number of state $\vert \mu \rangle $)
$$g(\mu,\omega)=\dfrac{1-n_\mu}{\omega-\epsilon_\mu+i\eta}+\dfrac{n_\mu}{\omega-\epsilon_\mu-i\eta},$$
or perhaps something like this but with $\epsilon_\mu \rightarrow \epsilon_\mu+\sum_{\beta} t_{\mu,\beta}$
2) Second method:  A direct inversion $(\omega \ \text{Id}-h)^{-1}$ is not well-defined because of all the issues with the poles. We need to insert convergent factors. $((\omega+i\eta) \ \text{Id}-h)^{-1}$ and $((\omega-i\eta) \ \text{Id}-h)^{-1}$ are the retarded and advanced Green function respectively, and in virtue of the identity $\frac{1}{x\pm iy}= P.V \frac{1}{x} \mp i\pi \delta(x)$ as $y \rightarrow 0,$ it is clear that by adding up the two of them we should get theGreen functions:
$$g(\mu,\omega)=((\omega+i\eta) \ \text{Id}-h)^{-1}(\mu,\mu)+((\omega-i\eta) \ \text{Id}-h)^{-1}(\mu,\mu).$$
But, in this case, for the simplest case possible, a diagonal hamiltonian, where all the $t_{\alpha,\beta}$ are zero, so that $\hat{h}=\sum_{\alpha}\epsilon_\alpha \hat{n}_\alpha,$ the formula above trivially yields the result:
 $$ g(\mu,\omega)= \dfrac{1}{\omega-\epsilon_\mu+i\eta}+\dfrac{1}{\omega-\epsilon_\mu-i\eta},$$
with no $n_\mu,1-n_\mu$ in the numerators at all, which makes little sense. -- -- -- - - - ---------------
So, I have two arguments and they are giving different results. Also, I would like to understand how to actually complete in a convincing way the reasonings I presented, specially in the first method. 
 A: It is not true that $(\omega \text{Id}-h)^{-1}$ is the causal Green function. Define as above $g(\mu,t-t^\prime)= \frac{-i}{\hbar}\langle T[\hat{c}_\mu(t)\hat{c}_\mu^\dagger(t^{\prime)}]\rangle.$ Suppose $\vert \phi_\lambda \rangle$ form a (orthonormal )basis of eigenstates of the one-electron hamiltonian, $\hat{h}\vert \phi_\lambda \rangle = e_\lambda \vert \phi_\lambda \rangle.$ Define an operator $\hat{\xi}$ such that the elements of the basis $\vert \phi_\lambda \rangle$ are its eigenfunctions, with eigenvalues 1 or -1 if $\vert \phi_\lambda \rangle$ is, respectively, an occupied or an unoccupied orbital (that is, if it appears or not in the Slater determinant forming the Ground state)
Then what holds is that the Fourier transform of this quantity is:
$$g(\mu,\omega)=\langle \mu \vert \hat{R}(\omega) \vert \mu \rangle,$$
where $\hat{R}(\omega)=(\omega \text{Id}-\hat{h}+i\eta\hat{\xi})^{-1},$ where the limit $\eta \rightarrow 0$ is implicit. That this is true can be seen by comparison with the formula obtained with the "many-body approach" in the first method, $\sum_ {m} \dfrac{\langle 0 \vert c_\mu \vert m \rangle \langle m \vert c_\mu^\dagger \vert 0 \rangle}{\omega-(\omega_m-\omega_0)+i\eta}+\sum_{m}\dfrac{\langle m \vert c_\mu \vert 0 \rangle \langle 0 \vert c_\mu^\dagger \vert m \rangle}{\omega-(\omega_0-\omega_m)-i\eta}$ (the "intuitive" reasoning after this formula in the first method is wrong unless $\vert \chi_\mu \rangle$ is itself one of the eigenstates of the Hamiltonian!).  
The function $\hat{R}(\omega)$ can thus be represented as 
$$\hat{R}(\omega)=\sum_{\lambda} \dfrac{\vert \phi_\lambda \rangle \langle \phi_\lambda \vert}{\omega-\omega_\lambda+i\eta\theta(\omega_F-\omega_\lambda)},$$
where $\omega_F$ is the Fermi level and $\theta(x)$ is a function taking the value 1 if $x>0$ and the value $-1$ for $x<0.$
Each element $\vert \phi_\lambda \rangle $ of the basis has an expansion in terms of the our basis orbitals:
$$\vert \phi_\lambda \rangle = \sum_{\mu} c_\mu (\lambda)\vert \chi_\mu \rangle.$$
Then, we have, splitting the occupied and the unoccupied part:
$$g(\mu,\omega)=\sum_{\lambda \ \text{occupied}} \dfrac{\vert c_\mu(\lambda) \vert ^2}{\omega-\omega_\lambda+i\eta}+\sum_{\lambda \ \text{unoccupied}} \dfrac{\vert c_\mu(\lambda) \vert ^2}{\omega-\omega_\lambda-i\eta}.$$
Notice also that $\sum_{\lambda \ \text{occupied}} \vert c_\mu(\lambda )\vert^2 = n_\mu $ and $\sum_{\lambda \ \text{unoccupied}} \vert c_\mu(\lambda )\vert^2 = 1-n_\mu. $  If $\vert \chi_\mu \rangle$ is one of the eigenstates $\vert \phi_\lambda \rangle$ is clear that we recover the well-known result $\dfrac{n_\mu}{\omega-\omega_\mu+i\eta}+\dfrac{1-n_\mu}{\omega-\omega_\mu-i\eta}$
