Reproduce: Schmidt et al., Phys. Rev. B 82, 100507(R) (2010)

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [2]:
from qutip import *

Introduction

This notebook reproduces some of the results in S. Schmidt et al., Phys. Rev. B 82, 100507(R) (2010) using QuTiP.

Problem parameters

Here we use units where $\hbar = 1$:

In [3]:
N = 21              # number of cavity fock states
nn = N-1

wc = ones(2) * 1.00 * 2 * pi        # cavity frequency
wa = ones(2) * 1.00 * 2 * pi         # atom frequency
J  = 0.01 * 2 * pi            # coupling strength: cavity-cavity

# localization threshold is 2.8 * sqrt(nn) * J
gc = 2.8 * sqrt(nn) * J

g  = ones(2) * gc * 2.0   # coupling strength: atom-cavity
kappa = ones(2) * 0.01 * J/(2*pi)      # cavity dissipation rate
gamma = ones(2) * 0.01 * J/(2*pi)      # atom dissipation rate

n_th = 0.0         # avg number of thermal bath excitation

T = 10 / (J/(2*pi))
tlist = linspace(0, T, 1001)

Setup the operators, the Hamiltonian and initial state

In [4]:
# start with a fock state witn nn photons in cavity #1
psi0 = tensor(basis(N,nn), basis(2,0), basis(N,0), basis(2,0))
#psi0 = tensor(coherent(N,sqrt(nn)), basis(2,0), basis(N,0), basis(2,0))
In [5]:
# cavity annihilation operator
C = [tensor(destroy(N), qeye(2), qeye(N), qeye(2)), tensor(qeye(N), qeye(2), destroy(N), qeye(2))]

# atomic annihilation operator
A = [tensor(qeye(N), destroy(2), qeye(N), qeye(2)), tensor(qeye(N), qeye(2), qeye(N), destroy(2))]

# number operators
NA = array([a.dag() * a for a in A])
NC = array([c.dag() * c for c in C])

# polariton operators
P = NA + NC
In [6]:
# Jaynes-Cumming Hamiltonian for system 0
H0c = wc[0] * C[0].dag() * C[0] 
H0a = wa[0] * A[0].dag() * A[0]
H0int = g[0] * (A[0].dag() * C[0] + A[0] * C[0].dag())
H0 = H0c + H0a + H0int

# Jaynes-Cumming Hamiltonian for system 1
H1c = wc[1] * C[1].dag() * C[1]
H1a = wa[1] * A[1].dag() * A[1] 
H1int = g[1] * (A[1].dag() * C[1] + A[1] * C[1].dag())
H1 = H1c + H1a + H1int

# cavity-cavity interaction
Hint = J * (C[0].dag() * C[1] + C[0] * C[1].dag())

H = H0 + H1 + Hint
In [7]:
H
Out[7]:
$\text{Quantum object: dims = [[21, 2, 21, 2], [21, 2, 21, 2]], shape = [1764, 1764], type = oper, isHerm = True}\\[1em]\begin{pmatrix}0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 6.283 & 1.574 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 1.574 & 6.283 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 12.566 & 2.225 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 2.225 & 12.566 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 251.327 & 6.859 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 6.859 & 251.327 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 257.611 & 7.037 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 7.037 & 257.611 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 263.894\\\end{pmatrix}$
In [8]:
energy_level_diagram([H0c + H1c, H0a + H1a, H0int + H1int, Hint], 23);

Create a list of collapse operators that describe the dissipation

In [9]:
c_ops = []

# cavity relaxation
for i in [0,1]:
    rate = kappa[i] * (1 + n_th)
    if rate > 0.0:
        c_ops.append(sqrt(rate) * C[i])

# cavity excitation, if temperature > 0
for i in [0,1]:
    rate = kappa[i] * n_th
    if rate > 0.0:
        c_ops.append(sqrt(rate) * C[i].dag())

# atom relaxation
for i in [0,1]:
    rate = gamma[i]
    if rate > 0.0:
        c_ops.append(sqrt(rate) * A[i])
In [10]:
print("Number of collapse operators = %d" % len(c_ops))
Number of collapse operators = 4

Evolve the system

In [11]:
opts = Odeoptions()
output = mesolve(H, psi0, tlist, c_ops, list(hstack([NA,NC,P,NC[1]-NC[0],NC[1]+NC[0]])), options=opts)
#output = mcsolve_f90(H, psi0, tlist, c_ops, list(hstack([NA,NC,P,NC[1]-NC[0],NC[1]+NC[0]])), ntraj=100)

Visualize the results

In [12]:
M = 2
In [13]:
fig, axes = subplots(1, 2, sharex=True, figsize=(12,4))

for m in range(M):
    axes[0].plot(tlist, real(output.expect[m]), label="Atom #%d" % m)

axes[0].legend(loc=0)
axes[0].set_xlabel('Time')
axes[0].set_title('Atomic Occupation probability')

for m in range(M):
    axes[1].plot(tlist, real(output.expect[M+m]), label="Cavity #%d" % m)

axes[1].legend(loc=0)
axes[1].set_xlabel('Time')
axes[1].set_title('Cavity Occupation probability');

For $g > g_c$, the energy is localized in one cavity-atom system. See Fig. 2 in Schmidt et al.

In [14]:
fig, axes = subplots(M, 2, sharex=True, figsize=(12,2*M))

for m in range(M):
    axes[m,0].plot(tlist, real(output.expect[m]), 'b', label="Atom #%d" % m)
    axes[m,0].legend(loc=0)
    axes[m,0].set_ylim(0,1)

    axes[m,1].plot(tlist, real(output.expect[M+m]), 'g', label="Cavity #%d" % m)
    axes[m,1].legend(loc=0)

    if m == 0:
        axes[m,0].set_xlabel('Time')
        axes[m,1].set_xlabel('Time')
        axes[m,0].set_title('Atomic Occupation probability')
        axes[m,1].set_title('Cavity Occupation probability')
In [15]:
"""
Plot polaritons
"""
fig, axes = subplots(M, 1, sharex=True, sharey=True, figsize=(12,2*M))

for m in range(M):
    axes[m].plot(tlist, real(output.expect[2*M+m]), label="Polariton system #%d" % m)
    axes[m].legend(loc=0)

    if m == 0:
        axes[m].set_xlabel('Time')
        axes[m].set_title('Polariton Occupation probability')
In [16]:
"""
Plot photon imbalance
"""

n = output.expect[3*M] / output.expect[3*M+1]

fig, axes = subplots(1, 1, sharex=True, sharey=True, figsize=(12,2*M))

axes.plot(tlist * J, n, label="Photon imbalance")
axes.set_xlabel('Time')
axes.set_ylim(-1.1, 1.1)
axes.set_title('Photon imbalance');

Large photon imbalance ($\sim - 1$) shows that the energy is localized in one cavity-atom system.