Calculating the ground state energy in an electron gas (Jellium) $\renewcommand{\phi}{\varphi}
\newcommand{\vec}[1]{\textbf{#1}}$
This is pretty lengthy because I kept having ideas while writing this.

Consider an electron gas in a cube of side length $L$, $V=L^3$, with periodic boundary conditions. 
The Hamiltonian is given by $\hat{H} = \hat{H}_\text{el} + \hat{H}_\text{b} + \hat{H}_\text{el$-$b}$, and we make a few approximations and limits to arrive at its final form (write $\lvert \vec{k}\rvert \equiv k)$,
\begin{equation*}
\hat{H} = \frac{e^2}{a_0 r_s^2}
\left(\sum_{\vec{k},\alpha}\frac{k^2}{2}\hat{a}_{\vec{k},\alpha}^\dagger\hat{a}_{\vec{k},\alpha}^{\phantom{\dagger}}+\frac{r_s}{2V}\sum_{{\underset{\vec{q}\neq0}{\vec{k,p,q}}}}\sum_{\alpha_1,\alpha_2}\frac{4\pi}{q}
\hat{a}_{\vec{k+q},\alpha_1}^\dagger
\hat{a}_{\vec{p-q},\alpha_2}^\dagger
\hat{a}_{\vec{k},\alpha_1}^{\phantom{\dagger}}
\hat{a}_{\vec{p},\alpha_2}^{\phantom{\dagger}}
\right),
\end{equation*}
with fermionic creation operators $\hat{a}^\dagger$. Here, 
$$
r_0\equiv\sqrt[3]{\frac{3V}{4\pi N}},\quad a_0 \equiv \frac{\hbar^2}{me^2},\quad r_s\equiv\frac{r_0}{a_0}.
$$
Call the first term $\mathcal{A}$, the second one $\mathcal{B}$. We see that $\mathcal{B}$ is a perturbation in the high-density limit $r_s\rightarrow0$. Thus, the ground state energy of this can be written as $E = E^{(0)}+E^{(1)}+\ldots$ where $E^{(0)}=\langle F|\mathcal{A}|F\rangle$ for the ground state $\vert F\rangle$ of $\mathcal{A}$, which is the Fermi sea, momenta filled up to $p_F = \hbar k_F$.
We wish to compute both $\langle F|\mathcal{A}|F\rangle$ and $E^{(1)}=\langle F|\mathcal{B}|F\rangle$. First, we need to determine $k_F$. We do this by looking at the expectation value of the number operator (in out limit, we replace sums over $\vec{k}$ by integrals times a factor).
And now my confusion begins. In the following calculation, I have sort of a handwavy argument for going from $(\star)$ to $(\star\star)$, but I'm not a hundred percent sure.
\begin{align*}
N = \langle F|\hat{N}|F\rangle&=\frac{V}{(2\pi)^3}\sum_\alpha\int\mathrm{d}^3k\langle F|\hat{a}_{\vec{k},\alpha}^\dagger \hat{a}_{\vec{k},\alpha}^{\phantom{\dagger}} |F\rangle\\[1em]
&=\frac{V}{4\pi^3}\int\mathrm{d}^3k\langle F|\hat{a}_{\vec{k}}^\dagger \hat{a}_{\vec{k}}^{\phantom{\dagger}} |F\rangle \tag{$\star$}\\[1em]
&=\frac{V}{4\pi^3}\int\mathrm{d}^3k \,\Theta(k_F - k)\tag{$\star\star$} \\[1em]
&=\frac{Vk_F^3}{3\pi^2}
\end{align*}
Basically, we're saying we have the Fermi sea, so there are no particles of momenta above $k_F$, i.e. the integrand should vanish (we'd then say the integral is zero, because constants are either 0 or 1 in physics). So we need the step function accounting for that. 
But I actually don't understand what $\hat{a}_\vec{k}\vert F\rangle$ is. I don't see what happens when we annihilate this particular ground state. Well, it can't be zero, can it?
And now, while writing this, I realized that we must have
$$
\vert F\rangle = \sum_\alpha \int \text{d}^3k\, \Theta(k_F - k) \hat{a}_{\vec{k},\alpha}^\dagger |0\rangle,
$$
right? Because then we can easily calculate
\begin{align*}
\hat{a}_{\vec{k},\alpha}|F\rangle = \Theta(k_F-k)\vert 0\rangle\implies\langle F|\hat{a}_{\vec{k},\alpha}^\dagger \hat{a}_{\vec{k},\alpha}^{\phantom{\dagger}} |F\rangle=\lVert \hat{a}_{\vec{k},\alpha}|F\rangle\rVert^2=\Theta(k_F-k),
\end{align*}
and thus the transition from $(\star)$ to $(\star\star)$ is explained.
But this wasn't everything by far. As mentioned above, we want to compute the expectation values of $\mathcal{A}$ and $\mathcal{B}$ in the Fermi sea, using the fact that $k_F\approx 1.92\cdot r_0^{-1}$. The results should be
\begin{align*}
\tag{$R\mathcal{A}$}\frac{2.21}{2} \frac{N e^2}{a_0 r_s}\\[1em]
\tag{$R\mathcal{B}$}-\frac{0.916}{2}\frac{N e^2}{a_0 r_s}.
\end{align*}
My problem is now that I can't even derive the first one, $(R\mathcal{A})$. This is what I did:
\begin{align*}
E^{(0)} &=\frac{V}{8 \pi^3}\frac{e^2}{a_0 r_s^2}\frac{1}{2}\sum_\alpha\int \text{d}^3 k\, k^2 \Theta(k_F-k)\\[1em]
&= \frac{V}{2 \pi^2}\frac{e^2}{a_0 r_s^2}\int_0^{k_F} \text{d}k\,k^4\\[1em]
&= \frac{V}{2 \pi^2}\frac{e^2}{a_0 r_s^2}\frac{k_F^5}{5}\\[1em]
&= \frac{1.92^5}{10 \pi^2}\frac{V e^2}{a_0 r_s^2}r_0^{-5} = \frac{1.92^5}{10 \pi^2}\frac{V e^2}{a_0 r_s^2r_0^2}\frac{4\pi N}{3V}\\[1em]
&= \frac{2\cdot 1.92^5}{15 \pi}\frac{N e^2}{a_0 r_s^2r_0^2}\\[1em]
&\approx \frac{2.21}{2}\frac{N e^2}{a_0 r_s^2r_0^2}.
\end{align*}
This is obviously not the same as $(R\mathcal{A})$, since $r_s r_0^2 = {}^{r_0^3}/_{a_0}\neq 1$.
And if that's wrong I don't even dare try my hand at $(R\mathcal{B})$, so:
What did I do wrong/not understand correctly?
Thanks!
 A: I can see some mistakes :

"I realized that we must have $|F\rangle=\sum_\alpha\int\;\mathrm{d}\textbf{k}\,\Theta(k_F-k)\,\hat{a}^\dagger_{\textbf{k},\alpha}|0\rangle$"

This is wrong : the ground state of an non-interacting fermion gas (Fermi sea) is not some sort of superpositions on the individual $\textbf{k}$ states as you seem to suggest. It has to be a real many-body state because of the Pauli exclusion principle playing on the spin state $\alpha=\uparrow$ or $\downarrow$ :
$$
|F\rangle=\prod_{k<k_F}\hat{a}^\dagger_{\textbf{k},\uparrow}\hat{a}^\dagger_{\textbf{k},\downarrow}|0\rangle
$$

Calculating the density $N/V$

Now we can compute properly the mean value of the particle number $\hat{N}$ in the ground state :
$$
N=\langle F|\hat{N}|F\rangle=\sum_{\textbf{k},\alpha}\langle F|\hat{a}^\dagger_{\textbf{k},\alpha}\hat{a}_{\textbf{k},\alpha}|F\rangle
$$
The operator $\hat{a}^\dagger_{\textbf{k},\alpha}\hat{a}_{\textbf{k},\alpha}$ counts the number of electrons occupying the state $\{\textbf{k},\alpha\}$. By definition of the ground state $|F\rangle$, all the $\{\textbf{k},\alpha\}$ states are occupied by one electron for $|\textbf{k}|<k_F$. Therefor :
\begin{eqnarray*}
\hat{a}^\dagger_{\textbf{k},\alpha}\hat{a}_{\textbf{k},\alpha}|F\rangle&=&1\times|F\rangle\quad\text{if}\;|\textbf{k}|<k_F\\
&=&0\times|F\rangle\quad\text{otherwise}
\end{eqnarray*}
which can contracted as :
$$
\hat{a}^\dagger_{\textbf{k},\alpha}\hat{a}_{\textbf{k},\alpha}|F\rangle=\Theta(k_F-|\textbf{k}|)\,|F\rangle
$$
Frequently we use the following normalisation :
$$
\sum_\textbf{k}\rightarrow\frac{V}{(2\pi)^3}\int\mathrm{d}\textbf{k}
$$
Thus, we get :
$$
N=2\times\frac{V}{(2\pi)^3}\int\mathrm{d}\textbf{k}\;\Theta(k_F-|\textbf{k}|)\langle F|F\rangle
$$
Using $\langle F|F\rangle=1$ and $\int\mathrm{d}\textbf{k}\;\Theta(k_F-|\textbf{k}|)=4\pi\,k^3_F/3$ by expressing the $\textbf{k}$ variable in the spherical coordinates, we finally get :
$$
\frac{N}{V}=\frac{k^3_F}{3\pi^2}
$$

Calulating the kinetic energy term

I will show you how to get the first term and let the second term as an exercise for you :-) .
Forgetting the factors, what we have left to compute is :
$$
\langle F|\sum_{\textbf{k},\alpha}\,k^2\,\hat{a}^\dagger_{\textbf{k},\alpha}\hat{a}_{\textbf{k},\alpha}|F\rangle=2\times\frac{V}{(2\pi)^3}\int\mathrm{d}\textbf{k}\,k^2\,\Theta(k_F-|\textbf{k}|)
$$
Similary to the previous calculation, we get the $\textbf{k}$ variable in the spherical coordinates $(k,\theta,\phi)$, so that :
$$
\int\mathrm{d}\textbf{k}\,k^2\,\Theta(k_F-|\textbf{k}|)=\int_0^{k_F}\mathrm{d}k\,k^4\;\times\;\int_0^{2\pi}\mathrm{d}\phi\;\times\;\int_0^{\pi}\mathrm{d}\theta\,\sin\theta=\frac{4\pi\,k_F^5}{5}
$$
which should give you the right answer :
$$
\langle F|\sum_{\textbf{k},\alpha}\,\frac{\hbar^2 k^2}{2m}\,\hat{a}^\dagger_{\textbf{k},\alpha}\hat{a}_{\textbf{k},\alpha}|F\rangle=\frac{3}{5}N E_F
$$
