How to calculate the second-order pertubation in an electron gas? 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. 
 A: From $\langle \Phi_0|H_1|m\rangle\langle m|H_1|\Phi_0\rangle$, consider $\langle m|H_1|\Phi_0\rangle$. We must have $p<k_F, k<k_F$. Further more, we could also obtain $|\vec{k}+\vec{q}|>k_F, |\vec{p}-\vec{q}|>k_F$, otherwise $H_1$ only exchange two particles under Fermi sea keeping ground state, which means $|m\rangle=|\Phi_0\rangle$. So step function $\theta(k_F-k)\theta(k_F-p)\theta(|\vec{k}+\vec{q}|-k_F)\theta(|\vec{p}-\vec{q}|-k_F)$ is obvious now.
The total contribution to the second energy correction is
$$
\frac{(4\pi)^2e^4}{(2V)^2}\frac{2m}{\hbar^2}\sum_{p,k,q,\lambda_1,\lambda_2}\frac{1}{q^2}\frac{\theta(k_F-k)\theta(k_F-p)\theta(|\vec{k}+\vec{q}|-k_F)\theta(|\vec{p}-\vec{q}|-k_F)}{2q^2+2\vec{q}\cdot(\vec{k}-\vec{p})}\\
\times\sum_{p',k',q',\lambda'_1,\lambda'_2}\frac{1}{q'^2}\langle\Phi_0|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'}|(\vec{k}+\vec{q})\lambda_1,(\vec{p}-\vec{q})\lambda_2;p\lambda_2,k\lambda_1\rangle
$$
where the ket above is $|m\rangle$, on the left of semicolon is the particles, and right the hole. The order in the ket is important, and the annihilation or creation operators can only operate one by one. That's say if I want the $a_{\vec{k}',\lambda_1'}$ to annihilate the particle $(\vec{p}-\vec{q})\lambda_2$ skipping over the particle $(\vec{k}+\vec{q})\lambda_1$, there is a factor $-1$.
Now consider the last factor. Obviously, we at least have one scheme to restore ground state, $\vec{k}'=\vec{k}+\vec{q}$, $\vec{p}'=\vec{p}-\vec{q}$, $\vec{p}'-\vec{q}'=\vec{p}$, $\vec{k}'+\vec{q}'=\vec{k}$, and from which one can see $\vec{q}'=-\vec{q}$. Besides, in this pairing scheme, the spin is arbitrary, $(\lambda_1,\lambda_2)=(++),(--),(+-),(-+)$, but when $(\lambda_1,\lambda_2)$ is selected, $(\lambda'_1,\lambda'_2)$ is also fixed. So the summation of spin index is now removed and there is a factor $4$. This contribution is (direct term)
$$
\frac{(4\pi)^2e^4m}{\hbar^2}V\int\frac{\text{d}^3(\vec{p},\vec{q},\vec{k})}{(2\pi)^9}\frac{1}{q^4}\frac{\theta(k_F-k)\theta(k_F-p)\theta(|\vec{k}+\vec{q}|-k_F)\theta(|\vec{p}-\vec{q}|-k_F)}{q^2+\vec{q}\cdot(\vec{k}-\vec{p})}
$$
We also have another pairing scheme, $\vec{k}'=\vec{k}+\vec{q}$, $\vec{p}'=\vec{p}-\vec{q}$, $\vec{p}'-\vec{q}'=\vec{k}$, $\vec{k}'+\vec{q}'=\vec{p}$, which gives $\vec{q}'=\vec{p}-\vec{k}-\vec{q}$. However, due to the exclusion principle, the spin $(\lambda_1,\lambda_2)$ could only be $(++),(--)$, otherwise the same state will be create twice which leads to zero. Besides, one operator "skips", there is a factor $-1$. So this contribution is (exchange term)
$$
-\frac{(4\pi)^2e^4m}{2\hbar^2}V\int\frac{\text{d}^3(\vec{p},\vec{q},\vec{k})}{(2\pi)^9}\frac{1}{q^2}\frac{1}{(\vec{p}-\vec{k}-\vec{q})^2}\frac{\theta(k_F-k)\theta(k_F-p)\theta(|\vec{k}+\vec{q}|-k_F)\theta(|\vec{p}-\vec{q}|-k_F)}{q^2+\vec{q}\cdot(\vec{k}-\vec{p})}
$$
In fact, there are other two pairing schemes, for direct part is $\vec{p}'=\vec{k}+\vec{q}$, $\vec{k}'=\vec{p}-\vec{q}$, $\vec{p}'-\vec{q}'=\vec{k}$, $\vec{k}'+\vec{q}'=\vec{p}$. This contributes the same [with a factor $(-1)^2$]. As a result, the final total contribution should $\times 2$.
Now adopt the variables substitution $\vec{k}\rightarrow-k_F\vec{k}$, $\vec{q}\rightarrow-k_F\vec{q}$, $\vec{p}\rightarrow+k_F\vec{p}$, and using $k_F^3=\frac{3\pi^2N}{V}$ and $a_0=\frac{\hbar^2}{me^2}$. The final result is consistent with the equations given by the book (see question description).
There may be some spelling\grammar mistakes, sorry.
