===== Reduced transition probability =====

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

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

==== Reduced matrix element of the electric operator for wave functions in a spherically symmetric potential  ====
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}



==== For wave functions in a spherically symmetric potential ====
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 [[https://dlmf.nist.gov/|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 [[https://dlmf.nist.gov/18.8|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 [[http://functions.wolfram.com/Polynomials/LaguerreL3/21/ShowAll.html|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|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}

==== Special cases ====

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}

=== $n_1=n_2=n$, $l_1=l_2=l$ ===

\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}

=== $n_1=n_2=n$, $l_1=l_2=l$, $\lambda=2$ ===

\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, ..., 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}

=== $l_1+\lambda=l_2$ ===

\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, ..., 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|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}
([[https://ja.wikipedia.org/wiki/%E3%83%9D%E3%83%83%E3%83%9B%E3%83%8F%E3%83%9E%E3%83%BC%E8%A8%98%E5%8F%B7|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}

==== $n_1=2,\ n_2=1,\ l_1=0,\ l_2=1,\ \lambda=1$ ====

From wolfram alpha,
<code>
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]
</code>

==== Backup ====


From [[http://functions.wolfram.com/Polynomials/LaguerreL3/21/ShowAll.html|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}

==== Appendix ====
=== Integral of three spherical harmonics ===
== Formula ==
Based on [[http://functions.wolfram.com/Polynomials/SphericalHarmonicY/21/ShowAll.html|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}

=== Notations ===
== Spherical harmonics and Clebsch-Gordan coefficients  ==
|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$ |

== References ==
  * NIST DLMF [[https://dlmf.nist.gov/14.30|Spherical and Spheroidal Harmonics]], [[http://dlmf.nist.gov/34.1| 3j,6j,9j Symbols - §34.1 Special Notation]]
  * Wolfram Research (functions.wolfram.com) [[http://functions.wolfram.com/Polynomials/SphericalHarmonicY/21/ShowAll.html|SphericalHarmonicY - Integration]]
  * I. Hamamoto [[http://link.aps.org/doi/10.1103/PhysRevC.69.041306|PRC69(2004)041306(R)]], [[http://link.aps.org/doi/10.1103/PhysRevC.81.021304|PRC81(2010)021304]]
  * Particle Data Group [[http://pdg.lbl.gov/2018/reviews/rpp2018-rev-clebsch-gordan-coefs.pdf|Clebsch-Gordan Coefficients, Spherical Harmonics, and d Functions]]
  * Wikipedia [[https://en.wikipedia.org/wiki/Spherical_harmonics|Spherical harmonics]], [[https://en.wikipedia.org/wiki/Clebsch%E2%80%93Gordan_coefficients|Clebsc–Gordan coefficients]]

=== Calculation notes ===

== Reduced matrix element for spin less particle ==
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}

== Definition of the reduced transition probability ==

** 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 ([[http://functions.wolfram.com/HypergeometricFunctions/ClebschGordan/23/ShowAll.html|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 [[http://pdg.lbl.gov/2018/reviews/rpp2018-rev-clebsch-gordan-coefs.pdf|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}

== Reduced matrix element of the electric operator ==
**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 ([[http://functions.wolfram.com/Polynomials/SphericalHarmonicY/21/ShowAll.html|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 ([[http://functions.wolfram.com/HypergeometricFunctions/ClebschGordan/23/ShowAll.html|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}