  The numerical analysis of stochastic time fractional evolution equations
presents considerable challenges due to the limited regularity of the model
caused by the nonlocal operator and the presence of noise.
  The existing time-stepping methods exhibit a significantly low order
convergence rate. In this work, we introduce a smoothing technique and develop
the novel high-order schemes for solving the linear stochastic fractional
evolution equations driven by integrated additive noise. Our approach involves
regularizing the additive noise through an $m$-fold integral-differential
calculus, and discretizing the equation using the $k$-step BDF convolution
quadrature. This novel method, which we refer to as the ID$m$-BDF$k$ method, is
able to achieve higher-order convergence in solving the stochastic models. Our
theoretical analysis reveals that the convergence rate of the ID$2$-BDF2 method
is $O(\tau^{\alpha + \gamma -1/2})$ for $1< \alpha + \gamma \leq 5/2$, and
$O(\tau^{2})$ for $5/2< \alpha + \gamma <3$, where $\alpha \in (1, 2)$ and
$\gamma \in (0, 1)$ denote the time fractional order and the order of the
integrated noise, respectively. Furthermore, this convergence rate could be
improved to $O(\tau^{\alpha + \gamma -1/2})$ for any $\alpha \in (1, 2)$ and
$\gamma \in (0, 1)$, if we employ the ID$3$-BDF3 method. The argument could be
easily extended to the subdiffusion model with $\alpha \in (0, 1)$. Numerical
examples are provided to support and complement the theoretical findings.
