  In this work, we calculate the sub-leading power contributions to the
radiative leptonic $D\to\gamma \,\ell \,\nu$ decay. For the first time, we
provide the analytic expressions of next-to-leading power contributions and the
error estimation associated with the power expansion of ${\cal O}(\Lambda_{\rm
QCD}/m_c)$. In our calculation, we adopt two different models of the $D$-meson
distribution amplitudes $\phi_{D,\rm I}^+$ and $\phi_{D,\rm II}^+$. Within the
framework of the QCD factorization as well as the dispersion relation, we
evaluate the soft contribution up to the next-to-leading logarithmic accuracy,
and the higher-twist contribution from the two-particle and three-particle
distribution amplitudes is also considered. Finally, we find that all the
sub-leading power contributions are significant at $\lambda_D(\mu_0)=354\,\rm
MeV$, and the next-to-leading power contributions will lead to $143\%$ in
$\phi_{D,\rm I}^+$ and $120\%$ in $\phi_{D,\rm II}^+$ corrections to leading
power vector form factors with $E_{\gamma}=0.5\,\rm GeV$. As the corrections
from the higher-twist and local sub-leading power contributions will be
enhanced with the growing of the inverse moment, it is difficult to extract an
appropriate inverse moment of the $D$-meson distribution amplitude. The
predicted branching fractions are $(1.88_{-0.29}^{+0.36})\times10^{-5}$ for
$\phi_{D,\rm I}^+$ and $(2.31_{-0.54}^{+0.65})\times10^{-5}$ for $\phi_{D,\rm
II}^+$, respectively.
