How to dress free Green functions with constant broadening? I want to find a way to dress free Keldysh Green functions with the simplest level broadening. But there seems to be some quite unexpected result.
Let's consider free Keldysh Green functions in textbooks for a single-level with time-translation symmetry
\begin{equation}\label{eq:simpleGs}
\begin{split}
    G^{\mathrm{R}(\mathrm{A})}(t) = \mp i \theta(\pm t)  e^{- i \varepsilon_0 t} &\longleftrightarrow \frac{1}{\varepsilon - \varepsilon_0 \pm  i  0} \\
    G^\mathrm{K}(t) = - i (1-2n_F)  e^{- i \varepsilon_0 t} &\longleftrightarrow -2\pi i (1-2n_F)\delta(\varepsilon - \varepsilon_0) \\
    G^<(t) =  i  n_F  e^{- i \varepsilon_0 t} &\longleftrightarrow 2\pi i  n_F \delta(\varepsilon - \varepsilon_0),
\end{split}
\end{equation}
where the right hand side shows the energy space expression and $n_F$ is the Fermi distribution valued at $\varepsilon_0$.
My expectation of the simplest possible level broadening is just a constant $\gamma$ as follows
\begin{equation}
\begin{split}
    \tilde{G}^{\mathrm{R}(\mathrm{A})}(t) = \mp i \theta(\pm t)   e^{- i \varepsilon_0 t}  e^{\mp \gamma t} &\longleftrightarrow  \frac{1}{\varepsilon - \varepsilon_0 \pm  i  \gamma} \\
    \tilde{G}^\mathrm{K}(t) = - i   (1-2n_F)  e^{- i \varepsilon_0 t}  e^{-\gamma |t|} &\longleftrightarrow -2 i   (1-2n_F) \frac{\gamma}{(\varepsilon-\varepsilon_0)^2+\gamma^2} \\
    \tilde{G}^<(t) =  i  n_F  e^{- i \varepsilon_0 t}  e^{-\gamma |t|} &\longleftrightarrow 2 i  n_F \frac{\gamma}{(\varepsilon-\varepsilon_0)^2+\gamma^2},
\end{split}\tag{1}\label{dressedGs}
\end{equation}
which are consistent with the general relation
$$G^< = \frac{1}{2}(G^\mathrm{K}-G^\mathrm{R}+G^\mathrm{A}).\tag{2}\label{G^<}$$
$\tilde{G}^{\mathrm{R}(\mathrm{A})}$ looks very natural as we just replace the infinitesimal imaginary number by a finite linewidth $\gamma$. The other two, $\tilde{G}^\mathrm{K}$ and $\tilde{G}^<$, are not really weird if one recalls the approximated $\delta$-function $\delta(\varepsilon-\varepsilon_0)=\lim_{\gamma\rightarrow0}\frac{1}{\pi}\frac{\gamma}{(\varepsilon-\varepsilon_0)^2+\gamma^2}$.
In order to get this, I tried adding the self-energy $\Sigma^<(\varepsilon)= i \gamma$ and hence $\Sigma^<(t) = \frac{i\gamma}{2\pi} \int d\omega e^{-i\omega(t)} = i\gamma\delta(t)$.
With the general relation $\Sigma^{\mathrm{R}(\mathrm{A})} = \mp \Sigma^<$ valid when we only have $\Sigma^<$, we can apply this to the Dyson equation
$$\tilde{G}^{\mathrm{R}(\mathrm{A})}(\varepsilon) = [1 - G^{\mathrm{R}(\mathrm{A})}\Sigma^{\mathrm{R}(\mathrm{A})}]^{-1} G^{\mathrm{R}(\mathrm{A})},$$
which exactly gives our expected $\tilde{G}^{\mathrm{R}(\mathrm{A})}$ in Eq. \eqref{dressedGs}.
Then let's apply the textbook formula [see Eq. (39) in this note]
$$\tilde{G}^< = (1+\tilde{G}^\mathrm{R}\Sigma^\mathrm{R}) G^< (1+\Sigma^\mathrm{A}\tilde{G}^\mathrm{A}) + \tilde{G}^\mathrm{R} \Sigma^< \tilde{G}^\mathrm{A}\tag{3}\label{G^<formula}$$ and Eq. \eqref{G^<}.
To my surprise, it gives the $\tilde{G}^\mathrm{K}$ and $\tilde{G}^<$ expressions in Eq. \eqref{dressedGs} with $n_F=\frac{1}{2}$, which is quite weird, i.e., the $n_F$ dependence is completely lost and actually $\tilde{G}^\mathrm{K}=0$. (BTW, it's straightforward to see this result. Here the whole first part before the last '+' in Eq. \eqref{G^<formula} vanishes as I've checked, which is indeed what happens for steady states as mentioned below Eq. (39) in that note. Then the remaining part $\tilde{G}^\mathrm{R} \Sigma^< \tilde{G}^\mathrm{A}$ obviously is missing $n_F$.)
Am I doing anything wrong here? Or how should one correctly get the simple level broadening to those Green functions?
 A: The simplest case where such a constant broadening arises naturally is a level coupled to a Fermi sea, in a broadband limit (i.e., with the constant density in the Fermi sea), with Hamiltonian:
$$
H=\epsilon_0 d^\dagger d + \sum_k\epsilon_k c_k^\dagger c_k + \sum_k\left(Vd^\dagger c_k + V^*c_k^\dagger d\right),
$$
The broadening in this case is given by $$\gamma = 2\pi|V|^2\rho,$$ where $$\rho(\epsilon) = \sum_k\delta(\epsilon - \epsilon_k)$$ is the density of states, which is assumed constant (the broad-band limit).
The self-energies are then given by
$$
\Sigma^{r,a} = \mp\frac{i\Gamma}{2},\\
\Sigma^<(\omega)= i\Gamma n(\omega),\\
\Sigma^>(\omega)= -i\Gamma \left[1-n(\omega)\right],
$$
which satisfy the basic relation
$$
\Sigma^>-\Sigma^<=\Sigma^r-\Sigma^a.
$$
Caution: I am relying on my memory here, so some factors may be incorrect. The problem is solved in some details in Jauho, Meir, and Wingreen's paper, but I am not necessarily following their notations. Another useful reference is the Rammer&Smith's review... everything else is math :)
Derivations
We separate the Hamilronian into a non-interacting part and the "interaction" (tuneling):
$$
H = H_0 + H_T,\\
H_0=\epsilon_0 d^\dagger d + \sum_k\epsilon_k c_k^\dagger c_k,\\
H_T = \sum_k\left(Vd^\dagger c_k + V^*c_k^\dagger d\right).
$$
The Keldysh Green's function for the level is:
$$
G_{dd}(t,t') = \left\langle T_K\left[Sd(t)d^\dagger(t')\right]\right\rangle,
$$
where $T_K$ is ordering along the Keldysh contour, whereas the scattering matrix is given by
$$
S=T_K\exp\left[-i\int_{C_K}d\tau H_T(\tau)\right],
$$
and all the operators are in the interaction representation (i.e. their evolution is determined by $H_0$).
We can now use Feynmann-Dyson expansion or another method to write down the Feynmann-Dyson equation as:
$$
G_{dd}(t,t') = g_{dd}(t,t') + \int_{C_K}d\tau_1d\tau_2g_{dd}(t,\tau_1)\Sigma(\tau_1,\tau_2)G_{dd}(\tau_2,t'),\\
\Sigma(\tau_1,\tau_2)=|V|^2\sum_k g_k(\tau_1,\tau_2).
$$
The non-interacting Green's functions are given by
$$
g_{dd}(t,t')=\left\langle T_K\left[d(t)d^\dagger(t')\right]\right\rangle,\\
g_{k,k'}(t,t')=\delta_{k,k'}g_k(t,t')=\delta_{k,k'}\left\langle T_K\left[c_k(t)c_k^\dagger(t')\right]\right\rangle.
$$
This allows us to write right away:
$$
\Sigma^{r,a}(\omega)=|V|^2\sum_{k}\frac{1}{\omega -\epsilon_k\pm i\eta}
=\rho|V|^2\int d\epsilon \left[ \mathcal{P}\frac{1}{\omega-\epsilon} \mp\pi\delta(\omega-\epsilon)\right] = \mp\frac{i\Gamma}{2},\\
\Sigma^<(\Omega)=|V|^2\sum_k i2\pi n(\epsilon_k)\delta(\omega-\epsilon_k)
= 2\pi i\rho|V|^2 \int d\epsilon n(\epsilon)\delta(\omega-\epsilon)
=i\Gamma n(\omega),\\
\Sigma^<(\Omega)=-|V|^2\sum_k i2\pi \left[1-n(\epsilon_k)\right]\delta(\omega-\epsilon_k)
= -2\pi i\rho|V|^2 \int d\epsilon \left[1-n(\epsilon)\right]\delta(\omega-\epsilon)
=-i\Gamma \left[1-n(\omega)\right].
$$
Using the Langreth algebra and switching to Fourier space the Dyson equation can be written as:
$$
G_{dd}^{r,a}(\omega)=g_{dd}^{r,a}(\omega) + g_{dd}^{r,a}(\omega)\Sigma^{r,a}(\omega)G_{dd}^{r,a}(\omega),\\
G_{dd}^{>,<}(\omega)=g_{dd}^{>,<}(\omega) + g_{dd}^{>,<}(\omega)\Sigma^{a}(\omega)G_{dd}^{a}(\omega) 
+ g_{dd}^{r}(\omega)\Sigma^{>,<}(\omega)G_{dd}^{a}(\omega)
 + g_{dd}^{r}(\omega)\Sigma^{r}(\omega)G_{dd}^{>,<}(\omega).
$$
The first equation can be recast as
$$
\left\{\left[g_{dd}^{r,a}(\omega)\right]^{-1} - \Sigma^{r,a}(\omega)\right\}G_{dd}^{r,a}(\omega)=\left[G_{dd}^{r,a}(\omega)\right]^{-1}G_{dd}^{r,a}(\omega)=1\\
\Longrightarrow G_{dd}^{r,a}(\omega) = \frac{1}{\left[g_{dd}^{r,a}(\omega)\right]^{-1} - \Sigma^{r,a}(\omega)} = \frac{1}{\omega-\epsilon_d\pm\frac{i\Gamma}{2}}
$$
The equation for the greater and the lesser functions can be written as (taking into account that $\left[g_{dd}^{r}(\omega)\right]^{-1}g_{dd}^{>,<}(\omega)=0$):
$$
\left\{\left[g_{dd}^{r}(\omega)\right]^{-1} - \Sigma^{r}(\omega)\right\}G_{dd}^{>,<}(\omega)
=\left[G_{dd}^{r}(\omega)\right]^{-1}G_{dd}^{>,<}(\omega)=\Sigma^{>,<}(\omega)G_{dd}^{a}(\omega)\\
\Longrightarrow 
G_{dd}^{>,<}(\omega) = G_{dd}^{r}(\omega)\Sigma^{>,<}(\omega)G_{dd}^{a}(\omega) = \frac{\Sigma^{>,<}(\omega)}{(\omega-\epsilon_d)^2+\frac{\Gamma^2}{4}}
$$
