  Using a total of $9.0~\mathrm{fb}^{-1}$ of $e^+e^-$ collision data with
center-of-mass energies between 4.15 and 4.30 GeV collected by the BESIII
detector, we search for the processes $e^+e^-\to \gamma X(3872)$ with
$X(3872)\to\pi^0\chi_{cJ}$ for $J=0,1,2$. We report the first observation of
$X(3872)\to\pi^{0}\chi_{c1}$, a new decay mode of the $X(3872)$, with a
statistical significance of more than 5$\sigma$. Normalizing to the previously
established process $e^+e^-\to \gamma X(3872)$ with
$X(3872)\to\pi^+\pi^-J/\psi$, we find ${\cal B}(X(3872) \to \pi^0
\chi_{c1})/{\cal B}(X(3872) \to \pi^+\pi^- J/\psi) =
0.88^{+0.33}_{-0.27}\pm0.10$, where the first error is statistical and the
second is systematic. We set 90% confidence level upper limits on the
corresponding ratios for the decays to $\pi^0\chi_{c0}$ and $\pi^0\chi_{c2}$ of
19 and 1.1, respectively.
