I'm trying to compute the 2nd order correction to the energy spectrum of a 1D quantum harmonic oscillator when a perturbation of the form $\gamma\,\hat{x}^k$ (with $\gamma\ll1$) is added to the Hamiltonian. This is
$$E_n^{(2)}\equiv\langle \psi_n^{(0)} | \,\gamma\hat{x}^k\, | \psi_n^{(1)} \rangle = \sum_{n'\neq n} \frac{{\left|\, \langle \psi_{n'}^{(0)} | \,\gamma\hat{x}^k\, | \psi_n^{(0)} \rangle\,\right|}^{\,2}}{E_n^{(0)}-E_{n'}^{(0)}} = \left(\frac{\hbar}{2m\omega}\right)^k \sum_{n'\neq n} \frac{{ \gamma^2 \left|\, \langle n' | \,(\hat{a}^\dagger+\hat{a})^k\, | n\rangle \,\right|}^{\,2}}{\hbar\omega\,(n-n')}\;.$$
where $\,\hat{x}=\sqrt{\frac{\hbar}{2mw}}\,(\,\hat{a}^\dagger + \hat{a}\,)\;$, $\;E_n^{(0)}=\hbar\omega(n+1/2)\;$ and $\;| \psi_n^{(0)} \rangle\equiv|n\rangle$.
Now, the question is how to compute this term:
$$\alpha_{n'n}\equiv \langle n' | \,(\hat{a}^\dagger+\hat{a})^k\, | n\rangle =\left(\frac{\langle \,0\,|\hat{a}^{n'} (\hat{a}^\dagger+\hat{a})^k {\hat{a}^\dagger}^n\, |\,0\, \rangle}{\sqrt{n'\,!\; n!}}\right)\;.$$
I think it can be "easily" obtained uaing Wick's theorem but I'm not sure about how to apply it. I'm struggling with the part of counting all of the possible ways to fully-contract the operators for any $k$. Do you know any shortcut or some hint to help me compute it?. In fact, is it even manageable?.