  In this paper, we propose a new approach for the time-discretization of the
incompressible stochastic Stokes equations with multiplicative noise. Our new
strategy is based on the classical Milstein method from stochastic differential
equations. We use the energy method for its error analysis and show a strong
convergence order of at most $1$ for both velocity and pressure approximations.
The proof is based on a new H\"older continuity estimate of the velocity
solution. While the errors of the velocity approximation are estimated in the
standard $L^2$- and $H^1$-norms, the pressure errors are carefully analyzed in
a special norm because of the low regularity of the pressure solution. In
addition, a new interpretation of the pressure solution, which is very useful
in computation, is also introduced. Numerical experiments are also provided to
validate the error estimates and their sharpness.
