  Motivated by very recent measurement of the radiative decays of
$\Upsilon(1S)$ to $\chi_{c1}$ at Belle, we use the nonrelativistic QCD
factorization theory and calculate the branching fractions of the radiative
decays of bottomonium into S-wave and P-wave charmonium, i.e. $\Upsilon(nS) \to
\eta_c(nS)+\gamma$ and $\Upsilon(nS)\to \chi_{cJ}+\gamma$. We systematically
studied the branching fractions of the radiative decays of bottomonium into
charmonium. Compared to the previous calculation, we obtained the analytical
expression for the decay widths and considered the color-octet contributions.
For $\Upsilon(nS) \to \eta_c(nS)+\gamma$, the relativistic corrections are also
obtained. Through the calculation, the theoretical prediction for
$\Upsilon(1S)\to\chi_{c1}+\gamma$ is still smaller than the recent Belle
measurement. Further theoretical work and experimental analysis are necessary
to understand the $\chi_{c1}$ production mechanism in upsilon decays.
