method to solve generic high dimensional ODE systems (or PDE systems combined with some particle / semiclassical methods, as in [@LuZhou2016; @CaiLu2017]). The overall idea of the algorithm contains two elements:

1.  a series expansion of the solution of the system based on a time-dependent perturbation theory, *i.e.*, from the ODE point of view, a repeated back-substitution in the integral form of the system based on Duhamel’s principle;

2.  a Monte Carlo method to stochastically evaluate the series expansion based on an efficient representation of part of the Hilbert space.

In this work, we will apply the above framework to develop methods for quantum many-body spin dynamics. For a particular system, to make the algorithm efficient, it is crucial to identify a suitable “small term” to be used in the series expansion from time-dependent perturbation theory. In the setting of spin dynamics, we will discuss two approaches: 1) one is based on a decomposition of the Hamiltonian into terms commuting with $\sigma_3$ and those not commuting (*i.e.*, diagonal and off-diagonal terms in the $Z$-basis of the spins); 2) the other approach is based on a splitting of the Hamiltonian into single-body term and two-body interactions.

Over the years, many numerical 