I'm having trouble with trying to find the expansion coefficients of a superposition of a Gaussian wave packet.
First I'm decomposing a Gaussian wave packet $$\psi(\textbf{r},0) = \frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{(\textbf{r} - \textbf{r}_0)^2}{4\sigma^2} + i\textbf{k}_0 \cdot \textbf{r}\right]$$ into a linear superposition of eigenfunctions of the hydrogen atom, where the initial position and momentum is chosen as $\textbf{r}_0$ and $\textbf{p}_0 = \hbar \textbf{k}_0$, and also a spatial width $\sigma$ such that the initial wave packet $\psi(\textbf{r},0)$ can to a sufficient approximation be superimposed by bound-state eigenfunctions $\varphi_\textit{nlm}(\textbf{r})$, $$\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})$$ Having solved the Schrodinger Equation of a hydrogen atom with a Hamiltonian containing the Coulomb potential, the eigenfunctions were obtained to be a product of the Radial Wave solution and Spherical Harmonics $$\varphi_\textit{nlm}(\textbf{r}) = R_{nl}(r) Y_{lm}(\theta, \phi)$$ where $$R_{nl}(r) = \frac{1}{a^{3/2}}\frac{2}{n^2}\sqrt{\frac{(n-l-1)!}{(n+1)!}}\left(\frac{2r}{na}\right)^l e^{-\frac{r}{na}} L^{2l+1}_{n-l-1}\left(\frac{2r}{na}\right)$$ where $L^{2l+1}_{n-l-1}\left(\frac{2r}{na}\right)$ is the Laguerre Polynomial which has the form of $$L^{k}_{p}\left(x\right) = \sum_{k}^{p} (-1)^s \binom {p+k}{p-s} \frac{x^s}{s!}$$ and $$Y_{lm}(\theta, \phi) = (-1)^m \sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}} P^{m}_{l} (\cos\theta) e^{im\phi}$$ with the Legendre Polynomial contained within the Associated Legendre Functions $$P^{m}_{l} (x) = (1 - x^2)^{\frac{m}{2}} \frac{d^m}{dx^m}(P_{l}(x))$$ $$P_{l}(x) = \frac{1}{2^l l! }\frac{d^l}{dx^l}[(x^2 - 1)^l]$$
What I ultimately would like to find is the closed form expansion coefficients $b_{nlm}$ of the eigenfunction superposition, now from my knowledge of finding the expansion coefficients of a single-indexed superposition of a wavefunction, the same method can be applied on this superposition which involves 3 different indices, such that $$ \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*} $$
My attempt: To avoid doing the scalar product directly and having to compute the whole integral, I was thinking of using the decomposition of stationary harmonic plane wave into partial waves, as such $$e^{\textbf{k}\cdot\textbf{r}} = e^{ikz} = e^{ikr\cos\theta} = \sum^{l=0}_{\infty}(2l+1)i^l j_l(kr) P_l(\cos\theta)$$ where $i$ is the imaginary number, $\textbf{k}$ and $\textbf{r}$ are the wave and position vectors, $j_l$ are the spherical Bessel functions, $P_l$ are Legendre Polynomials. But according to Wikipedia (https://en.wikipedia.org/wiki/Plane_wave_expansion), it can also be rewritten as $$e^{\textbf{k}\cdot\textbf{r}} = 4\pi \sum_{l=0}^{\infty}\sum_{m=-l}^{l} i^l j_l(kr)Y^m_l(\hat{\textbf{k}})Y^{m*}_l(\hat{\textbf{r}})$$
In theory, this should now mean that I can rewrite the scalar product in terms of these sum of spherical waves, so that $$ \begin{equation*} \begin{split} b_{nlm} &= \langle\varphi_\textit{nlm}(\textbf{r})|\psi(\textbf{r},0)\rangle\\ &= \langle R_{nl}(r) Y_{lm}(\theta, \phi)|\frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{(\textbf{r} - \textbf{r}_0)^2}{4\sigma^2}\right]\cdot 4\pi\sum_{l=0}^{\infty}\sum_{m=-l}^{l} i^l j_l(kr)Y^m_l(\hat{\textbf{k}})Y^{m*}_l(\hat{\textbf{r}})\rangle \end{split} \end{equation*} $$
and so this is where I'm at. I'm not sure how to go from here. I hope I've clarified my question well.