single_qubit_resonator.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jaseung Ku, 07/01/2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run this once"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "code_folding": [
     11
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n"
     ]
    }
   ],
   "source": [
    "%matplotlib \n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from qutip import *\n",
    "\n",
    "def compute_transmon_f01(Ej, Ecm, ng=0, N=10):\n",
    "    \"\"\"\n",
    "    Compute f01 qubit frequency of transmon\n",
    "    Args:\n",
    "        Ej, Ec : in GHz\n",
    "        ng: unitless gate charge\n",
    "        N: number of charge basis = 2N + 1\n",
    "    Return:\n",
    "        f01 qubit frequency\n",
    "    \"\"\"\n",
    "    H = Qobj(np.diag(4 * Ec * (np.arange(-N,N+1)-ng)**2) \n",
    "             + 0.5 * Ej * (np.diag(-np.ones(2*N), 1) + np.diag(-np.ones(2*N), -1)))\n",
    "    eigenvalues = H.eigenenergies()\n",
    "    f01 = eigenvalues[1] - eigenvalues[0]\n",
    "    \n",
    "    return f01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-level qubit + resonator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def compute_energies_two_level_qubit(wc, wa, g, N=10, use_rwa=False):\n",
    "    \"\"\"\n",
    "    Args:\n",
    "        wc : cavity frequency in GHz\n",
    "        wa : qubit frequency in GHz\n",
    "        N: order of fock state\n",
    "    Return:\n",
    "        tuple of eigenvalues and eigenvectors\n",
    "    \"\"\"\n",
    "    a  = tensor(destroy(N), qeye(2))\n",
    "    sm = tensor(qeye(N), destroy(2))\n",
    "    nc = a.dag() * a\n",
    "    na = sm.dag() * sm\n",
    "    use_rwa = False\n",
    "\n",
    "    if use_rwa: \n",
    "        H = wc * nc + wa * na + g * (a.dag() * sm + a * sm.dag())\n",
    "    else:\n",
    "        H = wc * nc + wa * na + g * (a.dag() + a) * (sm + sm.dag())\n",
    "\n",
    "    # find the groundstate of the composite system\n",
    "    evals, ekets = H.eigenstates()\n",
    "\n",
    "    return evals, ekets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10\n",
    "wa = 5\n",
    "g = 0.1\n",
    "wc_list = np.linspace(4,6,101)\n",
    "\n",
    "eigenvals = np.array([compute_energies_two_level_qubit(wc, wa, g, N=N)[0][0:3] for wc in wc_list])\n",
    "E_10 = eigenvals[:,1] - eigenvals[:,0]\n",
    "E_20 = eigenvals[:,2] - eigenvals[:,0]\n",
    "diff = E_20 - E_10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "split = 0.200 GHz\n"
     ]
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(8,6))\n",
    "ax.plot(wc_list, E_10, wc_list, E_20)\n",
    "ax.plot(wc_list, diff)\n",
    "ax.set_xlabel(r'Cavity Frequency (GHz)', fontsize=18)\n",
    "ax.set_ylabel(r'Frequency (GHz)', fontsize=18)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid()\n",
    "ax.set_title(f'Two-level qubit + resonator\\n Qubit Freq = {wa} GHz, g = {g} GHz', fontsize=18)\n",
    "print(f'split = {np.min(diff):.3f} GHz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transmon + resonator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def compute_energies_transmon_resonator(Ec, Ej, g01, wc, N=10, ng=0, M_res=3):\n",
    "    \"\"\"\n",
    "    Solve Hamiltoninian of transmon + resonator to compute eigenvalues/eigenvectors\n",
    "    Args:\n",
    "        Ec, Ej: Charging energy, Josephson energy in GHz\n",
    "        g01 : 0-1 coupling strength in GHz, See (3.3) and (3.4) in Koch's paper\n",
    "        wc: Cavity frequency in GHz\n",
    "        N : number of change basis = 2N + 1\n",
    "        M_res : order of fock state for cavity\n",
    "    Return:\n",
    "        tuple of eigenvalues and eigenvectors\n",
    "    \"\"\"\n",
    "    \n",
    "    a = destroy(M_res)\n",
    "    \n",
    "    H_tr = tensor(Qobj(np.diag(4 * Ec * (np.arange(-N,N+1)-ng)**2) \\\n",
    "        + 0.5 * Ej * (np.diag(-np.ones(2*N), 1) + np.diag(-np.ones(2*N), -1))), qeye(M_res))\n",
    "    \n",
    "    H_r = tensor(qeye(2*N+1), wc * a.dag() * a)\n",
    "    H_int = g01/(np.power(Ej/Ec/8, 1/4)/np.sqrt(2)) *  tensor( Qobj(np.diag(np.arange(-N,N+1))), a.dag() + a)\n",
    "    \n",
    "    H = H_tr + H_r + H_int\n",
    "    \n",
    "    evals, ekets = H.eigenstates()\n",
    "\n",
    "    return evals, ekets\n",
    "\n",
    "def calculate_numerical_mat_element_01(Ec, Ej, N=10, ng=0):\n",
    "    \"\"\" calculate <0|n|1>\"\"\"\n",
    "    H = Qobj(np.diag(4 * Ec * (np.arange(-N,N+1)-ng)**2) \n",
    "                 + 0.5 * Ej * (np.diag(-np.ones(2*N), 1) + np.diag(-np.ones(2*N), -1)))\n",
    "    eigenvals, eigenvecs = H.eigenstates()\n",
    "    n_op = Qobj(np.diag((np.arange(-N,N+1)-ng))) \n",
    "    matrix_element_01 = n_op.matrix_element(eigenvecs[0], eigenvecs[1])\n",
    "    return matrix_element_01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "split = 0.097 GHz\n"
     ]
    }
   ],
   "source": [
    "# plot energy diff vs wc\n",
    "N, M_res = 10, 3\n",
    "Ec, Ej = 0.35, 19.5\n",
    "g01 = 0.05 \n",
    "\n",
    "wc_list = compute_transmon_f01(Ej, Ec) + np.linspace(-g01*20 ,g01*20, 101)\n",
    "\n",
    "eigenvals = np.array([compute_energies_transmon_resonator(Ec, Ej, g01, wc, N=N, ng=0, M_res=M_res)[0][0:3] for wc in wc_list])\n",
    "E_10 = eigenvals[:,1] - eigenvals[:,0]  # energy difference between 0 and 1 \n",
    "E_20 = eigenvals[:,2] - eigenvals[:,0]  # energy difference between 0 and 2 \n",
    "diff = E_20 - E_10\n",
    "\n",
    "fig, ax = plt.subplots(1,1, figsize=(8,6))\n",
    "ax.plot(wc_list, E_10, wc_list, E_20)\n",
    "ax.plot(wc_list, diff, label='diff')\n",
    "ax.set_xlabel(r'Cavity Frequency (GHz)', fontsize=18)\n",
    "ax.set_ylabel(r'Frequency (GHz)', fontsize=18)\n",
    "ax.tick_params(axis='x', labelsize=18)\n",
    "ax.tick_params(axis='y', labelsize=18)\n",
    "ax.grid()\n",
    "ax.legend(fontsize=18)\n",
    "ax.set_title(f'Transmon + resonator\\nEc = {Ec}, Ej = {Ej} GHz, f01 = {compute_transmon_f01(Ej, Ec):.3f} GHz, g01 = {g01} GHz', fontsize=18)\n",
    "print(f'split = {np.min(diff):.3f} GHz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Transmon + cavity\\nEc = 0.35, Ej = 19.5 GHz, Cavity Freq = f01 = 7.020 GHz')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# g01 vs split comparion\n",
    "N, M_res = 10, 3\n",
    "Ec, Ej = 0.35, 19.5 # Ej = 15, 19.5, 25  --> qubit freq = 6.11, 7, 8 GHz\n",
    "g01_list = np.linspace(0.005, 0.05, 11)\n",
    "wc = compute_transmon_f01(Ej, Ec)\n",
    "\n",
    "diff_min_list = []\n",
    "for g in g01_list:\n",
    "    f01 = compute_transmon_f01(Ej, Ec)\n",
    "    wc_list = f01 + np.linspace(-g*2, g*2, 101)\n",
    "    eigenvals = np.array([compute_energies_transmon_resonator(Ec, Ej, g, wc, N=N, ng=0, M_res=M_res)[0][0:3] for wc in wc_list])\n",
    "    E_10 = eigenvals[:,1] - eigenvals[:,0]\n",
    "    E_20 = eigenvals[:,2] - eigenvals[:,0]\n",
    "    diff = E_20 - E_10\n",
    "    diff_min_list.append(np.min(diff))\n",
    "\n",
    "diff_min_list = np.array(diff_min_list)\n",
    "g01_list = g01_list/ (np.power(Ej/Ec/8, 1/4)/np.sqrt(2))\\\n",
    "        *np.absolute(calculate_numerical_mat_element_01(Ec,Ej))\n",
    "\n",
    "ratio_split_2g = np.round(diff_min_list/(2*g01_list), decimals=3)\n",
    "avg_ratio_split_2g = np.mean(ratio_split_2g)\n",
    "\n",
    "fig1, ax1 = plt.subplots(1,1,figsize=(8,6))\n",
    "ax1.plot(g01_list,ratio_split_2g ,'r-o' )\n",
    "ax1.grid()\n",
    "ax1.tick_params(axis='x', labelsize=18)\n",
    "ax1.tick_params(axis='y', labelsize=18)\n",
    "ax1.set_xlabel(r'g (GHz)', fontsize=18)\n",
    "ax1.legend([r'split/($2g_{01}$)'], fontsize=20)\n",
    "ax1.set_ylim([avg_ratio_split_2g * 0.98, avg_ratio_split_2g*1.02])\n",
    "ax1.set_title(f'Transmon + cavity\\nEc = {Ec}, Ej = {Ej} GHz, Cavity Freq = f01 = {compute_transmon_f01(Ej, Ec):.3f} GHz', fontsize=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f01 = 6.924 GHz\n",
      "<0|n|1> = 1.141258156873071\n"
     ]
    }
   ],
   "source": [
    "Ec, Ej = 0.35,19\n",
    "print(f'f01 = {compute_transmon_f01(Ej, Ec):.3f} GHz')\n",
    "print(f'<0|n|1> = {np.power((Ej/Ec/8), 1/4)/np.sqrt(2)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find f01 and resonator frequency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def compute_transmon_resonator_frequency(Ec, Ej, g01, wc, N=10, ng=0, M_res=3):\n",
    "    \"\"\"\n",
    "    Solve Hamiltoninian of transmon + resonator to compute eigenvalues/eigenvectors\n",
    "    Args:\n",
    "        Ec, Ej: Charging energy, Josephson energy in GHz\n",
    "        g01 : 0-1 coupling strength in GHz, See (3.3) and (3.4) in Koch's paper\n",
    "        wc: Cavity frequency in GHz\n",
    "        N : number of change basis = 2N + 1\n",
    "        M_res : order of fock state for cavity\n",
    "    Return:\n",
    "        tuple of qubit f01 and cavity frequency in GHz\n",
    "    \"\"\"\n",
    "    \n",
    "    eigenvals = compute_energies_transmon_resonator(Ec, Ej, g01, wc, N=N, ng=0, M_res=M_res)[0][0:3]\n",
    "    E_10 = eigenvals[1] - eigenvals[0]\n",
    "    E_20 = eigenvals[2] - eigenvals[0]\n",
    "        \n",
    "    return E_10, E_20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6.960129899714568, 7.059634798210018)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ec, Ej = 0.35, 19.5\n",
    "wc = 7 \n",
    "g01 = 0.05\n",
    "compute_transmon_resonator_frequency(Ec, Ej, g01, wc) # one corresponds to qubit f01 and the other to cavity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare $\\langle 0|\\hat{n}|1\\rangle$ and $(E_j/E_c/8)^{1/4}/\\sqrt(2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f01 = 7.020 GHz\n",
      "<0|n|1> = -1.119\n",
      "(Ej/Ec/8)^(1/4)/sqrt(2) = 1.149\n"
     ]
    }
   ],
   "source": [
    "Ec, Ej = 0.35, 19.5\n",
    "print(f'f01 = {compute_transmon_f01(Ej,Ec):.3f} GHz')\n",
    "print(f'<0|n|1> = {calculate_numerical_mat_element_01(Ec,Ej).real:.3f}')\n",
    "print(f'(Ej/Ec/8)^(1/4)/sqrt(2) = {np.power(Ej/Ec/8, 1/4)/np.sqrt(2):.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}