[proofplan]
We first transfer $L^2(\mathbb{R}^3)$ to polar coordinates, identifying it unitarily with the Bochner space $L^2((0,\infty), r^2\,d\mathcal{L}^1(r); L^2(S^2,\sigma))$. In that space, each radial slice is an element of $L^2(S^2,\sigma)$, so it can be expanded in the [orthonormal basis](/page/Orthonormal%20Basis) of spherical harmonics. [Parseval's identity](/theorems/434) on $S^2$, integrated over the radial variable, gives the norm identity and proves that the coefficient functions lie in the radial $L^2$ space. Orthogonality of the spherical harmonics gives orthogonality of the angular momentum sectors, and completeness gives the direct sum decomposition.
[/proofplan]
[step:Transfer $L^2(\mathbb{R}^3)$ to polar coordinates]
Let
\begin{align*}
P:(0,\infty)\times S^2 \to \mathbb{R}^3\setminus\{0\}
\end{align*}
be the polar coordinate map defined by $P(r,\omega) := r\omega$. Since $\mathcal{L}^3(\{0\}) = 0$, changing a function on $\{0\}$ does not change its class in $L^2(\mathbb{R}^3,\mathcal{L}^3)$.
Define the measure $\mu$ on $(0,\infty)$ by
\begin{align*}
d\mu(r) := r^2\,d\mathcal{L}^1(r).
\end{align*}
By the polar-coordinate formula for [Lebesgue measure](/page/Lebesgue%20Measure) (citing a result not yet in the wiki: polar-coordinate formula for Lebesgue measure), for every non-negative measurable function $f:\mathbb{R}^3 \to [0,\infty]$,
\begin{align*}
\int_{\mathbb{R}^3} f(x)\,d\mathcal{L}^3(x) = \int_0^\infty \int_{S^2} f(r\omega)\,d\sigma(\omega)\,r^2\,d\mathcal{L}^1(r).
\end{align*}
Applying this identity to $f = |\psi|^2$ shows that the map
\begin{align*}
U:L^2(\mathbb{R}^3,\mathcal{L}^3) \to L^2((0,\infty),\mu;L^2(S^2,\sigma))
\end{align*}
defined by
\begin{align*}
(U\psi)(r)(\omega) := \psi(r\omega)
\end{align*}
is an isometry. The same polar-coordinate formula also shows that every element of the Bochner space corresponds to an $L^2$ function on $\mathbb{R}^3$ after composition with $P^{-1}$ outside the null set $\{0\}$, so $U$ is unitary.
[guided]
The first point is that the origin is harmless in $L^2(\mathbb{R}^3,\mathcal{L}^3)$ because $\mathcal{L}^3(\{0\}) = 0$. Therefore a function on $\mathbb{R}^3$ is determined, as an $L^2$ class, by its values on $\mathbb{R}^3\setminus\{0\}$, where polar coordinates are available.
We introduce the polar coordinate map
\begin{align*}
P:(0,\infty)\times S^2 \to \mathbb{R}^3\setminus\{0\}
\end{align*}
by
\begin{align*}
P(r,\omega) := r\omega.
\end{align*}
We also define the radial measure $\mu$ on $(0,\infty)$ by
\begin{align*}
d\mu(r) := r^2\,d\mathcal{L}^1(r).
\end{align*}
The factor $r^2$ is exactly the Jacobian factor for polar coordinates in dimension $3$.
The polar-coordinate formula for Lebesgue measure says that for every non-negative measurable function $f:\mathbb{R}^3 \to [0,\infty]$,
\begin{align*}
\int_{\mathbb{R}^3} f(x)\,d\mathcal{L}^3(x) = \int_0^\infty \int_{S^2} f(r\omega)\,d\sigma(\omega)\,r^2\,d\mathcal{L}^1(r).
\end{align*}
This is the measure-theoretic statement behind the informal notation $d x = r^2\,dr\,d\omega$.
Now let $\psi \in L^2(\mathbb{R}^3,\mathcal{L}^3)$. Define
\begin{align*}
U:L^2(\mathbb{R}^3,\mathcal{L}^3) \to L^2((0,\infty),\mu;L^2(S^2,\sigma))
\end{align*}
by
\begin{align*}
(U\psi)(r)(\omega) := \psi(r\omega).
\end{align*}
Applying the polar-coordinate formula to $f = |\psi|^2$ gives
\begin{align*}
\|\psi\|_{L^2(\mathbb{R}^3)}^2 = \int_0^\infty \int_{S^2} |\psi(r\omega)|^2\,d\sigma(\omega)\,r^2\,d\mathcal{L}^1(r).
\end{align*}
The right-hand side is exactly
\begin{align*}
\|U\psi\|_{L^2((0,\infty),\mu;L^2(S^2,\sigma))}^2.
\end{align*}
Thus $U$ preserves the $L^2$ norm.
Finally, $U$ is onto because the polar coordinate map $P$ parametrizes $\mathbb{R}^3\setminus\{0\}$, and the missing point has $\mathcal{L}^3$-measure zero. Given a Bochner-square-integrable function $F:(0,\infty)\to L^2(S^2,\sigma)$, the formula $x=r\omega \mapsto F(r)(\omega)$ defines an $L^2(\mathbb{R}^3,\mathcal{L}^3)$ function after choosing representatives, and the polar-coordinate formula gives its square integrability. Hence $U$ is a unitary identification.
[/guided]
[/step]
[step:Expand each angular slice in spherical harmonics]
Let $F := U\psi$. For each admissible pair $(l,m)$, define the coefficient function
\begin{align*}
R_{l,m}:(0,\infty) \to \mathbb{C}
\end{align*}
by
\begin{align*}
R_{l,m}(r) := \int_{S^2} F(r)(\omega)\overline{Y_{l,m}(\omega)}\,d\sigma(\omega)
\end{align*}
for every $r$ for which $F(r) \in L^2(S^2,\sigma)$ is represented by a measurable angular function.
Since $(Y_{l,m})$ is an orthonormal basis of $L^2(S^2,\sigma)$, Parseval's identity in $L^2(S^2,\sigma)$ gives, for $\mu$-almost every $r \in (0,\infty)$,
\begin{align*}
\|F(r)\|_{L^2(S^2)}^2 = \sum_{l=0}^{\infty}\sum_{m=-l}^{l} |R_{l,m}(r)|^2.
\end{align*}
Also, for $\mu$-almost every $r$,
\begin{align*}
F(r) = \sum_{l=0}^{\infty}\sum_{m=-l}^{l} R_{l,m}(r)Y_{l,m}
\end{align*}
with convergence in $L^2(S^2,\sigma)$.
Integrating Parseval's identity with respect to $\mu$ gives
\begin{align*}
\int_0^\infty \|F(r)\|_{L^2(S^2)}^2\,r^2\,d\mathcal{L}^1(r) = \sum_{l=0}^{\infty}\sum_{m=-l}^{l}\int_0^\infty |R_{l,m}(r)|^2 r^2\,d\mathcal{L}^1(r).
\end{align*}
The left-hand side is finite because $F \in L^2((0,\infty),\mu;L^2(S^2,\sigma))$. Therefore each $R_{l,m}$ belongs to $L^2((0,\infty),r^2\,d\mathcal{L}^1(r))$.
[guided]
The unitary map $U$ has turned $\psi$ into a function $F := U\psi$ whose input is the radius $r$ and whose value $F(r)$ is an angular $L^2$ function on $S^2$. For each fixed admissible pair $(l,m)$, we extract the coefficient of $F(r)$ in the $Y_{l,m}$ direction.
Define
\begin{align*}
R_{l,m}:(0,\infty) \to \mathbb{C}
\end{align*}
by
\begin{align*}
R_{l,m}(r) := \int_{S^2} F(r)(\omega)\overline{Y_{l,m}(\omega)}\,d\sigma(\omega).
\end{align*}
This integral is meaningful for $\mu$-almost every $r$, because $F(r) \in L^2(S^2,\sigma)$ and $Y_{l,m} \in L^2(S^2,\sigma)$; the [Cauchy-Schwarz inequality](/theorems/432) on $L^2(S^2,\sigma)$ gives finite absolute value for the [inner product](/page/Inner%20Product).
Now we use the defining property of spherical harmonics: the family $(Y_{l,m})_{l \in \mathbb{N}_0,\,-l \le m \le l}$ is an orthonormal basis of $L^2(S^2,\sigma)$. Hence Parseval's identity in the [Hilbert space](/page/Hilbert%20Space) $L^2(S^2,\sigma)$ gives, for $\mu$-almost every radius $r$,
\begin{align*}
\|F(r)\|_{L^2(S^2)}^2 = \sum_{l=0}^{\infty}\sum_{m=-l}^{l} |R_{l,m}(r)|^2.
\end{align*}
The same Hilbert space expansion gives
\begin{align*}
F(r) = \sum_{l=0}^{\infty}\sum_{m=-l}^{l} R_{l,m}(r)Y_{l,m}
\end{align*}
with convergence in $L^2(S^2,\sigma)$ for $\mu$-almost every $r$.
The reason we integrate Parseval's identity in $r$ is that the theorem is not only an angular statement; it is a statement in the full space $L^2(\mathbb{R}^3)$. Integrating with respect to the radial measure $\mu$, we obtain
\begin{align*}
\int_0^\infty \|F(r)\|_{L^2(S^2)}^2\,r^2\,d\mathcal{L}^1(r) = \sum_{l=0}^{\infty}\sum_{m=-l}^{l}\int_0^\infty |R_{l,m}(r)|^2 r^2\,d\mathcal{L}^1(r).
\end{align*}
The left-hand side is finite because $F$ lies in the Bochner space $L^2((0,\infty),\mu;L^2(S^2,\sigma))$. Therefore each summand on the right is finite, which proves
\begin{align*}
R_{l,m} \in L^2((0,\infty),r^2\,d\mathcal{L}^1(r)).
\end{align*}
[/guided]
[/step]
[step:Recover convergence in $L^2(\mathbb{R}^3)$]
For $N \in \mathbb{N}_0$, define
\begin{align*}
F_N:(0,\infty) \to L^2(S^2,\sigma)
\end{align*}
by
\begin{align*}
F_N(r) := \sum_{l=0}^{N}\sum_{m=-l}^{l} R_{l,m}(r)Y_{l,m}.
\end{align*}
For $\mu$-almost every $r$, the partial sums $F_N(r)$ converge to $F(r)$ in $L^2(S^2,\sigma)$. Moreover, by Parseval's identity for the finite [orthogonal projection](/theorems/437) onto $\operatorname{span}\{Y_{l,m}:0 \le l \le N,\,-l \le m \le l\}$,
\begin{align*}
\|F(r)-F_N(r)\|_{L^2(S^2)}^2 = \|F(r)\|_{L^2(S^2)}^2 - \sum_{l=0}^{N}\sum_{m=-l}^{l}|R_{l,m}(r)|^2.
\end{align*}
The right-hand side decreases pointwise to $0$ and is bounded above by $\|F(r)\|_{L^2(S^2)}^2$, which is $\mu$-integrable. By dominated convergence,
\begin{align*}
\lim_{N\to\infty}\int_0^\infty \|F(r)-F_N(r)\|_{L^2(S^2)}^2\,r^2\,d\mathcal{L}^1(r) = 0.
\end{align*}
Thus $F_N \to F$ in $L^2((0,\infty),\mu;L^2(S^2,\sigma))$. Applying the inverse unitary $U^{-1}$ gives
\begin{align*}
\psi(r\omega) = \sum_{l=0}^{\infty}\sum_{m=-l}^{l} R_{l,m}(r)Y_{l,m}(\omega)
\end{align*}
with convergence in $L^2(\mathbb{R}^3,\mathcal{L}^3)$.
[/step]
[step:Identify the orthogonal angular momentum sectors]
For each admissible pair $(l,m)$, let
\begin{align*}
\mathcal{K}_{l,m} := \{F:(0,\infty)\to L^2(S^2,\sigma) : F(r)=R(r)Y_{l,m}\text{ for some }R\in L^2((0,\infty),\mu)\}
\end{align*}
inside $L^2((0,\infty),\mu;L^2(S^2,\sigma))$.
If $(l,m)\ne(l',m')$, $R \in L^2((0,\infty),\mu)$, and $S \in L^2((0,\infty),\mu)$, then orthonormality of spherical harmonics gives
\begin{align*}
\int_{S^2} R(r)Y_{l,m}(\omega)\overline{S(r)Y_{l',m'}(\omega)}\,d\sigma(\omega) = R(r)\overline{S(r)}\,0 = 0
\end{align*}
for $\mu$-almost every $r$. Integrating in $r$ yields
\begin{align*}
(RY_{l,m},SY_{l',m'})_{L^2((0,\infty),\mu;L^2(S^2))} = 0.
\end{align*}
Thus the subspaces $\mathcal{K}_{l,m}$ are pairwise orthogonal.
The expansion from the previous step shows that the closed linear span of the $\mathcal{K}_{l,m}$ is all of $L^2((0,\infty),\mu;L^2(S^2,\sigma))$. Therefore
\begin{align*}
L^2((0,\infty),\mu;L^2(S^2,\sigma)) = \bigoplus_{l=0}^{\infty}\bigoplus_{m=-l}^{l}\mathcal{K}_{l,m}.
\end{align*}
Applying $U^{-1}$ identifies $\mathcal{K}_{l,m}$ with $\mathcal{H}_{l,m}$ and gives
\begin{align*}
L^2(\mathbb{R}^3,\mathcal{L}^3) = \bigoplus_{l=0}^{\infty}\bigoplus_{m=-l}^{l}\mathcal{H}_{l,m}.
\end{align*}
[/step]
[step:Prove uniqueness and the norm identity]
Suppose
\begin{align*}
\psi(r\omega)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l} R_{l,m}(r)Y_{l,m}(\omega)
\end{align*}
in $L^2(\mathbb{R}^3,\mathcal{L}^3)$. Applying $U$ and then taking the $L^2(S^2,\sigma)$ inner product with $Y_{l,m}$ for $\mu$-almost every $r$ gives
\begin{align*}
R_{l,m}(r)=\int_{S^2} (U\psi)(r)(\omega)\overline{Y_{l,m}(\omega)}\,d\sigma(\omega).
\end{align*}
Thus the coefficient functions are unique as elements of $L^2((0,\infty),r^2\,d\mathcal{L}^1(r))$.
Finally, since $U$ is unitary and Parseval's identity holds on almost every angular slice,
\begin{align*}
\|\psi\|_{L^2(\mathbb{R}^3)}^2 = \|U\psi\|_{L^2((0,\infty),\mu;L^2(S^2))}^2.
\end{align*}
Using the integrated [Parseval identity](/theorems/248) from the angular expansion step gives
\begin{align*}
\|\psi\|_{L^2(\mathbb{R}^3)}^2 = \sum_{l=0}^{\infty}\sum_{m=-l}^{l}\int_0^\infty |R_{l,m}(r)|^2 r^2\,d\mathcal{L}^1(r).
\end{align*}
This is the asserted Hilbert space direct sum decomposition and the asserted expansion formula.
[/step]