Number conservation in imaginary time evolution It is clear that if we perform dynamics of the system with hamiltonian commuting with total particle number, this quantity will be an integral of a motion. Is it the case for imaginary time evolution?
To be precise, my question is:
assume we want to find a ground state of a system with a hamiltonian commuting with total number of particles. To do so, we will perform an imaginary time evolution. Our initial state has a definite number of particles. Will the ground obtained during the process have the same number of particles?
UPDATE: let us evolve a state with an operator $e^{-\hat{H}t}$. Then average number of particles evolves according to:
\begin{equation}
N=\frac{\left<\psi\right|e^{-\hat{H}t}\hat{N}e^{-\hat{H}t}\left|\psi\right>}{\left<\psi\right|e^{-2\hat{H}t}\left|\psi\right>}
\end{equation}
Then let us take a derivative by $t$:
\begin{equation}
\dot{N}=-\left<\{\hat{N},\hat{H}\}\right>+2\left<\hat{N}\right>\left<\hat{H}\right>
\end{equation}
Now, using that property that hamiltonian commutes with the particle number, we get
\begin{equation}
\dot{N}=-2\left<\hat{N}\hat{H}\right>+2\left<\hat{N}\right>\left<\hat{H}\right>
\end{equation}
The question is: does the equality
\begin{equation}
\left<\hat{N}\hat{H}\right>=\left<\hat{N}\right>\left<\hat{H}\right>
\end{equation}
neccesarily holds? 
This is the case when $\psi$ is an eigenvector for both energy and particle number operator, for example.
 A: On the track you're on, $⟨AB⟩$ can indeed differ from $⟨A⟩⟨B⟩$ even if $[A,B]=0$. For a trivial counterexample, use $A=B=x$ on $\mathcal H=L_2(\mathbb R)$, for which in general
$$
⟨AB⟩=⟨x^2⟩>⟨x⟩^2=⟨A⟩⟨B⟩
$$
unless the state is a delta function at the origin.
On the other hand, one expects any integral of motion to be conserved on both real and imaginary time evolution. This follows from the fact that
$$
n(t)=⟨\psi|e^{iHt}\hat Ne^{-iHt}|\psi⟩
$$
should be an analytical function around the real axis, and if it is constant on the real axis it must be constant everywhere.
More generally, if $[\hat N,H]=0$ then you expect any evolution under $H$, under both real and imaginary time, to preserve the eigenstates of $N$. Thus, if $\hat N|\psi⟩=N|\psi⟩$ then
$$
\hat Ne^{iH\tau}|\psi⟩=Ne^{iH\tau}|\psi⟩
$$
for all $\tau\in\mathbb C$. The time evolution under $\tau$ will then also preserve the expectation value of $\hat N$.

There is a problem, on the other hand, if you want to look at the expectation value of $\hat N$ in a superposition of shared eigenstates of $\hat N$ and $\hat H$. To see why, take $|\phi_1⟩$ and $|\phi_2⟩$ to be such shared eigenstates with $\hat N|\phi_i⟩=N_i|\phi_i⟩$ and $\hat H|\phi_i⟩=E_i|\phi_i⟩$, with $N_1\neq N_2$ and $E_1<E_2$. Then, if you propagate along the positive imaginary time axis, you get
\begin{align}
|\psi(it)⟩
&=e^{-iH(it)}\frac1{\sqrt{2}}\left(|\phi_1⟩+|\phi_2⟩\right)
\\&=\frac{e^{E_1t}}{\sqrt{2}}|\phi_1⟩+\frac{e^{E_2t}}{\sqrt{2}}|\phi_2⟩
\end{align}
and therefore
\begin{align}
⟨\psi(it)|\hat N|\psi(it)⟩
&=\frac{e^{2E_1t}}{2}N_1+\frac{e^{2E_2t}}{{2}}N_2
\end{align}
so that
\begin{align}
⟨\hat N⟩:=
\frac{⟨\psi(it)|\hat N|\psi(it)⟩}{⟨\psi(it)|\psi(it)⟩}
=\frac{
  e^{2E_1t}N_1+e^{2E_2t}N_2
  }{
  e^{2E_1t}+e^{2E_2t}
  }
=\frac{
  N_1+e^{-2(E_1-E_2)t}N_2
  }{
  1+e^{-2(E_1-E_2)t}
  }.
\end{align}
In particular $⟨\hat N⟩$ changes from $(N_1+N_2)/2$ at time $t=0$ to $N_1$ in the limit $it\to i\infty$. (Note, here, that the heuristic above is not violated, since $⟨\hat N⟩(t)=⟨\psi(t)|\hat N|\psi(t)⟩$ is oscillatory - and not constant - for real times.)

This is sort of your case, but the differences only make the non-conservation useful for you.

we want to find a ground state of a system with a hamiltonian commuting with total number of particles.

The above guarantees that (i) the ground state produced by the procedure will be an eigenstate of $\hat N$, (ii) if you specify at the beginning a particle number that you want in your ground state (i.e. by choosing your initial guess in that eigenspace, which is easy to do (and in fact it is very hard not to do)) then the ground state produced by the procedure will retain that eigenstate, and (iii) if for some reason you mess up and you specify a seed with components in more than one $\hat N$ eigenspace, the imaginary-time procedure will automatically select for you the correct ground state, $\hat N$ eigenstate and all.
