the number of trajectories, which is of course typical for Monte Carlo algorithms.

To further estimate ${\mathbb{E}}_{\Xi} \braket{\Phi_{\Xi}(t)}$, we use the definition of the trajectory to get $$\label{eq:Phi_norm}
\braket{\Phi_{\Xi}(t)} \leqslant
  \exp \left( 2\int_0^t \int_{\Omega} \lambda(s,\omega)
    \,\mathrm{d}\mu_{\omega} \,\mathrm{d}s \right)
  \left( \prod_{m=1}^M \| {\widetilde{A}}(t_m, \omega_m) \| \right)^2
  \braket{\Psi(0)},$$ where $\| \cdot \|$ is the operator norm, and we have used the fact that $U(\cdot, \cdot)$ is a unitary operator. If the intensity function $\lambda(t,\omega)$ is chosen such that there exist constants $\Lambda$ and ${\widetilde{\alpha}}$ satisfying $$\label{eq:Lam_al}
|\lambda(t,\omega)| \leqslant \Lambda \quad \text{and}
  \quad \|{\widetilde{A}}(t,\omega)\| \leqslant {\widetilde{\alpha}}, \qquad
  \forall t \in \mathbb{R}^+, \quad \forall \omega \in \Omega,$$ we can apply to get $$\label{eq:estimation}
{\mathbb{E}}_{\Xi} \braket{\Phi_{\Xi}(t)} \leqslant
  \mathrm{e}^{\Lambda t} \braket{\Psi(0)} \sum_{M=0}^{+\infty}
    \frac{t^M}{M!} \Lambda^M {\widetilde{\alpha}}^{2M}
  = \mathrm{e}^{\Lambda(1+{\widetilde{\alpha}}^2)t} \braket{\Psi(0)}.$$ Here one sees that the variance grows exponentially with respect to $t$.

Note that the equation becomes an equality if every operator $A(t,\omega)$ is a unitary operator multiplied by a positive constant $\alpha$. This condition holds for if we have either $\gamma_{jk}(t) \equiv 0$ or $\gamma_{jk}(t) \equiv \gamma$ for any $j,k$, and in this case, we have $\alpha = \norm{A(t, \omega)} = N_c \gamma$, where $N_c$ is the number of the pairs of coupling spins, i.e. the number 