Lecture 13: Resonance flourescence¶

Author: J. R. Johansson (robert@riken.jp), http://dml.riken.jp/~rob/

The latest version of this IPython notebook lecture is available at http://github.com/jrjohansson/qutip-lectures.

The other notebooks in this lecture series are indexed at http://jrjohansson.github.com.

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

In [2]:

from qutip import *

Introduction¶

$\displaystyle H_L = -\frac{\Omega}{2}(\sigma_+ + \sigma_-)$

$\displaystyle \frac{d}{dt}\rho = -i[H_L, \rho] + \gamma_0(N+1)\left(\sigma_-\rho(t)\sigma_+ - \frac{1}{2}\sigma_+\sigma_-\rho(t) - \frac{1}{2}\rho(t)\sigma_+\sigma_-\right) + \gamma_0 N \left(\sigma_+\rho(t)\sigma_- - \frac{1}{2}\sigma_-\sigma_+\rho(t) - \frac{1}{2}\rho(t)\sigma_-\sigma_+\right)$

Problem definition in QuTiP¶

In [3]:

Omega = 1.0 * 2 * pi

In [4]:

gamma0 = 0.05
w_th = 0.0 
N = n_thermal(Omega, w_th)

In [5]:

def system_spec(Omega, gamma0, N):
    HL = -0.5 * Omega * (sigmap() + sigmam())
    c_ops = [sqrt(gamma0 * (N + 1)) * sigmam(), sqrt(gamma0 * N) * sigmap()]
    return HL, c_ops

In [6]:

HL, c_ops = system_spec(Omega, gamma0, N)

In [7]:

e_ops = [sigmax(), sigmay(), sigmaz(), sigmam(), sigmap(), num(2)]

In [8]:

psi0 = basis(2, 0)

In [9]:

tlist = np.linspace(0, 20/(2*pi), 200)
result = mesolve(HL, psi0, tlist, c_ops, e_ops)

In [10]:

fig, axes = plt.subplots(2, 1, figsize=(12, 6), sharex=True)

axes[0].plot(result.times, result.expect[0], 'r', label=r'$\langle\sigma_x\rangle$')
axes[0].plot(result.times, result.expect[1], 'g', label=r'$\langle\sigma_y\rangle$')
axes[0].plot(result.times, result.expect[2], 'b', label=r'$\langle\sigma_z\rangle$')
axes[0].legend()
axes[0].set_ylim(-1, 1);


axes[1].plot(result.times, result.expect[5], 'b', label=r'$P_e$')

#axes[1].set_ylabel(r'$\langle\sigma_z\rangle$', fontsize=16)
axes[1].set_xlabel("time", fontsize=16)
axes[1].legend()
axes[1].set_ylim(0, 1);

In [11]:

fig, ax = plt.subplots(1, 1, figsize=(12, 6), sharex=True)

for idx, gamma0 in enumerate([0.1 * Omega, 0.5 * Omega, 1.0 * Omega]):

    HL, c_ops = system_spec(Omega, gamma0, N)
    result = mesolve(HL, psi0, tlist, c_ops, e_ops)

    ax.plot(result.times, result.expect[5], 'b', label=r'$\langle\sigma_z\rangle$')

ax.set_ylim(0, 1);

In [12]:

fig, ax = plt.subplots(1, 1, figsize=(12, 6), sharex=True)

for idx, gamma0 in enumerate([0.1 * Omega, 0.5 * Omega, 1.0 * Omega]):

    HL, c_ops = system_spec(Omega, gamma0, N)
    result = mesolve(HL, psi0, tlist, c_ops, e_ops)

    ax.plot(result.times, imag(result.expect[4]), label=r'im $\langle\sigma_+\rangle$')

ax.set_ylim(-.5, 0.5);

In [13]:

fig, axes = plt.subplots(1, 2, figsize=(12, 6))

taulist = np.linspace(0, 100, 10000)

for idx, gamma0 in enumerate([2 * Omega, 0.5 * Omega, 0.25 * Omega]):

    HL, c_ops = system_spec(Omega, gamma0, N)
    corr_vec = correlation_2op_1t(HL, None, taulist, c_ops, sigmap(), sigmam())
    w, S = spectrum_correlation_fft(taulist, corr_vec)

    axes[0].plot(taulist, corr_vec, label=r'$<\sigma_+(\tau)\sigma_-(0)>$')
    axes[1].plot(-w / (gamma0), S, 'b', label=r'$S(\omega)$')
    axes[1].plot( w / (gamma0), S, 'b', label=r'$S(\omega)$')

axes[0].set_xlim(0, 10)
axes[1].set_xlim(-5, 5);

/usr/lib/python3/dist-packages/numpy/core/numeric.py:460: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Software versions¶

In [14]:

from qutip.ipynbtools import version_table; version_table()

Out[14]:

Software	Version
Cython	0.20.1post0
OS	posix [linux]
Numpy	1.8.1
QuTiP	3.0.0.dev-5a88aa8
Python	3.4.1 (default, Jun 9 2014, 17:34:49) [GCC 4.8.3]
matplotlib	1.3.1
IPython	2.0.0
SciPy	0.13.3
Thu Jun 26 14:11:39 2014 JST