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 $$$$\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}$$$$ 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{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}$$$$ 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^ 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?

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}}$$

• Thank you for your answer. Is there an $\epsilon_k$ missing in the second term of $H$? It seems to give the desired result. But could you please elaborate on how to get these self-energies from $H$? Commented Jan 6, 2021 at 12:53
• Indeed, I have missed $\epsilon_k$. I will try to sketch the derivation when I have a bit more time. Commented Jan 6, 2021 at 14:43
• I have added the derivations. Commented Jan 8, 2021 at 10:26
• This is really a nice answer! Only one question left. Your last expression $\frac{\Sigma^{>,<}(\omega)}{(\omega-\epsilon_d)^2+\frac{\Gamma^2}{4}}$ has a complex form transformed to time. To get what I wrote in Eq.(1), one needs $\frac{\Sigma^{>,<}(\epsilon_d)}{(\omega-\epsilon_d)^2+\frac{\Gamma^2}{4}}$ instead. Any way to justify this? Commented Jan 9, 2021 at 0:33
• OK, we are on the same page and I don't object to any of the math. So you think the guess is actually wrong for $G^<$ and we just don't have a simple enough real-time form? I presume adding broadening like this is computationally convenient common practice to avoid singular $\delta$-functions or so. This is a bit surprising to me. Commented Jan 11, 2021 at 12:52