(Thermal) Average of quartic operators of free particles

I am having problems to compute explicitly averages of quartic terms, like $langle a_{mathbf{p}}^dagger a_{mathbf{q}} a_{mathbf{r}}^dagger a_{mathbf{s}}rangle$. The subscripts are linear momenta.

The system has $N$ free particles (let us consider fermions) in a volume $V$, whose Hamiltonian is given by
$$H = sum_mathbf{p} E_{mathbf{p}} a_{mathbf{p}}^dagger a_{mathbf{p}}. $$

As the particles are free, I am considering that the distribution of momentum is $f(mathbf{p})$, such that

$$int d mathbf{p}f(mathbf{p}) = 1.$$

I saw in Statistical Methods for Quantum Optics 1 (chap 1 – Carmichael) that thermal averages of $langle a_{mathbf{p}}^dagger a_{mathbf{q}}rangle = langle a_{mathbf{p}}^dagger a_{mathbf{p}}rangle delta_{mathbf{p},mathbf{q}}$. I understand that for different momenta $langle a_{mathbf{p}}^dagger a_{mathbf{q}}rangle = langle{a_{mathbf{p}}^dagger}rangle langle{a_{mathbf{q}}}rangle = 0$ in this type of context. It does make sense. Even so, I did not find an explicit derivation of this result and consequently the average for quartic terms are more complicated for me.

Finally, in the thermal state the density operator is a product and I can’t see how it works for a general distribution of free particles. How can I proceed?