QuTiP example: Bloch-Redfield Master Equation

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import numpy as np
In [3]:
from qutip import *

Two-level system

In [4]:
delta = 0.0 * 2 * np.pi
epsilon = 0.5 * 2 * np.pi
gamma = 0.25
times = np.linspace(0, 10, 100)
In [5]:
H = delta/2 * sigmax() + epsilon/2 * sigmaz()
H
Out[5]:
Quantum object: dims = [[2], [2]], shape = [2, 2], type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.571 & 0.0\\0.0 & -1.571\\\end{array}\right)\end{equation*}
In [6]:
psi0 = (2 * basis(2, 0) + basis(2, 1)).unit()
In [7]:
c_ops = [np.sqrt(gamma) * sigmam()]
a_ops = [sigmax()]
In [8]:
e_ops = [sigmax(), sigmay(), sigmaz()]
In [9]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [10]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, spectra_cb=[lambda w : gamma * (w > 0)])
In [11]:
plot_expectation_values([result_me, result_brme]);
/home/rob/py-envs/py3-stable/lib/python3.4/site-packages/numpy/core/numeric.py:462: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
In [12]:
b = Bloch()
b.add_points(result_me.expect, meth='l')
b.add_points(result_brme.expect, meth='l')
b.make_sphere()

Harmonic oscillator

In [13]:
N = 10

w0 = 1.0 * 2 * np.pi
g = 0.05 * w0
kappa = 0.15

times = np.linspace(0, 25, 1000)
In [14]:
a = destroy(N)
In [15]:
H = w0 * a.dag() * a + g * (a + a.dag())
In [16]:
# start in a superposition state
psi0 = ket2dm((basis(N, 4) + basis(N, 2) + basis(N,0)).unit())
In [17]:
c_ops = [np.sqrt(kappa) * a]
a_ops = [a + a.dag()]
In [18]:
e_ops = [a.dag() * a, a + a.dag()]

Zero temperature

In [19]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [20]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, spectra_cb=[lambda w : kappa * (w > 0)])
In [21]:
plot_expectation_values([result_me, result_brme]);

Finite temperature

In [22]:
times = np.linspace(0, 25, 250)
In [23]:
n_th = 1.5
c_ops = [np.sqrt(kappa * (n_th + 1)) * a, np.sqrt(kappa * n_th) * a.dag()]
In [24]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [25]:
w_th = w0/np.log(1 + 1/n_th)

def S_w(w):
    if w >= 0:
        return (n_th + 1) * kappa
    else:
        return (n_th + 1) * kappa * np.exp(w / w_th)
In [26]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, [S_w])
In [27]:
plot_expectation_values([result_me, result_brme]);

Storing states instead of expectation values

In [28]:
result_me = mesolve(H, psi0, times, c_ops, [])
In [29]:
result_brme = brmesolve(H, psi0, times, a_ops, [], [S_w])
In [30]:
n_me = expect(a.dag() * a, result_me.states)
In [31]:
n_brme = expect(a.dag() * a, result_brme.states)
In [32]:
fig, ax = plt.subplots()

ax.plot(times, n_me, label='me')
ax.plot(times, n_brme, label='brme')
ax.legend()
ax.set_xlabel("t");

Atom-Cavity

In [33]:
N = 10
a = tensor(destroy(N), identity(2))
sm = tensor(identity(N), destroy(2))
psi0 = ket2dm(tensor(basis(N, 1), basis(2, 0)))
a_ops = [(a + a.dag())]
e_ops = [a.dag() * a, sm.dag() * sm]

Weak coupling

In [34]:
w0 = 1.0 * 2 * np.pi
g = 0.05 * 2 * np.pi
kappa = 0.05
times = np.linspace(0, 5 * 2 * np.pi / g, 1000)
In [35]:
c_ops = [np.sqrt(kappa) * a]
H = w0 * a.dag() * a + w0 * sm.dag() * sm + g * (a + a.dag()) * (sm + sm.dag())
In [36]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [37]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, spectra_cb=[lambda w : kappa*(w > 0)])
In [38]:
plot_expectation_values([result_me, result_brme]);

In the weak coupling regime there is no significant difference between the Lindblad master equation and the Bloch-Redfield master equation.

Strong coupling

In [39]:
w0 = 1.0 * 2 * np.pi
g = 0.75 * 2 * np.pi
kappa = 0.05
times = np.linspace(0, 5 * 2 * np.pi / g, 1000)
In [40]:
c_ops = [np.sqrt(kappa) * a]
H = w0 * a.dag() * a + w0 * sm.dag() * sm + g * (a + a.dag()) * (sm + sm.dag())
In [41]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [42]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, spectra_cb=[lambda w : kappa*(w > 0)])
In [43]:
plot_expectation_values([result_me, result_brme]);

In the strong coupling regime there are some corrections to the Lindblad master equation that is due to the fact system eigenstates are hybridized states with both atomic and cavity contributions.

Versions

In [44]:
from qutip.ipynbtools import version_table

version_table()
Out[44]:
SoftwareVersion
SciPy0.14.1
Cython0.21.2
IPython2.3.1
Python3.4.0 (default, Apr 11 2014, 13:05:11) [GCC 4.8.2]
QuTiP3.1.0
Numpy1.9.1
matplotlib1.4.2
OSposix [linux]
Tue Jan 13 13:09:27 2015 JST