by applications in nuclear magnetic resonance (NMR) [@Ernst:90], where the nuclear spins react to magnetic fields. In such applications, for the coefficient of the single body term $m(t) = (m_1(t), m_2(t), m_3^{(k)})^{\top}$, $(m_1(t), m_2(t))$ is prescribed as a control field (such that the control magnetic field is only in $(x,y)$ direction) and the time independent $m_3^{(k)}$ is understood as an energy splitting of the $\ket{\uparrow}$ and $\ket{\downarrow}$ states of the $k$-th spin. The $\gamma_{jk}$ terms account for dipole-dipole interactions between the nuclear spins, and the magnitude of which decays very fast as the distance between nuclei increases and is usually quite small compared to the energy splitting and the external fields. The value of these coefficients might contain some uncertainty due to experimental imperfectness. One potential application of our method is robust control of NMR via pulse design (see e.g., [@FrieseckeHennekeKunisch] and references therein) when the spin-spin interactions are taken into account, which we will leave for future works.

The high dimensionality of the system naturally calls for Monte Carlo type methods. The motivation of the algorithm proposed in this work comes from a surface hopping method recently developed by us in [@CaiLu2017], which can be viewed as a stochastic 