Hamiltonians, where $m^{(k)}(t) = (m_1^{(k)}(t), m_2^{(k)}(t), m_3^{(k)}(t))^{\top}$ is the magnetic field acting on the $k$th spin. For the two-body interaction in the Hamiltonian , we have used the notations $$\begin{aligned}
  & \sigma^{(j)} \cdot \sigma^{(k)} = \sigma_1^{(j)} \sigma_1^{(k)} 
    + \sigma_2^{(j)} \sigma_2^{(k)} + \sigma_3^{(j)} \sigma_3^{(k)}, \\
  & \sigma_i^{(j)} \sigma_i^{(k)} = {\mathrm{Id}}^{\otimes (j-1)} \otimes \sigma_i \otimes {\mathrm{Id}}^{\otimes (k-j-1)} \otimes \sigma_i \otimes {\mathrm{Id}}^{\otimes (N-k)}.\end{aligned}$$ Thus $\sigma^{(j)} \cdot \sigma^{(k)}$ counts for a Heisenberg type interaction between the $j$-th and $k$-th spins with $\gamma_{jk}(t)$ being the interaction strength or coupling intensity.

This paper concerns numerical algorithm for the time evolution of the system: The many-body wave function $\ket{\Psi} \in \mathcal{H}$ is govern by the Schrödinger equation $$\label{eq:Psi_eq} \frac{\mathrm{d}}{\mathrm{d} t}
  \ket{\Psi} = -{\mathrm{i}}H \ket{\Psi}.$$ While is a linear ODE system, solving the system directly is impractical even for dozens of spins, as the size of the system $\dim \mathcal{H} = 2^N$ grows exponentially as the number of spins increases. In fact, even representing a particular state $\ket{\Psi}$ is challenging: for $N = 40$, the size of the vector is greater than $1$ trillion; not mentioning the computational cost involved in evaluation the matrix-vector product $H\ket{\Psi}$.

While the Hamiltonian is rather general, our algorithm development is mainly motivated 