Canonical / Grand-Canonical average annihilation operator Does anyone knows a simple way to understand why the average value of the creation (or annihilation) operator should be equal to zero in the Canonical Ensemble? Why instead if I'm dealing with a Grand-Canonical ensemble the same averages can be different from zero? 
I'm asking this because in the Bogoliubov approximation (in the case of weakly interacting Bose gas) basically we set
$a \sim \sqrt(N_0)$
so as far as I understood the average value of a should be different from 0 and the computation should be done in the Grand-Canonical ensemble.
 A: I'm unconvinced that the ensemble average value of an annihilation operator doesn't vanish in the grand canonical ensemble.  Here's my argument for why it does vanish:
Let $H$ be the Hamiltonian and $\mathscr N = \sum_{i=0}^\infty a_i^\dagger a_i$ be the number operator.  Let
$$
  |\mathbf n \rangle = |n_0, n_1, \dots\rangle
$$
be the occupation number basis which satisfies
\begin{align}
  H|\mathbf n \rangle &= E_{\mathbf n}|\mathbf n\rangle,\qquad  E_\mathbf n =     \sum_{i=0}^\infty \epsilon_in_i \\
  \mathscr N|\mathbf n \rangle &= N_{\mathbf n}|\mathbf n\rangle,\qquad  N_\mathbf n = \sum_{i=0}^\infty n_i
\end{align}
and
\begin{align}
  a_i|\mathbf n\rangle &= \sqrt{n_i}|n_0, \dots, n_{i-1}, n_i-1,n_{i+1}, \dots\rangle \\
a^\dagger_i|\mathbf n\rangle &= \sqrt{n_i+1}|n_0, \dots, n_{i-1}, n_i+1,n_{i+1}, \dots\rangle
\end{align}
The ensemble average of any annihilation operator $a_i$ is
\begin{align}
  \frac{1}{Z}\mathrm{tr}(e^{-\beta(H-\mu \mathscr N)} a_i) 
&= \frac{1}{Z}\sum_{\mathbf n}\langle \mathbf n| e^{-\beta(H-\mu \mathscr N)}a_i|\mathbf n\rangle\\
&= \frac{1}{Z}\sum_{\mathbf n}e^{-\beta(E_\mathbf n-\mu N_\mathbf n)}\langle \mathbf n |n_0, \dots, n_{i-1}, n_i-1,n_{i+1}, \dots\rangle\\
&= 0
\end{align}
where the last equality follows from the orthogonality of the basis vectors $|\mathbf n\rangle$.
Let me know if you find an error in this.
