I wish to solve this integral $$b_{lm}(k) = \frac{1}{2(\hbar)^{9/4}(2\pi)^{5/2}\sqrt{\sigma_{px} \sigma_{py} \sigma_{pz}}} \int_{\theta_k = 0}^{\pi}\int_{\varphi_k = 0}^{2\pi} i^l \text{exp}\left[ - \frac{1}{(2\hbar)^2}\left(\frac{(k_z - k_{z0})^2}{\sigma_{pz}^2} + \frac{(k_y - k_{y0})^2}{\sigma_{py}^2} + \frac{(k_x - k_{x0})^2}{\sigma_{px}^2}\right)\right]Y^m_l(\theta_k, \varphi_k) \sin\theta_k d\phi_k d\theta_k$$
which I've obtained from decomposition of a Gaussian into the eigenfunctions of Hydrogen atom. $$\psi(\textbf{r},0) = \sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l} b_{\textit{nlm}}\varphi_\textit{nlm}(\textbf{r})$$ by performing the scalar product to obtain the expansion coefficient, as such. $$\begin{equation*} \begin{split} \langle\varphi_\textit{nlm}(\textbf{r})|\psi(\textbf{r},0)\rangle &= \langle\varphi_\textit{nlm}(\textbf{r})|\sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l} b_{\textit{nlm}}\varphi_\textit{nlm}(\textbf{r})\rangle\\ &= b_{\textit{nlm}}\sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l}\langle\varphi_\textit{n'l'm'}(\textbf{r})| \varphi_\textit{nlm}(\textbf{r})\rangle\\ &= b_{\textit{nlm}}\sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l}\delta_{n'n}\delta_{l'l}\delta_{m'm}\\ &= b_{nlm} \end{split} \end{equation*}$$ Now using two different way of expressing the Gaussian $\psi(\textbf{r},0)$ allowed me to end up with 2 different terms of $b_{nlm}$ which I was able to compare. $$b_{nlm} = \frac{2}{\pi}\int^{\infty}_{k=0}\int^{\infty}_{r=0} R^*_{nl}(\textbf{r}) b_{lm}(k) j_l(kr) r^2 k^2 dr dk <=> N \int^{\pi}_{\theta_k=0}\int^{2\pi}_{\varphi_k=0}\int^{\infty}_{k=0}\int^{\infty} _{r=0} R^*_{nl}(\textbf{r}) i^l j_l(kr) \eta(\textbf{k}-\textbf{k}_0)Y^m_l(\theta_k, \varphi_k) r^2 k^2 \sin\theta_k dr dk d\varphi_k d\theta_k$$
My approach
Now I thought if I set the initial conditions $k_{x0} = k_{y0} = k_{z_0} = 0$ and set all the $\sigma_{pz} = \sigma_{py} = \sigma_{pz} = \sigma$ we see that the integral in the beginning $$b_{lm}(k) = \frac{1}{2(\hbar)^{9/4}(2\pi)^{5/2}\sqrt{\sigma_{px} \sigma_{py} \sigma_{pz}}} \int_{\theta_k = 0}^{\pi}\int_{\varphi_k = 0}^{2\pi} i^l \text{exp}\left[ - \frac{|\textbf{k}|^2}{(2\sigma\hbar)^2}\right]Y^m_l(\theta_k, \varphi_k) \sin\theta_k d\phi_k d\theta_k$$ Now that means I'll be able to move the Gaussian out of the integral $$b_{lm}(k) = \frac{1}{2(\hbar)^{9/4}(2\pi)^{5/2}\sqrt{\sigma_{px} \sigma_{py} \sigma_{pz}}} \text{exp}\left[ - \frac{k^2}{(2\sigma\hbar)^2}\right]\int_{\theta_k = 0}^{\pi}\int_{\varphi_k = 0}^{2\pi} i^l Y^m_l(\theta_k, \varphi_k) \sin\theta_k d\phi_k d\theta_k$$ and then I can perform the rest of the integral and all is fine. But the problem is my project requires for me such that $k_{y0} \neq 0$ and $k_{x0} = k_{z0} = 0$ And for that to happen I'll have to change my Gaussian into spherical polar coordinate using $$k_x = k\cos\theta_k\sin\varphi_k$$ $$k_y = k\sin\theta_k\sin\varphi_k$$ $$k_z = k\cos\varphi$$ $$k^2 = k_x^2 + k_y^2 + k_z^2$$ which makes the previous manoeuvre not possible. From what my supervisor said he told me to find a form of Spherical Harmonics that's not defined at the origin? So I'm assuming that there's a form of Spherical Harmonics which I could start at a coordinates same as the Gaussian then I should be able to solve the integral.