Solution (source code)

= Solution

Expand the <Grassmann fields> in fermionic <Matsubara frequencies> $\omega_l=(2l+1)\pi/\beta$. For $U=0$, the <Grassmann Gaussian integral> factorizes into thermal <fermionic one-level thermal determinants>:
$$
Z_0=\prod_{k,\sigma}(1+e^{-\beta\xi_k})=\prod_k(1+e^{-\beta\xi_k})^2.
$$
Each factor counts the empty and occupied states of one spin orbital. Equivalently, the regularized <Matsubara sum> gives the logarithm of that factor; the raw infinite determinant requires <thermal normalization of a fermion determinant>. Therefore
$$
\boxed{\Omega_0=-\frac1\beta\ln Z_0=-\frac2\beta\sum_k\ln(1+e^{-\beta(\epsilon_k-\mu)}).}
$$
\b[The printed $F$ is the grand potential $\Omega$], because the trace contains $H-\mu\hat N$. The <Helmholtz free energy> at fixed mean <particle number> is $F_{\rm H}=\Omega+\mu\langle\hat N\rangle$.