013_nonmarkov-transfer-tensor-method.md

jupyter
jupytext kernelspec text_representation extension format_name format_version jupytext_version .md markdown 1.3 1.14.5 display_name language name Python 3 (ipykernel) python python3

QuTiP Example: The Transfer Tensor Method for Non-Markovian Open Quantum Systems

Arne L. Grimsmo
Université de Sherbrooke
[email protected] $\newcommand{\ket}[1]{\left|#1\right\rangle}$ $\newcommand{\bra}[1]{\left\langle#1\right|}$

Introduction

The "Transfer Tensor Method" was introduced by Cerrillo and Cao in Phys. Rev. Lett 112, 110401 (2014) (arXiv link), as a general method for evolving non-Markovian open quantum systems.

The method takes as input a set of dynamical maps $\mathcal{E}_k$, such that

$$ \rho(t_k) = \mathcal{E}_k \rho(0) $$

for an intial set of times $t_k$. This set of dynamical maps could be the result of experimental process tomography of they could be precomputed through some other (typically costly) method. The idea is that based on knowledge of these maps, one can try to exptrapolate the, in general non-Markovian, time-evolution to larger times, $t_n > t_k$. The method assumes that there is no explicit time-dependence in the total system-bath Hamiltonian.

Preamble

Imports

import numpy as np
import qutip as qt
import qutip.solver.nonmarkov.transfertensor as ttm
from qutip.ipynbtools import version_table

Plotting Support

%matplotlib inline
import matplotlib.pyplot as plt

Jaynes-Cummings model, with the cavity as a non-Markovian bath

As a simple example, we consider the Jaynes-Cummings mode, and the non-Markovian dynamics of the qubit when the cavity is traced out. In this example, the dynamical maps $\mathcal{E}_k$ are the reduced time-propagators for the qubit, after evolving and tracing out the cavity, i.e.

$$ \mathcal{E}k \rho = {\rm tr}{\rm cav} \left[ {\rm e}^{\mathcal{L} t_k} \rho \otimes \rho_{0,{\rm cav}} \right], $$

where $\mathcal{L}$ is the Lindbladian for the dissipative JC model (defined below) and $\rho_{0,{\rm cav}}$ is the initial state of the cavity.

Problem setup

kappa = 1.0  # cavity decay rate
wc = 0.0 * kappa  # cavity frequency
wa = 0.0 * kappa  # qubit frequency
g = 10.0 * kappa  # coupling strength
N = 3  # size of cavity basis

# intial state
psi0c = qt.basis(N, 0)
rho0c = qt.ket2dm(psi0c)
rho0a = qt.ket2dm(qt.basis(2, 0))
rho0 = qt.tensor(rho0a, rho0c)
rho0avec = qt.operator_to_vector(rho0a)

# identity superoperator
Id = qt.tensor(qt.qeye(2), qt.qeye(N))
E0 = qt.sprepost(Id, Id)

# partial trace over the cavity, reprsented as a superoperator
ptracesuper = qt.tensor_contract(E0, (1, 3))

# intial state of the cavity, represented as a superoperator
superrho0cav = qt.sprepost(
    qt.tensor(qt.qeye(2), psi0c), qt.tensor(qt.qeye(2), psi0c.dag())
)

# operators
a = qt.tensor(qt.qeye(2), qt.destroy(N))
sm = qt.tensor(qt.sigmam(), qt.qeye(N))
sz = qt.tensor(qt.sigmaz(), qt.qeye(N))

# Hamiltonian
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [np.sqrt(kappa) * a]

Exact timepropagators to learn from

The function dynmap generates an exact timepropagator for the qubit $\mathcal{E}_{k}$ for a time $t_k$.

prop = qt.Propagator(qt.liouvillian(H, c_ops))


def dynmap(t):
    # reduced dynamical map for the qubit at time t
    return ptracesuper @ prop(t) @ superrho0cav


dynmap(1.0)

Exact time evolution using standard mesolve method

exacttimes = np.arange(0, 5, 0.02)
exactsol = qt.mesolve(H, rho0, exacttimes, c_ops, [sz])

Approximate solution using the Transfer Tensor Method for different learning times

times = np.arange(0, 5, 0.05)  # total extrapolation time
ttmsols = []
maxlearningtimes = [0.5, 2.0]  # maximal learning times
for T in maxlearningtimes:
    learningtimes = np.arange(0, T, 0.05)
    # generate exact dynamical maps to learn from
    learningmaps = [
        dynmap(t) for t in learningtimes
    ]
    # extrapolate using TTM
    ttmsols.append(ttm.ttmsolve(learningmaps, rho0a, times))

Visualize results

fig, ax = plt.subplots(figsize=(10, 7))
ax.plot(exactsol.times, exactsol.expect[0], "-b", linewidth=3.0)
style = ["og", "or"]
for i, ttmsol in enumerate(ttmsols):
    ax.plot(
        ttmsol.times,
        qt.expect(qt.sigmaz(), ttmsol.states),
        style[i],
        linewidth=1.5,
    )
ax.legend(["exact", str(maxlearningtimes[0]), str(maxlearningtimes[1])])
ax.set_xlabel(r"$\kappa t$", fontsize=20)
ax.set_ylabel(r"$\sigma_z$", fontsize=20)

Discussion

The figure above illustrates how the transfer tensor method needs a sufficiently long set of learning times to get good results. The green dots show results for learning times $t_k=0,0.1,\dots,0.5$, which is clearly not sufficient. The red dots show results for $t_k=0,0.1,\dots,2.0$, which gives results that are in very good agreement with the exact solution.

Epilouge

version_table()