Suppose I have a many-body system with creation/annihilation operators $\hat{c}^\dagger_n$, $\hat{c}_n$, and Hamiltonian: $$\hat{H}=\sum_n E_n \hat{c}^\dagger_n\hat{c}_n$$ If I wanted to write down the density matrix for a thermal state in the grand canonical ensemble with inverse temperature $\beta$ and chemical potential $\mu$, I could simply use: $$ \hat{\rho}_{\text{GC}}=\frac{\exp(-\beta(\hat{H}-\mu\sum_n \hat{c}^\dagger_n\hat{c}_n))}{\text{Tr}\left\{\exp(-\beta(\hat{H}-\mu\sum_n \hat{c}^\dagger_n\hat{c}_n))\right\}} $$ where the mean particle number is fixed by $\mu$, but there will of course be some number variance.
I am a little confused on how one would do this in the case of the canonical ensemble, for a many-body system. According to Wikipedia, the canonical ensemble corresponds to the density matrix: $$\hat{\rho}_{\text{C}}=\exp(\beta(F-\hat{H}))$$ where $F$ is chosen to normalise the density matrix, requiring $e^{-\beta F}=\text{Tr}\left\{e^{-\beta \hat{H}}\right\}$. Is this for a one-body system? If I were to choose some arbitrary total particle number $N$, in the case of the grand canonical ensemble I could just tune $\mu$ to fix the average to it. But there doesn't seem to be any parameter to control particle number in the given expression for the canonical ensemble. It isn't clear to me how I would obtain the density matrix for a canonical ensemble with $N$ particles, and the density matrix ought to be different for different $N$. Is there a generalisation for many-body systems that I don't know about, or am I missing something subtle?