The Quantum Jump Method¶
The quantum jump, or Monte Carlo wave function (MCWF) method [13, 14] provides an alternative to solving the Lindblad equation of an open quantum system. It relies on propagating stochastic trajectories \(\ket{\Psi_n(t)}\) in Hilbert space. For a sufficient number of trajectories, the density matrix can be constructed as
More importantly, expectation values can be calculated as the “ensemble average”
with
The advantage of the MCWF method is that the trajectories can be calculated independently in parallel. Also, additional trajectories can be added at a later time to improve precision. Using MPI, the MCWF method can scale to system sizes where a full density-matrix propagation is no longer be feasible.
Algorithm¶
The propagation of a single trajectory proceeds as follows [15, 16]
Define an effective Hamiltonian
\[H_{\text{eff}} = \Op{H} - \frac{i\hbar}{2} \sum_i \Op{L}_i^{\dagger} \Op{L}_i\]operators, and \(\Op{H}\) is the Hamiltonian for the unitary evolution. The decay term that is being added to the Hamiltonian causes a reduction in the norm of the propagated wave function, and is used to track the probability of possible quantum jumps.
Draw a random number \(r \in [0, 1)\), propagate until \(\Braket{\Psi(t)}{\Psi(t)} = r\). The time at which this occurs is the jump time \(t_j\). It is important determine \(t_j\) as exactly as possible.
Apply an instantaneous quantum jump \(\Ket{\Psi(t_j)} \rightarrow \Op{L}_n \Ket{\Psi(t_j)}\) using one of the Lindblad operators. The Lindblad operator \(\Op{L}_n\) that is used is randomly chosen from the full set of Lindblad operators, with a relative probability of \(\langle \Psi(t_j) \vert \Op{L}_n^{\dagger} \Op{L}_n \vert \Psi(t_j) \rangle\).
After the jump, normalize \(\ket{\Psi(t_j)}\), draw a new random number \(r \in [0, 1)\) and continue the propagation.
At any point (and especially at final time \(T\)), we may normalize the wave function \(\Psi(t)\) and draw a new random number \(r \in [0, 1)\), even if no quantum jump has occurred. We should do this whenever we need the propagated state \(\ket{\Psi(t)}\) directly, for example to calculate expectation values. We may even do this at every time step (or rather, plotting step). However, in this case we must still identify the exact jump time \(t_j\) within a particular time step, and apply the jump instantaneously.
Setting up an MCWF propagation in QDYN¶
The MCWF method is set up entirely through the config file. The config
file must define all the Lindblad operators and ensure that they are
stored as-is, instead of being converted into a super-operator
(conv_to_superop = F). Secondly, all observables of interest should
be defined. Lastly, and most importantly, the flag use_mcwf = T must
be set in the prop section. For example, the following config file
would be appropriate to simulate spontaneous decay in a two-level
system:
tgrid: t_start = 0, t_stop = 10_iu, nt = 100
prop: method = exact, use_mcwf = T
ham:
* type = matrix, n_surf = 2, filename = H0.dat, specrad_method = diag, &
sparsity_model = banded, op_type = potential
dissipator:
* type = lindblad_ops, sparsity_model = banded, conv_to_superop = F, &
filename = L1.dat, add_to_H_jump = banded
psi:
* type = file, filename = psi0.dat
observables:
* outfile = pop.dat, exp_unit = unitless, is_real = T, time_unit = iu, &
column_label = <P_1>, type = matrix, n_surf = 2, sparsity_model = banded, &
filename = O1.dat
The file L1.dat contains the jump operator
\(\sqrt{\gamma}\Ket{0}\Bra{1}\). Note the setting add_to_H_jump =
banded. If given, it instructs QDYN to automatically add the decay term to
the Hamiltonian described in the ham section to obtain
\(\Op{H}_{\text{eff}}\). The value of add_to_H_jump
defined the sparsity model for the extra term (it should take the same value
for all Lindblad operators). If not given, one may alternatively set up the
effective Hamiltonian “manually” in the ham section. This may allow for
greater efficiency. For example, the file H0.dat could be set up to
directly contain the decay term, as its imaginary part. Obviously, the
effective Hamiltonian must be consistent with the defined Lindblad operators,
or mayhem will ensue. Forgetting to set add_to_H_jump
(and not setting up the effective Hamiltonian manually) has no catastrophic
consequences: the evolution will simply be unitary, as if use_mcwf = F.
Internally, the use_mcwf flag is taken into account by the
prop() routine, which implements the algorithm outlined above. Since the
effective Hamiltonian is non-Hermitian, and also because determining the exact
jump time requires the ability to propagate with a variable time step, the
cheby propagation method cannot be used in conjunction with MCWF. Note that
QDYN includes the qdyn_prop_traj utility, which may be used to propagate of
an arbitrary number of trajectories, and obtain ensemble averages for
observables, from an arbitrary given config file. When QDYN is compile with MPI
support (./configure --use-mpi), the qdyn_prop_traj program can
calculate the different trajectories in parallel.
Estimating the statistical error¶
For a reasonable number of trajectories (“reasonable” being generally still greater than 1000), any calculated expectation values have a non-negligible error, see chapter 7 of Ref. [14].
The statistical error of an expectation value \(\Avg{\Op{O}(t)}\) obtained from \(N\) trajectories can be estimated as
where the statistical variance is
The variance typically converges quickly and is approximately independent of \(N\), so that the overall statistical error \(\delta\Avg{\Op{O}(t)}\) scales as \(1/\sqrt{N}\). The statistical error (or a multiple of it) should be included as error bars in any plot of the averaged expectation value \(\Avg{\Op{O}(t)}\).
The averaged statistical squares
\(\frac{1}{N} \sum\limits_{n=1}^{N} \left(\Avg{\Op{O}(t)}_n\right)^2\)
can be obtained by giving the --write-statistical-squares to the
qdyn_prop_traj utility. In this case, for every expectation value
output file, an additional file is written with the same name but an
appendix .statsq that contains
\(\sqrt{\frac{1}{N}\sum\limits_{n=1}^{N} \left(\Avg{\Op{O}(t)}_n\right)^2}\);
the square-root is solely so that the columns are in the same units as
the original expectation values.
Alternatively, an upper bound for the statistical variance can be obtained from the quantum variance
The expectation values of the squares of the observables,
\(\Avg{\Op{O}^2(t)}\), must be configured in the observables
section of the config file, e.g.
observables: is_real = T, outfile = exp_values.dat
* type = r, column_label='<r>', square='<r^2>'
* type = mom, column_label='<p>', square='<p^2>'
In some cases (when the system is in an eigenstate of \(\Op{O}\) at all times) the quantum variance and the statistical variance match exactly.