I am trying to calculate the time evolution of the operator \begin{equation} h(k)=\sum_k b_k^{\dagger}b_k\, . \end{equation} Therefore, I go to the Heisenberg picture $$ h(k ,t) \equiv e^{\frac{i}{\hbar}Ht}\,\left( \sum_k b_k^{\dagger}b_k\right)\, e^{-\frac{i}{\hbar}Ht} \, , $$ where $H$ is the Bose-Hubbard Hamiltonian in $k-$space $$ H = \sum_k (\epsilon_k-\mu) b_k^{\dagger}b_k +\frac{g}{2V}\sum_{k,p,q}b_{k+q}^{\dagger}b_{p-q}^{\dagger}b_k b_p\, , $$ and the $b-$operators fulfil the bosonic commutation relation $[b_k,b_p^{\dagger}]=\delta_{k,p}$. Now, I want to evaluate the above time dependent operator $h(k,t)$ by using the BCH formula $$ e^XYe^{-X}=\sum_{m=0}^{\infty}\frac{1}{m!}[X,Y]_m\qquad\text{with}\quad [X,Y]_0=Y \, \, \text{and}\, \, [X,Y]_m=[X,[X,Y]_{m-1}]\, . $$ For $m=1$, I calculated the following commutator relation $$ [H,h_{KE}(k)] = \\ = \frac{g}{2V}\sum_{k,p,q,u} b_{k +q}^{\dagger}b_{p -q}^{\dagger}b_k b_u\delta_{p ,u}-b_{u}^{\dagger}b_{p -q}^{\dagger}b_k b_p\delta_{u ,k +q}+ b_{k +q}^{\dagger}b_{p -q}^{\dagger}b_p b_u\delta_{k ,u}-b_{u}^{\dagger}b_{k +q}^{\dagger}b_k b_p\delta_{u ,p -q}\, . $$ As you can see, there are Kronecker deltas in every term. I am not sure how to evaluate them. My first attempt: Separate the terms since the summation $\Sigma$ is linear $$\begin{align} &[H,h_{KE}(k)] = &\\ &= \frac{g}{2V}\left(\sum_{k,p,q,u} b_{k +q}^{\dagger}b_{p -q}^{\dagger}b_k b_u\delta_{p ,u}-\sum_{k,p,q,u}b_{u}^{\dagger}b_{p -q}^{\dagger}b_k b_p\delta_{u ,k +q}+ \\ +\sum_{k,p,q,u}b_{k +q}^{\dagger}b_{p -q}^{\dagger}b_p b_u\delta_{k ,u}-\sum_{k,p,q,u}b_{u}^{\dagger}b_{k +q}^{\dagger}b_k b_p\delta_{u ,p -q}\right)\, . \end{align}$$ Now, we can consider every Kronecker delta separately and set in the first term $p=u$, in the second $u=k+q$, in the third $k=u$ and in the fourth $u=p-q$. However, this renders $0$. I felt like that my approach is to naive and I actually would not expect $0$ to be the correct answer. So as a second attempt, I tried to get rid of index summations in the Kronecker deltas. So in the 2nd term, I set $k+q=n$ and in the fourth term, I set $p-q=m$ , and in both cases $q$ was replaced by writing it in terms of n and m, respectively. Then I executed the Kronzucker deltas and set $n=q$ in the second term and $m=q$ in the fourth term. $$\begin{align} &[H,h_{KE}(k)] =\\ &= \frac{g}{2V}\left(\sum_{k,p,q} b_{k +q}^{\dagger}b_{p -q}^{\dagger}b_k b_p-\sum_{k,p,q}b_{q}^{\dagger}b_{p -q+k}^{\dagger}b_k b_p+\\+\sum_{k,p,q} b_{k +q}^{\dagger}b_{p -q}^{\dagger}b_k b_p-\sum_{k,p,q}b_{q}^{\dagger}b_{k+p-q}^{\dagger}b_k b_p\right) \end{align}$$
My Question is: Is this a correct way of calculating this expression? I was hoping that the Kronzucker deltas will simplify my expression a little more, after all I have to calculate higher order terms $m$ and it will get too messy... Thank you in advance
Ilias