This problem is from the book Quantum theory of many particles systems by Fetter & Walecka (1971), exercise 1.4.
Problem discription:
The Hamiltonian could be divided into $$H_0=\sum_\limits{\vec{k},\lambda}\frac{\hbar^2 k^2}{2m}a^\dagger_{\vec{k},\lambda}a_{\vec{k},\lambda}$$ and the perturbation $$H_1=\frac{e^2}{2V}\sum_{\vec{p},\vec{k},\vec{q}\neq0}\sum_{\lambda_1,\lambda_2}\frac{4\pi}{q^2}a^\dagger_{\vec{k}+\vec{q},\lambda_1}a^\dagger_{\vec{p}-\vec{q},\lambda_2}a_{\vec{p},\lambda_2}a_{\vec{k},\lambda_1}$$ where $V$ is the volumn of the electron gas (and $N$ is the particle number, see below).
Using the standerd second order perturbation theory, the second order correction is $$E^{(2)}=\sum_{n\neq0}\frac{\langle0|H_1|n\rangle\langle n|H_1|0\rangle}{E_0-E_n}$$ where $|0\rangle$ is the ground state (Fermi sea) and $|n\rangle$ excited state.
The result is given by $$E^{(2)}=\frac{Ne^2}{2a_0}(\epsilon_2^r+\epsilon_2^b)$$ where the direct term is $$\epsilon_2^r=-\frac{3}{8\pi^5}\int\text{d}q\frac{1}{q^4}\int_{|\vec{p}+\vec{q}|>1}\text{d}k^3\int_{|\vec{p}+\vec{q}|>1}\text{d}^3p\frac{\theta(1-k)\theta(1-p)}{q^2+\vec{q}\cdot(\vec{k}+\vec{p})}$$ and the exchange term is $$\epsilon_2^b=\frac{3}{16\pi^5}\int\text{d}q\frac{1}{q^2}\int_{|\vec{p}+\vec{q}|>1}\text{d}k^3\int_{|\vec{p}+\vec{q}|>1}\text{d}^3p\frac{\theta(1-k)\theta(1-p)}{(\vec{q}+\vec{k}+\vec{p})^2[q^2+\vec{q}\cdot(\vec{k}+\vec{p})]}$$
I tried:
I think now that $\langle0|H_1$ is a bra with two electrons excited from the ground state, $|n\rangle$ must be two electrons excited state like $a^\dagger_{\vec{k}'+\vec{q}',\lambda_1'}a^\dagger_{\vec{p}'-\vec{q}',\lambda_2'}a_{\vec{p}',\lambda_2'}a_{\vec{k}',\lambda_1'}|0\rangle$ otherwise the factor $\langle0|H_1|n\rangle$ would be zero. For a fixed $|n\rangle=a^\dagger_{\vec{k}'+\vec{q}',\lambda_1'}a^\dagger_{\vec{p}'-\vec{q}',\lambda_2'}a_{\vec{p}',\lambda_2'}a_{\vec{k}',\lambda_1'}|0\rangle$, the summation $\langle0|\sum'_{q,k,p}\sum_{\lambda_1,\lambda_2}\frac{4\pi}{q^2}a^\dagger_{\vec{k}+\vec{q},\lambda_1}a^\dagger_{\vec{p}-\vec{q},\lambda_2}a_{\vec{p},\lambda_2}a_{\vec{k},\lambda_1}|n\rangle$ only contributes terms in which annihilation and creation operators are "paired". There are two situations, $(\lambda_1\rightarrow\lambda_1',\lambda_2\rightarrow\lambda_2')$ and $(\lambda_1\rightarrow\lambda_2',\lambda_2\rightarrow\lambda_1')$. In both situation the denominator is $1/q'^2$. The factor $\langle n|H_1|0\rangle$ is similar. Lastly consider $1/(E_0-E_n)$ and the summation is over $(p',k',q')$. However, I could not obtain $\epsilon_2^b$ at all.
PS. I found this result is given by Gell-mann in 1957 : Phys. Rev. 106, 364 (1957), but only results are listed without derivation. I suspect that this problem is not difficult, but I couldn't solve it.