Observation of the decay $X(3872) to pi^0 chi_{c1}(1P)$

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)topi^0chi_{cJ}$ for $J=0,1,2$. We report the first observation of $X(3872)topi^{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)topi^+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^0chi_{c0}$ and $pi^0chi_{c2}$ of 19 and 1.1, respectively.