03: Jaynes-Cummings Model and the Rotating Wave Approximation#

In this tutorial we will construct the Jaynes-Cummings Hamiltonian (with and without the RWA) and see how the system evolves under the Schrodinger equation (that is, without dissipation) . We will use this to investigate the limits of the RWA in the JCM.

The Jaynes-Cumming model is the simplest possible model of quantum mechanical light-matter interaction, describing a single two-level atom interacting with a single electromagnetic cavity mode. The full Hamiltonian for the system (in dipole interaction form) is given by

\[H = H_{\rm atom} + H_{\rm cavity} + H_{\rm interact}\]

The atom Hamiltonian we use in this case is

\[\frac{1}{2} \hbar \omega_{a} \sigma_z\]

where \(\omega_a\) is the system frequency.

Note that the Hamiltonian for the atom may take numerous forms. Any Hermitean operator on a two-level state is possible, but it is useful to nomalize the operator so that the difference between its eigenvalues is \(1\) so that \(\omega_a\) has consistent units.

The cavity Hamiltonian is given by

\[H_{\rm cavity} = \hbar \omega_c a^\dagger a\]

where \(\omega_c\) is \(\omega_a\) is the frequencies of the cavity and \(a\) and \(a^\dagger\) are the annihilation and creation operators of the cavity respectively.

The interaction Hamiltonian is given by

\[H_{\rm interact} = \hbar g(a^\dagger + a)(\sigma_- + \sigma_+)\]

or with the rotating-wave approximation

\[H_{\rm interact-RWA} = \hbar g(a^\dagger\sigma_- + a\sigma_+)\]

where \(\sigma_-\) and \(\sigma_+\) are the annihilation and creation operators for the atom respectively.

Note that in this notebook we will work in units where \(\hbar=1\).

Tasks#

Imports#

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

Helper functions#

def display_eigenstates(op):
    """ Display the eigenvalues and eigenstates of an operator. """
    evals, evecs = op.eigenstates()
    print("Eigenvalues:", evals)
    print()
    print("Eigenstates")
    print("===========")
    for v in evecs:
        display(v)
def jcm_h(wc, wa, g, N, atom):
    """ Construct the Jaynes-Cummings Hamiltonian (non-RWA). """
    a = qutip.tensor(qutip.destroy(N), qutip.qeye(2))
    sm = qutip.tensor(qutip.qeye(N), qutip.sigmam())
    atom = qutip.tensor(qutip.qeye(N), atom)
    
    H = wc * a.dag() * a + wa * atom + g * (a.dag() + a) * (sm + sm.dag())
    return H
def jcm_rwa_h(wc, wa, g, N, atom):
    """ Construct the Jaynes-Cummings Hamiltonian (RWA). """
    a = qutip.tensor(qutip.destroy(N), qutip.qeye(2))
    sm = qutip.tensor(qutip.qeye(N), qutip.sigmam())
    atom = qutip.tensor(qutip.qeye(N), atom)

    H = wc * a.dag() * a + wa * atom + g * (a.dag() * sm + a * sm.dag())
    return H

Construct the Hamiltonian#

  • add variables for the atom and cavity parameters

  • create the operators for the JCM Hamiltonian

  • combine into the JCM Hamiltonian (no RWA)

  • look at the energy eigenvalues and eigenstates of the Hamiltonian

Here are some example parameter values to start with:

\( \omega_c = 2 \pi \\ \omega_a = 2 \pi \\ g = 0.05 \cdot 2 \pi \\ \)

# system parameters
wc = 1.0 * 2 * np.pi  # cavity frequency
wa = 1.0 * 2 * np.pi  # atom frequency
g = 0.05 * 2 * np.pi # 0.05 * 2 * np.pi  # coupling strength
N = 15  # number of cavity fock states
# operators for the JCM
# this is the annihilation operator for the cavity (note it acts trivally on the atom) 
a = qutip.tensor(qutip.destroy(N), qutip.qeye(2))
# note the creation operator for the cavity is given by a.dag() - that is a "dagger"

# this is annihilation operator for the atom (note it acts trivally on the cavity)
sm = qutip.tensor(qutip.qeye(N), qutip.sigmam())
# hamiltonian of atom
H_atom = 0.5 * qutip.sigmaz()
# hamiltonian (non-rwa)
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() + a) * (sm + sm.dag())
display_eigenstates(H)
Eigenvalues: [-7.85889344e-03  5.96126089e+00  6.58938206e+00  1.21144878e+01
  1.30025059e+01  1.82980462e+01  1.93852982e+01  2.44973215e+01
  2.57523736e+01  3.07066509e+01  3.21093947e+01  3.69231295e+01
  3.84592664e+01  4.31450310e+01  4.48037151e+01  4.93712306e+01
  5.11438657e+01  5.56009463e+01  5.74804999e+01  6.18336091e+01
  6.38141870e+01  6.80687891e+01  7.01453568e+01  7.43061538e+01
  7.64743463e+01  8.05458680e+01  8.28020375e+01  8.68443187e+01
  8.91878173e+01  9.43576629e+01]

Eigenstates
===========
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\-1.000\\0.025\\0.0\\0.0\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.698\\0.0\\0.0\\-0.715\\0.025\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.716\\0.0\\0.0\\0.698\\-0.025\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\-0.018\\-0.694\\0.0\\0.0\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.017\\0.719\\0.0\\0.0\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}6.687\times10^{-04}\\0.0\\0.0\\0.026\\0.691\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}5.872\times10^{-04}\\0.0\\0.0\\0.024\\0.722\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\-1.502\times10^{-05}\\-0.001\\0.0\\0.0\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\1.230\times10^{-05}\\0.001\\0.0\\0.0\\\vdots\\1.015\times10^{-12}\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}3.916\times10^{-07}\\0.0\\0.0\\3.044\times10^{-05}\\0.002\\\vdots\\0.0\\0.0\\-1.111\times10^{-12}\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}2.960\times10^{-07}\\0.0\\0.0\\2.434\times10^{-05}\\0.001\\\vdots\\0.0\\0.0\\1.675\times10^{-12}\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\7.562\times10^{-09}\\8.888\times10^{-07}\\0.0\\0.0\\\vdots\\-3.258\times10^{-09}\\7.235\times10^{-11}\\0.0\\0.0\\-1.666\times10^{-12}\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\-5.343\times10^{-09}\\-6.541\times10^{-07}\\0.0\\0.0\\\vdots\\-4.807\times10^{-09}\\1.100\times10^{-10}\\0.0\\0.0\\-2.613\times10^{-12}\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}1.603\times10^{-10}\\0.0\\0.0\\1.881\times10^{-08}\\1.561\times10^{-06}\\\vdots\\0.0\\0.0\\-4.379\times10^{-09}\\1.007\times10^{-10}\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}-1.054\times10^{-10}\\0.0\\0.0\\-1.293\times10^{-08}\\-1.121\times10^{-06}\\\vdots\\0.0\\0.0\\-6.665\times10^{-09}\\1.585\times10^{-10}\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\-2.744\times10^{-12}\\-4.312\times10^{-10}\\0.0\\0.0\\\vdots\\6.387\times10^{-06}\\-1.876\times10^{-07}\\0.0\\0.0\\5.715\times10^{-09}\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\1.692\times10^{-12}\\2.754\times10^{-10}\\0.0\\0.0\\\vdots\\9.109\times10^{-06}\\-2.805\times10^{-07}\\0.0\\0.0\\8.954\times10^{-09}\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\7.889\times10^{-12}\\8.757\times10^{-10}\\\vdots\\0.0\\0.0\\-7.605\times10^{-06}\\2.313\times10^{-07}\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\-4.723\times10^{-12}\\-5.443\times10^{-10}\\\vdots\\0.0\\0.0\\-1.108\times10^{-05}\\3.543\times10^{-07}\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\-0.005\\1.980\times10^{-04}\\0.0\\0.0\\-8.907\times10^{-06}\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\0.006\\-2.722\times10^{-04}\\0.0\\0.0\\1.325\times10^{-05}\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\0.0\\0.0\\0.005\\-2.214\times10^{-04}\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\0.0\\0.0\\-0.006\\3.091\times10^{-04}\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\-0.734\\0.061\\0.0\\0.0\\-0.005\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\-0.673\\0.067\\0.0\\0.0\\-0.007\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\0.0\\0.0\\-0.735\\0.063\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\0.0\\0.0\\-0.671\\0.069\\0.0\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\-0.069\\-0.688\\0.0\\0.0\\0.722\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\-0.059\\-0.720\\0.0\\0.0\\-0.692\\\end{matrix}\right)$\end{split}\]
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}0.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\0.0\\0.0\\-0.093\\-0.996\\0.0\\\end{matrix}\right)$\end{split}\]

Solve the Schrodinger equation#

  • create the initial state of the system (use the state with no photons and the spin system in its excited state)

  • evolve the system for some time, saving the result.

If you need to remind yourself of how sesolve, remember that you can type qutip.sesolve? into a notebook cell to bring up the documentation.

# initial state
psi0 = qutip.basis([N, 2], [0, 0])  # start with an excited atom
psi0
\[\begin{split}Quantum object: dims = [[15, 2], [1, 1]], shape = (30, 1), type = ket $ \\ \left(\begin{matrix}1.0\\0.0\\0.0\\0.0\\0.0\\\vdots\\0.0\\0.0\\0.0\\0.0\\0.0\\\end{matrix}\right)$\end{split}\]
# Solve using sesolve
tlist = np.linspace(0, 25, 101)
result = qutip.sesolve(H, psi0, tlist)

Visualise the evolution#

  • create expectation operators for observing the state of the system.

  • add these to sesolve

  • plot the expectation values together on a set of axes

  • change the value of g and see how it affects the period

Two good expectation operators to use are the projectors on the light and matter sub-systems.

# Operators to determine the expectation values of:
eop_a = a.dag() * a  # light
eop_sm = sm.dag() * sm  # matter
# Solve using sesolve
tlist = np.linspace(0, 25, 101)
result = qutip.sesolve(H, psi0, tlist, e_ops=[eop_a, eop_sm])
plt.plot(tlist, result.expect[0], label="Light")
plt.plot(tlist, result.expect[1], label="matter")
plt.legend();
../_images/03-qutip-jcm-i_23_0.png
# vary g
H = jcm_h(wc, wa, 0.1 * 2 * np.pi, N, H_atom)

tlist = np.linspace(0, 25, 101)

result = qutip.sesolve(H, psi0, tlist, e_ops=[eop_a, eop_sm])

plt.plot(tlist, result.expect[0], label="Light")
plt.plot(tlist, result.expect[1], label="matter")
plt.legend();
../_images/03-qutip-jcm-i_24_0.png

Compare the RWA and non-RWA#

  • construct another Hamiltonian that uses the RWA

  • evolve the system under the RWA Hamiltonian.

  • add the results to the plot

  • experiment with parameters to determine where the RWA non-RWA diverges

# Construct the RWA Hamiltonian
H = jcm_h(wc, wa, g, N, H_atom)
H_RWA = jcm_rwa_h(wc, wa, g, N, H_atom)
# Hmm. The results look the same.
tlist = np.linspace(0, 25, 101)

result = qutip.sesolve(H, psi0, tlist, e_ops=[eop_a, eop_sm])
plt.plot(tlist, result.expect[0], label="Light")
plt.plot(tlist, result.expect[1], label="matter")

result_rwa = qutip.sesolve(H_RWA, psi0, tlist, e_ops=[eop_a, eop_sm])
plt.plot(tlist, result_rwa.expect[0], label="Light (RWA)")
plt.plot(tlist, result_rwa.expect[1], label="matter (RWA)")

plt.legend();
../_images/03-qutip-jcm-i_27_0.png
# Try with different frequencies w_c and w_a:
tlist = np.linspace(0, 25, 101)
f = 0.9

H = jcm_h(wc, f * wa, g, N, H_atom)
H_RWA = jcm_rwa_h(wc, f * wa, g, N, H_atom)

result = qutip.sesolve(H, psi0, tlist, e_ops=[eop_a, eop_sm])
plt.plot(tlist, result.expect[0], label="Light")
plt.plot(tlist, result.expect[1], label="matter")

result_rwa = qutip.sesolve(H_RWA, psi0, tlist, e_ops=[eop_a, eop_sm])
plt.plot(tlist, result_rwa.expect[0], label="Light (RWA)")
plt.plot(tlist, result_rwa.expect[1], label="matter (RWA)")

plt.legend();
../_images/03-qutip-jcm-i_28_0.png