{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2513e46f",
   "metadata": {},
   "source": [
    "# Circuit Quantum Electrodynamics (cQED) and the Jaynes-Cummings Interaction Model Using the HCPB Class and QuTip"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55d414af",
   "metadata": {},
   "source": [
    "> 💡 **Using this tutorial without the Qt GUI**\n",
    "> \n",
    "> 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(...)`.\n",
    "> \n",
    "> See [1.4 Headless quick view](../1-Overview/1.4-Headless-quick-view-%28no-Qt-GUI%29.ipynb) for a complete runnable walkthrough and [`docs/headless-usage.rst`](../../../docs/headless-usage.rst) for the full reference."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28a5580c",
   "metadata": {},
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c98b86b2",
   "metadata": {},
   "source": [
    "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. \n",
    "\n",
    "The basic form of the interaction Hamiltonian has the form:\n",
    "\n",
    "$H = H_{transmon} + H_{CPW} + H_{interaction}$ \n",
    "\n",
    "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: \n",
    "\n",
    "$H = H_{transmon} + \\omega_{r} a^{\\dagger} a + gN(a+a^{\\dagger})  $\n",
    "\n",
    "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.\n",
    "\n",
    "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. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40317e0a",
   "metadata": {},
   "source": [
    "## Helpful Background Material"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05727eae",
   "metadata": {},
   "source": [
    "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: \n",
    "\n",
    "https://scqubits.readthedocs.io/en/latest/guide/ipynb/hilbertspace.html#\n",
    "\n",
    "https://nbviewer.ipython.org/github/jrjohansson/qutip-lectures/blob/master/Lecture-0-Introduction-to-QuTiP.ipynb\n",
    "\n",
    "https://nbviewer.ipython.org/github/jrjohansson/qutip-lectures/blob/master/Lecture-1-Jaynes-Cumming-model.ipynb\n",
    "\n",
    "https://nbviewer.ipython.org/github/jrjohansson/qutip-lectures/blob/master/Lecture-3B-Jaynes-Cumming-model-with-ultrastrong-coupling.ipynb\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0dfbf53",
   "metadata": {},
   "source": [
    "## Table of Contents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12d3371d",
   "metadata": {},
   "source": [
    "Here's what we'll work through in this tutorial notebook:\n",
    "\n",
    "*  Transmon Qubit Hamiltonian (Using HCPB Class: anharmonic)\n",
    "*  CPW Resonator Hamiltonian (Using QuTiP: harmonic)\n",
    "*  Composite Hamiltonian w/ Tensor Products\n",
    "*  Time Evolution: Defining the Initial State(s) \n",
    "*  Time Evolution: Defining Occupation Operators\n",
    "*  Time Evolution: Steady-State Evolution \n",
    "*  Time Evolution: Time Evolution w/ CPW Photon Loss (Vacuum-Rabi Oscillations) \n",
    "*  How to Set Up Another Example "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15fe9e34",
   "metadata": {},
   "source": [
    "Let's start by loading the relevant modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d6a3b13",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%config IPCompleter.greedy = True\n",
    "%matplotlib inline\n",
    "%config Completer.use_jedi = False\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "from qiskit_metal.analyses.hamiltonian.transmon_charge_basis import Hcpb\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from qutip import *\n",
    "from math import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "079831a1",
   "metadata": {},
   "source": [
    "## Defining the Transmon Qubit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af405be7",
   "metadata": {},
   "source": [
    "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: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1b3bfba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This contructs the HCPB Hamiltonian for the qubit\n",
    "H1 = Hcpb(nlevels=20, Ej=20.0, Ec=0.4, ng=0.00)  # Hamiltonian definition\n",
    "\n",
    "# convert to qutip matrices\n",
    "Ham = H1.h0_to_qutip(4)  # Hamiltonian matrix for transmon\n",
    "Num = H1.n_to_qutip(4)  # Number operator matrix for transmon\n",
    "\n",
    "# This calculates the 0->1 transition frequency of the qubit\n",
    "f01 = H1.fij(0, 1)\n",
    "f01"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad6a5b88",
   "metadata": {},
   "source": [
    "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: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e301267b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Ham"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dbec7bd",
   "metadata": {},
   "source": [
    "## Defining the CPW Resonator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0bb308e",
   "metadata": {},
   "source": [
    "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.) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bd9ef64",
   "metadata": {},
   "outputs": [],
   "source": [
    "wc = 7.5  # frequency of the CPW resonator cavity\n",
    "N = 5  # number of Fock states in the resonator cavity\n",
    "a = destroy(5)  # lowering operator\n",
    "Ham_CPW = wc * a.dag() * a  # resonator Hamiltonian\n",
    "Ham_CPW"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0cc2657",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d256274c",
   "metadata": {},
   "source": [
    "## Tensor Products for the Composite Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "118c3283",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "This turns both the 4x4 transmon Hamiltonian matrix as well as the 5x5 CPW Hamiltonian matrix into 20x20 matrices. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e0efd11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tensor Product for Transmon Hamiltonian (N=number of Fock states in the CPW)\n",
    "Ham_t = tensor(qeye(N), Ham)\n",
    "\n",
    "# Tensor Product for Transmon Charge Number Operator  (N=number of Fock states in the CPW)\n",
    "Num_t = tensor(qeye(N), Num)\n",
    "\n",
    "# Tensor Product for CPW resonator creation/annihilation operators  (4=number of Fock states in the transmon)\n",
    "a = tensor(destroy(N), qeye(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e3e90da",
   "metadata": {},
   "source": [
    "Now we have all the components we need to construct the Jaynes-Cummings Hamiltonian:\n",
    "\n",
    "$H = H_{transmon} + \\omega_{r} a^{\\dagger} a + gN(a+a^{\\dagger})  $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed7d3e57",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = 0.1  # transmon-cavity coupling\n",
    "\n",
    "Ham_JC = Ham_t + wc * a.dag() * a + g * Num_t * (a + a.dag())\n",
    "Ham_JC"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "570933be",
   "metadata": {},
   "source": [
    "## Time Evolution: Defining the Initial State"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe5f760f",
   "metadata": {},
   "source": [
    "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):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b477e9d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This defines the initial state of the cavity; recall that we set N=5\n",
    "# (N,n) = (Number of Fock states in the cavity, initial state of the cavity)\n",
    "Psi0_cavity = basis(N, 0)  # start with ground state of cavity\n",
    "\n",
    "# This defines the initial state of the transmon qubit\n",
    "# Ham.eigenstates()[1][x] returns the xth state\n",
    "Psi_qubit = Ham.eigenstates()[1][1]  # first excited state\n",
    "\n",
    "# Tensor Product for Initial State of Composite System\n",
    "psi0 = tensor(Psi0_cavity, Psi_qubit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed41b078",
   "metadata": {},
   "source": [
    "## Time Evolution: Defining Occupation Operators "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e49b88e3",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "First let's define the operators for both the transmon and the CPW:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "911005cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create occupation operators for the transmon qubit (using HCPB class)\n",
    "s0 = Ham.eigenstates()[1][0]  # n=0 state\n",
    "n0 = s0 * s0.dag()  # occupation operator for n=0\n",
    "s1 = Ham.eigenstates()[1][1]  # n=1 state\n",
    "n1 = s1 * s1.dag()  # occupation operator for n=1\n",
    "s2 = Ham.eigenstates()[1][2]  # n=2 state\n",
    "n2 = s2 * s2.dag()  # occupation operator for n=2\n",
    "s3 = Ham.eigenstates()[1][3]  # n=3 state\n",
    "n3 = s3 * s3.dag()  # occupation operator for n=3\n",
    "\n",
    "# Create occupation operators for the CPW resonator (using QuTiP)\n",
    "c0 = basis(N, 0)  # n=0 state\n",
    "n0_cavity = c0 * c0.dag()  # occupation operator for n=0\n",
    "c1 = basis(N, 1)  # n=1 state\n",
    "n1_cavity = c1 * c1.dag()  # occupation operator for n=1\n",
    "c2 = basis(N, 2)  # n=2 state\n",
    "n2_cavity = c2 * c2.dag()  # occupation operator for n=2\n",
    "c3 = basis(N, 3)  # n=3 state\n",
    "n3_cavity = c3 * c3.dag()  # occupation operator for n=3\n",
    "c4 = basis(N, 4)  # n=4 state\n",
    "n4_cavity = c4 * c4.dag()  # occupation operator for n=4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdcb7977",
   "metadata": {},
   "source": [
    "Now let's construct the corresponding tensor products:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "068e93de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# occupation probability matrices for the qubit\n",
    "n0_qubit = tensor(qeye(N), n0)\n",
    "n1_qubit = tensor(qeye(N), n1)\n",
    "n2_qubit = tensor(qeye(N), n2)\n",
    "n3_qubit = tensor(qeye(N), n3)\n",
    "\n",
    "# occupation probability matrices for the CPW resonator\n",
    "n0_resonator = tensor(n0_cavity, qeye(4))\n",
    "n1_resonator = tensor(n1_cavity, qeye(4))\n",
    "n2_resonator = tensor(n2_cavity, qeye(4))\n",
    "n3_resonator = tensor(n3_cavity, qeye(4))\n",
    "n4_resonator = tensor(n4_cavity, qeye(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afeb2ccd",
   "metadata": {},
   "source": [
    "## Time Evolution: Steady State (No Dissipation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d833db2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# time will vary from t=0 to t=50s with 1001 increments\n",
    "tlist = np.linspace(0, 50, 1001)\n",
    "\n",
    "# initial list of collapse operators (empty here)\n",
    "c_ops = []\n",
    "\n",
    "# solve master equation and for output calculate the occupation numbers\n",
    "output = mesolve(\n",
    "    Ham_JC,\n",
    "    psi0,\n",
    "    tlist,\n",
    "    c_ops,\n",
    "    [n0_resonator, n1_resonator, n2_resonator, n0_qubit, n1_qubit, n2_qubit],\n",
    ")\n",
    "\n",
    "# tabulate the occupation probabilities for resonator and qubit\n",
    "occ0_resonator = output.expect[0]\n",
    "occ1_resonator = output.expect[1]\n",
    "occ2_resonator = output.expect[2]\n",
    "occ0_qubit = output.expect[3]\n",
    "occ1_qubit = output.expect[4]\n",
    "occ2_qubit = output.expect[5]\n",
    "\n",
    "fig, axes = plt.subplots(1, 1, figsize=(10, 5))\n",
    "\n",
    "plt.figure(1)\n",
    "plt.subplot(121)\n",
    "plt.title(\"Resonator\")\n",
    "plt.plot(tlist, occ0_resonator, \"b\", label=\"n=0\")\n",
    "plt.plot(tlist, occ1_resonator, \"r\", label=\"n=1\")\n",
    "plt.plot(tlist, occ2_resonator, \"g\", label=\"n=2\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.legend()\n",
    "plt.title(\"Qubit\")\n",
    "plt.plot(tlist, occ0_qubit, \"b\", label=\"n=0\")\n",
    "plt.plot(tlist, occ1_qubit, \"r\", label=\"n=1\")\n",
    "plt.plot(tlist, occ2_qubit, \"g\", label=\"n=2\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ac381e5",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31a61264",
   "metadata": {},
   "source": [
    "## Time Evolution: Observing Vacuum-Rabi Oscillations by Adding Dissipation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ad01f16",
   "metadata": {},
   "source": [
    "Here we'll add some photon dissipation in the CPW resonator and then evolve the system in time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e74b8d74",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = 0.1  # transmon-cavity coupling; set to be strong\n",
    "\n",
    "# rebuild Hamiltonian matrix of composite system\n",
    "Ham_JC = Ham_t + wc * a.dag() * a + g * Num_t * (a + a.dag())\n",
    "\n",
    "# time will vary from t=0 to t=100s\n",
    "tlist = np.linspace(0, 100, 1001)\n",
    "\n",
    "# photon relaxation rates of the transmon and the CPW cavity\n",
    "kappa = 0.05  # 0.005 # cavity photon relaxation rate\n",
    "n_th_a = 0.0  # avg number of thermal bath excitations\n",
    "\n",
    "# initial list of collapse operators\n",
    "c_ops = []\n",
    "\n",
    "# cavity relaxation\n",
    "rate = kappa * (1 + n_th_a)\n",
    "if rate > 0.0:\n",
    "    c_ops.append(sqrt(rate) * a)\n",
    "\n",
    "# solve master equation and for output calculate the occupation numbers\n",
    "output = mesolve(\n",
    "    Ham_JC,\n",
    "    psi0,\n",
    "    tlist,\n",
    "    c_ops,\n",
    "    [n0_resonator, n1_resonator, n2_resonator, n0_qubit, n1_qubit, n2_qubit],\n",
    ")\n",
    "\n",
    "# tabulate the occupation probabilities for resonator and qubit\n",
    "occ0_resonator = output.expect[0]\n",
    "occ1_resonator = output.expect[1]\n",
    "occ2_resonator = output.expect[2]\n",
    "occ0_qubit = output.expect[3]\n",
    "occ1_qubit = output.expect[4]\n",
    "occ2_qubit = output.expect[5]\n",
    "\n",
    "fig, axes = plt.subplots(1, 1, figsize=(10, 5))\n",
    "\n",
    "plt.figure(1)\n",
    "plt.subplot(121)\n",
    "plt.title(\"Resonator\")\n",
    "plt.plot(tlist, occ0_resonator, \"b\", label=\"n=0\")\n",
    "plt.plot(tlist, occ1_resonator, \"r\", label=\"n=1\")\n",
    "plt.plot(tlist, occ2_resonator, \"g\", label=\"n=2\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.legend()\n",
    "plt.title(\"Qubit\")\n",
    "plt.plot(tlist, occ0_qubit, \"b\", label=\"n=0\")\n",
    "plt.plot(tlist, occ1_qubit, \"r\", label=\"n=1\")\n",
    "plt.plot(tlist, occ2_qubit, \"g\", label=\"n=2\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d7e7793",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac082863",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cb2839e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 1, figsize=(10, 5))\n",
    "\n",
    "plt.figure(1)\n",
    "plt.subplot(121)\n",
    "plt.title(\"n=0\")\n",
    "plt.plot(tlist, occ0_resonator, \"b\", label=\"resonator\")\n",
    "plt.plot(tlist, occ0_qubit, \"r\", label=\"qubit\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.legend()\n",
    "plt.title(\"n=1\")\n",
    "plt.plot(tlist, occ1_resonator, \"b\", label=\"resonator\")\n",
    "plt.plot(tlist, occ1_qubit, \"r\", label=\"qubit\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c3b1735",
   "metadata": {},
   "source": [
    "## Another Example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e05e29e3",
   "metadata": {},
   "source": [
    "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). \n",
    "\n",
    "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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6c45023",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This defines the initial state of the cavity; recall that we set N=5\n",
    "# (N,n) = (Number of Fock states in the cavity, initial state of the cavity)\n",
    "Psi0_cavity = basis(N, 2)  # start with ground state of cavity\n",
    "\n",
    "# This defines the initial state of the transmon qubit\n",
    "# Ham.eigenstates()[1][x] returns the xth state\n",
    "Psi_qubit = Ham.eigenstates()[1][0]  # first excited state\n",
    "\n",
    "# Tensor Product for Initial State of Composite System\n",
    "psi0 = tensor(Psi0_cavity, Psi_qubit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74eeba2a",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1713ecc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# time will vary from t=0 to t=50s with 1001 increments\n",
    "tlist = np.linspace(0, 50, 1001)\n",
    "\n",
    "# initial list of collapse operators (empty here)\n",
    "c_ops = []\n",
    "\n",
    "# solve master equation and for output calculate the occupation numbers\n",
    "output = mesolve(\n",
    "    Ham_JC,\n",
    "    psi0,\n",
    "    tlist,\n",
    "    c_ops,\n",
    "    [n0_resonator, n1_resonator, n2_resonator, n0_qubit, n1_qubit, n2_qubit],\n",
    ")\n",
    "\n",
    "# tabulate the occupation probabilities for resonator and qubit\n",
    "occ0_resonator = output.expect[0]\n",
    "occ1_resonator = output.expect[1]\n",
    "occ2_resonator = output.expect[2]\n",
    "occ0_qubit = output.expect[3]\n",
    "occ1_qubit = output.expect[4]\n",
    "occ2_qubit = output.expect[5]\n",
    "\n",
    "fig, axes = plt.subplots(1, 1, figsize=(10, 5))\n",
    "\n",
    "plt.figure(1)\n",
    "plt.subplot(121)\n",
    "plt.title(\"Resonator\")\n",
    "plt.plot(tlist, occ0_resonator, \"b\", label=\"n=0\")\n",
    "plt.plot(tlist, occ1_resonator, \"r\", label=\"n=1\")\n",
    "plt.plot(tlist, occ2_resonator, \"g\", label=\"n=2\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.legend()\n",
    "plt.title(\"Qubit\")\n",
    "plt.plot(tlist, occ0_qubit, \"b\", label=\"n=0\")\n",
    "plt.plot(tlist, occ1_qubit, \"r\", label=\"n=1\")\n",
    "plt.plot(tlist, occ2_qubit, \"g\", label=\"n=2\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44eabfb6",
   "metadata": {},
   "source": [
    "And here's the time evolution with photon loss included in the CPW resonator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f40c95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# time will vary from t=0 to t=100s\n",
    "tlist = np.linspace(0, 100, 1001)\n",
    "\n",
    "# photon relaxation rates of the transmon and the CPW cavity\n",
    "kappa = 0.05  # 0.005 # cavity photon relaxation rate\n",
    "n_th_a = 0.0  # avg number of thermal bath excitations\n",
    "\n",
    "# initial list of collapse operators\n",
    "c_ops = []\n",
    "\n",
    "# cavity relaxation\n",
    "rate = kappa * (1 + n_th_a)\n",
    "if rate > 0.0:\n",
    "    c_ops.append(sqrt(rate) * a)\n",
    "\n",
    "# solve master equation and for output calculate the occupation numbers\n",
    "output = mesolve(\n",
    "    Ham_JC,\n",
    "    psi0,\n",
    "    tlist,\n",
    "    c_ops,\n",
    "    [n0_resonator, n1_resonator, n2_resonator, n0_qubit, n1_qubit, n2_qubit],\n",
    ")\n",
    "\n",
    "# tabulate the occupation probabilities for resonator and qubit\n",
    "occ0_resonator = output.expect[0]\n",
    "occ1_resonator = output.expect[1]\n",
    "occ2_resonator = output.expect[2]\n",
    "occ0_qubit = output.expect[3]\n",
    "occ1_qubit = output.expect[4]\n",
    "occ2_qubit = output.expect[5]\n",
    "\n",
    "fig, axes = plt.subplots(1, 1, figsize=(10, 5))\n",
    "\n",
    "plt.figure(1)\n",
    "plt.subplot(121)\n",
    "plt.title(\"Resonator\")\n",
    "plt.plot(tlist, occ0_resonator, \"b\", label=\"n=0\")\n",
    "plt.plot(tlist, occ1_resonator, \"r\", label=\"n=1\")\n",
    "plt.plot(tlist, occ2_resonator, \"g\", label=\"n=2\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.legend()\n",
    "plt.title(\"Qubit\")\n",
    "plt.plot(tlist, occ0_qubit, \"b\", label=\"n=0\")\n",
    "plt.plot(tlist, occ1_qubit, \"r\", label=\"n=1\")\n",
    "plt.plot(tlist, occ2_qubit, \"g\", label=\"n=2\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Occupation Probability\")\n",
    "plt.legend(title=\"Energy Levels\", loc=\"upper right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a41d9a4",
   "metadata": {},
   "source": [
    "Given a long enough time, we see that both the CPW and the qubit eventually settle into the n=0 states. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b54a77e7",
   "metadata": {},
   "source": [
    "## Toward Multi-Qubit Gates..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e469260f",
   "metadata": {},
   "source": [
    "With the above methodology, you can perform analysis for more complicated but realistic systems like two-qubit CR or iSWAP gates!!!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "metal_qamp_v2",
   "language": "python",
   "name": "metal_qamp_v2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}