  Using a data sample corresponding to an integrated luminosity of 11.3 $\rm
fb^{-1}$ collected at center-of-mass energies from $4.23$ to $4.70$~GeV with
the BESIII detector, we observe the process $e^+e^-\to\pi^0\pi^0\psi_2(3823)$
for the first time with a statistical significance of $6.0\sigma$. The ratio of
average cross sections for $e^+e^-\to\pi^0\pi^0\psi_2(3823)$ over
$\pi^+\pi^-\psi_2(3823)$ is determined to be
$\mathcal{R}=\frac{\sigma[e^+e^-\to\pi^0\pi^0\psi_2(3823)]}{\sigma[e^+e^-\to\pi^+\pi^-\psi_2(3823)]}=0.57\pm
0.14\pm0.05$, which is consistent with expectations from isospin symmetry. Here
and below, the first uncertainties are statistical and the second are
systematic. The mass of the $\psi_2(3823)$ is measured to be
$M[\psi_2(3823)]=3824.5\pm 2.4\pm 1.0$~MeV/$c^2$. Due to the limited
statistics, the intrinsic width of the $\psi_2(3823)$ is measured as an upper
limit of $18.8$~MeV at $90\%$ confidence level.
