Using Symmetry to Identify When Electron-Phonon Interaction Energy Corrections Are Zero General question:
Under what circumstances are terms involving the expectation value of products of fermionic/bosonic number operators zero? How can one use symmetry to determine when this will be the case?
Specific  Example:
I'm reading through chapter 6, section 6.2 of Taylor's "A Quantum Approach to Condensed Matter Physics", on Electron-Phonon Interactions. The section applies second-order perturbation theory to calculate the effect of the (assumed weak) electron-phonon interactions on the total energy of the system, as described by the Fröhlich Hamiltonian
\begin{equation}
\mathcal{H} = \sum_{\mathbf{k}}\mathcal{E}_{k}\,c^{\dagger}_{\mathbf{k}}\,c_{\mathbf{k}}+\sum_{\mathbf{q}}\hbar\omega_{\mathbf{q}}\,a^{\dagger}_{\mathbf{q}}\,a_{\mathbf{q}} + \sum_{\mathbf{k},\mathbf{k'}}M_{\mathbf{k}\mathbf{k'}}\,(a^{\dagger}_{-\mathbf{q}}+a_{\mathbf{q}})\,c^{\dagger}_{\mathbf{k}}\,c_{\mathbf{k}}.
\end{equation}
In Eq. (6.2.2), the contribution of the second-order terms is found to be
\begin{equation}
\mathcal{E}_{2}=\sum_{\mathbf{k},\mathbf{k'}}|M_{\mathbf{k}\mathbf{k'}}|^{2}\langle n_{\mathbf{k}}(1-n_{\mathbf{k'}}) \rangle\left(\frac{\langle n_{\mathbf{-q}}\rangle}{\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}}+\hbar\omega_{\mathbf{-q}}}+\frac{\langle n_{\mathbf{q}}+1\rangle}{\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}}-\hbar\omega_{\mathbf{-q}}}\right),
\end{equation}
where $\langle n_{\mathbf{k}}\rangle,\langle n_{\mathbf{k'}}\rangle$ are electron occupation numbers while $\langle n_{\mathbf{q}}\rangle,\langle n_{\mathbf{-q}}\rangle$ are phonon occupation numbers. I understood the steps up to this point. To reach the next equation, the author writes both terms in the parentheses with a common denominator and multiplies in the $\langle (1-n_{\mathbf{k'}}) \rangle$ term, obtaining
\begin{equation}
\mathcal{E}_{2}=\sum_{\mathbf{k},\mathbf{k'}}|M_{\mathbf{k}\mathbf{k'}}|^{2}\langle n_{\mathbf{k}} \rangle\left[\frac{2(\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}})\langle n_{\mathbf{q}}\rangle}{(\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}})^2-(\hbar\omega_{\mathbf{q}})^2}+\frac{1-\langle n_{\mathbf{k'}}\rangle}{\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}}-\hbar\omega_{\mathbf{q}}}\right].
\end{equation}
However, to obtain this equation, one is required to drop a term in $\langle n_{\mathbf{k}}n_{\mathbf{k'}}n_{\mathbf{q}}\rangle$. It is claimed that this term goes to zero "by symmetry", but what symmetry is being applied in this case? Does it pertain to the parity of the many-body wave function under particle permutations?
 A: The equation you have
\begin{equation}
\mathcal{E}_{2}=\sum_{\mathbf{k},\mathbf{k'}}|M_{\mathbf{k}\mathbf{k'}}|^{2}\langle n_{\mathbf{k}}(1-n_{\mathbf{k'}}) \rangle\left(\frac{\langle n_{\mathbf{-q}}\rangle}{\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}}+\hbar\omega_{\mathbf{-q}}}+\frac{\langle n_{\mathbf{q}}+1\rangle}{\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}}-\hbar\omega_{\mathbf{-q}}}\right) \tag{1}
\end{equation}
by simple calculations, is equal to
\begin{equation}
\mathcal{E}_{2}=\sum_{\mathbf{k},\mathbf{k'}}|M_{\mathbf{k}\mathbf{k'}}|^{2}\big(\langle n_{\mathbf{k}} \rangle-\langle n_{\mathbf{k}}n_{\mathbf{k'}} \rangle\big)\left(\frac{2(\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}})\langle n_{\mathbf{q}}\rangle}{(\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}})^2-(\hbar\omega_{\mathbf{q}})^2}+\frac{1}{\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}}-\hbar\omega_{\mathbf{q}}}\right) \tag{2}
\end{equation}
Eq. (2) can be simplified to your final result if
\begin{equation}
\mathcal{E}_3=\sum_{\mathbf{k},\mathbf{k'}}|M_{\mathbf{k}\mathbf{k'}}|^{2}\langle n_{\mathbf{k}}n_{\mathbf{k'}} \rangle\left(\frac{2(\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}})\langle n_{\mathbf{q}}\rangle}{(\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}})^2-(\hbar\omega_{\mathbf{q}})^2}\right)=0 \tag{3}
\end{equation}
This is true by the symmetry that we can exchange $\mathbf{k}$ and $\mathbf{k'}$, having
\begin{equation}
\begin{split}
\mathcal{E}_3 & = \sum_{\mathbf{k},\mathbf{k'}}|M_{\mathbf{k}\mathbf{k'}}|^{2}\langle n_{\mathbf{k}}n_{\mathbf{k'}} \rangle\left(\frac{2(\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}})\langle n_{\mathbf{q}}\rangle}{(\mathcal{E}_{\mathbf{k}}-\mathcal{E}_{\mathbf{k'}})^2-(\hbar\omega_{\mathbf{q}})^2}\right) \\
& = \sum_{\mathbf{k'},\mathbf{k}}|M_{\mathbf{k'}\mathbf{k}}|^{2}\langle n_{\mathbf{k'}}n_{\mathbf{k}} \rangle\left(\frac{2(\mathcal{E}_{\mathbf{k'}}-\mathcal{E}_{\mathbf{k}})\langle n_{\mathbf{q}}\rangle}{(\mathcal{E}_{\mathbf{k'}}-\mathcal{E}_{\mathbf{k}})^2-(\hbar\omega_{\mathbf{q}})^2}\right)=-\mathcal{E}_3
\end{split}
\tag{4}
\end{equation}
Eq. (4) encertains $\mathcal{E}_3=0$, and the term containing $\langle n_{\mathbf{k}}n_{\mathbf{k'}}n_{\mathbf{q}} \rangle$ can be dropped off, which answer your question.
