I have figured this out myself... The crucial point is that the non-equal time commutation relation $[a(t), a(t')]=0$ for free Bosonic fields (the same for $a^\dagger$).
(In the context of high energy physics, people usually use field operators $\phi(x)$ which involves with both $a_k$ and $a_k^\dagger$. This is because the dispersion relation is given by $E^2=p^2+m^2$. However, in our case, we will consider the non-relativistic limit $E=\frac{p^2}{2m}$ thus the field operator is just $\psi(x)=a(x)$ which is only related to $a_k$ but not $a_k^\dagger$.)
Now comes to the Wick's theorem. Things are perfectly nice when we are dealing with the vacuum state. However, here we are considering a ground state with $N$ particles. Of course, one may pull them out as $\vert\Phi_0\rangle=\frac{(a_0^\dagger)^N}{\sqrt{N!}}\vert vac\rangle$, but that is too many operator to contract! Let's use Wick Theorem here:
$$
\begin{eqnarray}
\hat{T}[a^\dagger(x,t)a(x,t)a^\dagger(0,0)a(0,0)]&=&:a^\dagger(x,t)a(x,t)a^\dagger(0,0)a(0,0):\\
&+&:a^\dagger(x,t)a(0,0):\langle vac\vert \hat{T}a(x,t)a^\dagger(0,0)\vert vac\rangle\\
&+&:a(x,t)a^\dagger(0,0):\langle vac\vert \hat{T}a(0,0)a^\dagger(x,t)\vert vac\rangle\\
&+&\langle vac \vert \hat{T}a^\dagger(x,t)a(0,0)\vert vac\rangle\langle vac\vert\hat{T}a(x,t)a^\dagger(0,0)\vert vac\rangle.
\end{eqnarray}
$$
(Terms containing $\langle vac\vert\hat{T}a^\dagger(x,t)a(x,t)\vert vac\rangle$ and $\langle vac\vert\hat{T}a^\dagger(0,0)a(0,0)\vert vac\rangle$ vanishes).
The last term containing $\theta(t)\theta(-t)$ so that it vanishes as well. So we are left with three terms.
The second and the third are easy to compute. I won't do it here (it is the first question under this problem in the book). What was confusing me is the first term:
How do I normal ordering it?
Should I first evolve $a(0,0)$ to $a(0,t)$ to make them at the same time then do the normal ordering? This will be too much compicated. Or should I just put the $a^\dagger$'s to the left of $a$'s ? If so, shall it be $a^\dagger(x,t)a^\dagger(0,0)a(x,t)a(0,0)$ or $a^\dagger(0,0)a^\dagger(x,t)a(x,t)a(0,0)$ or $a^\dagger(x,t)a^\dagger(0,0)a(0,0)a(x,t)$, or ...? The point is, I didn't realise it that $[a(x,t),a(0,0)]=0$ , since the Bosons are free, so that the order doesn't actually matter.
It is easy to see in $k$-space,
$$
[H,a_k]=-\frac{k^2}{2m}a_k\Longrightarrow e^{iHt}a_k e^{-iHt}=\exp(-i\frac{k^2}{2m}t)a_k,
$$
that the evolution of free Bosonic operator only gives rise to a phase factor so that
$$
[a_k(t),a(0)]=0
$$
for all $t$. (That is, non-equal-time commutation relation is trivial as the equal-time one. This is only true if the field is free.) And then performing a Fourier transformation shows that $[a(x,t),a(0,0)]=0$ as well. I myself think this is not that trivial but few books talk about it. They just ignore all the time label when doing normal-ordering. A good reference about this is David Tong's lecture notes on QFT, Page 37. Moreover, non-equal-time commutation relation for $a$ and $a^\dagger$ is just a $c$-number rather than an operator for a free field. (We don't need this point here)
I think this means that the Wich's theorem is only a pertubative tool, which is not exact if interaction is turned on.
Now evaluate the operators in $\vert \Phi_0\rangle$, yields
$$
\begin{eqnarray}
\langle\Phi_0\vert\hat{T}\rho(x,t)\rho(0,0)\vert \Phi_0\rangle&=&\langle\Phi_0\vert (a(x,t)a(0,0))^\dagger (a(x,t)a(0,0))\vert \Phi_0\rangle\\
&+&\langle \Phi_0\vert a^\dagger(x,t)a(0,0)\vert \Phi_0\rangle\langle vac\vert\hat{T}a(x,t)a^\dagger(0,0)\vert vac\rangle\\
&+&\langle \Phi_0\vert a(x,t)a^\dagger(0,0)\vert \Phi_0\rangle\langle vac\vert\hat{T}a^\dagger(x,t)a(0,0)\vert vac\rangle.
\end{eqnarray}
$$
The first term can be evaluated by some trick: $(a(x)a(0))^\dagger$ can be somehow regarded as a pair creation operator so that $(aa)^\dagger(aa)$ is a pair-number operator which count how many pairs are in the state. Now that we have got $N$ Bosons in $\vert \Phi_0\rangle$, for each Boson there are $N-1$ to pair thus the number of pairs is $N(N-1)$ and the pair density is given by $\frac{N(N-1)}{V^2}$.
And then we end up with the same result:
$$
\begin{eqnarray}
\langle\Phi_0\vert\hat{T}\rho(x,t)\rho(0,0)\vert \Phi_0\rangle&=&
\frac{N(N-1)}{V^2}\\
&+&\frac{N}{V}\int\frac{\mathrm{d}k}{2\pi}[\theta(t)\exp(ikx-it\frac{k^2}{2m})\\
&+&\theta(-t)\exp(-ikx+it\frac{k^2}{2m})].
\end{eqnarray}
$$