than the control field, the growing rate of $\eta(t)$ in is larger than that in , which indicates larger variance in the numerical solution with $Z$-basis.

Analysis of the sampling variance
=================================

In this section, we discuss the evolution of the numerical error for the above algorithm. Suppose $N_{\mathrm{traj}}$ trajectories are used in the simulation. Then the numerical solution is $$\label{eq:Psi_num}
\ket{\Psi_{\mathrm{num}}(t)} = \frac{1}{N_{\mathrm{traj}}}
  \sum_{i=1}^{N_{\mathrm{traj}}} \ket{\Phi_{\Xi_i}(t)},$$ where $\Xi_i$ is the realization of the marked point process corresponding to the $i$-th trajectory, and $\Xi_i$ and $\Xi_j$ are independent of each other if $i \neq j$. Define $\ket{\Psi_{\mathrm{err}}(t)}$ as the difference between the numerical solution and the exact solution: $$\label{eq:Psi_err}
\ket{\Psi^{\mathrm{err}}(t)} =
  \ket{\Psi^{\mathrm{num}}(t)} - \ket{\Psi(t)}.$$ Apparently ${\mathbb{E}}\ket{\Psi^{\mathrm{err}}(t)} = 0$. Therefore the variance of $\ket{\Psi^{\mathrm{err}}(t)}$ is $$\label{eq:sigma}
\sigma^2 = {\mathbb{E}}\braket{\Psi^{\mathrm{err}}(t)},$$ which will be estimated below.

Inserting and into and using the fact that ${\mathbb{E}}_{\Xi_i} \ket{\Phi_{\Xi_i}(t)} = \ket{\Psi(t)}$ for any $i$, we have $$\sigma^2 = \frac{1}{N_{\mathrm{traj}}^2}
  \sum_{i=1}^{N_{\mathrm{traj}}} \sum_{j=1}^{N_{\mathrm{traj}}}
    {\mathbb{E}}_{\Xi_i,\Xi_j} \braket{\Phi_{\Xi_i}(t)}{\Phi_{\Xi_j}(t)} - \braket{\Psi(t)}.$$ Since all trajectories are independent, in the above sum, the terms with $i
\neq j$ can be directly evaluated, which yields $$\sigma^2 = \frac{1}{N_{\mathrm{traj}}} {\mathbb{E}}_{\Xi} \braket{\Phi_{\Xi}(t)}
  - \frac{1}{N_{\mathrm{traj}}} \braket{\Psi(t)}.$$ This equation shows that the the variance is proportional to the inverse of 