In order to demonstrate the stochastic trajectory method let us consider a simple differential equation \[ \frac{\mathrm{d}}{\mathrm{d}t}x(t)=-\gamma x(t) \] that describes exponential decay of the variable \(x\). The solution of this equation is \[ x(t)=x(0)\mathrm{e}^{-\gamma t}\,. \] After a short time interval \(\Delta t\) (\(\Delta t\ll\gamma^{-1}\)) we can approximate the variable \(x\) as \[ x(t+\Delta t)=x(t)-\gamma x(t)\Delta t =x(t)(1-\gamma\Delta t)+0\cdot\gamma\Delta t\,. \] This equation can be interpreted in the following way: during the time interval \(\Delta t\) two possibilities can occur. Either the variable \(x\) stays the same, \(x(t+\Delta t)=x(t)\), with the probability \(p_{\mathrm{same}}=1-\gamma\Delta t\) or \(x\) becomes \(0\) with the probability \(p_{\mathrm{zero}}=\gamma\Delta t\). This interpretation describes a stochastic trajectory \(\tilde{x}(t)\). Differential equation can be replaced by a stochastic process, with the solution of the equation recovered by averaging the trajectories \(\tilde{x}(t)\). Numerically the trajectory \(\tilde{x}(t)\) can be obtained in the following way: Replace the time \(t\) with discrete values \(t_n\) having the step \(\Delta t\). For each time \(t_n\) generate a random number \(r_n\) uniformly distributed in the interval \([0,1]\). The next value \(\tilde{x}(t_{n+1})\) is determined according to the rule \[ \tilde{x}(t_{n+1})= \begin{cases} \tilde{x}(t_n)\,, & r_n > p_{\mathrm{zero}}\,,\\ 0\,, & r_n < p_{\mathrm{zero}}\,. \end{cases} \]
Stochastic simulation methods provide an alternative to the description of a quantum system using a density matrix. A particular stochastic method, the theory of quantum trajectories, has been developed by many authors [1–8]. This method describes the evolution of the quantum system in time by means of quantum trajectories for the states of the system subjected to random quantum jumps. The quantum trajectory is calculated by integrating the time-dependent Schrödinger equation using a non-Hermitian effective Hamiltonian. Incoherent processes such as spontaneous emission are incorporated as random quantum jumps that cause a collapse of the wave function to a single state [6]. The results for the density matrix are obtained by repeating the stochastic simulations several times and calculating the average. Using stochastic methods one can examine the behavior of individual trajectories, thus such methods can provide a description of an experiment on a single system in a more direct way.
We start by assuming that the Markovian approximation for the description of time evolution is valid. The dynamics of the quantum system is described by a master equation \[ \frac{\partial}{\partial t}\rho(t)=\mathcal{M}\rho\,, \] where \(\mathcal{M}\) is the superoperator describing the time evolution and \(\rho\) is the density matrix of the quantum system. Let us consider the superoperator \(\mathcal{M}\) that can be separated into two parts: \[ \mathcal{M}=\mathcal{L}+\mathcal{J}\,. \] The part \(\mathcal{J}\) is interpreted as describing quantum jumps, and \(\mathcal{L}\) describes the jump-free evolution. After a short time interval \(\Delta t\) the density matrix is \[ \rho(t+\Delta t)=\rho(t)+\mathcal{L}\rho(t)\Delta t+\mathcal{J}\rho(t)\Delta t\,. \] Since the master equation should preserve the trace of the density matrix, we have the equality \[ \mathrm{Tr}\{\mathcal{L}\rho(t)\}+\mathrm{Tr}\{\mathcal{J}\rho(t)\}=0\,. \] Thus we can write \[ \rho(t+\Delta t)=\frac{\rho+\mathcal{L}\rho(t)\Delta t}{1+\mathrm{Tr}\{\mathcal{L}\rho(t)\}\Delta t}(1-\mathrm{Tr}\{\mathcal{J}\rho(t)\}\Delta t) +\frac{\mathcal{J}\rho(t)}{\mathrm{Tr}\{\mathcal{J}\rho(t)\}}\mathrm{Tr}\{\mathcal{J}\rho(t)\}\Delta t\,. \] This equation can be interpreted in the following way: during the time interval \(\Delta t\) two possibilities can occur. Either after time \(\Delta t\) the density matrix is equal to conditional density matrix \[ \rho_{\mathrm{jump}}(t+\Delta t)=\frac{\mathcal{J}\rho(t)}{\mathrm{Tr}\{\mathcal{J}\rho(t)\}} \] with the probability \[ p_{\mathrm{jump}}(t)=\mathrm{Tr}\{\mathcal{J}\rho(t)\}\Delta t \] or to the density matrix \[ \rho_{\text{no-jump}}(t+\Delta t)=\frac{\rho+\mathcal{L}\rho(t)\Delta t}{1+\mathrm{Tr}\{\mathcal{L}\rho(t)\}\Delta t} \] with the probability \(1-p_{\mathrm{jump}}(t)\). This interpretation describes a stochastic trajectory \(\tilde{\rho}(t)\). The master equation can be replaced by a stochastic process, with the solution of the equation recovered by averaging the trajectories \(\tilde{\rho}(t)\).
Let us further assume that the superoperators \(\mathcal{L}\) and \(\mathcal{J}\) have the form \[ \mathcal{L}\rho=\frac{1}{i\hbar}(H_{\mathrm{eff}}\rho-\rho H_{\mathrm{eff}}^{\dagger}) \] and \[ \mathcal{J}\rho=C\rho C^{\dagger}\,. \] The operators \(H_{\mathrm{eff}}\) and \(C\) are non-Hermitian in general. If the superoperators \(\mathcal{L}\) and \(\mathcal{J}\) have this form and the density matrix at the time \(t\) factorizes as \(\rho(t)=|\Psi(t)\rangle\langle\Psi(t)|\) then after time interval \(\Delta t\) the density matrices \(\rho_{\mathrm{jump}}(t+\Delta t)\) and \(\rho_{\text{no-jump}}(t+\Delta t)\) factorize also. Therefore, the equation for the density matrix can be replaced by the corresponding equations for the state vectors. The state vector after time \(\Delta t\) in which a jump is recorded is given by \[ |\Psi_{\mathrm{jump}}(t+\Delta t)\rangle= \frac{1}{\sqrt{\langle\Psi(t)|C^{\dagger}C|\Psi(t)\rangle}}C|\Psi(t)\rangle\,. \] The probability of a jump occurring in the time interval \(\Delta t\) is \[ p_{\mathrm{jump}}(t)=\langle\Psi(t)|C^{\dagger}C|\Psi(t)\rangle\Delta t\,. \] If no jump occurs, the state vector evolves according to the non-Hermitian Hamiltonian \(H_{\mathrm{eff}}\), \[ |\Psi_{\text{no-jump}}(t+\Delta t)\rangle =\frac{1}{\sqrt{\langle\Psi(t)|\left(1+\frac{\mathrm{i}}{\hbar}(H_{\mathrm{eff}}^{\dagger}-H_{\mathrm{eff}})\Delta t\right)|\Psi(t)\rangle}} \left(1-\frac{\mathrm{i}}{\hbar}H_{\mathrm{eff}}\Delta t\right)|\Psi(t)\rangle\,. \]
Numerical simulation takes place over discrete time with time step \(\Delta t\). When the wave function \(|\Psi(t_n)\rangle\) is given, the wave function \(|\Psi(t_{n+1})\rangle\) is determined by the following algorithm [6]: