  Using $4.48 \times 10^{8}$ $\psi(3686)$ events collected with the BESIII
detector, we search for the decays $\chi_{cJ} \rightarrow \mu^{+}\mu^{-}J/\psi$
through the radiative decays $\psi(3686) \rightarrow \gamma\chi_{cJ}$, where
$J=0,1,2$. The decays $\chi_{c1,2} \rightarrow \mu^{+}\mu^{-}J/\psi$ are
observed, and the corresponding branching fractions are measured to be
$\mathcal{B}(\chi_{c1} \rightarrow \mu^{+}\mu^{-}J/\psi) = (2.51 \pm 0.18 \pm
0.20)\times10^{-4}$ and $\mathcal{B}(\chi_{c2} \rightarrow
\mu^{+}\mu^{-}J/\psi) = (2.33 \pm 0.18 \pm 0.29)\times10^{-4}$, where the first
uncertainty is statistical and the second one systematic. No significant
$\chi_{c0} \rightarrow \mu^{+}\mu^{-}J/\psi$ decay is observed, and the upper
limit on the branching fraction is determined to be $2.0\times10^{-5}$ at 90%
confidence level. Also, we present a study of di-muon invariant mass dependent
transition form factor for the decays $\chi_{c1,2} \rightarrow
\mu^{+}\mu^{-}J/\psi$.
