Author: J.R. Johansson, robert@riken.jp
Latest version of this ipython notebook lecture is available at: http://github.com/jrjohansson/qutip-lectures
# setup the matplotlib graphics library and configure it to show
# figures inline in the notebook
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
# make qutip available in the rest of the notebook
from qutip import *
Consider a single atom coupled to a single cavity mode. If there atom excitation rate exceeds the relaxation rate, a population inversion can occur in the atom, and if coupled to the cavity the atom can then act as a photon pump on the cavity.
The coherent dynamics in this model is described by the Hamiltonian
where $\omega_0$ is the cavity energy splitting, $\omega_a$ is the atom energy splitting and $g$ is the atom-cavity interaction strength.
In addition to the coherent dynamics the following incoherent processes are also present: 1) $\kappa$ relaxation and thermal excitations of the cavity, $\Gamma$ atomic excitation rate (pumping process).
The Lindblad master equation for the model is:
$+ \kappa (1 + n_{\rm th}) \left(a\rho a^\dagger - \frac{1}{2}a^\dagger a\rho - \frac{1}{2}\rho a^\dagger a\right)$
$+ \kappa n_{\rm th} \left(a^\dagger\rho a - \frac{1}{2}a a^\dagger \rho - \frac{1}{2}\rho a a^\dagger\right)$
in units where $\hbar = 1$.
References:
w0 = 1.0 * 2 * pi # cavity frequency
wa = 1.0 * 2 * pi # atom frequency
g = 0.05 * 2 * pi # coupling strength
kappa = 0.04 # cavity dissipation rate
gamma = 0.00 # atom dissipation rate
Gamma = 0.35 # atom pump rate
N = 50 # number of cavity fock states
n_th_a = 0.0 # avg number of thermal bath excitation
tlist = linspace(0, 150, 101)
# intial state
psi0 = tensor(basis(N,0), basis(2,0)) # start without excitations
# operators
a = tensor(destroy(N), qeye(2))
sm = tensor(qeye(N), destroy(2))
sx = tensor(qeye(N), sigmax())
# Hamiltonian
H = w0 * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() + a) * sx
H
# collapse operators
c_ops = []
rate = kappa * (1 + n_th_a)
if rate > 0.0:
c_ops.append(sqrt(rate) * a)
rate = kappa * n_th_a
if rate > 0.0:
c_ops.append(sqrt(rate) * a.dag())
rate = gamma
if rate > 0.0:
c_ops.append(sqrt(rate) * sm)
rate = Gamma
if rate > 0.0:
c_ops.append(sqrt(rate) * sm.dag())
Here we evolve the system with the Lindblad master equation solver, and we request that the expectation values of the operators $a^\dagger a$ and $\sigma_+\sigma_-$ are returned by the solver by passing the list [a.dag()*a, sm.dag()*sm]
as the fifth argument to the solver.
opt = Odeoptions(nsteps=2000) # allow extra time-steps
output = mesolve(H, psi0, tlist, c_ops, [a.dag() * a, sm.dag() * sm], options=opt)
Here we plot the excitation probabilities of the cavity and the atom (these expectation values were calculated by the mesolve
above).
n_c = output.expect[0]
n_a = output.expect[1]
fig, axes = subplots(1, 1, sharex=True, figsize=(12,8))
axes.plot(tlist, n_c, label="Cavity")
axes.plot(tlist, n_a, label="Atom excited state")
axes.set_xlim(0, 150)
axes.legend(loc=0)
axes.set_xlabel('Time')
axes.set_ylabel('Occupation probability')
<matplotlib.text.Text at 0x2e08e50>
rho_ss = steadystate(H, c_ops)
fig, axes = subplots(1, 2, figsize=(16,6))
xvec = linspace(-5,5,200)
rho_cavity = ptrace(rho_ss, 0)
W = wigner(rho_cavity, xvec, xvec)
wlim = abs(W).max()
axes[1].contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-wlim,wlim), cmap=get_cmap('RdBu'))
axes[1].set_xlabel(r'Im $\alpha$', fontsize=18)
axes[1].set_ylabel(r'Re $\alpha$', fontsize=18)
axes[0].bar(arange(0, N), real(rho_cavity.diag()), color="blue", alpha=0.6)
axes[0].set_ylim(0, 1)
axes[0].set_xlim(0, N)
axes[0].set_xlabel('Fock number', fontsize=18)
axes[0].set_ylabel('Occupation probability', fontsize=18)
<matplotlib.text.Text at 0x3273b50>
opt = Odeoptions(nsteps=5000) # allow extra time-steps
output = mesolve(H, psi0, linspace(0, 25, 5), c_ops, [], options=opt)
rho_ss_sublist = output.states
xvec = linspace(-5,5,200)
fig, axes = subplots(2, len(rho_ss_sublist), figsize=(3*len(rho_ss_sublist), 6))
for idx, rho_ss in enumerate(rho_ss_sublist):
# trace out the cavity density matrix
rho_ss_cavity = ptrace(rho_ss, 0)
# calculate its wigner function
W = wigner(rho_ss_cavity, xvec, xvec)
# plot its wigner function
wlim = abs(W).max()
axes[0,idx].contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-wlim,wlim), cmap=get_cmap('RdBu'))
# plot its fock-state distribution
axes[1,idx].bar(arange(0, N), real(rho_ss_cavity.diag()), color="blue", alpha=0.8)
axes[1,idx].set_ylim(0, 1)
axes[1,idx].set_xlim(0, 15)
References:
def calulcate_avg_photons(N, Gamma):
# collapse operators
c_ops = []
rate = kappa * (1 + n_th_a)
if rate > 0.0:
c_ops.append(sqrt(rate) * a)
rate = kappa * n_th_a
if rate > 0.0:
c_ops.append(sqrt(rate) * a.dag())
rate = gamma
if rate > 0.0:
c_ops.append(sqrt(rate) * sm)
rate = Gamma
if rate > 0.0:
c_ops.append(sqrt(rate) * sm.dag())
# Ground state and steady state for the Hamiltonian: H = H0 + g * H1
rho_ss = steadystate(H, c_ops)
n_cavity = expect(a.dag() * a, rho_ss)
return n_cavity
Gamma_max = 2 * (4*g**2) / kappa
Gamma_vec = linspace(0.0, Gamma_max, 50)
n_avg_vec = [calulcate_avg_photons(N, Gamma) for Gamma in Gamma_vec]
fig, axes = subplots(1, 1, figsize=(12,6))
axes.plot(Gamma_vec * kappa / (4*g**2), n_avg_vec, color="blue", alpha=0.6, label="numerical")
axes.set_xlabel(r'$\Gamma\kappa/(4g^2)$', fontsize=18)
axes.set_ylabel(r'Occupation probability $\langle n \rangle$', fontsize=18)
axes.legend()
<matplotlib.legend.Legend at 0x616dc10>
Here we see that lasing is suppressed for $\Gamma\kappa/(4g^2) > 1$.
Let's look at the fock-state distribution at $\Gamma\kappa/(4g^2) = 0.5$ (lasing regime) and $\Gamma\kappa/(4g^2) = 1.5$ (suppressed regime):
Gamma = 0.5 * (4*g**2) / kappa
c_ops = [sqrt(kappa * (1 + n_th_a)) * a, sqrt(kappa * n_th_a) * a.dag(), sqrt(gamma) * sm, sqrt(Gamma) * sm.dag()]
rho_ss = steadystate(H, c_ops)
fig, axes = subplots(1, 2, figsize=(16,6))
xvec = linspace(-10,10,200)
rho_cavity = ptrace(rho_ss, 0)
W = wigner(rho_cavity, xvec, xvec)
wlim = abs(W).max()
axes[1].contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-wlim,wlim), cmap=get_cmap('RdBu'))
axes[0].bar(arange(0, N), real(rho_cavity.diag()), color="blue", alpha=0.6)
axes[0].set_ylim(0, 1)
axes[0].set_xlim(0, N)
(0, 50)
Gamma = 1.5 * (4*g**2) / kappa
c_ops = [sqrt(kappa * (1 + n_th_a)) * a, sqrt(kappa * n_th_a) * a.dag(), sqrt(gamma) * sm, sqrt(Gamma) * sm.dag()]
rho_ss = steadystate(H, c_ops)
fig, axes = subplots(1, 2, figsize=(16,6))
xvec = linspace(-10,10,200)
rho_cavity = ptrace(rho_ss, 0)
W = wigner(rho_cavity, xvec, xvec)
wlim = abs(W).max()
axes[1].contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-wlim,wlim), cmap=get_cmap('RdBu'))
axes[0].bar(arange(0, N), real(rho_cavity.diag()), color="blue", alpha=0.6)
axes[0].set_ylim(0, 1)
axes[0].set_xlim(0, N)
(0, 50)
Too large pumping rate $\Gamma$ kills the lasing process: reversed threshold.