Standing Spherical Wave normalization after decomposition in Plane Waves  
I'm stuck in the following problem, which I tried to outline in the title.
For starting, consider a spherical cavity with its oscillation modes (Spherical Bessel functions and Spherical Harmonics). Let's suppose that this cavity can communicate with a square cavity, with its proper oscillation modes (plane waves). It's axiomatic within quantum statistics that a wavefunction with defined energy corresponds to a single state (apart from spin-degeneracy considerations).However, it's not clear (for me), how this counting is maintained when a particle passes from a system to another. In other words, if we take somehow a boson from the Cartesian Box and throw it in the Spherical Box, could it be that the boson "pulverizes" into X bosons,  X being different from 1? The problem should be possible to solve, because Spherical Waves and Plane Waves are both complete basis for functions over $R^3$. However, things seem not to be so straightforward, since Spherical Waves can be adequately normalized, while Plane Waves gives infinity after a volume integration. 
Looking at the development of a spherical wave [Stratton, p. 417] :
$$\begin{equation} 
\label{eq:1}
J_l (kr) Y_l^m (\theta,\phi) = \frac{{(  -i)}^l}{(4\pi)} \int\limits_0^{2\pi} \int\limits_0^\pi e^{ikr\cos(\gamma)} Y_l^m (\alpha,\beta)\sin(\alpha)d\alpha d\beta 
\end{equation}$$
Being $\gamma$ the angle between the direction $(\theta,\phi)$ and the direction $(\alpha,\beta)$ , such that:
$$
\cos(\gamma) = \sin(\theta)\sin(\alpha)\cos(\phi - \beta) + \cos(\theta)\cos(\alpha)
$$
My idea was: multiply the first eq. for its complex conjugate, change the integration order, isolate the term which gives the modulus of the exponential, and take the rest, after integration over $\theta, \phi$ (not shown) as the measure of X:
$$
J_l (kr) J_l (kr) Y_l^m (\theta,\phi) \overline{{Y_l^m}}(\theta,\phi) = \frac{{( - 1)}^l}{(16 \pi^2 )} \int\limits_0^{\pi} \int\limits_0^{2\pi} \int\limits_0^{\pi} \int\limits_0^{2\pi} e^{ikr\cos(\gamma)} e^{ikr\cos(\gamma')} Y_l^m (\alpha,\beta) \overline{{Y_l^m}}(\alpha',\beta')\sin(\alpha')d\alpha'd\beta'\sin(\alpha) d\alpha d\beta
$$
In order to do this, the angles  $\gamma$ and $\gamma'$ (or their cosines) have to be equal (describing in fact a cone of axis ($\theta, \phi$)). This is equivalent to say that the four integration angles $\alpha, \beta, \alpha', \beta'$  are not independent, and one of them ($\beta$ seems the best choice) can be eliminated using the relation $\cos(\gamma)=\cos(\gamma')$. The calculations, unluckily, are quite cumbersome, with a plethora of trig and inverse functions nested together, so I was wondering if someone can see a better way, maybe by an appropriate change of variables. Thank you in advance.
Update: It could be that complicate angles formulas can be avoided with the following reasoning:
 the angles  $\gamma$ and $\gamma'$ must be equal in order to "kill" the oscillations of the exponential term, but also the angles $\beta$ and $\beta'$ have to be equal for the same schema applied to the exponentials in ${Y_l^m}(\alpha,\beta)\overline{{Y_l^m}}(\alpha',\beta')$, so that at the end $\alpha=\alpha'$ and $\beta=\beta'$, and the integral reduces to:
$$
J_l (kr) J_l (kr) Y_l^m (\theta,\phi) \overline{{Y_l^m}}(\theta,\phi) = \frac{{( - 1)}^l}{(16 \pi^2 )} \int\limits_0^{\pi}  P_l^m (\cos({\alpha}))^2 \sin(\alpha) d\alpha 
$$
This object can be taken as a constant approximation of $J_l (kr) J_l (kr) Y_l^m (\theta,\phi) \overline{{Y_l^m}}(\theta,\phi)$ on $R^3$, being the identity of the two objects non rigorously shown.
Update:
About the proposal given in answer n°1, that is , to use formula
$$
e^{ikr\cos(\theta)}=\sum\limits_{l=1}^\infty{(2l+1)(-i)^lJ_l(kr)P_l(\cos(\theta))}
$$
Once taken the complex conjugate and multiplied, one gets:
$$
1 = \sum\limits_{l=1}^\infty{(2l+1)^2 J_l(kr)^2 P_l(\cos(\theta))^2}
$$
Where I made use of the orthogonality relation for Bessel functions to reduce the double sum. Now, once integrated over $R^3$, but in practice one ends with $\infty= \infty$, which is surely correct, but not very useful.
 A: The question is a little unclear. I am interpreting it as follows:
Given a quantum particle in with wave function
\begin{equation}
\psi_{k \ell m}(\vec{x}) = N j_{\ell}(k r) Y_{\ell m}(\theta,\phi),
\end{equation}
where N is a normalization constant, how can I express the wave function in a plane wave basis? In other words what is the function $\phi_{k\ell m}(\vec{p})$ in the following 
\begin{equation}
\psi_{k \ell m}(\vec{x}) = \int \frac{d^3 p}{(2\pi)^3} e^{i \vec{p}\cdot \vec{x}} \phi_{k\ell m} (\vec{p}).
\end{equation}
The answer is that you need to multiply by the complex conjugate of the new basis function, aka $e^{-i \vec{p}\cdot \vec{x}}$,iand integrate over $\vec{x}$. In bra ket notation this is $\phi_{k\ell m}(\vec{p})=\langle \vec{p} | k \ell m\rangle$. 
This amounts to computing the Fourier tranfrom of $\psi_{k \ell m}(\vec{x})$
\begin{equation}
\phi_{k\ell m}(\vec{p}) = \int d^3 x e^{-i\vec{p}\cdot \vec{x}} \psi_{k\ell m} (\vec{x}).
\end{equation}
In fact this problem has a known solution. Google 'plane wave expansion in spherical Bessel functions'. You should be able to find the answer to this intergal, or at least an expansion of $e^{i \vec{p}\cdot\vec{x}}$ in terms of spherical Bessel functions and spherical harmonics (or Legendre polynomials), which will make the integral easy to do using orthogonality of the $j_{\ell}$ and $Y_{\ell m}$. (If you get stuck on this integral and you've tried looking for find a good resource, including stack exchange, this might be worth a separate question, perhaps on math.SE). 
Finally some comments. You say a state can only be occupied by one boson. In fact a state can be occupied by as many bosons as you like, it is fermions that are forbidden from occupying the same state by the Pauli exclusion principle.  Also you don't need the word 'oscillatory' in 'oscillatory state'. Lastly no particles get 'pulverized'--unless relativistic effects are large. What can happen though is that a single particle's wave function can be a superposition over many basis states. 

Update 12/13/2016
Based on discussions below I think I have a clearer idea of the original question. Apologies for not getting it earlier.
Let's go back to the wavefunction, except choosing $N=k \sqrt{2/\pi}$ so that it is properly normalized:
\begin{equation}
\psi_{k \ell m}(\vec{x}) = \sqrt{\frac{2}{\pi}}\ k \ j_{\ell}(k r) Y_{\ell m}(\theta,\phi).
\end{equation}
This function obeys the normalization condition
\begin{equation}
\int_0^\infty dk \sum_{\ell = 0}^\infty \sum_{m=-\ell}^\ell \psi^\star_{k \ell m}(\vec{x}) \psi_{k \ell m}(\vec{x}') = \delta^{(3)}(\vec{x}-\vec{x}')
\end{equation}
To see this we will use the following identities:
\begin{equation}
\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell Y^\star_{\ell m}(\theta',\phi') Y_{\ell m}(\theta \phi) = \frac{1}{\sin \theta} \delta(\theta-\theta')\delta(\phi-\phi')
\end{equation}
and
\begin{equation}
\int_0^\infty dk k^2 j_\ell (kr) j_\ell (k r') = \frac{\pi}{2 r^2} \delta(r-r')
\end{equation}
Putting everything together we have
\begin{eqnarray}
\int_0^\infty dk \sum_{\ell = 0}^\infty \sum_{m=-\ell}^\ell \psi^\star_{k \ell m}(\vec{x}) \psi_{k \ell m}(\vec{x}') &=& \frac{2}{\pi} \sum_{\ell =-\infty}^\infty \sum_{m=-\ell}^{\ell}  Y^\star_{\ell m}(\theta',\phi') Y_{\ell m}(\theta, \phi) \left(\int_0^\infty dk k^2 j_\ell (kr) j_\ell (k r')\right) \\
&=& \frac{1}{ r^2 \sin(\theta)} \delta(r-r')\delta(\theta-\theta')\delta(\phi-\phi') \\
&=& \delta^{(3)}(\vec{x}-\vec{x}')
\end{eqnarray}
To see that the last line follows from the line before, note that
\begin{equation}
1 = \int d^3 x \delta^{(3)}(\vec{x}-\vec{x}') = \int_0^\infty dr \int_0^\pi d\theta \int_0^{2\pi} d \phi r^2 \sin \theta \delta^{(3)} (\vec{x}-\vec{x}')
\end{equation}
which is consistent with the decomposition of $\delta^{(3)} (\vec{x}-\vec{x}')$ used above.
Note that there is also the 'backwards' normalization condition (which follows from a similar argument)
\begin{equation}
\int d^3 x\ \psi_{k \ell m}(\vec{x})^\star \psi_{k' \ell' m'}(\vec{x}) = \delta(k-k') \delta_{\ell \ell'}\delta_{mm'}
\end{equation}
(In bra ket language, what we've shown schematically is that $\sum_A |A \rangle \langle A | = \sum_B |B \rangle \langle B | = 1$, where $A$ is the $k,\ell,m$ basis and $B$ is the $\vec{x}$ basis).
One last point is how to see that the normalization condition still holds even after expressing $\psi_{k \ell m}(\vec{x})$ to the plane wave basis. In fact this is pretty easy based on what we already have:
\begin{eqnarray}
\delta(k-k') \delta_{\ell \ell'}\delta_{mm'} &=& \int d^3 x\ \psi_{k \ell m}(\vec{x})^\star \psi_{k' \ell' m'}(\vec{x}) \\
&=& \int \frac{d^3 p}{(2\pi)^3} \frac{d^3 p'}{(2\pi)^3}d^3 x\ e^{i (\vec{p} - \vec{p}')\cdot \vec{x}} \phi_{k\ell m}(\vec{p}')^\star \phi_{k' \ell' m'}(\vec{p}) \\
&=& \int \frac{d^3 p}{(2\pi)^3}\phi_{k\ell m}(\vec{p})^\star \phi_{k' \ell' m'}(\vec{p})
\end{eqnarray}
You could also in principle verify this by working out the $\phi_{k \ell m}(\vec{p})$ explicitly using the plane wave expansion, plugging the explicit formulas into the above expression, and show that it all works. But, that's sort of like showing $1+1=2$ by adding and then subtracting some ugly number to both sides.
Let me know if there anything is unclear.
