  Anomalous diffusion in the presence or absence of an external force field is
often modelled in terms of the fractional evolution equations, which can
involve the hyper-singular source term.
  For this case, conventional time stepping methods may exhibit a severe order
reduction. Although a second-order numerical algorithm is provided for the
subdiffusion model with a simple hyper-singular source term $t^{\mu}$,
$-2<\mu<-1$ in [arXiv:2207.08447], the convergence analysis remain to be
proved.
  To fill in these gaps, we present a simple and robust smoothing method for
the hyper-singular source term, where the Hadamard finite-part integral is
introduced. This method is based on the smoothing/ID$m$-BDF$k$ method proposed
by the authors [Shi and Chen, SIAM J. Numer. Anal., to appear] for subdiffusion
equation with a weakly singular source term.
  We prove that the $k$th-order convergence rate can be restored for the
diffusion-wave case $\gamma \in (1,2)$ and sketch the proof for the
subdiffusion case $\gamma \in (0,1)$, even if the source term is hyper-singular
and the initial data is not compatible.
  Numerical experiments are provided to confirm the theoretical results.
