Reduced transition probability $B(i\rightarrow f)$
\begin{align} B(i\rightarrow f) \equiv \sum_{\mu M_2} \left| \left\langle J_f M_f \left| \mathcal{O}(\pi\lambda)_\mu \right| J_i M_i \right\rangle\right|^2. \end{align} By the Wigner-Eckart theorem, this equation becomes \begin{align} B(i\rightarrow f) = \frac{1}{2J_i+1} \left| \left\langle J_f \left|\left| \mathcal{O}(\pi\lambda) \right|\right| J_i \right\rangle\right|^2, \end{align} where $\mathcal{O}(\pi\lambda)$ is a one-body operator and a sum over the operators for the individual nucleons $k$ \begin{align} \mathcal{O}(\pi\lambda) = \sum_k O(\pi\lambda, k). \end{align} The electro transition operator is given by: \begin{align} O(E\lambda) = r^\lambda Y^\lambda_\mu(\hat{r})e_q e, \end{align} where $Y^\lambda_\mu(\hat{\boldsymbol{r}})(=Y^\lambda_\mu(\theta, \varphi))$ is a spherical harmonics, $\hat{\boldsymbol{r}} $ is a unit vector, and $e$ is elementary charge. The $e_q$ are the electric charges for the proton ($q=p$) and neutron ($q=n$) in units of $e$. For example, for free nucleon, $e_p = 1$ for proton and $e_n=0$ for neutron (Brown's lecture note).
From Bohr & Mottelson Vol. I, Eq. (3C-34),
\begin{align} B(E\lambda; j_1 \rightarrow j_2) = e^2 \frac{2\lambda+1}{4\pi} \left\langle j_2 \left| r^\lambda \right| j_1 \right\rangle^2 \left\langle j_1 \frac{1}{2} \lambda 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2. \end{align}
For the radial part of wave function $R_{nl}(r)$ in a spherically symmetric potential
\begin{align} R_{nl}(r)&=N_{nl}r^{l}e^{-\frac{1}{2}{2\nu}r^2}L_{n-1}^{l+\frac{1}{2}}({2\nu}r^2)\\ N_{nl} &= \sqrt{\frac{2(2\nu)^{l+3/2}(n-1)!}{\Gamma(n+l+1/2)}}\\ \nu &= \frac{m\omega}{2\hbar}\\ N &= 2(n-1)+l\\ E_{nl}&=\left(2(n-1)+l+\frac{3}{2}\right)\hbar\omega=\left(N+\frac{3}{2}\right)\hbar\omega\\ n &\in\mathbb{N}, \end{align}
$\langle j_2 \left| r^\lambda \right| j_1 \rangle$ is
\begin{align} \langle j_2 \left| r^\lambda \right| j_1 \rangle =&\ (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{\lambda}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}}\\ &\times\Gamma\left(\frac{l_1+l_2+\lambda+3}{2}\right) \sum_{k=0}^{\min(n_1-1,n_2-1)}\binom{\frac{-l_1+l_2+\lambda}{2}}{n_1-k-1}\binom{\frac{l_1-l_2+\lambda}{2}}{n_2-k-1}\binom{k+\frac{l_1+l_2+\lambda+1}{2}}{k}. \end{align}
Therefore,
\begin{align} B(E\lambda; j_1 \rightarrow j_2) =&\ e^2 \frac{2\lambda+1}{4\pi}\left\langle j_1 \frac{1}{2} \lambda 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2 \left[(-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{\lambda}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}}\right.\\ & \left.\times\Gamma\left(\frac{l_1+l_2+\lambda+3}{2}\right) \sum_{k=0}^{\min(n_1-1,n_2-1)}\binom{\frac{-l_1+l_2+\lambda}{2}}{n_1-k-1}\binom{\frac{l_1-l_2+\lambda}{2}}{n_2-k-1}\binom{k+\frac{l_1+l_2+\lambda+1}{2}}{k}\right]^2\\ =&\ e^2 \frac{2\lambda+1}{4\pi}\left\langle j_1 \frac{1}{2} \lambda 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2 \left(\frac{\hbar}{m\omega}\right)^{\lambda}\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}\\ & \times\left[\Gamma\left(\frac{l_1+l_2+\lambda+3}{2}\right) \sum_{k=0}^{\min(n_1-1,n_2-1)}\binom{\frac{-l_1+l_2+\lambda}{2}}{n_1-k-1}\binom{\frac{l_1-l_2+\lambda}{2}}{n_2-k-1}\binom{k+\frac{l_1+l_2+\lambda+1}{2}}{k}\right]^2\\ \end{align}
If $l_1+\lambda=l_2$,
\begin{align} \langle{j_2}\left|r^\lambda\right|{j_1}\rangle = (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!\Gamma(n_2+l_2+1/2)}{(n_2-1)!\Gamma(n_1+l_1+1/2)}} \binom{l_2-l_1}{n_1-n_2}. \end{align}
Therefore,
\begin{align} B(E\lambda; j_1 \rightarrow j_2) =&\ e^2 \frac{2\lambda+1}{4\pi}\left\langle j_1 \frac{1}{2} \lambda 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2 \left[(-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!\Gamma(n_2+l_2+1/2)}{(n_2-1)!\Gamma(n_1+l_1+1/2)}} \binom{l_2-l_1}{n_1-n_2}\right]^2\\ =&\ e^2 \frac{2\lambda+1}{4\pi}\left\langle j_1 \frac{1}{2} \lambda 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2 \left(\frac{\hbar}{m\omega}\right)^{\lambda}\frac{(n_1-1)!\Gamma(n_2+l_2+1/2)}{(n_2-1)!\Gamma(n_1+l_1+1/2)} \binom{\lambda}{n_1-n_2}^2. \end{align}
If $l_1+\lambda=l_2$ and $\lambda=1$,
\begin{align} B(E\lambda; j_1 \rightarrow j_2) =&\ e^2 \frac{2+1}{4\pi}\left\langle j_1 \frac{1}{2} 1 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2 \frac{\hbar}{m\omega}\frac{(n_1-1)!\Gamma(n_2+l_2+1/2)}{(n_2-1)!\Gamma(n_1+l_1+1/2)} \binom{1}{n_1-n_2}^2\\ =&\ e^2 \frac{3}{4\pi}\left\langle j_1 \frac{1}{2} 1 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2 \frac{\hbar}{m\omega}\frac{(n_1-1)!\Gamma(n_2+l_2+1/2)}{(n_2-1)!\Gamma(n_1+l_1+1/2)}. \end{align}
The radial part of the Schrödinger equation is
\begin{align} & \left[-\frac{\hbar^2}{2m}\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right)+\frac{l(l+1)\hbar^2}{2mr^2}+V(r)\right]R(r)=ER(r) \\ & \Leftrightarrow \left[-\frac{\hbar^2}{2m}\frac{1}{r}\left(\frac{d^2}{dr^2}\right)r+\frac{l(l+1)\hbar^2}{2mr^2}+V(r)\right]R(r)=ER(r)\\ & \Leftrightarrow -\frac{\hbar^2}{2m}\frac{1}{r}\frac{d^2}{dr^2}\left(rR(r)\right)+\frac{l(l+1)\hbar^2}{2mr^2}R_{nl}(r)+V(r)R(r)=ER(r)\\ & \Leftrightarrow -\frac{\hbar^2}{2m}\frac{d^2}{dr^2}\left(rR(r)\right)+\frac{l(l+1)\hbar^2}{2mr^2}\left(rR(r)\right)+V(r)\left(rR(r)\right)=E\left(rR(r)\right). \end{align}
By using $u(r)=rR(r)$, \begin{align} & -\frac{\hbar^2}{2m}\frac{d^2}{dr^2}u(r)+\frac{l(l+1)\hbar^2}{2mr^2}u(r)+V(r)u(r)=Euleft(r)\\ & \Leftrightarrow \left[-\frac{\hbar^2}{2m}\frac{d^2}{dr^2}+\frac{l(l+1)\hbar^2}{2mr^2}+V(r)\right]u(r)=Eu(r). \end{align}
For the harmonic oscillator potential
\begin{align} V(r) = \frac{1}{2}m\omega^2{r^2}, \end{align}
the radial part of the Schrödinger equation becomes
\begin{align} \left[-\frac{\hbar^2}{2m}\frac{d^2}{dr^2}+\frac{l(l+1)\hbar^2}{2mr^2}+\frac{1}{2}m\omega^2r^2\right]u(r)=Eu(r). \end{align}
How should we solve this differential equation? We can find the solution in NIST Digital Library of Mathematical Functions (DLMF). If we think of a very general equation
\begin{align} A(x)f''(x) + B(x)f'(x) + C(x) f(x) + \lambda_n f(x) = 0, \end{align}
the solutions are summarized in 18 Orthogonal Polynomials - Classical Orthogonal Polynomials - §18.8 Differential Equations. In this table, we can choose the case of $A(x) = 1$, $B(x) = 0$, $C(x) = -x^2+\left(\frac{1}{4} -\alpha^2\right) x^{-2}$, $\lambda_n = 4n+2\alpha+2$, namely,
\begin{align} f''(x) + \left(-x^2+\frac{\frac{1}{4} -\alpha^2}{x^2}\right)f(x) + \left(4n+2\alpha+2\right) f(x) = 0. \end{align} the solution is \begin{align} f(x)=e^{-\frac{1}{2}x^2}x^{\alpha+\frac{1}{2}}L_{n}^{(\alpha)}(x^2), \end{align}
where $L_{n}^{(\alpha)}(x)$ is a generalized Laguerre function. (Here this is not a associated generalized Laguerre function.) So in order to compare this differential equation with the Schrödinger equation, the Schrödinger equation should be rearranged. At first, by multiplying the both sides of the equation with $-\frac{2}{\hbar\omega}$, it becomes
\begin{align} &\left[\frac{\hbar}{m\omega}\frac{d^2}{dr^2}-\frac{\hbar}{m\omega}\frac{l(l+1)}{r^2}-\frac{m\omega}{\hbar}r^2\right]u(r)=-\frac{2E}{\hbar\omega}u(r)\\ &\Leftrightarrow \left[\frac{\hbar}{m\omega}\frac{d^2}{dr^2}-\frac{\hbar}{m\omega}\frac{l(l+1)}{r^2}-\frac{m\omega}{\hbar}r^2+\frac{2E}{\hbar\omega}\right]u(r)=0 \end{align}
Using $q = \sqrt{\frac{m\omega}{\hbar}}r$ and $f(q)=u\left(\sqrt{\frac{\hbar}{m\omega}}q\right)=u(r)$, the equation becomes
\begin{align} &\left[\frac{d^2}{dq^2}-\frac{l(l+1)}{q^2}-q^2+\frac{2E}{\hbar\omega}\right]f(q)=0\\ &\Leftrightarrow \frac{d^2}{dq^2}f(q)+\left[-q^2-\frac{l(l+1)}{q^2}\right]f(q)+\frac{2E}{\hbar\omega}f(q)=0.\\ &\Leftrightarrow \frac{d^2}{dq^2}f(q)+\left[-q^2+\frac{\frac{1}{4}-\left(l+\frac{1}{2}\right)^2}{q^2}\right]f(q)+\frac{2E}{\hbar\omega}f(q)=0. \end{align}
Using $\frac{2E}{\hbar\omega}=4(n-1)+2\left(l+\frac{1}{2}\right)+2$, namely, $E=\left(2(n-1)+l+\frac{3}{2}\right)\hbar\omega$,
\begin{align} \frac{d^2}{dq^2}f(q)+\left[-q^2+\frac{\frac{1}{4}-\left(l+\frac{1}{2}\right)^2}{q^2}\right]f(q)+\left[4(n-1)+2\left(l+\frac{1}{2}\right)+2\right]f(q)=0. \end{align}
When $n$ is an integer, this equation has solutions. In other words, when the energy $E$ is equal to specific values $E_{nl}=\left(2(n-1)+l+\frac{3}{2}\right)\hbar\omega$ for $n$ of integers, this equation has solutions. By comparing this equation with the above equation, the solution is
\begin{align} u(r)=f(q)=e^{-\frac{1}{2}q^2}q^{l+1}L_{n-1}^{l+\frac{1}{2}}(q^2). \end{align}
Using $\nu = \frac{m\omega}{2\hbar}$ and $q=\sqrt{2\nu}r$,
\begin{align} u(r)=f\left(\sqrt{2\nu}r\right)=e^{-\frac{1}{2}{2\nu}r^2}(2\nu)^{\frac{l+1}{2}}r^{l+1}L_{n-1}^{l+\frac{1}{2}}({2\nu}r^2). \end{align}
By remembering $u(r)=rR(r)$,
\begin{align} R(r)=u(r)/r=e^{-\frac{1}{2}{2\nu}r^2}(2\nu)^{\frac{l+1}{2}}r^{l}L_{n-1}^{l+\frac{1}{2}}({2\nu}r^2). \end{align}
Using a normalization coefficient $N_{nl}$, we can define the solutions $R_{nl}(r)$ like
\begin{align} R(r)=R_{nl}(r)=N_{nl}r^{l}e^{-\frac{1}{2}{2\nu}r^2}L_{n-1}^{l+\frac{1}{2}}({2\nu}r^2). \end{align}
The normalization condition is
\begin{align} \int_0^\infty r^2|R(r)|^2dr=1. \end{align}
Therefore,
\begin{align} \int_0^{\infty}r^2|R(r)|^2dr &= \int_0^{\infty}r^2\left|N_{nl}r^{l}e^{-\frac{1}{2}{2\nu}r^2}L_{n-1}^{l+\frac{1}{2}}({2\nu}r^2)\right|^2dr\\ &= \left|N_{nl}\right|^2\int_0^{\infty}r^{2l+2}e^{-{2\nu}r^2}\left[L_{n-1}^{l+\frac{1}{2}}({2\nu}r^2)\right]^2dr. \end{align}
Using $t={2\nu}r^2$, $dt=2({2\nu})rdr$, and $\frac{1}{2}(2\nu)^{-1/2}t^{-1/2}dt=dr$ \begin{align} \int_0^{\infty}r^2|R(r)|^2dr &= \left|N_{nl}\right|^2\int_0^{\infty}\left(\frac{t}{2\nu}\right)^{l+1}e^{-t}\left[L_{n-1}^{l+\frac{1}{2}}(t)\right]^2\frac{1}{2}(2\nu)^{-1/2}t^{-1/2}dt\\ &= \frac{\left|N_{nl}\right|^2}{2}(2\nu)^{-l-3/2}\int_0^{\infty}t^{l+1/2}e^{-t}\left[L_{n-1}^{l+\frac{1}{2}}(t)\right]^2dt. \end{align}
From Generalized Laguerre polynomials: Integration,
\begin{align} \int_0^{\infty}t^{\lambda}e^{-t}L_{m}^{(\lambda)}(t)L_{n}^{(\lambda)}(t)dt = \frac{\Gamma(n+\lambda+1)\delta_{mn}}{n!}. \end{align}
By using this formula,
\begin{align} \int_0^{\infty}r^2|R(r)|^2dr = \frac{\left|N_{nl}\right|^2}{2}(2\nu)^{-l-3/2}\frac{\Gamma(n+l+1/2)}{(n-1)!}=1. \end{align}
Therefore,
\begin{align} N_{nl} = \sqrt{\frac{2(2\nu)^{l+3/2}(n-1)!}{\Gamma(n+l+1/2)}}. \end{align}
In summary, the radial part of the Schrödinger equation with a spherically symmetric potential is
\begin{align} \left[-\frac{\hbar^2}{2m}\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d}{dr}\right)+\frac{l(l+1)\hbar^2}{2mr^2}+\frac{1}{2}m\omega^2r^2\right]R(r)=ER(r). \end{align}
The solutions are
\begin{align} R(r)&=R_{nl}(r)=N_{nl}r^{l}e^{-\frac{1}{2}{2\nu}r^2}L_{n-1}^{l+\frac{1}{2}}({2\nu}r^2)\\ N_{nl} &= \sqrt{\frac{2(2\nu)^{l+3/2}(n-1)!}{\Gamma(n+l+1/2)}}\\ \nu &= \frac{m\omega}{2\hbar}\\ N &= 2(n-1)+l\\ E &= E_{nl}=\left(2(n-1)+l+\frac{3}{2}\right)\hbar\omega=\left(N+\frac{3}{2}\right)\hbar\omega\\ n &\in\mathbb{N}. \end{align}
By using the result, one can calculate the matrix element
\begin{align} \langle j_2 \left| r^\lambda \right| j_1 \rangle & = \int_0^{\infty}r^{\lambda}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r)r^2dr\\ & = \int_0^{\infty}r^{\lambda+2}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r)dr \end{align} In case of the spherically symmetric potential, $\mathscr{R}_{nlj}(r)=R_{nl}(r)$, and \begin{align} \langle j_2 \left| r^\lambda \right| j_1 \rangle = \int_0^{\infty}r^{\lambda+2}R_{n_2l_2}(r)R_{n_1l_1}(r)dr. \end{align}
By substituting the above solutions to this formula,
\begin{align} \langle j_2 \left| r^\lambda \right| j_1 \rangle &= \int_0^{\infty}r^{\lambda+2}\left[N_{n_2l_2}r^{l_2}e^{-\frac{1}{2}{2\nu}r^2}L_{n_2-1}^{l_2+\frac{1}{2}}({2\nu}r^2)\right]\left[N_{n_1l_1}r^{l_1}e^{-\frac{1}{2}{2\nu}r^2}L_{n_1-1}^{l_1+\frac{1}{2}}({2\nu}r^2)\right]dr\\ &= N_{n_1l_1}N_{n_2l_2}\int_0^{\infty}r^{l_1+l_2+\lambda+2}e^{-{2\nu}r^2}L_{n_1-1}^{l_1+\frac{1}{2}}({2\nu}r^2)L_{n_2-1}^{l_2+\frac{1}{2}}({2\nu}r^2)dr\\ \end{align}
Using $t={2\nu}r^2$ (i.e., $r=(t/({2\nu}))^{1/2}$), $dt=2({2\nu})rdr$, and $\frac{1}{2}(2\nu)^{-1/2}t^{-1/2}dt=dr$,
\begin{align} \langle j_2 \left| r^\lambda \right| j_1 \rangle &= N_{n_1l_1}N_{n_2l_2}\int_0^{\infty}\left(\frac{t}{2\nu}\right)^{\frac{l_1+l_2+\lambda}{2}+1}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)\frac{1}{2}(2\nu)^{-1/2}t^{-1/2}dt\\ &= \frac{1}{2}N_{n_1l_1}N_{n_2l_2}({2\nu})^{-\frac{l_1+l_2+\lambda+3}{2}}\int_0^{\infty}t^{\frac{l_1+l_2+\lambda+1}{2}}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)dt\\ &= \frac{1}{2}\sqrt{\frac{2(2\nu)^{l_1+3/2}(n_1-1)!}{\Gamma(n_1+l_1+1/2)}}\sqrt{\frac{2(2\nu)^{l_2+3/2}(n_2-1)!}{\Gamma(n_2+l_2+1/2)}}({2\nu})^{-\frac{l_1+l_2+\lambda+3}{2}} \int_0^{\infty}t^{\frac{l_1+l_2+\lambda+1}{2}}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)dt\\ &= \sqrt{\frac{(n_1-1)!(n_2-1)!}{({2\nu})^\lambda\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \int_0^{\infty}t^{\frac{l_1+l_2+\lambda+1}{2}}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)dt\\ &= \left(\frac{\hbar}{m\omega}\right)^{\frac{\lambda}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \int_0^{\infty}t^{\frac{l_1+l_2+\lambda+1}{2}}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)dt. \end{align}
From Eq. (10) of https://doi.org/10.1016/S0893-9659(03)90106-6,
\begin{align} &\int_0^{\infty}x^{\mu}e^{-x}L_{m}^{\alpha}(x)L_{n}^{\beta}(x)dx\\ &=(-1)^{m+n}\Gamma(\mu+1)\sum_{k=0}^{\min(m,n)}\binom{\mu-\alpha}{m-k}\binom{\mu-\beta}{n-k}\binom{k+\mu}{k}\\ &(\mathfrak{R}(\mu)>-1;\ m,n\in\mathbb{N}_0). \end{align}
Using this formula,
\begin{align} &\int_0^{\infty}t^{\frac{l_1+l_2+\lambda+1}{2}}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)dt\\ &=(-1)^{(n_1-1)+(n_2-1)}\Gamma\left(\frac{l_1+l_2+\lambda+1}{2}+1\right) \sum_{k=0}^{\min(n_1-1,n_2-1)}\binom{\frac{l_1+l_2+\lambda+1}{2}-\left(l_1+\frac{1}{2}\right)}{(n_1-1)-k}\binom{\frac{l_1+l_2+\lambda+1}{2}-\left(l_2+\frac{1}{2}\right)}{(n_2-1)-k}\binom{k+\frac{l_1+l_2+\lambda+1}{2}}{k}\\ &=(-1)^{n_1+n_2}\Gamma\left(\frac{l_1+l_2+\lambda+3}{2}\right) \sum_{k=0}^{\min(n_1-1,n_2-1)}\binom{\frac{-l_1+l_2+\lambda}{2}}{n_1-k-1}\binom{\frac{l_1-l_2+\lambda}{2}}{n_2-k-1}\binom{k+\frac{l_1+l_2+\lambda+1}{2}}{k}. \end{align}
Therefore,
\begin{align} \langle j_2 \left| r^\lambda \right| j_1 \rangle &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{\lambda}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(\frac{l_1+l_2+\lambda+3}{2}\right) \sum_{k=0}^{\min(n_1-1,n_2-1)}\binom{\frac{-l_1+l_2+\lambda}{2}}{n_1-k-1}\binom{\frac{l_1-l_2+\lambda}{2}}{n_2-k-1}\binom{k+\frac{l_1+l_2+\lambda+1}{2}}{k}\\\end{align}
For wave functions $|{j}\rangle=|nl\rangle$ in a potential of a 3D isotropic harmonic oscillator, the matrix element $\langle{n_2l_2}\left|r^\lambda\right|{n_1l_1}\rangle$ is
\begin{align} \langle{j_2}\left|r^\lambda\right|{j_1}\rangle &=\langle{n_2l_2}\left|r^\lambda\right|{n_1l_1}\rangle\\ &= \left(\frac{\hbar}{m\omega}\right)^{\frac{\lambda}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \int_0^{\infty}t^{\frac{l_1+l_2+\lambda+1}{2}}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)dt\\ &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{\lambda}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(\frac{l_1+l_2+\lambda+3}{2}\right) \sum_{k=0}^{\min(n_1-1,n_2-1)}\binom{\frac{-l_1+l_2+\lambda}{2}}{n_1-k-1}\binom{\frac{l_1-l_2+\lambda}{2}}{n_2-k-1}\binom{k+\frac{l_1+l_2+\lambda+1}{2}}{k}. \end{align}
\begin{align} \langle{nl}\left|r^\lambda\right|{nl}\rangle &= \left(\frac{\hbar}{m\omega}\right)^{\frac{\lambda}{2}}\frac{(n-1)!}{\Gamma(n+l+1/2)} \int_0^{\infty}t^{l+\frac{\lambda+1}{2}}e^{-t}\left[L_{n-1}^{l+\frac{1}{2}}(t)\right]^2dt\\ &= \left(\frac{\hbar}{m\omega}\right)^{\frac{\lambda}{2}}\frac{(n-1)!}{\Gamma(n+l+1/2)} \Gamma\left(l+\frac{\lambda+3}{2}\right)\sum_{k=0}^{n-1}\binom{\frac{\lambda}{2}}{n-k-1}^2\binom{k+l+\frac{\lambda+1}{2}}{k}. \end{align}
\begin{align} \langle{nl}\left|r^2\right|{nl}\rangle &= \frac{\hbar}{m\omega}\frac{(n-1)!}{\Gamma(n+l+1/2)} \int_0^{\infty}t^{l+\frac{3}{2}}e^{-t}\left[L_{n-1}^{l+\frac{1}{2}}(t)\right]^2dt\\ &= \frac{\hbar}{m\omega}\frac{(n-1)!}{\Gamma(n+l+1/2)} \Gamma\left(l+\frac{5}{2}\right)\sum_{k=0}^{n-1}\binom{1}{n-k-1}^2\binom{k+l+\frac{3}{2}}{k}. \end{align}
For $n=1$,
\begin{align} \langle{nl}\left|r^2\right|{nl}\rangle &= \frac{\hbar}{m\omega}\frac{0!}{\Gamma(1+l+1/2)} \Gamma\left(l+\frac{5}{2}\right)\sum_{k=0}^{0}\binom{1}{1-k-1}^2\binom{k+l+\frac{3}{2}}{k}\\ &= \frac{\hbar}{m\omega}\frac{\Gamma(l+5/2)}{\Gamma(l+3/2)}\binom{1}{0}^2\binom{0+l+\frac{3}{2}}{0}\\ &= \frac{\hbar}{m\omega}\left(l+\frac{3}{2}\right) \end{align}
In case of $n{\gt}1$, for $k=0, 1, \ldots, n-1$, $\binom{1}{n-k-1}$ has non-zero value only for $k=n-2$ and $k=n-1$. Therefore,
\begin{align} \langle{nl}\left|r^2\right|{nl}\rangle &= \frac{\hbar}{m\omega}\frac{(n-1)!}{\Gamma(n+l+1/2)} \Gamma\left(l+\frac{5}{2}\right) \sum_{k=n-2}^{n-1}\binom{1}{n-k-1}^2\binom{k+l+\frac{3}{2}}{k}\\ &= \frac{\hbar}{m\omega}\frac{(n-1)!}{\Gamma(n+l+1/2)} \Gamma\left(l+\frac{5}{2}\right) \left[\binom{1}{n-(n-2)-1}^2\binom{(n-2)+l+\frac{3}{2}}{n-2}+\binom{1}{n-(n-1)-1}^2\binom{(n-1)+l+\frac{3}{2}}{n-1}\right]\\ &= \frac{\hbar}{m\omega}\frac{(n-1)!}{\Gamma(n+l+1/2)} \Gamma\left(l+\frac{5}{2}\right) \left[\binom{1}{1}^2\binom{n+l-\frac{1}{2}}{n-2}+\binom{1}{0}^2\binom{n+l+\frac{1}{2}}{n-1}\right]\\ &= \frac{\hbar}{m\omega}\frac{(n-1)!}{\Gamma(n+l+1/2)} \Gamma\left(l+\frac{5}{2}\right) \left[\binom{n+l-\frac{1}{2}}{n-2}+\binom{n+l+\frac{1}{2}}{n-1}\right]\\ &= \frac{\hbar}{m\omega}\frac{(n-1)!}{\Gamma(n+l+1/2)} \Gamma\left(l+\frac{5}{2}\right) \left[\frac{\Gamma(n+l+1/2)}{\Gamma(n-1)\Gamma(l+5/2)}+\frac{\Gamma(n+l+3/2)}{\Gamma(n)\Gamma(l+5/2)}\right]\\ &= \frac{\hbar}{m\omega}\frac{\Gamma(n)}{\Gamma(n+l+1/2)} \Gamma\left(l+\frac{5}{2}\right) \left[\frac{(n-1)\Gamma(n+l+1/2)}{\Gamma(n)\Gamma(l+5/2)}+\frac{(n+l+1/2)\Gamma(n+l+1/2)}{\Gamma(n)\Gamma(l+5/2)}\right]\\ &= \frac{\hbar}{m\omega} \left[(n-1)+(n+l+1/2)\right]\\ &= \frac{\hbar}{m\omega}\left[2(n-1)+l+\frac{3}{2}\right]. \end{align}
Therefore, for $n\in\mathbb{N}$,
\begin{align} \langle{nl}\left|r^2\right|{nl}\rangle &= \frac{\hbar}{m\omega}\left[2(n-1)+l+\frac{3}{2}\right]\\ &= \frac{\hbar}{m\omega}\left(N+\frac{3}{2}\right). \end{align}
\begin{align} &\langle{n_2l_2}\left|r^\lambda\right|{n_1l_1}\rangle\\ &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(l_2+\frac{3}{2}\right) \sum_{k=0}^{\min(n_1-1,n_2-1)}\binom{l_2-l_1}{n_1-k-1}\binom{0}{n_2-k-1}\binom{k+l_2+\frac{1}{2}}{k}. \end{align}
$\binom{0}{n_2-k-1}$ has a non-zero value only for $k=n_2-1$. Therefore, in case of $n_1 < n_2$, $\binom{0}{n_2-k-1}$ is equal to zero for $k=0, 1, \ldots, n_1-1$, and thus
\begin{align} &\langle{n_2l_2}\left|r^\lambda\right|{n_1l_1}\rangle\\ &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(l_2+\frac{3}{2}\right) \sum_{k=0}^{n_1-1}\binom{l_2-l_1}{n_1-k-1}\binom{0}{n_2-k-1}\binom{k+l_2+\frac{1}{2}}{k}\\ &=0. \end{align}
For $n_1 \ge n_2$, $\binom{0}{n_2-k-1}$ is equal to zero for $k=n_2-1$, and thus
\begin{align} &\langle{n_2l_2}\left|r^\lambda\right|{n_1l_1}\rangle\\ &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(l_2+\frac{3}{2}\right) \sum_{k=0}^{n_2-1}\binom{l_2-l_1}{n_1-k-1}\binom{0}{n_2-k-1}\binom{k+l_2+\frac{1}{2}}{k}\\ &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(l_2+\frac{3}{2}\right) \binom{l_2-l_1}{n_1-(n_2-1)-1}\binom{0}{n_2-(n_2-1)-1}\binom{(n_2-1)+l_2+\frac{1}{2}}{n_2-1}\\ &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(l_2+\frac{3}{2}\right) \binom{l_2-l_1}{n_1-n_2}\binom{0}{0}\binom{n_2+l_2-\frac{1}{2}}{n_2-1}\\ &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(l_2+\frac{3}{2}\right) \binom{l_2-l_1}{n_1-n_2}\frac{\Gamma(n_2+l_2+1/2)}{\Gamma(n_2)\Gamma(l_2+3/2)}\\ &= (-1)^{n_1+n_2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!\Gamma(n_2+l_2+1/2)}{(n_2-1)!\Gamma(n_1+l_1+1/2)}} \binom{l_2-l_1}{n_1-n_2}. \end{align}
Another derivation is the follows.
\begin{align} \langle{n_2l_2}\left|r^\lambda\right|{n_1l_1}\rangle = \left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \int_0^{\infty}t^{l_2+\frac{1}{2}}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)dt. \end{align}
From Eq. (18) of https://doi.org/10.1016/S0893-9659(03)90106-6,
\begin{align} &\int_0^{\infty}x^{\beta}e^{-{\sigma}x}L_{m}^{\alpha}(x)L_{n}^{\beta}(x)dx =\frac{\Gamma(\beta+1)}{\sigma^{\beta+1}}\binom{n+\beta}{n}\frac{(\alpha-\beta)_{m-n}}{(m-n)!}\\ &(\mathfrak{R}(\beta)>-1;\ \mathfrak{R}(\sigma)>0;\ m{\ge}n{\ge}0\ (m,n\in\mathbb{N}_0)). \end{align}
From http://functions.wolfram.com/GammaBetaErf/Pochhammer/introductions/FactorialBinomials/ShowAll.html,
\begin{align} \binom{n}{k}=\frac{(-1)^k(-n)_k}{k!}. \end{align} (Wikipedia - Pochhammer symbol (in Japanese) seems to be wrong.)
Therefore,
\begin{align} &\int_0^{\infty}x^{\beta}e^{-{\sigma}x}L_{m}^{\alpha}(x)L_{n}^{\beta}(x)dx =(-1)^{m+n}\frac{\Gamma(\beta+1)}{\sigma^{\beta+1}}\binom{n+\beta}{n}\binom{\alpha-\beta}{m-n}\\ &(\mathfrak{R}(\beta)>-1;\ \mathfrak{R}(\sigma)>0;\ m{\ge}n{\ge}0\ (m,n\in\mathbb{N}_0)). \end{align}
Using this formula,
\begin{align} \int_0^{\infty}t^{l_2+\frac{1}{2}}e^{-t}L_{n_1-1}^{l_1+\frac{1}{2}}(t)L_{n_2-1}^{l_2+\frac{1}{2}}(t)dt &=(-1)^{n1+n2}\Gamma(l_2+3/2)\binom{n_2+l_2-1/2}{n_2-1}\binom{l_2-l_1}{n_1-n_2}, \end{align}
and
\begin{align} &\langle{n_2l_2}\left|r^\lambda\right|{n_1l_1}\rangle\\ &= (-1)^{n1+n2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!(n_2-1)!}{\Gamma(n_1+l_1+1/2)\Gamma(n_2+l_2+1/2)}} \Gamma\left(l_2+\frac{3}{2}\right) \binom{n_2+l_2-\frac{1}{2}}{n_2-1}\binom{l_2-l_1}{n_1-n_2}\\ &= (-1)^{n1+n2}\left(\frac{\hbar}{m\omega}\right)^{\frac{l_2-l_1}{2}}\sqrt{\frac{(n_1-1)!\Gamma(n_2+l_2+1/2)}{(n_2-1)!\Gamma(n_1+l_1+1/2)}} \binom{l_2-l_1}{n_1-n_2}. \end{align}
From wolfram alpha,
n1 = 2, n2 = 1, l1 = 0, l2=1, lambda = 1, Integrate[x^((l1+l2+lamda+1)/2)*Exp[-x]*LaguerreL[n1-1,l1+1/2,x]*LaguerreL[n2-1,l2+1/2,x],{x,0,Infinity}]
n1 = 2, n2 = 1, l1 = 0, l2=1, Gamma[l2+3/2]*Binomial[n2+l2-1/2,n2-1]*Binomial[l2-l1,n1-n2]
From Generalized Laguerre polynomials: Integration,
\begin{align} \int_0^{\infty}t^{\alpha-1}e^{-pt}L_{m}^{\lambda}(pt)L_{n}^{\beta}(pt)dt = \frac{p^{-\alpha}\Gamma(\alpha)\Gamma(n-\alpha+\beta+1)\Gamma(m+\lambda+1)}{m!n!\Gamma(1-\alpha+\beta)\Gamma(\lambda+1)}{}_3F_2(-m,\alpha,\alpha-\beta;-n+\alpha-\beta,\lambda+1;1)\\ \Rightarrow \int_0^{\infty}t^{\alpha-1}e^{-t}L_{m}^{\lambda}(t)L_{n}^{\beta}(t)dt = \frac{\Gamma(\alpha)\Gamma(n-\alpha+\beta+1)\Gamma(m+\lambda+1)}{m!n!\Gamma(1-\alpha+\beta)\Gamma(\lambda+1)}{}_3F_2(-m,\alpha,\alpha-\beta;-n+\alpha-\beta,\lambda+1;1). \end{align}
Based on SphericalHarmonicY - Integration, \begin{align} \int_0^{\pi}\int_0^{2\pi} \sin(\vartheta) Y_{n_1}^{m_1}(\vartheta,\phi) Y_{n_2}^{m_2}(\vartheta,\phi) \overline{Y_{n_3}^{m_3}(\vartheta,\varphi)} d\varphi d\vartheta = \sqrt{\frac{(2n_1+1)(2n_2+1)}{4\pi(2n_3+1)}} \langle n_1n_200|n_1n_2n_3 0 \rangle \langle n_1n_2m_1m_2|n_1n_2n_3m_3 \rangle. \end{align}
Therefore,
\begin{align} \left\langle \ell m \left| Y_\lambda^\mu(\theta,\phi) \right| \ell_0 m_0 \right\rangle &= \int_0^{2\pi}\int_0^{\pi} Y_\ell^{m*}(\theta,\phi) Y_\lambda^\mu(\theta,\phi) Y_{\ell_0}^{m_0}(\theta,\phi) \sin\theta d\theta d\phi\\ &= \sqrt{\frac{(2\ell_0+1)(2\lambda+1)}{4\pi(2\ell+1)}} \langle \ell_0\lambda 00|\ell_0\lambda \ell 0 \rangle \langle \ell_0\lambda m_0 \mu|\ell_0\lambda \ell m \rangle. \end{align}
| Reference | Spherical harmonics | Clebsch-Gordan coefficients |
| NIST Digital Library of Mathematical Functions (DLMF) | $Y_{\ell,m}(\theta,\phi)$ | $(j_1m_1j_2m_2|j_1j_2j_3m_3)$ |
| Wolfram Research (functions.wolfram.com) | $Y_n^m(\vartheta,\varphi)$ | $\langle n_1n_2m_1m_2|n_1n_2n_3m_3 \rangle$ |
| Bohr & Mottelson | $Y_{L M}$ | $\langle j_1m_1j_2m_2|JM\rangle$ |
| I. Hamamoto | $Y_{\ell m}(\hat{r})$ | $C(j_1, j_2, J; m_1, m_2, M)$ |
| Particle Data Group | $Y_\ell^m$ | $\langle j_1j_2m_1m_2|j_1j_2JM \rangle$ |
| Wikipedia on Mar. 4, 2020 | $Y_\ell^m(\theta,\varphi)$ | $\langle j_1m_1j_2m_2|JM\rangle$ |
By the Wigner–Eckart theorem, \begin{align} \left\langle \ell m \left| Y_\lambda^\mu(\theta,\phi) \right| \ell_0 m_0 \right\rangle = \langle \ell_0\lambda m_0 \mu|\ell_0\lambda \ell m \rangle \frac{\left\langle \ell \left|\left| Y_\lambda(\theta,\phi) \right|\right| \ell_0 \right\rangle}{\sqrt{2\ell+1}}. \end{align} By using \begin{align} \left\langle \ell m \left| Y_\lambda^\mu(\theta,\phi) \right| \ell_0 m_0 \right\rangle = \sqrt{\frac{(2\ell_0+1)(2\lambda+1)}{4\pi(2\ell+1)}} \langle \ell_0\lambda 00|\ell_0\lambda \ell 0 \rangle \langle \ell_0\lambda m_0 \mu|\ell_0\lambda \ell m \rangle, \end{align}
one can get
\begin{align} \left\langle \ell \left|\left| Y_\lambda(\theta,\phi) \right|\right| \ell_0 \right\rangle = \sqrt{\frac{(2\ell_0+1)(2\lambda+1)}{4\pi}} \langle \ell_0\lambda 00|\ell_0\lambda \ell 0 \rangle. \end{align}
Goal
From Bohr & Mottelson Vol.1 (3C-17),
\begin{align} &B(E(M)\lambda;I_1{\rightarrow}I_1)\equiv\\ & \sum_{\mu M_2} \left| \left\langle{I_1}{M_1}\left|\mathscr{M}(E(M)\lambda,\mu)\right|{I_2}{M_2}\right\rangle\right|^2 = \frac{1}{2I_1+1} \left| \left\langle{I_1}\left|\left|\mathscr{M}(E(M)\lambda)\right|\right|{I_2}\right\rangle\right|^2. \end{align}
Derivation (including mistakes! Under construction!)
By the Wigner–Eckart theorem, \begin{align} \left\langle{I_1}{M_1}\left|\mathscr{M}(E(M)\lambda,\mu)\right|{I_2}{M_2}\right\rangle = \left\langle{I_1}{M_1}\lambda \mu|{I_2}{M_2}\right\rangle \frac{\left\langle{I_1}\left|\left|\mathscr{M}(E(M)\lambda)\right|\right|{I_2}\right\rangle}{\sqrt{2I_2+1}}, \end{align} one can get \begin{align} \sum_{\mu M_2}\left|\left\langle{I_1}{M_1}\left|\mathscr{M}(E(M)\lambda,\mu)\right|{I_2}{M_2}\right\rangle\right|^2 &= \sum_{\mu{M_2}}\left|\left\langle{I_1}{M_1}\lambda \mu|{I_2}{M_2}\right\rangle \frac{\left\langle{I_1}\left|\left|\mathscr{M}(E(M)\lambda)\right|\right|{I_2}\right\rangle}{\sqrt{2I_2+1}}\right|^2\\ &= \frac{1}{2I_2+1}\left|\left\langle{I_1}\left|\left|\mathscr{M}(E(M)\lambda)\right|\right|{I_2}\right\rangle\right|^2 \sum_{\mu{M_2}}\left|\left\langle{I_1}{M_1}\lambda\mu|{I_2}{M_2}\right\rangle\right|^2. \end{align} By the formula of summation of two Clebsh-Gordan coefficients (Clebsch-Gordan coefficients: Summation - Wolfram Functions), \begin{align} \sum_{m_1=-j_1}^{j_1}\sum_{m=-j}^{j} \left\langle{j_1}{j_2}{m_1}{m_2}|{j_1}{j_2}{j}{m}\right\rangle \left\langle{j_1}{j_2'}{m_1}{m_2'}|{j_1}{j_2'}{j}{m}\right\rangle &=\frac{2j+1}{2j_2+1}\delta_{j_2j_2'}\delta_{m_2m_2'}. \end{align} From Clebsch-Gordan Coefficients, Spherical Harmonics, and d Functions, \begin{align} \left\langle{j_1}{j_2}{m_1}{m_2}|{j_1}{j_2}JM\right\rangle = (-1)^{J-j_1-j_2}\left\langle{j_2}{j_1}{m_2}{m_1}|{j_2}{j_1}JM\right\rangle, \end{align}
one can get
\begin{align} \sum_{m_1=-j_1}^{j_1}\sum_{m=-j}^{j} \left\langle{j_1}{j_2}{m_1}{m_2}|{j_1}{j_2}{j}{m}\right\rangle \left\langle{j_1}{j_2'}{m_1}{m_2'}|{j_1}{j_2'}{j}{m}\right\rangle &=\sum_{m_1=-j_1}^{j_1}\sum_{m=-j}^{j} \left[(-1)^{j-j_1-j_2}\left\langle{j_2}{j_1}{m_2}{m_1}|{j_2}{j_1}jm\right\rangle\right]\left[(-1)^{j-j_1-j_2'}\left\langle{j_2'}{j_1}{m_2}{m_1}|{j_2'}{j_1}jm\right\rangle\right]\\ &=(-1)^{2j-2j_1-j_2-j_2'}\sum_{m_1=-j_1}^{j_1}\sum_{m=-j}^{j} \left\langle{j_2}{j_1}{m_2}{m_1}|{j_2}{j_1}jm\right\rangle \left\langle{j_2'}{j_1}{m_2}{m_1}|{j_2'}{j_1}jm\right\rangle\\ &=\frac{2j+1}{2j_2+1}\delta_{j_2j_2'}\delta_{m_2m_2'}. \end{align}
Therefore,
\begin{align} \sum_{m_1=-j_1}^{j_1}\sum_{m=-j}^{j} \left\langle{j_2}{j_1}{m_2}{m_1}|{j_2}{j_1}jm\right\rangle \left\langle{j_2'}{j_1}{m_2}{m_1}|{j_2'}{j_1}jm\right\rangle &=(-1)^{2(j_1-j)}\frac{2j+1}{2j_2+1}\delta_{j_2j_2'}\delta_{m_2m_2'}. \end{align}
Since $j_1-j$ is a integer,
\begin{align} \sum_{m_1=-j_1}^{j_1}\sum_{m=-j}^{j} \left\langle{j_2}{j_1}{m_2}{m_1}|{j_2}{j_1}jm\right\rangle \left\langle{j_2'}{j_1}{m_2}{m_1}|{j_2'}{j_1}jm\right\rangle &=\frac{2j+1}{2j_2+1}\delta_{j_2j_2'}\delta_{m_2m_2'}. \end{align}
By using this formula,
\begin{align} \sum_{\mu{M_2}}\left|\left\langle{I_1}{M_1}\lambda \mu|{I_2}{M_2}\right\rangle\right|^2 = \frac{2I_2+1}{2I_1+1}, \end{align} where we replaced $j_1\rightarrow\lambda$, $m_1\rightarrow\mu$, $j_2\rightarrow{I_1}$, $m_2\rightarrow{M_1}$, $j\rightarrow{I_2}$, $m\rightarrow{M_2}$. Therefore, \begin{align} \sum_{\mu M_2} \left| \left\langle{I_1}{M_1}\left|\mathscr{M}(E(M)\lambda,\mu)\right|{I_2}{M_2}\right\rangle\right|^2 = \frac{1}{2I_1+1} \left| \left\langle{I_1}\left|\left|\mathscr{M}(E(M)\lambda)\right|\right|{I_2}\right\rangle\right|^2. \end{align}
Goal
Bohr & Mottelson Vol. I, Eq. (3C-34)
\begin{align} B(E\lambda; j_1 \rightarrow j_2) = e^2 \frac{2\lambda+1}{4\pi} \left\langle j_2 \left| r^\lambda \right| j_1 \right\rangle^2 \left\langle j_1 \frac{1}{2} \lambda 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2 \end{align}
Derivation
The single particle wave function [Bohr & Mottelson Vol. I, Eq. (3A-1)] is
\begin{align} \psi_{nljm} =\mathscr{R}_{nlj}(r) \sum_{m_l,m_s} \left\langle \left.lm_l\frac{1}{2}m_s \right| jm \right\rangle i^lY_l^{m_l}(\theta,\phi)\chi_{m_s}. \end{align}
Therefore,
\begin{align} &\left\langle\left.j_2m_2\right|i^\lambda\mathscr{M}(E\lambda,\mu)\left|j_1m_1\right.\right\rangle\\ &=\left\langle\left.n_2l_2j_2m_2\right|i^{\lambda}er^{\lambda}Y_\lambda^\mu(\theta,\phi)\left|n_1l_1j_1m_1\right.\right\rangle\\ &= e\int_0^{2\pi} \int_0^\pi \int_0^\infty \psi_{n_2l_2j_2m_2}^* i^{\lambda}r^{\lambda}Y_\lambda^\mu(\theta,\phi) \psi_{n_1l_1j_1m_1} r^2\sin\theta drd\theta d\phi\\ &= e\int_0^{2\pi} \int_0^\pi \int_0^\infty\left[\mathscr{R}_{n_2l_2j_2}(r) \sum_{m_{l2},m_{s2}} \left\langle \left.l_2m_{l2}\frac{1}{2}m_{s2} \right| j_2m_2 \right\rangle i^{l_2}Y_{l_2}^{m_{l2}}(\theta,\phi)\chi_{m_{s2}}\right]^* i^{\lambda}r^{\lambda}Y_\lambda^\mu(\theta,\phi) \left[\mathscr{R}_{n_1l_1j_1}(r)\sum_{m_{l1},m_{s1}} \left\langle \left.l_1m_{l1}\frac{1}{2}m_{s1} \right| j_1m_1 \right\rangle i^{l_1}Y_{l_1}^{m_{l1}}(\theta,\phi)\chi_{m_{s1}}\right] r^2\sin\theta drd\theta d\phi\\ &= e\int_0^{2\pi} \int_0^\pi \int_0^\infty\left[\mathscr{R}_{n_2l_2j_2}^*(r) \sum_{m_{l2},m_{s2}} \left\langle \left.l_2m_{l2}\frac{1}{2}m_{s2} \right| j_2m_2 \right\rangle i^{-l_2}Y_{l_2}^{m_{l2}*}(\theta,\phi)\chi_{m_{s2}}^*\right] i^{\lambda}r^{\lambda}Y_\lambda^\mu(\theta,\phi) \left[\mathscr{R}_{n_1l_1j_1}(r)\sum_{m_{l1},m_{s1}} \left\langle \left.l_1m_{l1}\frac{1}{2}m_{s1} \right| j_1m_1 \right\rangle i^{l_1}Y_{l_1}^{m_{l1}}(\theta,\phi)\chi_{m_{s1}}\right] r^2\sin\theta drd\theta d\phi\\ &= ei^{-l_2}i^{\lambda}i^{l_1} \int_0^\infty r^{\lambda}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r) r^2dr \sum_{m_{l1},m_{s1},m_{l2},m_{s2}} \left\langle\left.{l_1}m_{l1}\frac{1}{2}m_{s1}\right|j_1m_1\right\rangle \left\langle \left.l_2m_{l2}\frac{1}{2}m_{s2}\right|j_2m_2\right\rangle \chi_{m_{s2}}^*\chi_{m_{s1}} \int_0^{2\pi} \int_0^\pi Y_{l_2}^{m_{l2}*}(\theta,\phi)Y_\lambda^\mu(\theta,\phi)Y_{l1}^{m_{l1}}(\theta,\phi)\sin\theta d\theta d\phi\\ &= ei^{l_1-l_2+\lambda}\int_0^\infty r^{\lambda}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r)r^2dr \sum_{m_{l1},m_{s},m_{l2}} \left\langle \left.l_1m_{l1}\frac{1}{2}m_{s}\right|j_1m_1\right\rangle \left\langle\left.l_2m_{l2}\frac{1}{2}m_{s}\right|j_2m_2\right\rangle \int_0^{2\pi}\int_0^{\pi}Y_{l_2}^{m_{l2}*}(\theta,\phi) Y_\lambda^\mu(\theta,\phi)Y_{l_1}^{m_{l1}}(\theta,\phi)\sin{\theta}d{\theta}d\phi. \end{align}
By using the following formula (SphericalHarmonicY - Integration) \begin{align} \int_0^{\pi}\int_0^{2\pi} \sin(\vartheta) Y_{n_1}^{m_1}(\vartheta,\phi) Y_{n_2}^{m_2}(\vartheta,\phi) \overline{Y_{n_3}^{m_3}(\vartheta,\varphi)} d\varphi d\vartheta = \sqrt{\frac{(2n_1+1)(2n_2+1)}{4\pi(2n_3+1)}} \langle n_1n_200|n_1n_2n_3 0 \rangle \langle n_1n_2m_1m_2|n_1n_2n_3m_3 \rangle, \end{align}
one can get
\begin{align} &\left\langle\left.n_2l_2j_2m_2\right|i^{\lambda}er^{\lambda}Y_\lambda^\mu(\theta,\phi)\left|n_1l_1j_1m_1\right.\right\rangle\\ &= ei^{l_1-l_2+\lambda}\int_0^\infty r^{\lambda}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r)r^2dr \sum_{m_{l1},m_{s},m_{l2}} \left\langle \left.l_1m_{l1}\frac{1}{2}m_{s}\right|j_1m_1\right\rangle \left\langle \left.l_2m_{l2}\frac{1}{2}m_{s}\right|j_2m_2\right\rangle \sqrt{\frac{(2l_1+1)(2{\lambda}+1)}{4\pi(2l_2+1)}}{\langle}l_10{\lambda}0|l_20{\rangle}{\langle}l_1m_{l1}{\lambda}{\mu}|l_2m_{l2}{\rangle}\\ &= ei^{l_1-l_2+\lambda}\sqrt{\frac{(2l_1+1)(2{\lambda}+1)}{4\pi(2l_2+1)}}{\langle}l_10{\lambda}0|l_20{\rangle} \int_0^{\infty}r^{\lambda}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r)r^2dr \sum_{m_{l1},m_{s},m_{l2}} \left\langle\left.l_1m_{l1}\frac{1}{2}m_{s}\right|j_1m_1\right\rangle \left\langle\left.l_2m_{l2}\frac{1}{2}m_{s}\right|j_2m_2\right\rangle {\langle}l_1m_{l1}{\lambda}{\mu}|l_2m_{l2}{\rangle}\\ \end{align}
By using the following formula (ClebschGordan - Summation - Involving three Clebsch Gordan coefficients)
\begin{align} &\sum_{m_1=-j_1}^{j_1}\sum_{m_2=-j_2}^{j_2}\sum_{m_6=-j_6}^{j_6} \left\langle\left.j_1j_2m_1m_2\right|j_1j_2j_3m_3\right\rangle \left\langle\left.j_6j_2m_6m_2\right|j_6j_2j_4m_4\right\rangle \left\langle\left.j_1j_5m_1m_5\right|j_1j_5j_6m_6\right\rangle\\ &=(-1)^{j_2+j_3+j_5+j_6}\sqrt{2j_3+1}\sqrt{2j_6+1}\left\langle\left.j_3j_5m_3m_5\right|j_3j_5j_4m_4\right\rangle \left\{\begin{array}{ccc} j_1 & j_2 & j_3 \\ j_4 & j_5 & j_6 \end{array}\right\}, \end{align}
one can get
\begin{align} &\sum_{m_{l1},m_{s},m_{l2}} \left\langle\left.l_1m_{l1}\frac{1}{2}m_{s}\right|j_1m_1\right\rangle \left\langle\left.l_2m_{l2}\frac{1}{2}m_{s}\right|j_2m_2\right\rangle {\langle}l_1m_{l1}{\lambda}{\mu}|l_2m_{l2}{\rangle}\\ &=(-1)^{1/2+j_1+\lambda+l_2}\sqrt{2j_1+1}\sqrt{2l_2+1}\left\langle\left.j_1m_1\lambda\mu\right|j_2m_2\right\rangle \left\{\begin{array}{ccc} l_1 & 1/2 & j_1 \\ j_2 & \lambda & l_2 \end{array}\right\}, \end{align}
where we replaced $j_1 \rightarrow l_1$, $j_2 \rightarrow 1/2$, $j_3 \rightarrow j_1$, $j_4 \rightarrow j_2$, $j_5 \rightarrow \lambda$, $j_6 \rightarrow l_2$, $m_1 \rightarrow m_{l1}$, $m_2 \rightarrow m_s$, $m_3 \rightarrow m_1$, $m_4 \rightarrow m_2$, $m_5 \rightarrow \mu$, and $m_6 \rightarrow m_{l2}$. Then,
\begin{align} &\left\langle\left.n_2l_2j_2m_2\right|i^{\lambda}er^{\lambda}Y_\lambda^\mu(\theta,\phi)\left|n_1l_1j_1m_1\right.\right\rangle\\ &= e(-1)^{1/2+j_1+\lambda+l_2}i^{l_1-l_2+\lambda}\sqrt{\frac{(2l_1+1)(2j_1+1)(2{\lambda}+1)}{4\pi}}{\langle}l_10{\lambda}0|l_20{\rangle}\left\langle\left.j_1m_1\lambda\mu\right|j_2m_2\right\rangle \left\{\begin{array}{ccc} l_1 & 1/2 & j_1 \\ j_2 & \lambda & l_2 \end{array}\right\} \int_0^{\infty}r^{\lambda}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r)r^2dr. \end{align}
From the formula [Bohr & Mottelson Vol. I, Eq. (3A-16a)], for even values of $l_1+\lambda-l_2$, \begin{align} \left\{\begin{array}{ccc} l_1 & 1/2 & j_1 \\ j_2 & \lambda & l_2 \end{array}\right\} = (-1)^{l_2+1/2+j_2}(2l_1+1)^{-1/2}(2j_2+1)^{-1/2}\frac{\left\langle\left.j_1\frac{1}{2}{\lambda}0\right|j_2\frac{1}{2}\right\rangle}{{\langle}l_10{\lambda}0|l_20{\rangle}}. \end{align}
Since $l_2+1/2+j_2$ is a integer, \begin{align} (-1)^{l_2+1/2+j_2} = (-1)^{-l_2-1/2-j_2}, \end{align}
and
\begin{align} \left\{\begin{array}{ccc} l_1 & 1/2 & j_1 \\ j_2 & \lambda & l_2 \end{array}\right\} = (-1)^{-l_2-1/2-j_2}(2l_1+1)^{-1/2}(2j_2+1)^{-1/2}\frac{\left\langle\left.j_1\frac{1}{2}{\lambda}0\right|j_2\frac{1}{2}\right\rangle}{{\langle}l_10{\lambda}0|l_20{\rangle}}. \end{align}
Then,
\begin{align} &\left\langle\left.n_2l_2j_2m_2\right|i^{\lambda}er^{\lambda}Y_\lambda^\mu(\theta,\phi)\left|n_1l_1j_1m_1\right.\right\rangle\\ &= e(-1)^{j_1-j_2+\lambda}i^{l_1-l_2+\lambda}\sqrt{\frac{(2j_1+1)(2{\lambda}+1)}{4\pi(2j_2+1)}}\left\langle\left.j_1m_1\lambda\mu\right|j_2m_2\right\rangle \left\langle\left.j_1\frac{1}{2}{\lambda}0\right|j_2\frac{1}{2}\right\rangle \int_0^{\infty}r^{\lambda}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r)r^2dr. \end{align}
By the Wigner–Eckart theorem [Bohr & Mottelson Vol. I, Eq. (1A-60)] \begin{align} \left\langle\left.n_2l_2j_2m_2\right|i^{\lambda}r^{\lambda}Y_\lambda^\mu(\theta,\phi)\left|n_1l_1j_1m_1\right.\right\rangle = \left\langle\left.j_1m_1\lambda\mu\right|j_2m_2\right\rangle \frac{\left\langle\left.\left.n_2l_2j_2\right|\right|i^{\lambda}r^{\lambda}Y_\lambda(\theta,\phi)\left|\left|n_1l_1j_1\right.\right.\right\rangle}{\sqrt{2j_2+1}}, \end{align}
one can get
\begin{align} \left\langle\left.\left.n_2l_2j_2\right|\right|i^{\lambda}er^{\lambda}Y_\lambda(\theta,\phi)\left|\left|n_1l_1j_1\right.\right.\right\rangle &= e(-1)^{j_1-j_2+\lambda}i^{l_1-l_2+\lambda}\sqrt{\frac{(2j_1+1)(2{\lambda}+1)}{4\pi}} \int_0^{\infty}r^{\lambda}\mathscr{R}_{n_2l_2j_2}^*(r)\mathscr{R}_{n_1l_1j_1}(r)r^2dr \left\langle\left.j_1\frac{1}{2}{\lambda}0\right|j_2\frac{1}{2}\right\rangle\\ &= e(-1)^{j_1-j_2+\lambda}i^{l_1-l_2+\lambda}\sqrt{\frac{(2j_1+1)(2{\lambda}+1)}{4\pi}} \left\langle j_2\left|r^\lambda\right|j_1\right\rangle \left\langle\left.j_1\frac{1}{2}{\lambda}0\right|j_2\frac{1}{2}\right\rangle. \end{align}
As a result,
\begin{align} \left\langle\left.\left.j_2\right|\right|i^\lambda\mathscr{M}(E\lambda)\left|\left|j_1\right.\right.\right\rangle = e(-1)^{j_1-j_2+\lambda}i^{l_1-l_2+\lambda}\sqrt{\frac{(2j_1+1)(2{\lambda}+1)}{4\pi}} \left\langle j_2\left|r^\lambda\right|j_1\right\rangle \left\langle\left.j_1\frac{1}{2}{\lambda}0\right|j_2\frac{1}{2}\right\rangle. \end{align} This is written as Eq. (3C-33) of Bohr & Mottelson Vol. I. By substituting this formula into the following formula [Bohr & Mottelson Vol. I, Eq. (3C-17)] \begin{align} B(E\lambda; j_1 \rightarrow j_2) = \frac{1}{2j_1+1}\left| \left\langle j_2 \left|\left| i^\lambda\mathscr{M}(E\lambda) \right|\right| j_1 \right\rangle\right|^2, \end{align} one can get \begin{align} B(E\lambda; j_1 \rightarrow j_2) = e^2 \frac{2\lambda+1}{4\pi} \left\langle j_2 \left| r^\lambda \right| j_1 \right\rangle^2 \left\langle j_1 \frac{1}{2} \lambda 0 \left| j_2 \frac{1}{2}\right.\right\rangle^2. \end{align}