Grand canonical hamiltonian and thermodynamic ensembles In many-body theory (and quantum field theory I suppose) we often work in the grand canonical ensemble, where the number of particles in the system is only fixed on average. The density operator used to compute expectation values is 
$$
\rho = \frac{e^{-\beta (H - \mu N)}}{Z_{G}} = \frac{e^{-\beta K}}{Z_{G}}
$$
where $K$ is the so-called grand canonical hamiltonian.
My problem is when we substitute $H$ for $K$ in the expression of the evolution operator
$
U(t) = e^{-iHt} \rightarrow e^{-iKt} 
$, which is done most of the time because it simplifies calculations. It seems equivalent to saying that Schrodinger equation is invariant under the change $H \rightarrow K$.
The justifications I've seen so far, are based on the fact that since the original hamiltonian conserves the number of particles and thus commutes with the operator $N$, this replacement is just a displacement in energy and does not essentially changes the dynamics of the system. I'm not really convinced by this because the replacement with $K$ is not equivalent with adding a simple constant to the hamiltonian. 
Additionally, it seems that if this argument is true, then in general one would be justified to construct a new hamiltonian $H' = H + \mathcal{O}$ to describe the dynamics of a system, as long as $[H, \mathcal{O}]=0$.
Given that, my questions are as follow :


*

*Is the above statement that you can replace a hamiltonian $H' = H + \mathcal{O}$ when $[H, \mathcal{O}]=0$ true? If it'sn't, is there anything special with the case $\mathcal{O} = \mu N$, or are there some caveats?

*Thermodynamics ensembles are often defined by density operators on the form $\rho=e^{-\beta S}$ with $S$ that is a linear combination of $H$ and various conserved operators $c_j \mathcal{O}_j$. Could we get rid of those ensembles by just including the operators $c_j \mathcal{O}_j$ in the hamiltonian in the first place and finding $c_j$ such that $\langle \mathcal{O_j} \rangle$ is the value we want? $c_j$ would have the physical interpretation of a classical field that couples with $\mathcal{O}_j$.

 A: Actually, replacing $H \to H - \mu N$ does affect the dynamics substantially, but it does so in a rather trivial way. This is precisely because $H$ and $N$ commute. If two operators $A$ and $B$ commute, and only then, we can factor an exponential as $e^{A + B} = e^{A} e^{B}$. This allows us to conclude that $e^{-iHt} = e^{iN\mu t} e^{-iKt}$. That is, we are free to evolve with $K$ as long as we remember that at the end of the day we have to apply $e^{iN\mu t}$ to recover the "true" time evolution. This is generally pretty easy since $N$ is a very simple operator. Another way to think about it is that we are working in a "rotating frame" -- we study the evolution not of the actual state of the system $|\psi(t)\rangle$, but of the "rotated" state $e^{iN\mu t} |\psi(t)\rangle$.
The answer above, which claimed that the differences are insignificant because $\sqrt{ \langle (N - \bar{N})^2 \rangle }$ is small, is not correct.  In the thermodynamic limit, this quantity still goes to infinity. Its ratio with $\bar{N}$ goes to zero, but it would have to go to zero in absolute value to not affect the dynamics, and this is not the case. In fact, one can verify that, for example, the expectation value (in a bosonic system) of a boson creation operator $\langle a^{\dagger} \rangle$ rotates at a different frequency under $K$ than under $H$ (specifically, the difference in frequency is $\mu$). This is precisely accounted for by the extra rotation by $e^{iN \mu t}$ discussed above.
A: The Schrödinger equation is not invariant under the change $H\rightarrow K$. However, in the grand canonical ensemble with some values of inverse temperature $\beta$ and chemical potential $\mu$, the fluctuation of the number of particles $N$ is of the order of $\sqrt{\langle N \rangle }$, where $\langle N \rangle $ is the average number of particles in this ensemble. This fluctuation is much less than $\langle N\rangle$, so  the states of this ensemble most probably have a number of particles that is very close to $ \langle N \rangle $, and $K=H-\mu N$ gives pretty much the same evolution as $K^\prime=H-\mu \langle N\rangle$ for this specific grand canonical ensemble. And $K^\prime$ is indeed just $H$ with a constant shift.  
A: The discussion seems to be finished long ago, but I'll put in my 5 cents anyway.
This is a very good question that I have been trying to answer in order to understand Matsubara's Green function formalism. Surprisingly, this question is not addressed in the textbooks that I am aware of.
The condition $[\hat{\mathcal{O}}, \hat H] = 0$ is not sufficient to be able to replace the time evolution operator with $e^{i(\hat H + \hat{\mathcal{O}})t}$ and rotate the state to $e^{-i\hat{\mathcal{O}}t}|\psi\rangle$. The special thing about the particle number operator $\hat N$ is that it is used in the second quantization picture and acts on the states represented by occupation numbers, $|n_1, n_2, \ldots\rangle$. Then, the rotated state $e^{-i\mu\hat Nt}|n_1, n_2, \ldots\rangle = e^{-i\mu (n_1 + n_2 + \ldots)t}|n_1, n_2, \ldots\rangle$ differs from the original state by only a phase factor. Because the phase of a wave function is not observable, the rotated state $e^{-i\mu\hat Nt}|n_1, n_2, \ldots\rangle$ is physically equivalent to the original state $|n_1, n_2, \ldots\rangle$. For an arbitrary operator $e^{i\hat{\mathcal{O}}t}$, this is not necessarily the case, because $|\psi\rangle$ may be a superposition of many eigenstates of $\hat{\mathcal{O}}$; if you rotate each of those contributions with $e^{-i\hat{\mathcal{O}}t}$, they will acquire different phase factors.
The term $-\mu\hat N$ is added to the Hamiltonian, thereby turning it into the grand canonical Hamiltonian, in order to be able to develop the Green's function method at finite temperatures. It is nice to have the same Hamiltonian in the statistical operator and in the time evolution operator, so objects like $G(t) = -i\,\text{Tr}(e^{-\beta(\hat H - \mu\hat N)}[e^{i(\hat H - \mu\hat N)t}\,\hat A\,e^{-i(\hat H - \mu\hat N)t}, \hat B])$ could be manipulated more easily.
There is a further assumption inherent in the use of the same Hamiltonian in the canonical (or grand canonical) statistical operator, $e^{-\beta\hat{H}}$, and in the time evolution operator, $e^{i\hat{H}t}$. Heisenberg-picture operators evolve according to $e^{i\hat H t}\hat Ae^{-i\hat H t}$ only if the system is isolated. But (grand) canonical distribution implies that the system is not isolated; rather, it is coupled to a much larger heat bath, with which our system exchanges energy (and particles in the grand canonical case). So, strictly speaking, we should include the Hamiltonian that describes the system-bath coupling into the operator: $\hat A(t) = e^{i(\hat H + \hat H_{SB})t}\hat Ae^{-i(\hat H + H_{SB})t}$. Neglecting the coupling Hamiltonian is an approximation. The system-bath coupling must be weak.
