• A. Bohr and B. R. Mottelson, Nuclear Structure (Benjamin, Reading, MA, 1975), Vol. I.
• B. A. Brown, Lecture Notes in Nuclear Structure Physics, 2011.

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).

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}

Therefore,

\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 &=\Gamma(l_2+3/2)\binom{n_2+l_2-1/2}{n_2-1}\frac{(l_1-l_2)_{n_1-n_2}}{(n_1-n_2)!}\\ &=\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\\ &= \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}\\ &= \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,l1+1/2,x]*LaguerreL[n2,l2+1,x],{x,0,Infinity}]

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}

• NIST DLMF Spherical and Spheroidal Harmonics, 3j,6j,9j Symbols - §34.1 Special Notation
• Wolfram Research (functions.wolfram.com) SphericalHarmonicY - Integration
• I. Hamamoto PRC69(2004)041306(R), PRC81(2010)021304
• Particle Data Group Clebsch-Gordan Coefficients, Spherical Harmonics, and d Functions
• Wikipedia Spherical harmonics, Clebsc–Gordan coefficients

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}