Appendices¶

Description of the MCWF method¶

The MCWF method cite{carmichael87,dalibard92,dum92,molmer93} aims at the simulation of open quantum systems based on a stochastic (“Monte Carlo”) trajectory. In terms of dimensionality, this is certainly a huge advantage as compared to solving the Master equation directly. On the other hand, stochasticity requires us to run many trajectories, but the method provides an optimal sampling of the ensemble density operator so that the relative error is inversely proportional to the number of trajectories.

The optimal sampling is achieved by evolving the state vector in two steps, one deterministic and one stochastic (quantum jump). Suppose that the Master equation of the system is of the form

(1) $\dot\rho=\frac1{i\hbar}\comm{H}\rho+\Liou\rho\equiv\frac1{i\hbar}\comm{H}\rho+\sum_m\lp J_m\rho J_m^\dag-\frac12\comm{J_m^\dag J_m}{\rho}_+\rp\equiv\frac1{i\hbar}\comm{\HnH}\rho+\sum_mJ_m\rho J_m^\dag,$

the usual form in quantum optics, and, in fact, the most general (so-called Lindblad) form. The non-Hermitian Hamiltonian is defined as

(2) $\HnH=H-\frac{i\hbar}2\sum_m J^\dag_m J_m.$

At time $t$ the system is in a state with normalised state vector $\ket{\Psi(t)}$ . To obtain the state vector at time $t+\delta t$ up to first order in $\delta t$ :

The state vector is evolved according to the nonunitary dynamics

(3) $i\hbar\frac{d\ket{\Psi}}{dt}=\HnH\ket{\Psi}$

to obtain (up to first order in $\delta t$ )

(4) $\ket{\Psi_{\text{nH}}(t+\delta t)}=\lp1-\frac{i\HnH\,\delta t}\hbar\rp \ket{\Psi(t)}.$

Since $\HnH$ is non-Hermitian, this new state vector is not normalised. The square of its norm reads

$\braket{\Psi_{\text{nH}}(t+\delta t)}{\Psi_{\text{nH}}(t+\delta t)}=\bra{\Psi(t)}\lp1+\frac{iH^\dag_{\text{nH}}\,\delta t}\hbar\rp\lp1-\frac{i\HnH\,\delta t}\hbar\rp\ket{\Psi(t)}\equiv 1-\delta p,$

where $\delta p$ reads

$\delta p&=\delta t\,\frac i\hbar \bra{\Psi(t)}\HnH-H^\dag_{\text{nH}}\ket{\Psi(t)}\equiv\sum_m\delta p_m,$

$\delta p_m&=\delta t\,\bra{\Psi(t)} J^\dag_m J_m\ket{\Psi(t)}\geq 0.$

Note that the timestep $\delta t$ should be small enough so that this first-order calculation be valid. In particular, we require that

(5) $\delta p\ll1.$
A possible quantum jump with total probability $\delta p$ . For the physical interpretation of such a jump see e.g. Refs.citep{dum92,molmer93}. We choose a random number $\epsilon$ between 0 and 1, and if $\delta p<\epsilon$ , which should mostly be the case, no jump occurs and for the new normalised state vector at $t+\delta t$ we take

(6) $\ket{\Psi(t+\delta t)}=\frac{\ket{\Psi_{\text{nH}}(t+\delta t)}}{\sqrt{1-\delta p}}.$

If $\epsilon<\delta p$ , on the other hand, a quantum jump occurs, and the new normalised state vector is chosen from among the different state vectors $J_m\ket{\Psi(t)}$ according to the probability distribution $\Pi_m=\delta p_m/\delta p$ :

$\ket{\Psi(t+\delta t)}=\sqrt{\delta t}\frac{J_m\ket{\Psi(t)}}{\sqrt{\delta p_m}}.$

Nonorthogonal basis sets¶

Formalism¶

We adopt the so-called covariant-contravariant formalism (see also cite{artacho91}), which in physics is primarily known from the theory of relativity. Assume we have a basis $\lbr\ket{i}\rbr_{i\in\mathbb{N}}$ , where the basis vectors are nonorthogonal, so that the metric tensor

$g_{ij}\equiv\braket{i}{j}$

is nondiagonal. The contravariant components of a state vector $\ket\Psi$ are then defined as the expansion coefficients

$\ket\Psi\equiv\Psi^i\ket{i},$

where we have adopted the convention that there is summation for indeces appearing twice. In this case one index has to be down, while the other one up, otherwise we have as it were a syntax error. This ensures that the result is independent of the choice of basis, as one would very well expect e.g. from the trace of an operator. The covariant components are the projections

$\Psi_i\equiv\braket{i}\Psi.$

The two are connected with the metric tensor:

$\Psi_i=g_{ij}\Psi^j,\quad\Psi^i=g^{ij}\Psi_j,$

where it is easy to show that $g^{ij}$ is the matrix inverse of the above defined $g_{ij}$ :

$g_{ik}g^{kj}\equiv g_i^j=\delta_{ij}.$

For the representation of operators assume $A$ is an operator, then

$\ket\Phi\equiv A\ket\Psi=\Psi^jA\ket{j},$

multiplying from the left by $\bra{i}$ yields

$\Phi_i=\bra iA\ket j\Psi^j\equiv A_{ij}\Psi^j.$

While the definition of the up-indeces matrix of the operator can be read from

$\Phi=\Phi^i\ket i\equiv A^{ij}\Psi_j\ket{i}=A^{ij}\braket{j}\Psi\ket{i}=\lp A^{ij}\ket{i}\bra{j}\rp\ket\Psi=A\ket\Psi.$

Hermitian conjugation¶

From the above it is easy to see that

$\lp A^\dag\rp_{ij}=\lp A_{ji}\rp^*,\;\text{and}\;\lp A^\dag\rp^{ij}=\lp A^{ji}\rp^*.$

The situation, however, is more involved with the mixed-index matrices since conjugation mixes the two kinds:

${\lp A^\dag\rp_i}^j=\lp {A^j}_i\rp^*,\;\text{and}\;{\lp A^\dag\rp^i}_j=\lp {A_j}^i\rp^*.$

Here we prove the second:

${A_j}^i=A_{jk}g^{ki}=\lp\lp A^\dag\rp_{kj}\rp^*g^{ki},$

so that

$\lp{A_j}^i\rp^*=\lp A^\dag\rp_{kj}\lp g^{ki}\rp^*=g^{ik}\lp A^\dag\rp_{kj}={\lp A^\dag\rp^i}_j,$

where we have used the Hermiticity of the metric tensor.

Implications¶

If one is to apply the above formalism for density matrices in nonorthogonal bases, one has to face some strange implications: The well-known properties of Hermiticity and unit trace are split: $\rho_{ij}$ and $\rho^{ij}$ will be Hermitian matrices, but their trace is irrelevant, in fact, it is not conserved during time evolution, while ${\rho_i}^j$ and ${\rho^i}_j$ have unit conserved trace.

This means that e.g. the quantumtrajectory::Master driver has to use either $\rho_{ij}$ or $\rho^{ij}$ (by convention, it uses the latter, which also means that quantumtrajectory::MCWF_Trajectory uses $\Psi^i$ ) to be able to calculate the evolution using the real part as in Eqref{eq:Master-MCWF} (otherwise, it is not the real part which appears there). On the other hand, for all normalizations one index has to be pulled, and this (at least in the case of density operators) has to be done in place to avoid excessive copying, which also means that the index has to be pulled back again afterwards.

The situation is even more delicate if one wants to calculate averages or jump probabilities for a subsystem. Then, as explained in Sref{sec:Jump} and Sref{sec:Displayed}, the partial density operator has to be used. Imagine that we want to calculate the expectation value of an operator $O$ acting on a given subsystem. Imagine that we are using an MCWF trajectory, so that we dispose of the expansion coefficients $\Psi^{iI}$ , where $i$ is the index (quantum number) for the given subsystem, while $I$ is a multi-index, fusing the (dummy) indeces of the other subsystems, which are all up. Now the matrix of partial density operator of the given subsystem reads

begin{equation} label{eq:NO_PartialTrace} rho^{ij}=lp{Psi^i}_Irp^*Psi^{jI}, end{equation}

and then, since we dispose of $O_{ji}$ (this is most easily obtained, since only scalar products have to be calculated), the desired expectation value is obtained as

begin{equation} avr O=rho^{ij}O_{ji}. end{equation}

This means that somehow all the indeces but the index of the given subsystem have to be pulled. Then, when we move to the next subsystem (this is what code{Composite}’s code{Display} does cf Sref{sec:Composite}) to calculate the expectation value of one of its operators, we have to pull up again its index, and pull down the index of the previous subsystem, and so on. This is tedious, and possibly also wasteful.

What happens in practice is that the code{MCWF_Trajectory} driver, when it comes to display, pulls emph{all indeces} by calling the function code{Pull} for the complete code{NonOrthogonal}, and passes both the covariant and contravariant vectors to code{Display}. The summation in Eqref{eq:NO_PartialTrace} is then performed by code{ElementDisplay}, while the code{Average} of the given subsystem works with ${\rho_i}^j$ (check definition!!!!!!) to calculate $\avr O={\rho_i}^j{O_j}^i$ , so that the additional pull has to be implemented here on $O$ . The same procedure is repeated with code{Jump}, since to obtain the jump probability, also an expectation value has to be calculated. The drivers code{Master} and code{EnsembleMCWF} working with density operators $\rho^{IJ}$ pull one (multi-)index in place to obtain ${\rho_I}^J$ , passes this for display, and afterwards pull the index back again.

A given subsystem’s code{Average} function can therefore receive four possible sets of parameters:

begin{enumerate} item for orthogonal basis and state-vector evolution only a single array of code{double}s, representing $\Psi_I=\Psi^I$ . item for nonorthogonal basis and state-vector evolution two arrays, representing $\Psi_I$ and $\Psi^I$ separately. item for orthogonal basis and density-operator evolution a single array, representing the matrix of $\rho$ . item for nonorthogonal basis and density-operator evolution a single array, representing ${\rho_I}^J$ . end{enumerate}

However, an actual implementor of such a function does not see the difference between these cases, since as explained in Chref{chap:Indexing}, they are hidden behind the code{Indexing} interface, an abstract class, which is then implemented to cover the four separate cases. In fact, if the given subsystem is represented in an orthogonal basis, it does not really have to know about the existence of nonorthogonal bases, the only thing which may come as a surprise is that the diagonal elements of ${\rho_I}^J$ are in general not real. This is why the diagonal function of code{Indexing} returns complex...

For reference we finally list the Master equation for $\rho^{ij}$ :

Some conventions¶

The Schr”odinger equation is: begin{equation} label{eq:SchEq} ihbarfrac{dketPsi}{dt}=HketPsi. end{equation} It may sound strange to call this a convention, but the position of $i$ emph{is} conventional in this equation. Of course, it has to be consistent with other postulates, in particular the canonical commutation relations, e.g. begin{equation} comm xp=ihbar. end{equation} The consistency of these two equations can be seen e.g. from Ehrenfest’s theorem.

newcommandHinter{H^{text{int}}} newcommandIAP[1]{{#1_{text{I}}}}

From Eqref{eq:SchEq} we can derive the necessary definitions for interaction picture. If the Hamiltonian is of the form begin{equation} Hequiv H^{(0)}+Hinter, end{equation} where $H^{(0)}$ is something we can solve, then the operator $U$ transforming between Schr”odinger picture and the interaction picture defined by $H^{(0)}$ begin{equation} ketPsiequiv U(t)ket{IAPPsi} end{equation} reads begin{equation} U(t)=e^{-frac ihbar H^{(0)} t}. end{equation} The Schr”odinger equation reads begin{equation} label{eq:SchEqIAP} ihbarfrac{dket{IAPPsi}}{dt}=U^{-1}(t)Hinter U(t)ket{IAPPsi} equivIAPHinter. end{equation} If $\mathcal O$ is an operator then in interaction picture it evolves as begin{equation} frac{dIAP{mathcal O}}{dt}= frac ihbarcomm{H^{(0)}}{IAP{mathcal O}}+ U^{-1}(t)frac{partialmathcal O}{partial t}U(t) end{equation}

For p look at how Sakurai does

Table Of Contents

Previous topic

Next topic

Appendices¶

Description of the MCWF method¶

Refinement of the method¶

An adaptive MCWF method¶

Exploiting interaction picture¶