  In this paper, the semileptonic decays $D_s^+ \to \pi^0\ell^+ \nu_\ell$ with
$\ell=(e,\mu)$ are investigated by using the light-cone sum rule approach.
Firstly, the neutral meson mixing scheme between $\pi^0$, $\eta$, $\eta^\prime$
and pseudoscalar gluonium $G$ is discussed in a unified way, which leads to the
direct connection between two different channels for $D_s^+\to
\pi^0\ell^+\nu_\ell$ and $D_s^+ \to \eta\ell^+\nu_\ell$ by the $\pi^0-\eta$
mixing angle. Then we calculated the $D_s\to \pi^0$ transition form factors
(TFFs) within QCD light-cone sum rule approach up to next-to-leading order
correction. At the large recoil point, we have
$f_+^{D_s^+\pi^0}(0)=0.0113_{-0.0019}^{+0.0024}$ and
$f_-^{D_s^+\pi^0}(0)=0.0020_{-0.0009}^{+0.0008}$. Furthermore, the TFFs are
extrapolated to the whole physical $q^2$-region by using the simplified
$z(q^2)$-series expansion. The behaviors of TFFs and related three angular
coefficient functions $a_{\theta_\ell}(q^2)$, $b_{\theta_\ell}(q^2)$ and
$c_{\theta_\ell}(q^2)$ are given. The differential decay widths for $D_s^+ \to
\pi^0\ell^+ \nu_\ell$ with respect to $q^2$ and $\cos\theta_\ell$ are
presented, and also lead to the branching fractions ${\cal B}(D_s^+\to \pi
^0e^+\nu_e) =2.60_{-0.51}^{+0.57}\times 10^{-5}$ and ${\cal B}(D_s^+\to \pi
^0\mu^+\nu _\mu )= 2.58_{-0.51}^{+0.56}\times 10^{-5}$. These results show well
agreement with the recent BESIII measurements and theoretical predictions. Then
the differential distributions and integrated predictions for three angular
observables, {\it i.e.} forward-backward asymmetries, $q^2$-differential flat
terms and lepton polarization asymmetries are given separately. Lastly, we
estimate the ratio for different decay channels ${\cal R}_{\pi
^0/\eta}^{\ell}=1.108_{-0.071}^{+0.039}\times 10^{-3}$.
