Canonical ensemble in many-body quantum mechanics

First of all, speaking about thermodynamical ensembles(microcanonical, canonical, grand canonical) only makes sense when speaking about many-body systems, since thermodynamics only makes sense for many-body-systems.

The density matrix for the canonical ensemble from Wikipedia is written in a confusing form. Since one needs $text{Tr}rho=1$, one can write $exp (beta F) = frac{1}{text{Tr}exp (-beta hat{H})}$ and $$hat{rho}_C = frac{exp (-beta hat{H})}{text{Tr}exp (-beta hat{H})},$$ in correspondence with the density matrix for the grand canonical ensemble.

The thermodyncamical ensembles(meaning their density operator or partition sum) can be derived from the principle of maximal entropy using the method of Lagrange multipliers. The difference between the ensembles is which quantities are allowed to vary: For the microcanonical ensemble, energy $E_i$ and particle number $N_i$ of every particle in the multi-particle-system is fixed, for the canonical ensemble, one allows the energy $E_i$ of individual particles to vary with fixed total energy $E$ and fixed $N_i$, and for the grand canonical ensemble one allows $E_i$ and $N_i$ to vary with fixed $E$ and $N$.

This corresponds to a closed system(no heat- and particle flux) being perfectly described by the microcanonical ensemble, a system that can exchange heat but no particles for the canonical ensemble and an open system(heat- and particle flux) for the grand canonical system. Since one mostly describes gases in thermodynamics and gases as open systems are naturally described by the grand canonical ensemble, one mostly uses this ensemble. Of course, one can also describe a system with the "wrong" ensemble, but this is unnatural and leads to clumsy calculations. One also has this correspondence in thermodynamics when one has to choose a suitable thermodynamical potential: $U$ for microcanonical ensemble, $F$ for canonical ensemble, $Omega$ for grand canonical ensemble.

With this preparations, I can answer your question. The canonical ensemble is used naturally to describe systems without particle flux, so one can just label the states(particles) and does not need to use the occupation number representation(the $N_n=c_n^dagger c_n$ stuff that one uses when particle flux is allowed in the grand canonical ensemble). If one likes the occupation number representation, one can also use it in the canonical ensemble. Since the canonical ensemble is constructed to not have changes in individual particle numbers, the $exp (-mu sum_n N_n )$ terms would just drop out in the density matrix.

The density matrix only depends on the total particle number $N$ through the trace since the trace sums over matrix elements for all($N$) states.