I see this question was asked several times before but I don't think any answer can explain the issue perfectly. I am studying many body theory and encounters finite temperature Green's function. At first glance, it seems to me that the correct time dependence of, let's say a greater green's function, would be
$$ \begin{align*} G^{>}\left(\vec{r},\sigma,t;\vec{r}',\sigma',t'\right)&=-i\left\langle c_{\vec{r},\sigma}(t),c_{\vec{r}',\sigma'}^{\dagger}(t')\right\rangle \\ &=-\frac{i}{Z}\mathrm{tr}\left[e^{-\beta\left(\hat{H}-\mu\hat{N}\right)}e^{i\hat{H}t}c_{\vec{r},\sigma}e^{-i\hat{H}t}e^{i\hat{H}t'}c_{\vec{r}',\sigma'}^{\dagger}e^{-i\hat{H}t'}\right] \end{align*} $$
However, some texts used the "grand canonical Hamiltonian" $H_{G}=H-\mu N$ for the operators. Some said it is a natural way to "define" the time evolution, some don't even say anything at the beginning and include $\mu$ in all the Hamiltonian, some say it is merely a parameter but not chemical potential.
Of course if $\left[\hat{H},\hat{N}\right]=0$, such a change in Hamiltonian just introduces an extra phase factor like $e^{i\mu\left(t-t'\right)}$ to the Green's function and the change can be reversed at the final result easily. BUT if $\left[\hat{H},\hat{N}\right]\neq 0$ like a BCS hamiltonian, diagonalizing $\hat{H}$ and diagonalizing $\hat{H}_{G}$ give two completely different set of "eigen-operators" like $$ \begin{align*} \hat{H}&=\sum_{i\sigma}E_{i\sigma}\hat{\gamma}^{\dagger}_{i\sigma}\hat{\gamma}_{i\sigma}\\ \hat{H}_{G}&=\sum_{i\sigma}\tilde{E}_{i\sigma}\hat{\tilde{\gamma}}^{\dagger}_{i\sigma}\hat{\tilde{\gamma}}_{i\sigma} \end{align*} $$ In this case, changing $H$ to $H_G$ in the green's function will not just produce a phase factor but instead you are really calculating something very different.
Is there some theoretical reason that requires the time evolution is generated by $H_G$ instead of $H$?
