Demonstrate how to pass liouvillians to the mesolve

This notebook demonstrates how the mesolve solver can be used in different ways to solve the same problem, by passing a Hamiltonian and list of collapse operators, or by manually creating the Liouvillian and pass it to the solver in place of the Hamiltonian or collapse operators. These methods can be used to solve master equations that are not on standard Linblad form.

In [1]:
%pylab inline
In [2]:
from qutip import *

Visualization of the mesolve results

In [3]:
def visualize(result):
    fig, ax = subplots(1, 1, figsize=(8,5))
    
    for e in result.expect:
        ax.plot(result.times, e)

    ax.set_xlabel('Time')
    ax.set_ylabel('Occupation probability')
    ax.set_ylim(0, 1)

Parameters

In [4]:
wc = 1.0  * 2 * pi  # cavity frequency
wa = 1.0  * 2 * pi  # atom frequency
g  = 0.05 * 2 * pi  # coupling strength
kappa = 0.005       # cavity dissipation rate
gamma = 0.05        # atom dissipation rate
N = 15              # number of cavity fock states
tlist = linspace(0,25,100)

Operators and initial states

In [5]:
psi0 = tensor(basis(N,0), basis(2,1))
a  = tensor(destroy(N), qeye(2))
sm = tensor(qeye(N), destroy(2))

Standard method

In [6]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())

c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

result = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm])

visualize(result)

Passing a liouvillian instead of an Hamiltonian

The Hamiltonian and the collapse operators are used to contruct a Liouvillian, which is passed to mesolve in place of the Hamiltonian (and with no collapse operators)

In [7]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = liouvillian(H, c_ops)

result = mesolve(L, psi0, tlist, [], [a.dag() * a, sm.dag() * sm])

visualize(result)

Pass liouvillian instead of collapse operators

Works with list-function and list-string time-dependencies too:

In [8]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = liouvillian(None, c_ops)

result = mesolve(H, psi0, tlist, [[L, '1.0']], [a.dag() * a, sm.dag() * sm])

visualize(result)
In [9]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = liouvillian(None, c_ops)

result = mesolve(H, psi0, tlist, [[L, lambda t, args: 1.0]], [a.dag() * a, sm.dag() * sm])

visualize(result)

Liouvillian instead of a list of collapse operators

In [10]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = liouvillian(None, c_ops)

result = mesolve(H, psi0, tlist, L, [a.dag() * a, sm.dag() * sm])

visualize(result)

Mix normal collapse operators and Liouvillians, passed as the collapse-operator argument

In [11]:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, liouvillian(None, [sqrt(gamma) * sm])]

result = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm])

visualize(result)

Make a custom non-standard Liouvillian

In [12]:
psi0 = tensor(basis(N,0), basis(2,1))    # start with an excited atom

a  = tensor(destroy(N), qeye(2))
sm = tensor(qeye(N), destroy(2))

H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
c_ops = [sqrt(kappa) * a, sqrt(gamma) * sm]

L = sum([spre(c)*spost(c.dag()) - 0.25 * spre(c.dag()*c) - 0.5 * spost(c.dag()*c) for c in c_ops])

result = mesolve(H, psi0, tlist, L, [a.dag() * a, sm.dag() * sm])

visualize(result)

Software versions

In [13]:
from qutip.ipynbtools import version_table
version_table()
Out[13]:
SoftwareVersion
Cython0.16
SciPy0.10.1
QuTiP2.2.0.dev-793af74
Python2.7.3 (default, Sep 26 2012, 21:51:14) [GCC 4.7.2]
IPython0.13
OSposix [linux2]
Numpy1.7.0.dev-3f45eaa
matplotlib1.3.x
Mon Feb 18 16:06:16 2013