Solution (source code)

= Solution

At $U=0$, use fermionic <Matsubara frequencies> $\omega_r=(2r+1)\pi/\beta$. The <fermionic Gaussian integral> has inverse propagator $-i\omega_r+\xi_{\mathbf k}$. With its thermal normalization, the determinant for one single-particle level is $1+e^{-\beta\xi}$, as also follows from its empty and occupied states. Two spin species therefore give
$$
\mathcal Z_0=\prod_{\mathbf k}[1+e^{-\beta(\epsilon_{\mathbf k}-\mu)}]^2,
\qquad
\boxed{\Omega_0=-\frac1\beta\ln\mathcal Z_0=-2k_BT\sum_{\mathbf k}\ln[1+e^{-\beta(\epsilon_{\mathbf k}-\mu)}].}
$$
The expression labeled <Helmholtz free energy> in the PDF is the <grand potential>, since the trace contains $H-\mu\hat N$. The thermodynamic <Helmholtz free energy> at a prescribed electron number $N_e$ is \b[$F=\Omega_0+\mu N_e$], with $\mu$ fixed by $N_e=2\sum_{\mathbf k}n_F(\epsilon_{\mathbf k})$. The Matsubara determinant requires a normalization prescription; a raw infinite sum of logarithms is not an ordinary convergent sum.