is strictly positive can be used, a better choice will however reduce the sampling variance of the algorithm. We will further discuss the choice of the intensity function in the next section.

The following identity is essential to our method: It turns the problem of calculating $\ket{\Psi(t)}$ into a sampling problem: $$\label{eq:expectation}
\ket{\Psi(t)} = {\mathbb{E}}_{\Xi} \ket{\Phi_{\Xi}(t)},$$ if for a given realization of the marked point process $\Xi$, the state $\ket{\Phi_{\Xi}(t)}$ is defined by $$\label{eq:traj_def}
\begin{split}
\ket{\Phi_{\Xi}(t)} = & \exp \left(
  \int_0^t\int_{\Omega} \lambda(s,\omega) \,\mathrm{d}\mu_{\omega} \,\mathrm{d}s
\right) U(t,t_M) {\widetilde{A}}(t_M,\omega_M) \times \\
& \qquad \times U(t_M, t_{M-1}) {\widetilde{A}}(t_{M-1}, \omega_{M-1}) \cdots
  U(t_2,t_1) {\widetilde{A}}(t_1,\omega_1) U(t_1,0) \ket{\Psi(0)},
\end{split}$$ where $M$ is the number of marks in $\Xi$ before time $t$, and we have used the short hands $$\label{eq:hat_A}
  {\widetilde{A}}(t_m, \omega_m) = -{\mathrm{i}}A(t_m, \omega_m) / \lambda(t_m,\omega_m),
  \qquad m = 1, \cdots, M.$$

The equality follows as for any function ${\mathcal{F}}(\Xi)$ depending only on the part of $\Xi$ in the time interval $[0,t)$, we have $$\label{eq:expectation_F}
\begin{split}
{\mathbb{E}}_{\Xi} {\mathcal{F}}(\Xi) &= \sum_{M=0}^{+\infty} \int_0^t \int_{\Omega}
  \int_0^{t_M} \int_{\Omega} \cdots \int_0^{t_2} \int_{\Omega}
\exp \left(
  -\int_0^t\int_{\Omega} \lambda(s,\omega) \,\mathrm{d}\mu_{\omega} \,\mathrm{d}s
\right) \times \\
& \qquad \times \left( \prod_{m=1}^M \lambda(t_m, \omega_m) \right) {\mathcal{F}}(\Xi)
  \,\mathrm{d}\mu_{\omega_1} \,\mathrm{d}t_1 \cdots
  \,\mathrm{d}\mu_{\omega_{M-1}} \,\mathrm{d}t_{M-1} \,\mathrm{d}\mu_{\omega_M}
  \,\mathrm{d}t_M.
\end{split}$$

An algorithm to evaluate 