C++QEDCore  v2 Milestone 10
a framework for simulating open quantum dynamics – core
Description of the Master-equation evolution

Table of Contents

The Master equation for the evolution of a density operator reads

\[\dot\rho=\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}\lp\HnH\,\rho-\rho\HnH^\dag\rp+\sum_mJ_m\rho{J_m^\dag}=2\Re\lbr\frac\HnH{i\hbar}\rho\rbr+\sum_mJ_m\lp{J_m\rho}\rp^\dag,\]

where, after the second equality, the non-Hermitian “Hamiltonian” is defined as (cf. Description of the MCWF method)

\[\HnH=H-\frac{i\hbar}2\sum_m J^\dag_m J_m.\]

Since both the Hamiltonian and Lindblad operators are expressed in a finite discrete basis, this is simply an ODE evolution.

Since each Hamiltonian and/or lossy physical system in the framework provides services to calculate the effect of $\HnH/(i\hbar)$ (cf. structure::Hamiltonian::addContribution) and the jumps $J_m$ (cf. structure::Liouvillean::actWithJ) on state vectors, these being indispensable for the Monte Carlo wave-function method, no separate services are needed to calculate Master-equation evolution for these systems.

The steps are as follows:

  1. Iterate through the columns of the (in general, multi-)matrix $\rho$ to get $\HnH\,\rho/(i\hbar)$.
  2. Calculate two times the real part of this matrix in place.
  3. For each $m$:
    1. Take a temporary copy of $\rho$ and iterate through its columns to get $J_m\rho$.
    2. Take the Hermitian conjugate of this matrix in place, and repeat step (a).
    3. Add the contribution obtained in steps (a)-(b)-(a) to the matrix obtained in step 2.

The drawback of this method is that it necessitates a temporary copy of the full density operator. Alternatively, systems may provide specialized code to calculate $J_m\rho J_m^\dag$ directly, hence circumventing steps (a)-(b)-(a).

Todo:
Implement this possibility (cf. also this tracker)

Limitations

Consider these notes, whose notation we adopt. It is easy to see that a density-operator evolution of the form

\[\dot\rho_T=2\Re\lbr\mathcal{A}^T\rho_T\rbr+\sum_mJ_m^T\rho_T\lp J_m^T\rp^\dag\]

can always be unravelled up to first order in $\delta t$ via an MCWF-like two-step stochastic process as:

  1. $\ket{\Psi_T(t+\delta t)}=\frac{1+\delta t\mathcal{A}^T}{\sqrt{1-\delta p}}\ket{\Psi_T(t)}$, with probability $1-\delta p$.
  2. $\ket{\Psi_T(t+\delta t)}=\frac{\sqrt{\delta t}J_m^T}{\sqrt{\delta p_m}}\ket{\Psi_T(t)}$, with probability $\delta p_m$ for either $m$, such that $\delta p=\sum_m\delta p_m$.

If we require that the second possibility do not alter the norm, as in the case of the MCWF method (the first may in the case when $\mathcal{A}^T$ is non-Hermitian), that is, if the norm may change only in time (jumps considered instantaneous), then we get the further requirement that $\delta p_m=\bra{\Psi_T(t)}\lp J_m^T\rp^\dag J_m^T\ket{\Psi_T(t)}$.

This shows that the Master driver will only work correctly if

\[J_m^T\rho\lp J_m^T\rp^\dag=J_m\rho J_m^\dag\quad\forall\rho,\;m,\]

because structure::Liouvillean::actWithJ calculates the latter, while quantumtrajectory::Master expects the former.

See also
structure::ExactCommon::applicableInMaster

An example

The requirement for the applicability of the quantumtrajectory::Master driver might seem formidable, but the situation is actually not so bad. First of all, the problem does not arise at all for systems that do not inherit from structure::Exact.

Consider a harmonic-oscillator mode in some rotating frame, as in this Section. Here, the transformation operator reads $T=e^{-z\,a^\dag a\,t}$, and the transformed jump operators: $a^T=e^{-zt},\;\lp a^T\rp^\dag=\lp a^\dag\rp^T=e^{-z^*t}$, which makes that in the case of a purely imaginary $z$, which is the case for a purely rotating (and not shrinking) frame, the requirement is fulfilled.