Finite temperature greens function in grand canonical ensemble

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)&=-ileftlangle c_{vec{r},sigma}(t),c_{vec{r}’,sigma’}^{dagger}(t’)rightrangle
&=-frac{i}{Z}mathrm{tr}left[e^{-betaleft(hat{H}-muhat{N}right)}e^{ihat{H}t}c_{vec{r},sigma}e^{-ihat{H}t}e^{ihat{H}t’}c_{vec{r}’,sigma’}^{dagger}e^{-ihat{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^{imuleft(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_{isigma}E_{isigma}hat{gamma}^{dagger}_{isigma}hat{gamma}_{isigma}
hat{H}_{G}&=sum_{isigma}tilde{E}_{isigma}hat{tilde{gamma}}^{dagger}_{isigma}hat{tilde{gamma}}_{isigma}
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$?