the 2-norm of numerical error $\ket{\Psi^{\mathrm{err}}(t)} = \ket{\Psi^{\mathrm{num}}(t)} -
\ket{\Psi^{\mathrm{ref}}(t)}$: $$\label{eq:error}
e_{N_{\mathrm{traj}}}(t) = \sqrt{\braket{\Psi^{\mathrm{err}}(t)}},$$ where $\ket{\Psi^{\mathrm{num}}(t)}$ is the numerical solution of QKMC with $N_{\mathrm{traj}}$ trajectories. Figure \[fig:error\] shows the reduction of the error as the number of trajectories increases, and the order of convergence is calculated in Figure \[fig:order\]. It is obvious that the numerical order is around $1/2$ for all $t$, which indicates that the expected convergence rate in the Monte Carlo method is achieved in our numerical test.

![Evolution of the numerical error $e_{N_{\mathrm{traj}}}$ (defined in ) for 4 spins. The $x$-axis is the time $t$, and the $y$-axis is the 2-norm error of the solution.[]{data-label="fig:error"}](spin4error.pdf){width=".6\textwidth"}

![Evolution of the numerical order for 4 spins. The $x$-axis is the time $t$, and the $y$-axis is the numerical order obtained as indicated in the legend.[]{data-label="fig:order"}](spin4order.pdf){width=".6\textwidth"}

High dimensional tests
----------------------

The cases with more spins are shown in Figure \[fig:10-40spins\]. Here we consider weaker coupling intensity $\gamma_0 = 0.01$, while the number of spins ranges from 10 to 40. In these cases, the method with $Z$-basis needs a huge number of trajectories to get meaningful results, which is out of our current computational capacity. However, by using $8,000,000$ trajectories, two complete 