Using 106~million $\psi(3686)$ events collected with the BESIII detector, we measure multipole amplitudes for the decay $\psi(3686) \rightarrow \gamma\chi_{c1,2}\to\gamma\gamma J/\psi$ beyond the dominant electric-dipole amplitudes. The normalized magnetic-quadrupole amplitude for $\psi(3686) \rightarrow \gamma\chi_{c1,2} \rightarrow \gamma\gamma J/\psi$ and the normalized electric-octupole amplitudes for $\psi(3686) \rightarrow\gamma\chi_{c2}$,~$\chi_{c2} \rightarrow\gamma J/\psi$ are determined. The M2 amplitudes for $\psi(3686) \rightarrow\gamma\chi_{c1}$ and $\chi_{c1,2} \rightarrow\gamma J/\psi$ are found to differ significantly from zero and are consistent with theoretical predictions. We also obtain the ratios of M2 contributions of $\psi(3686)$ and $J/\psi$ decays to $\chi_{c1,2}$, $b_{2}^{1}/b_{2}^{2} = 1.35\pm0.72$ and $a_{2}^{1}/a_{2}^{2} = 0.617\pm0.083$, which agree well with theoretical expectations. By considering the multipole contributions of $\chi_{c1,2}$, we measure the product branching fractions for the cascade decays $\psi(3686) \rightarrow\gamma\chi_{c0,1,2}\to\gamma\gamma J/\psi$ and search for the process $\eta_{c}(2S)\to\gamma J/\psi$ through $\psi(3686)\rightarrow\gamma\eta_{c}(2S)$. The product branching fraction for $\psi(3686) \rightarrow\gamma\chi_{c0}\to\gamma\gamma J/\psi$ is 3$\sigma$ larger than published measurements, while those of $\psi(3686) \rightarrow\gamma\chi_{c1,2}\to\gamma\gamma J/\psi$ are consistent. No significant signal for the decay $\psi(3686)\to\gamma \eta_c(2S)\to\gamma \gamma J/\psi$ is observed, and the upper limit of the product branching fraction at the 90\% confidence level is determined.