Matsubara Green function of anderson impurity model

I am currently having trouble computing the imaginary-time Green's functions of a model similar to the single-impurity anderson model.

The hamiltonian is given as:

\begin{equation} H = \Sigma_{\vec{k} } \epsilon_k c^{\dagger}_{\vec{k}} c_{\vec{k}} + V \Sigma_{\vec{k}} ( c^{\dagger}_{\vec{k}} d+c_{\vec{k}} d^{\dagger} ) + \epsilon_d d^\dagger d \end{equation}

where $c^{(\dagger)}_{\vec{k}}$ are the creation/annihilation operators for the fermionic bath (in momentum space) and $d^{(\dagger)}$ are the operators for the fermionic impurity and $V\in\mathbb{R}$. The model is taken from Elementary Introduction to the Hubbard Model (Exercise 45).

The task is to determine the Matsubara GF $G_d(i\omega_n)$ of the impurity.

Starting from the time-ordered($\mathcal{T}$) representation of the Green's function $G(t) = -\langle \mathcal{T}[ a(t) a^\dagger(0) ]\rangle$ (for now omitting the indices $k,d$) and differentiating w.r.t $t$ one obtains:

\begin{equation} \partial_t G(t) = -\delta(t)\langle \lbrace a(t),a^\dagger(0)\rbrace \rangle -\theta(t)\langle [H,a(t)]a^\dagger(0)\rangle + \theta(-t)\langle a^\dagger(0) [H,a(t)]\rangle \end{equation} where $a\in\lbrace c_{\vec{k}},d\rbrace$ for brevity. If one now takes the FT in Matsubara formalism as $G(i\omega_n) = \frac{1}{\beta}\int_0^\beta G(t)e^{i\omega_n t}dt$ with $\beta=(k_bT)^{-1}$ it follows (using $G(t+\beta) = -G(t)$ tand integration by parts) \begin{align*} \frac{-i\omega_n}{\beta}\int_0^\beta G(t)e^{i\omega_n t}dt &= -\frac{1}{\beta}\int_0^\beta \delta(t)\langle \lbrace a(t),a^\dagger(0)\rbrace\rangle e^{i\omega_n t}dt\\ &- \frac{1}{\beta}\int_0^\beta \theta(t)\langle [H,a(t)]a^\dagger(0)\rangle e^{i\omega_n t}dt\\ &+ \frac{1}{\beta}\int_0^\beta \theta(-t)\langle [H,a(t)]a^\dagger(0)\rangle e^{i\omega_n t}dt\\ &=\frac{-1}{\beta} - \frac{1}{\beta}\int_0^\beta \langle [H,a(t)]a^\dagger(0)\rangle e^{i\omega_n t}dt \end{align*} where we loose the $\theta(-t)$ term due to integration over $[0,\beta]$.

It is stated in the exercise that one should follow the above procedure for $G_d(t)$ and $G_{\vec{k}}(t)$ to obtain two equations in two unknowns and then solve the system to obtain: \begin{equation} G_d(i\omega_n) = - \left( i\omega_n - \epsilon_d -V^2\sum_{\vec{k}} \frac{1}{i\omega_n - \epsilon_{\vec{k}}} \right)^{-1} \end{equation}

The last term in the bracket looks a lot like the sum over the Green functions of the separate modes $G_{\vec{k}}(i\omega_n)=\frac{1}{i\omega_n - \epsilon_{\vec{k}} }$ which are obtained from $H = \sum_{\vec{k} } \epsilon_{ \vec{k} }c^{\dagger}_{ \vec{k}} c_{\vec{k}}$.

My problem is, that I do not see how to perform the last integral on the RHS of the FT. As far as I have been able to calculate the commutator of $H$ with either $d(t)$ or $c_{\vec{k}}(t)$ will always result in an additional term with the other operator, for example: \begin{equation} \frac{1}{\beta}\int_0^\beta \langle [H,d(t)]d^\dagger(0)\rangle e^{i\omega_n t}dt = \frac{\epsilon_d}{\beta}\int_0^\beta \langle d(t)d^{\dagger}(0)\rangle e^{i\omega_n t}dt + \frac{U}{\beta}\sum_{\vec{k}}\int_0^\beta \langle c_{\vec{k}}(t)d^{\dagger}(0)\rangle e^{i\omega_n t}dt \end{equation} If this is correct, then I am stuck in evaluating the expectation values which do not seem to want to yield to me.

My question thus:

Is there a reference where I could find a comprehensible path from the hamiltonian to the Green function of the impurity $G_d$? It would be even better if somebody could point me to a way how to evaluate the aforementioned FT integral (I think i am missing something somewhere).