Numerical Simulation of Stochastic Master Equation using Stochastic Schrödinger Equation (Wave Function Monte Carlo)

Question

Consider a time independent system coupled to a Markovian bath, the equation of motion for the density matrix of the system has to take the form

\begin{equation} \dot{\rho} = - i \left[H,\rho\right] - \sum_m \left(c_m^{\dagger}c_m\rho+ \rho c_m^{\dagger}c_m - 2 c_m \rho c_m^{\dagger}\right). \end{equation}

I will call this a Lindblad master equation (LME). This equation describes the full deterministic dynamics of the density matrix, including decoherence and dissipation.

Alternatively, one can follow the trajectory of a pure state by means of a stochastic Schrödinger equation (SSE). The LME can be obtained by taking the ensemble average of the SSE over noise realizations. Different types of noise in an SSE can correspond to the same LME. For instance, the SSE

\begin{equation} \textrm{d}\lvert \psi \rangle = \left[-i H \textrm{d}t - \sum_m \left(c_m^{\dagger}c_m - \langle c_m^{\dagger}c_m \rangle \right)\textrm{d}t +\sum_m \xi_m\left(\frac{c_m}{\sqrt{\langle c^{\dagger}_m c_m \rangle}}-1\right) \right]\lvert\psi\rangle \end{equation}

describes a possible unraveling of the above LME. Here the $\xi_m$ are binary random increments that satisfy $\xi_m\xi_n = \delta_{nm}\xi_m$ and $\langle\langle\xi_m \rangle\rangle = 2\langle c^{\dagger}_m c_m \rangle \textrm{d}t$, where the double brackets indicate ensemble average over noise realizations.

Numerical simulation of an ensemble of SSEs can be advantageous over simulation of an LME as the number of entries of the density matrix scales as $N^2$ where $N$ is the dimension of the Hilbert space.

$\textbf{Question}$: if the system is continuously monitored AND damped, the evolution of the density matrix will become conditioned on the (random) measurement record and correspondingly a stochastic term is added to the master equation (making it a stochastic master equation, i.e. SME) which encodes the back-action of the measurement on the system. Can an equivalent formulation in terms of an SSE can be found? Are there any caveats? It seems to me that this must be the case as generalized measurements and decoherence are closely related.

If this is the case, there will now be two types of stochastic terms in the SSE, describing the damping and the measurement respectively. Can the SME be simulated by averaging the SSE only with respect to the "damping noise"?

In particular, in Kurt Jacobs' book on quantum measurement theory in chapter 4.3.4 there is a section titled "Monte Carlo method for Stochastic Master Equations" which describes an algorithm that involves evolving the Schmidt coefficients of the density matrix as well. Why does this need to be done in the case of continuous measurement AND damping while it is not necessary when considering either of the two?

Answer 1

Before tackling the OP's questions, let us first quickly establish some conventions and notation. A stochastic master equation (SME) describing quantum-jump trajectories can be written in the form $$ \mathrm{d}\rho = \mathrm{d}t \mathcal{L}_0\rho - {\rm d}t\sum_m \mathcal{H}[ \tfrac{1}{2}M^\dagger_m M_m]\rho + \sum_m {\rm d}\mu_m \mathcal{G}[M_m]\rho. \qquad (1)$$ Here, $\rho$ is the quantum state conditioned on a given sequence of measurement outcomes and we defined the non-linear superoperators $\mathcal{H}[L]\rho = L\rho + \rho L^\dagger - {\rm Tr}[(L+L^\dagger)\rho]$ and $\mathcal{G}[L]\rho = L\rho L^\dagger /{\rm Tr}[L^\dagger L\rho] - \rho$, while $\mathcal{L}_0$ is any Lindblad generator of the form $$ \mathcal{L}_0 = \mathcal{H}[-{\rm i}H] + \sum_n\mathcal{D}[N_n],$$ with $H$ a Hamiltonian ($\hbar=1$) and $\mathcal{D}[L] \rho = L\rho L^\dagger - \tfrac{1}{2}\{L^\dagger L,\rho\}$ a dissipator with Lindblad operator $L$. Finally, we have the stochastic Poisson increments ${\rm d} \mu_m = {\rm d}\mu_m^2,$ which are binary random variables (either 0 or 1) with conditional expectation value $$ \mathbb{E}[{\rm d}\mu_m|\rho] = {\rm d}t \,{\rm Tr}[M^\dagger_m M_m\rho],$$ and which are statistically independent, thus satisfying ${\rm d} \mu_m{\rm d} \mu_n = \delta_{mn}{\rm d} \mu_m$.

The particular sequence of jumps ${\rm d}\mu_m$ gives rise to a measurement record $ \mathbf{\mu}_m(t)$ for each dissipation channel that is monitored (by direct detection, e.g. photon counting). For example, an experiment in which the fluorescence from a two-level atom with spontaneous emission rate $\Gamma$ is monitored with detection efficiency $\eta$ is described by a stochastic master equation with one monitored channel $M = \sqrt{\eta\Gamma}\sigma^-$ and one unmonitored channel $N = \sqrt{(1-\eta)\Gamma}\sigma^-$ (describing decoherence due to the undetected photons). The measurement record is the sequence of photodetector clicks obtained in a given realisation of the experiment.

Averaging the SME (1) over the measurement record according to the jump statistics above, one obtains a Lindblad equation for the ensemble-averaged state $\bar{\rho} = \mathbb{E}[\rho]$ given by $$ \frac{{\rm d}\bar{\rho}}{{\rm d}t} = \mathcal{L}_0 \bar{\rho} + \sum_m\mathcal{D}[M_m]\bar{\rho}.$$ For more mathematical details, see the textbook by Wiseman & Milburn.

The OP first asks whether it is possible to unravel the master equation (1) into pure-state trajectories in this way, and indeed the answer is affirmative. If we unravel the generator $\mathcal{L}_0$, we obtain a SME $${\rm d}\rho = {\rm d}t\mathcal{H}[-{\rm i}H]\rho - {\rm d}t\sum_m \mathcal{H}[ \tfrac{1}{2}M^\dagger_m M_m]\rho - {\rm d}t\sum_n \mathcal{H}[ \tfrac{1}{2}N^\dagger_n N_n]\rho+ \sum_m {\rm d}\mu_m \mathcal{G}[M_m]\rho + \sum_n {\rm d}\nu_n \mathcal{G}[N_n]\rho, $$ where the stochastic Poisson increments ${\rm d}\nu_n = {\rm d}\nu_n^2$ obey the statistics $$ \mathbb{E}[{\rm d}\nu_n|\rho] = {\rm d}t \,{\rm Tr}[N^\dagger_n N_n\rho],$$ and are statistically independent from each other and from the ${\rm d}\mu_m$.

This evolution preserves purity and can therefore be also written as a stochastic Schroedinger equation (SSE) for (conditional) pure states $$ {\rm d}|\psi\rangle = -{\rm i}H|\psi\rangle {\rm d}t - \tfrac{1}{2}\sum_l \left(L^\dagger_l L_l - \langle \psi\lvert L^\dagger_l L_l|\psi\rangle \right)|\psi\rangle {\rm d}t + \sum_l{\rm d}\lambda_l \left( \frac{L_l}{\sqrt{\langle\psi|L^\dagger_l L_l|\psi\rangle}} -1 \right)|\psi\rangle, \qquad (2)$$ where all the jump operators $L_l = \{M_m,N_n\}$ and stochastic increments ${\rm d}\lambda_l = \{{\rm d}\mu_m,{\rm d}\nu_n\}$ have been grouped together for brevity. The SSE can be interpreted as a description of an efficient continuous measurement process on the environment responsible for the dissipation channels $N_n$, in addition to the environment responsible for $M_m$. In our example above, it simply corresponds to a perfectly efficient fluorescence measurement on our two-level atom, in which no photon goes undetected.

Now, the OP asks if one can recover the first SME by averaging over the "measurement record" $\nu_n$ corresponding only to the jumps $N_n$ occurring in the dissipation channels. The answer is yes, formally speaking, without caveats, since the stochastic increments ${\rm d}\nu_n $ are independent from the ${\rm d}\mu_m$. However, in practice this would be a very challenging simulation. The reason for this is that the SME describes a different evolution for each possible realisation of the measurement record(s) $\mu_m$. In order to recover the SME (1) for a given $\mu_m$, one needs to average Eq. (2) over many trajectories in which exactly the same sequence of jumps ${\rm d}\mu_m$ occurred. However, it is extremely unlikely to obtain two trajectories with exactly (or near-exactly) the same sequence of jumps. Therefore I expect that a very (i.e. prohibitively) large number of trajectories of Eq. (2) will be needed in order to recover Eq. (1) for a given $\mu$.

Answer 2

I will try to answer to the question

if the system is continuously monitored AND damped, the evolution of the density matrix will become conditioned on the (random) measurement. Can an > equivalent formulation in terms of an SSE can be found? Are there any caveats?

First of all, you are not only talking about continuously monitoring the quantum system (which is what the environment is doing in the Lindblad Master equation) but you are talking about introducing some feedback mechanism conditioned on the result of your monitoring. When you do this, you end up with a Master equation that is of the same Lindblad form (with different loss super-operators), and not a Stochastic master equation. Let me sketch the reasoning in the following:

From SSE to Lindblad form and back

Let me start by saying that a master equation of the form $$ \dot{\rho} = \mathcal{L} \rho$$ encodes the average evolution of a system. In a general non-equilibrium Quantum System we can suppose that the mixedness of the density matrix is the result of an average of several possible pure states $|\psi_i\rangle$ each with it's probability $p_i$. The time-evolution of those pure states can then be unravelled into the stochastic schroedinger equation, which has the general form $$ d|\psi_i(t)\rangle = f(|\psi_i(t)\rangle, t)dt + g(|\psi_i(t)\rangle, t)dW(t) $$ where dW(t) is a Wiener-Process normally giving rise to a delta-correlated noise. This is formally equivalent to the second equation you wrote ( which I believe has a couple of typos).

If you define the pure-state density matrix $\rho_i = |\psi_i\rangle\langle\psi_i|$, then you can rewrite the SSE into a stochastic equation for this pure-state density matrix, obtaining $$ d\rho_i(t) = \mathcal{L}\rho_i dt + g(\rho_i(t), t)dW(t) $$

From this equation you can recover the Lindblad ME by simply summing over all possible pure states, as

• if you sum over all possible pure states you have $\rho=\sum_i\rho_i$;
• $\langle dW\rangle = \sum_i dW = 0$ because it is a Wiener process.

Active feedback

So, if you want to introduce feedback in your system, you usually do that by saying that you apply a certain operator $\hat{B}$ at time $t$ conditioned on something you measured at $t-\tau$, where $\tau$ represents the delay of the back action (and you then might want to see what happens if $\tau\rightarrow0$).

Essentially, you can show that if you add the feedback mechanism to your SSE you obtain an equation of the form:

$$ d\rho_i(t) = \mathcal{L}\rho_i dt + (g(\rho_i(t), t)+b(\rho_i(t), t))dW(t) + h(\rho_i(t), t)dW(t)dW(t-\tau) $$

where $\mathcal{L}$ is the Liouvillian and $g$ the stochastic effect of the losses on the system with no feedback, while b and h are due to the feedback. The proof to get to this formula is quite elaborate, but you should also find it in the book you mentioned.

Thanks to the rules of stochastic calculus $dW(t)dW(t-\tau)\propto dt$, so that you get a new Stochastic Equation for pure states

$$ d\rho_i(t) = (\mathcal{L}\rho_i + h(\rho_i(t), t))dt + g(\rho_i(t), t)dW(t)$$

which again, if you average across all pure states, will give you a Lindblad Master Equation with a modified Lindbladian, proving that adding feedback to a system is very similar to coupling it to a particular environment.

An active research area in Quantum System is exactly that: how can we engineer dissipation (the environment) so that it's effect is that of performing a feedback (or feedforward) onto the system.