Note

This page was generated from tut/4-Analysis/cQED-with-the-Jaynes-Cummings-Interaction-Model.ipynb.

Circuit Quantum Electrodynamics (cQED) and the Jaynes-Cummings Interaction Model Using the HCPB Class and QuTip

💡 Using this tutorial without the Qt GUI

This tutorial uses the desktop MetalGUI. To follow along on Colab, Binder, JupyterHub, or any environment where Qt isn’t available, replace any ``gui.rebuild()`` / ``gui.screenshot()`` call with ``qm.view(design)`` — it renders the design to a matplotlib Figure you can display inline or save with fig.savefig(...).

See 1.4 Headless quick view for a complete runnable walkthrough and `docs/headless-usage.rst <../../../docs/headless-usage.rst>`__ for the full reference.

Overview

This tutorial describes how to use the Cooper Pair Box Hamiltonian (HCPB) class in Qiskit Metal to model the interaction between a transmon qubit and a CPW resonator. We will build the so-called Jaynes-Cummings Hamiltonian and then perform some analysis.

The basic form of the interaction Hamiltonian has the form:

\(H = H_{transmon} + H_{CPW} + H_{interaction}\)

While we could model both the transmon and the CPW as harmonic oscillators, we’ll use the HCPB class in Qiskit Metal to model the transmon and we’ll model the CPW as just a regular harmonic oscillator using ladder operators defined in QuTiP. The resulting Jaynes-Cummings Hamiltonian will then have the form:

$H = H_{transmon} + \omega_{r} a^{\dagger} a + gN(a+a^{\dagger}) $

where \(\omega_{r}\) is the resonator frequency, \(a\) and \(a^{\dagger}\) are the annihilation and creation operators, respectively, for photons in the resonator cavity, g is the transmon-cavity coupling and \(N\) is the charge number operator of the transmon.

In this tutorial, we’ll construct the above Hamiltonian using the HCPB class in Qiskit Metal and QuTiP. Then we study how the system evolves in time using the QuTiP master equation solver. We’ll add some photon dissipation to the CPW resonator and show that impacts the dynamics.

Helpful Background Material

This tutorial draws heavily from existing tutorials in QuTiP as well as scqbuits, and you may find it helpful to review some of those tutorials to build up some additional background knowledge:

https://scqubits.readthedocs.io/en/latest/guide/ipynb/hilbertspace.html#

https://nbviewer.ipython.org/github/jrjohansson/qutip-lectures/blob/master/Lecture-0-Introduction-to-QuTiP.ipynb

https://nbviewer.ipython.org/github/jrjohansson/qutip-lectures/blob/master/Lecture-1-Jaynes-Cumming-model.ipynb

https://nbviewer.ipython.org/github/jrjohansson/qutip-lectures/blob/master/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.ipynb

Table of Contents

Here’s what we’ll work through in this tutorial notebook:

  • Transmon Qubit Hamiltonian (Using HCPB Class: anharmonic)

  • CPW Resonator Hamiltonian (Using QuTiP: harmonic)

  • Composite Hamiltonian w/ Tensor Products

  • Time Evolution: Defining the Initial State(s)

  • Time Evolution: Defining Occupation Operators

  • Time Evolution: Steady-State Evolution

  • Time Evolution: Time Evolution w/ CPW Photon Loss (Vacuum-Rabi Oscillations)

  • How to Set Up Another Example

Let’s start by loading the relevant modules:

[ ]:

%load_ext autoreload
%autoreload 2
%config IPCompleter.greedy = True
%matplotlib inline
%config Completer.use_jedi = False
%config InlineBackend.figure_format = 'svg'

from qiskit_metal.analyses.hamiltonian.transmon_charge_basis import Hcpb
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from qutip import *
from math import *

Defining the Transmon Qubit

We’ll use the HCPB (Cooper Pair Box Hamiltonian) class within qiskit metal to define a transmon by supplying the values of charging energy and Josephson energy:

[ ]:

# This contructs the HCPB Hamiltonian for the qubit
H1 = Hcpb(nlevels=20, Ej=20.0, Ec=0.4, ng=0.00)  # Hamiltonian definition

# convert to qutip matrices
Ham = H1.h0_to_qutip(4)  # Hamiltonian matrix for transmon
Num = H1.n_to_qutip(4)  # Number operator matrix for transmon

# This calculates the 0->1 transition frequency of the qubit
f01 = H1.fij(0, 1)
f01

This transmon has a transition frequency of about 7.6 GHz and we’ve used a wrapper around QuTiP to generate the diagonalized Hamiltonian matrix and the Number operator matrix for this transmon. Here’s what the Hamiltonian matrix looks like:

[ ]:

Ham

Defining the CPW Resonator

Rather than use the HCPB class for the CPW resonator, we’ll just use ladder operators in QuTiP. Note this means that the CPW resonator will be perfectly harmonic (unlike the transmon.)

[ ]:

wc = 7.5  # frequency of the CPW resonator cavity
N = 5  # number of Fock states in the resonator cavity
a = destroy(5)  # lowering operator
Ham_CPW = wc * a.dag() * a  # resonator Hamiltonian
Ham_CPW

Note that the number of Fock states in the transmon does not need to be the same as the number of Fock states in the CPW resonator. As we’ve defined them above, we have N=4 in the transmon and N=5 in the CPW.

Tensor Products for the Composite Hamiltonian

Before we can construct the composite Hamiltonian, we need to take the appropriate tensor products between the transmon basis and the CPW basis. We do this using the N-dimensional identity matrix given by qeye(N) in QuTiP.

This turns both the 4x4 transmon Hamiltonian matrix as well as the 5x5 CPW Hamiltonian matrix into 20x20 matrices.

[ ]:

# Tensor Product for Transmon Hamiltonian (N=number of Fock states in the CPW)
Ham_t = tensor(qeye(N), Ham)

# Tensor Product for Transmon Charge Number Operator  (N=number of Fock states in the CPW)
Num_t = tensor(qeye(N), Num)

# Tensor Product for CPW resonator creation/annihilation operators  (4=number of Fock states in the transmon)
a = tensor(destroy(N), qeye(4))

Now we have all the components we need to construct the Jaynes-Cummings Hamiltonian:

$H = H_{transmon} + \omega_{r} a^{\dagger} a + gN(a+a^{\dagger}) $

[ ]:

g = 0.1  # transmon-cavity coupling

Ham_JC = Ham_t + wc * a.dag() * a + g * Num_t * (a + a.dag())
Ham_JC

Time Evolution: Defining the Initial State

Just as we did for the composite Hamiltonian matrix, we need to take the tensor product of the initial transmon/CPW states. Let’s suppose that the CPW starts in the ground state (n=0) while the transmon starts in the first excited state (n=1):

[ ]:

# This defines the initial state of the cavity; recall that we set N=5
# (N,n) = (Number of Fock states in the cavity, initial state of the cavity)
Psi0_cavity = basis(N, 0)  # start with ground state of cavity

# This defines the initial state of the transmon qubit
# Ham.eigenstates()[1][x] returns the xth state
Psi_qubit = Ham.eigenstates()[1][1]  # first excited state

# Tensor Product for Initial State of Composite System
psi0 = tensor(Psi0_cavity, Psi_qubit)

Time Evolution: Defining Occupation Operators

Once we let the system evolve in time, it will be interesting to see how the occupation of the transmon and/or CPW changes. In order to do this, we’ll have to define the occupation (number) operators and construct their corresponding tensor products.

First let’s define the operators for both the transmon and the CPW:

[ ]:

# Create occupation operators for the transmon qubit (using HCPB class)
s0 = Ham.eigenstates()[1][0]  # n=0 state
n0 = s0 * s0.dag()  # occupation operator for n=0
s1 = Ham.eigenstates()[1][1]  # n=1 state
n1 = s1 * s1.dag()  # occupation operator for n=1
s2 = Ham.eigenstates()[1][2]  # n=2 state
n2 = s2 * s2.dag()  # occupation operator for n=2
s3 = Ham.eigenstates()[1][3]  # n=3 state
n3 = s3 * s3.dag()  # occupation operator for n=3

# Create occupation operators for the CPW resonator (using QuTiP)
c0 = basis(N, 0)  # n=0 state
n0_cavity = c0 * c0.dag()  # occupation operator for n=0
c1 = basis(N, 1)  # n=1 state
n1_cavity = c1 * c1.dag()  # occupation operator for n=1
c2 = basis(N, 2)  # n=2 state
n2_cavity = c2 * c2.dag()  # occupation operator for n=2
c3 = basis(N, 3)  # n=3 state
n3_cavity = c3 * c3.dag()  # occupation operator for n=3
c4 = basis(N, 4)  # n=4 state
n4_cavity = c4 * c4.dag()  # occupation operator for n=4

Now let’s construct the corresponding tensor products:

[ ]:

# occupation probability matrices for the qubit
n0_qubit = tensor(qeye(N), n0)
n1_qubit = tensor(qeye(N), n1)
n2_qubit = tensor(qeye(N), n2)
n3_qubit = tensor(qeye(N), n3)

# occupation probability matrices for the CPW resonator
n0_resonator = tensor(n0_cavity, qeye(4))
n1_resonator = tensor(n1_cavity, qeye(4))
n2_resonator = tensor(n2_cavity, qeye(4))
n3_resonator = tensor(n3_cavity, qeye(4))
n4_resonator = tensor(n4_cavity, qeye(4))

Time Evolution: Steady State (No Dissipation)

[ ]:

# time will vary from t=0 to t=50s with 1001 increments
tlist = np.linspace(0, 50, 1001)

# initial list of collapse operators (empty here)
c_ops = []

# solve master equation and for output calculate the occupation numbers
output = mesolve(
    Ham_JC,
    psi0,
    tlist,
    c_ops,
    [n0_resonator, n1_resonator, n2_resonator, n0_qubit, n1_qubit, n2_qubit],
)

# tabulate the occupation probabilities for resonator and qubit
occ0_resonator = output.expect[0]
occ1_resonator = output.expect[1]
occ2_resonator = output.expect[2]
occ0_qubit = output.expect[3]
occ1_qubit = output.expect[4]
occ2_qubit = output.expect[5]

fig, axes = plt.subplots(1, 1, figsize=(10, 5))

plt.figure(1)
plt.subplot(121)
plt.title("Resonator")
plt.plot(tlist, occ0_resonator, "b", label="n=0")
plt.plot(tlist, occ1_resonator, "r", label="n=1")
plt.plot(tlist, occ2_resonator, "g", label="n=2")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

plt.subplot(122)
plt.legend()
plt.title("Qubit")
plt.plot(tlist, occ0_qubit, "b", label="n=0")
plt.plot(tlist, occ1_qubit, "r", label="n=1")
plt.plot(tlist, occ2_qubit, "g", label="n=2")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

We see that the CPW resonator is in the n=0 state and the transmon qubit is in the n=1 state, just as we’ve defined. There is no dissipation in the system so average occupations do not change over time.

Time Evolution: Observing Vacuum-Rabi Oscillations by Adding Dissipation

Here we’ll add some photon dissipation in the CPW resonator and then evolve the system in time:

[ ]:

g = 0.1  # transmon-cavity coupling; set to be strong

# rebuild Hamiltonian matrix of composite system
Ham_JC = Ham_t + wc * a.dag() * a + g * Num_t * (a + a.dag())

# time will vary from t=0 to t=100s
tlist = np.linspace(0, 100, 1001)

# photon relaxation rates of the transmon and the CPW cavity
kappa = 0.05  # 0.005 # cavity photon relaxation rate
n_th_a = 0.0  # avg number of thermal bath excitations

# initial list of collapse operators
c_ops = []

# cavity relaxation
rate = kappa * (1 + n_th_a)
if rate > 0.0:
    c_ops.append(sqrt(rate) * a)

# solve master equation and for output calculate the occupation numbers
output = mesolve(
    Ham_JC,
    psi0,
    tlist,
    c_ops,
    [n0_resonator, n1_resonator, n2_resonator, n0_qubit, n1_qubit, n2_qubit],
)

# tabulate the occupation probabilities for resonator and qubit
occ0_resonator = output.expect[0]
occ1_resonator = output.expect[1]
occ2_resonator = output.expect[2]
occ0_qubit = output.expect[3]
occ1_qubit = output.expect[4]
occ2_qubit = output.expect[5]

fig, axes = plt.subplots(1, 1, figsize=(10, 5))

plt.figure(1)
plt.subplot(121)
plt.title("Resonator")
plt.plot(tlist, occ0_resonator, "b", label="n=0")
plt.plot(tlist, occ1_resonator, "r", label="n=1")
plt.plot(tlist, occ2_resonator, "g", label="n=2")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

plt.subplot(122)
plt.legend()
plt.title("Qubit")
plt.plot(tlist, occ0_qubit, "b", label="n=0")
plt.plot(tlist, occ1_qubit, "r", label="n=1")
plt.plot(tlist, occ2_qubit, "g", label="n=2")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

Here we see that the resonator starts in the ground (n=0) state, is briefly excited into the n=1 state and then gradually settles back into the ground state. Likewise, the qubit begins in the first excited state (n=1) and then gradually settles into the ground state.

Let’s look at the same data in a slightly different way. By plotting the n=0 probabilities and the n=1 probabilities together, we can see that the CPW and transmon excitations are always out of phase:

[ ]:

fig, axes = plt.subplots(1, 1, figsize=(10, 5))

plt.figure(1)
plt.subplot(121)
plt.title("n=0")
plt.plot(tlist, occ0_resonator, "b", label="resonator")
plt.plot(tlist, occ0_qubit, "r", label="qubit")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

plt.subplot(122)
plt.legend()
plt.title("n=1")
plt.plot(tlist, occ1_resonator, "b", label="resonator")
plt.plot(tlist, occ1_qubit, "r", label="qubit")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

Another Example

In the previous example, we considered the initial state where the CPW resonator was in the n=0 state and the qubit was in the first excited state (n=1).

Let’s see what happens if we consider the CPW resonator in the second excited state (n=2) and the transmon qubit in the ground state.

[ ]:

# This defines the initial state of the cavity; recall that we set N=5
# (N,n) = (Number of Fock states in the cavity, initial state of the cavity)
Psi0_cavity = basis(N, 2)  # start with ground state of cavity

# This defines the initial state of the transmon qubit
# Ham.eigenstates()[1][x] returns the xth state
Psi_qubit = Ham.eigenstates()[1][0]  # first excited state

# Tensor Product for Initial State of Composite System
psi0 = tensor(Psi0_cavity, Psi_qubit)

The definitions of the Hamiltonian and operators remain unchanged, so we can directly evolve the system in time they same way we did above. Here’s the steady-state time evoluation:

[ ]:

# time will vary from t=0 to t=50s with 1001 increments
tlist = np.linspace(0, 50, 1001)

# initial list of collapse operators (empty here)
c_ops = []

# solve master equation and for output calculate the occupation numbers
output = mesolve(
    Ham_JC,
    psi0,
    tlist,
    c_ops,
    [n0_resonator, n1_resonator, n2_resonator, n0_qubit, n1_qubit, n2_qubit],
)

# tabulate the occupation probabilities for resonator and qubit
occ0_resonator = output.expect[0]
occ1_resonator = output.expect[1]
occ2_resonator = output.expect[2]
occ0_qubit = output.expect[3]
occ1_qubit = output.expect[4]
occ2_qubit = output.expect[5]

fig, axes = plt.subplots(1, 1, figsize=(10, 5))

plt.figure(1)
plt.subplot(121)
plt.title("Resonator")
plt.plot(tlist, occ0_resonator, "b", label="n=0")
plt.plot(tlist, occ1_resonator, "r", label="n=1")
plt.plot(tlist, occ2_resonator, "g", label="n=2")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

plt.subplot(122)
plt.legend()
plt.title("Qubit")
plt.plot(tlist, occ0_qubit, "b", label="n=0")
plt.plot(tlist, occ1_qubit, "r", label="n=1")
plt.plot(tlist, occ2_qubit, "g", label="n=2")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

And here’s the time evolution with photon loss included in the CPW resonator:

[ ]:

# time will vary from t=0 to t=100s
tlist = np.linspace(0, 100, 1001)

# photon relaxation rates of the transmon and the CPW cavity
kappa = 0.05  # 0.005 # cavity photon relaxation rate
n_th_a = 0.0  # avg number of thermal bath excitations

# initial list of collapse operators
c_ops = []

# cavity relaxation
rate = kappa * (1 + n_th_a)
if rate > 0.0:
    c_ops.append(sqrt(rate) * a)

# solve master equation and for output calculate the occupation numbers
output = mesolve(
    Ham_JC,
    psi0,
    tlist,
    c_ops,
    [n0_resonator, n1_resonator, n2_resonator, n0_qubit, n1_qubit, n2_qubit],
)

# tabulate the occupation probabilities for resonator and qubit
occ0_resonator = output.expect[0]
occ1_resonator = output.expect[1]
occ2_resonator = output.expect[2]
occ0_qubit = output.expect[3]
occ1_qubit = output.expect[4]
occ2_qubit = output.expect[5]

fig, axes = plt.subplots(1, 1, figsize=(10, 5))

plt.figure(1)
plt.subplot(121)
plt.title("Resonator")
plt.plot(tlist, occ0_resonator, "b", label="n=0")
plt.plot(tlist, occ1_resonator, "r", label="n=1")
plt.plot(tlist, occ2_resonator, "g", label="n=2")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

plt.subplot(122)
plt.legend()
plt.title("Qubit")
plt.plot(tlist, occ0_qubit, "b", label="n=0")
plt.plot(tlist, occ1_qubit, "r", label="n=1")
plt.plot(tlist, occ2_qubit, "g", label="n=2")
plt.xlabel("Time")
plt.ylabel("Occupation Probability")
plt.legend(title="Energy Levels", loc="upper right")

Given a long enough time, we see that both the CPW and the qubit eventually settle into the n=0 states.

Toward Multi-Qubit Gates…

With the above methodology, you can perform analysis for more complicated but realistic systems like two-qubit CR or iSWAP gates!!!



For more information, review the Introduction to Quantum Computing and Quantum Hardware lectures below

  • Superconducting Qubits I: Quantizing a Harmonic Oscillator, Josephson Junctions Part 1
 Lecture Video Lecture Notes Lab
  • Superconducting Qubits I: Quantizing a Harmonic Oscillator, Josephson Junctions Part 2
 Lecture Video Lecture Notes Lab
  • Superconducting Qubits I: Quantizing a Harmonic Oscillator, Josephson Junctions Part 3
 Lecture Video Lecture Notes Lab
  • Superconducting Qubits II: Circuit Quantum Electrodynamics, Readout and Calibration Methods Part 1
 Lecture Video Lecture Notes Lab
  • Superconducting Qubits II: Circuit Quantum Electrodynamics, Readout and Calibration Methods Part 2
 Lecture Video Lecture Notes Lab
  • Superconducting Qubits II: Circuit Quantum Electrodynamics, Readout and Calibration Methods Part 3
 Lecture Video Lecture Notes Lab