QuTiPv5 Paper Example: Monte Carlo Solver for non-Markovian Baths¶

Authors: Maximilian Meyer-Mölleringhof (m.meyermoelleringhof@gmail.com), Paul Menczel (paul@menczel.net), Neill Lambert (nwlambert@gmail.com)

Introduction¶

When a quantum system experiences non-Markovian effects, it can no longer be described by the standard Lindblad formalism. To compute the time evolution of the density matrix, we can however use time-convolutionless (TCL) projection operators that lead to a differential equation in time-local form:

$\dot{\rho} (t) = - \dfrac{i}{\hbar} [H(t), \rho(t)] + \sum_n \gamma_n(t) \mathcal{D}_n(t) [\rho(t)]$

with

$\mathcal{D}_n[\rho(t)] = A_n \rho(t) A^\dagger_n - \dfrac{1}{2} [A^\dagger_n A_n \rho(t) + \rho(t) A_n^\dagger A_n]$.

These equations include the system Hamiltonian $H(t)$ and jump operators $A_n$. Contrary to a Lindblad equation, the coupling rates $\gamma_n(t)$ may be negative here.

In QuTiP v5, the non-Markovian Monte Carlo solver is introduced. It enables the mapping of the general master equation given above, to a Lindblad equation on the same Hilbert space. This is achieved by the introduction of the so called "influence martingale" which acts as trajectory weighting [1, 2, 3]:

$\mu (t) = \exp \left[ \alpha \int_0^t s(\tau) d \tau \right] \Pi_k \dfrac{\gamma_{n_k} (t_k)}{\Gamma_{n_k} (t_k)}$.

Here, the product runs over all jump operators on the trajectory with jump channels $n_k$ and jump times $t_k < t$. To finally arrive at the Lindblad form, the shift function

$s(t) = 2 \left| \min \{ 0, \gamma_1(t), \gamma_2(t), ... \} \right|$

is applied, such that the shifted rates $\Gamma_n (t) = \gamma_n(t) + s(t)$ are non-negative. We obtain the completely positive Lindblad equation

$\dot{\rho}'(t) = - \dfrac{i}{\hbar} [ H(t), \rho'(t) ] + \sum_n \Gamma(t) \mathcal{D}_n[\rho'(t)]$,

and $\rho'(t) = \mathbb{E} \{\ket{\psi(t)} \bra{\psi(t)}\}$ using the regular MCWF method. Here, $\ket{\psi (t)}$ are the generated trajectories and $\mathbb{E}$ is the ensemble average over the trajectories. This can furthermore be used to finally reconstruct the original states via $\rho(t) = \mathbb{E}\{\mu(t) \ket{\psi(t)} \bra{\psi(t)}\}$.

Note that, for this technique, the jump operators are required to fulfill the completeness relation $\sum_n A_n^\dagger A_n = \alpha \mathbb{1}$ for $\alpha > 0$. This condition is automatically taken care of by QuTiP's nm_mcsolve() function.

Damped Jaynes-Cumming Model¶

To illustrate the application of nm_mcsolve(), we want to look at the damped Jaynes-Cumming model. It describes a two-level atom coupled to a damped cavity mode. Such a model can be accurately model as a two-level system coupled to an environment with the power spectrum [4]:

$S(\omega) = \dfrac{\lambda \Gamma^2}{(\omega_0 - \Delta - \omega)^2 + \Gamma^2}$,

where $\lambda$ is the atom-cavity coupling strength, $\omega_0$ is the transition frequency of the atom, $\Delta$ the detuning of the cavity and $\Gamma$ the spectral width. At zero temperature and after performing the rotating wave approximation, the dynamics of the two-level atom can be described by the master equation

$\dot{\rho}(t) = \dfrac{A(t)}{2i\hbar} [ \sigma_+ \sigma_-, \rho(t) ] + \gamma (t) \mathcal{D}_- [\rho (t)]$,

with the state $\rho(t)$ in the interaction picture, the the ladder operators $\sigma_\pm$, the dissipator $\mathcal{D}_-$ for the Lindblad operator $\sigma_-$, and lastly $\gamma (t)$ and $A(t)$ being the real and imaginary parts of

$\gamma(t) + i A(t) = \dfrac{2 \lambda \Gamma \sinh (\delta t / 2)}{\delta \cosh (\delta t / 2) + (\Gamma - i \Delta) \sinh (\delta t/2)}$,

with $\delta = [(\Gamma - i \Delta)^2 - 2 \lambda \Gamma]^{1/2}$.

Comparison to other methods¶

We want to compare this computation with the standard mesolve(), the HEOM and the Bloch-Refield solver. For the mesolve() we can directly use the operators we have created for nm_mcsolve():

For the other methods we directly apply a spin-boson model and the free reservoir auto-correlation function

$C(t) = \dfrac{\lambda \Gamma}{2} e^{-i (\omega - \Delta) t - \lambda |t|}$,

which corresponds to the power spectrum defined above. We use the Hamiltonian $H = \omega_0 \sigma_+ \sigma_-$ in the Schrödinger picture and the coupling operator $Q = \sigma_+ + \sigma_-$. Here, we chose $\omega_0 \gg \Delta$ to ensure validity of the rotating wave approximation.

Firstly, the HEOM solver can directly be applied after separating the auto-correlation function into its real and imaginary parts:

Secondly, for the Bloch-Redfield solver we can directly use the power spectrum as input:

Finally, in order to compare these results with the nm_mesolve() method, we transform them to the interaction picture:

Plotting the Time Evolution¶

We choose to plot three comparisons to highlight the different solvers' computational results. In all plots, the gray areas show when $\gamma(t)$ is negative.

In these plots we notice two things. First, the result from the Bloch-Redfield equation deviates greatly from the others, showing us that non-Markovian effects are strongly influencing the dynamics. Second, in the grey areas - when $\gamma(t)$ is negative - the atom state restores coherence. In the last plot especially we see that during these times, the average influence martingale fluctuates, but is constant otherwise. Its deviation from unity tells us how well the simulation has converged.

References¶

[1] Donvil and Muratore-Ginanneschi. Nat Commun (2022).

[2] Donvil and Muratore-Ginanneschi. New J. Phys. (2023).

[3] Donvil and Muratore-Ginanneschi. Open Systems & Information Dynamics.

[4] Breuer and Petruccione The Theory of Open Quantum Systems.

[5] QuTiP 5: The Quantum Toolbox in Python

About¶

Testing¶