I'm currently doing some calculations which require evaluating various standard thermal expectation values in the canonical ensemble (both bosons and fermions). Now, in order to make my theoretical machinations easier, I am actually using the grand canonical ensemble, where the chemical potential acts as a Lagrange multiplier enforcing the constraint $\langle \hat{N} \rangle = N$, where $N/V$ is the fixed density of the physical system. The justification for this is that the relative fluctuations in $\langle \hat{N} \rangle$ should vanish in the thermodynamic limit, in which case I expect fixing the average number to be physically equivalent to fixing the number once and for all. (Also, this approach seems to be adopted by a several presumably trustworthy references, see for example Simons & Altland Section 6.3.) This intuition seems reasonable, but I wonder if matters may be more subtle than this argument implies.
Do thermal averages in the thermodynamic limit of the grand canonical and canonical ensembles coincide?
I'm hoping for either a more rigorous justification supporting this procedure, or examples where it can go horribly wrong. Pointers to appropriate references would also be much appreciated.