\cdots \,\mathrm{d}\mu_{\omega_{M-1}}
  \,\mathrm{d}t_{M-1} \,\mathrm{d}\mu_{\omega_M} \,\mathrm{d}t_M,
\end{split}$$ where the last step is just renaming the integration variables. Using dominated convergence, it is easy to see that if the operators $A(\omega)$ are uniformly bounded, the above series expansion converges absolutely. According to the properties B) and C), the integrand in the above equation is always a state in ${\mathcal{A}}$, and is therefore easy to represent. The challenge then lies in the high dimensional integration both in the stochastic space and in the time sequence, for which we turn to Monte Carlo method.

Before discussing the algorithm, let us note that the above can be extended to time-dependent Hamiltonian operators, with the assumption now changes to for any $\ket{\Psi} \in {\mathcal{A}}$ $${\mathbb{E}}_{\omega} A(t,\omega) \ket{\Psi} = H_{\mathrm{hard}}(t) \ket{\Psi}, \qquad \text{and}  \qquad A(t, \omega) \ket{\Psi} \in {\mathcal{A}},$$ and the semigroup generated by $H_{\mathrm{easy}}(t)$ preserves the state in ${\mathcal{A}}$ (it is also possible to consider the more general case that the set ${\mathcal{A}}$ depends on time). The expansion changes to $$\label{eq:Psi_expansion_td}
\begin{split}
  \ket{\Psi(t)} &= \sum_{M=0}^{+\infty} \int_0^t
  \int_{\Omega} \int_0^{t_M} \int_{\Omega} \cdots
  \int_0^{t_2} \int_{\Omega} (-{\mathrm{i}})^M \times \\
  & \qquad \times U(t,t_M) A(t_M,\omega_M)
  U(t_M,t_{M-1}) A(t_{M-1},\omega_{M-1}) \cdots \\
  & \qquad \times U(t_2,t_1) A(t_1,\omega_1) U(t_1,0)
  \ket{\Psi(0)} \,\mathrm{d}\mu_{\omega_1} \,\mathrm{d}t_1 \cdots
  \,\mathrm{d}\mu_{\omega_{M-1}} \,\mathrm{d}t_{M-1} 