Rabi_Oscillation_Python.ipynb

{
  "cells": [
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "2Ioy9wScEdFc"
      },
      "source": [
        "#Modelling Rabi Oscillations in Three Level Optical Systems"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "zX_M8CXJpew5"
      },
      "source": [
        "The following document outlines the theory behind and implements the procedure for modelling **Rabi oscillations** in the excitation probability of a **biexciton-exciton** scheme invoked in a **quantum dot** excited by a laser. This model will make use of **QuTiP**, an python SDK written for modelling open quantum systems.\n",
        "\n",
        "For a more thorough explaination of the theory outlined in this document and the decisions made to optimize this simulation, see this document."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "u8BwzWwwEXMd"
      },
      "outputs": [],
      "source": [
        "!pip install qutip"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "PdpRXW5arfIB"
      },
      "source": [
        "This project will take advantage of QuTiP's built in Master Equation Solver, which is equipt to solve a given **Lindblad master equation** of the form:\n",
        "\n",
        "$$\\dot{\\rho}(t) = -\\frac{i}{\\hbar}[H(t),\\rho(t)] + \\sum \\frac{1}{2}[2A_n\\rho(t)A_n^\\dagger - \\rho(t)A_n^\\dagger A_n - A_n^\\dagger A_n\\rho(t)]$$\n",
        "\n",
        "where $A_n$ are collapse operators and $\\rho(t)$ is the density matrix of the system at t = 0 (implying $\\rho$ is the time evolved system) [1].\n",
        "\n",
        "In order to see Rabi Oscillations in the aforemention excitation scheme, the system must be subject to laser pulses\n",
        " of varying intensity $\\Omega_{max}$. In short, it is expected this range of frequenies, when a constant pulse length is held, will invoke a **sinusoidal excitation probability** as $\\Omega_{max}$ increases, fluctuating periodically between maxima and minima.\n",
        "\n",
        "This system, when introduced to the system, will drive the transition between ground state $|g\\rangle$, exciton state $|e\\rangle$, and biexciton state $|b\\rangle$.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "aj5ErKoGEaEw"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import qutip as qt\n",
        "from qutip import Options, qutrit_basis, dag, fock, about, basis, destroy, mesolve, qeye, sigmaz, tensor\n",
        "\n",
        "#Parameters and Values\n",
        "#n = 50\n",
        "gamma_b = 0.5 #biexciton decay constant (THz)\n",
        "gamma_x = 1 #exiton decay constant (THz)\n",
        "globaldt = 0.15/gamma_b\n",
        "minpointsperpulse = 5.0\n",
        "tmin = 0.001\n",
        "tmax = 40*gamma_b\n",
        "\n",
        "#t_list = np.linspace(0.01, 40*gamma_b, n) #Lifetime of the simulation\n",
        "t_i = 3.5 #Centre of the Gaussian pulse (offset) (ns)\n",
        "sigma = 4*gamma_b #Pulse Length (ns) - multiply by gamma?\n",
        "\n",
        "if sigma/minpointsperpulse > globaldt:\n",
        "    t_list = np.linspace(tmin, tmax, int((tmax - tmin)/globaldt))\n",
        "    n = int((tmax - tmin)/globaldt)\n",
        "    # otherwise, increase resolution inside the pulse only.\n",
        "else:\n",
        "    tlist1 = np.linspace(tmin, sigma*t_i*2.0, \\\n",
        "                             int(minpointsperpulse*(t_i*2.0 - tmin)))[:-1]\n",
        "    tlist2 = np.linspace(sigma*t_i*2.0, tmax, \\\n",
        "                             int((tmax - sigma*t_i*2.0)/globaldt))\n",
        "    t_list = np.append(tlist1, tlist2)\n",
        "\n",
        "omega_max = np.linspace(0.001, np.pi, n) #Maximum intensity of the laser pulse\n",
        "\n",
        "#System states\n",
        "gnd = fock(3,0)\n",
        "ex = fock(3,1)\n",
        "biex = fock(3,2)"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "WRn1vDTC1nie"
      },
      "source": [
        "In order to accurately translate a pulse adhering to the above paramters, the **Rabi Frequency**  $\\Omega_R$ must be calculated and the pulse must be normalized. For this excitation scheme, a **Gaussian pulse** is applied, which allows for accurate intensity variation. From the supplementary material to Ref. [1],\n",
        "\n",
        "$$\\Omega_R = \\Omega_{max} e^{{-ln(2)(t-t_i)^2}\\over\\sigma^2}$$\n",
        "\n",
        "Thus, a distinct $\\Omega_R$ list must be calculated for each constant $\\Omega_{max}$ value to be tested.\n",
        "\n",
        "To normalize, a normalization factor can be calculated for each Gaussian pulse as the inverse of the area beneath the gaussian.\n",
        "\n",
        "$$F_{norm} = (\\int_{0}^{t_f}\\Omega_Rdt)^{-1}$$\n",
        "\n",
        "This factor, when applied to the respective gaussian, will result in a normalized pulse.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "i5GUof2EFB_r"
      },
      "outputs": [],
      "source": [
        "#Frequency\n",
        "omega_list = []\n",
        "final_pulse = []\n",
        "\n",
        "for freq in range(0,n):\n",
        "  omega = (np.exp(-(np.log(2)*((t_list-t_i)**2))/(sigma**2))) #Rabi frequency\n",
        "  omega_list.append(omega)\n",
        "  #final_pulse.append(omega)\n",
        "\n",
        "for i in range(0,n):\n",
        "  norm_val = 1/(np.trapz((omega_list[i]), t_list))\n",
        "  norm = norm_val*omega_max[i]*sigma\n",
        "  normalized = (omega_list[i]*norm)\n",
        "  final_pulse.append(normalized)"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "qzemehnHHzKI"
      },
      "source": [
        "To properly model the evolution of the system, four **Louvillian operators** must be defined. Two operators ($L_1$ and $L_2$) describe spontaneous photon emission leading to transitions between adjascent states and the other two ($L_3$ and $L_4$) describe dephasing between adjascent states [1].\n",
        "\n",
        "$$L_1 = {\\gamma_b \\over 2}D(|b\\rangle\\langle x|) $$\n",
        "$$L_2 = {\\gamma_x \\over 2}D(|x\\rangle\\langle g|) $$\n",
        "$$L_3 = {\\gamma_x^d \\over 2}D(|b\\rangle\\langle b| - |x\\rangle\\langle x|)$$\n",
        "$$L_4 = {\\gamma_b^d \\over 2}D(|x\\rangle\\langle x| - |g\\rangle\\langle g|)$$\n",
        "\n",
        "Where $\\gamma_b$ and $\\gamma_x$ are the spontaneous emission rates for the biexiton and exiton, and  $\\gamma_x^d$ and $\\gamma_b^d$ are the dephasing rates between of the biexiton-exiton and exiton-ground states, respectively [1].\n",
        "\n",
        "Also, $D = 2A_n\\rho (t)A_n^† - \\rho(t)A_n^†A_n - A_n^†A_n\\rho (t)$ and can be omitted from the operator definition due to the capabilities of the QuTiP master equation solver, which assumes operators will be input in this form.\n",
        "\n",
        "$\\gamma_x^d$ and $\\gamma_b^d$ can either be **time dependent** or **constant**. In the case of time dependence, these rates can be calculated in relation to the aforementioned guassian pulse.\n",
        "\n",
        "$$\\gamma_x^d = \\gamma_b^d = \\gamma_{I_0} \\Omega_{div}^{n_p}$$\n",
        "\n",
        "In this expression, $\\Omega_{div} = {\\Omega_R \\over \\Omega_{max}}$ and $\\gamma_{I_0}$ is the rate of amplitude dependent dephasing. $n_p$ is a scaling factor that dictates the geometric behaviour of the system. Assuming $n_p = 2$ allows for a linear relationship suitable for this application [1]."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "LHv72p-PFQyM"
      },
      "outputs": [],
      "source": [
        "#Collapse Operators\n",
        "#hard code decay rate\n",
        "gamma_I0 = 0.05 #Amplitude of intensity-dependent dephasing rate\n",
        "n_p = 2\n",
        "\n",
        "#Run 0 edge case\n",
        "L1 = (gamma_b/2)*(ex*biex.dag()) #biex-ex spontaneous decay\n",
        "L2 = (gamma_x/2)*(gnd*ex.dag()) #ex-gnd spontaneous decay\n",
        "L3 = []\n",
        "L4 = []\n",
        "\n",
        "#extrapolate gamma values?\n",
        "for i in (range(0,n)):\n",
        "  omega_working = final_pulse[i]\n",
        "  holder3 = []\n",
        "  holder4 = []\n",
        "  for j in (range(0,n)):\n",
        "    coeff = ((gamma_I0*((omega_working[j]/omega_max[i])**n_p))/2)\n",
        "    # coeff4 = 1/211\n",
        "    # coeff3 = 1/119\n",
        "    val3 = coeff*((biex*biex.dag())-(ex*ex.dag())) #biex-ex dephasing\n",
        "    val4 = coeff*((ex*ex.dag())-(gnd*gnd.dag())) #ex-gnd dephasing\n",
        "    holder3.append(val3)\n",
        "    holder4.append(val4)\n",
        "  L3.append(holder3)\n",
        "  L4.append(holder4)"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "OuzntJkRz5Mo"
      },
      "source": [
        "There are several ways in which the system Hamiltonian can be written depending on the frame and picture chosen. Optimal results found when the **rotating frame** is considered. The rotating frame assumes a frame of reference that rotates at the laser's frequency, thus suppressing rapidly rotating terms. With this in mind, the system can be described by the following Hamiltonian, as defined in Ref. [2]:\n",
        "$$H_{sys} = H_0 + H_E(t)$$\n",
        "where,\n",
        "$$H_0 = {E_b\\over 2} |i\\rangle \\langle i|$$\n",
        "and $$H_E(t) = {(\\mu E(t))^2\\over E_b}(|g\\rangle\\langle e|+|e\\rangle\\langle g|)$$\n",
        "\n",
        "such that $E(t)$ refers to the pulse driving the system, $\\mu$ refers to the coupling strength, and $E_b$ refers to the binding energy.\n",
        "\n",
        "To further simplify, choosing to implement this Hamiltonian in the **interaction picture** allows $H_0$ to be eliminated, since the intermediate step is not necessary in this picture to accurately represent the system.\n",
        "\n",
        "Once the Hamiltonian has been defined for each Gaussian pulse, the master equation solver can be implemented. This function accepts arguments as **mesolve(Hamiltonian, initial state, time interval, Collapse Operators, Evolution Operators, ...).**\n",
        "\n",
        "Once implemented, the expectation values for each Hamiltonian (and, thus, each state of the system) can be obtained via **output.expect[$x$]** where $x$ is the index corresponding to the desired set of values, dictated by the order of the evolution operators (ie. if the evolution operators are defined as [exciton, biexciton], the exciton expectation values can be obtained using output.expect[0]).\n",
        "\n",
        "**Note:** Since QuTiP assumes $\\hbar = 1$, $\\hbar$ can be omitted from the Hamiltonian definition."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lGetDbp5FUMP"
      },
      "outputs": [],
      "source": [
        "#Hamiltonian\n",
        "C = 1 #dipole coupling strength and binding energy\n",
        "deltax = 2*np.pi*0.335\n",
        "H_list = []\n",
        "for i in range(0,n):\n",
        "  H0 = qeye(3)\n",
        "  omega_working = final_pulse[i]\n",
        "  holder = []\n",
        "  #H_I = (C/2)*(biex*gnd.dag() + gnd*biex.dag()) #fischer\n",
        "  #for j in range(0,n):\n",
        "  H_I = ((1/2)*(gnd*ex.dag()+ ex*gnd.dag() + ex*biex.dag() + biex*ex.dag())) + deltax*(ex*ex.dag()) #huber\n",
        "    #holder.append(H_I)\n",
        "  #H_list.append(holder)\n",
        "  H_list.append([H_I,omega_working])\n",
        "\n",
        "#Evolution\n",
        "master_output_ex = []\n",
        "master_output_biex = []\n",
        "\n",
        "for val in range(0,n):\n",
        "  H_working = H_list[val]\n",
        "  L3_working = L3[val]\n",
        "  L4_working = L4[val]\n",
        "  holder_ex = []\n",
        "  holder_biex = []\n",
        "\n",
        "  output = mesolve(H_working, gnd, t_list, [L1,L2,L3_working,L4_working], [ex*ex.dag(),biex*biex.dag()])\n",
        "\n",
        "  master_output_ex.append(output.expect[0])\n",
        "  master_output_biex.append(output.expect[1])"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "ytcJL9VnOUmx"
      },
      "source": [
        "To quickly revisit the theory behind this model, it is expected that the maximum expectation value for each Gaussian pulse will vary periodically. As such, the area beneath the curves formed by these expectation values will also vary periodically. Thus, the system's Rabi oscillations can be obtained via scaled integration, where:\n",
        "\n",
        "$$P_b(t_f) = \\gamma_b \\int_0^{t_f} \\langle b| \\rho(t)|b\\rangle dt$$\n",
        "$$P_x(t_f) = \\gamma_x \\int_0^{t_f} \\langle x| \\rho(t)|x\\rangle dt$$\n",
        "\n",
        "denote the occupation probability of the biexciton and exiton, respectively [1]."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PL1yoEfCFWyX"
      },
      "outputs": [],
      "source": [
        "int_freq_ex = []\n",
        "int_freq_biex = []\n",
        "max_val_ex = []\n",
        "max_val_biex = []\n",
        "\n",
        "for vals in range(0,n):\n",
        "  integration_ex = gamma_x*np.trapz(master_output_ex[vals], t_list)\n",
        "  max_ex = max(master_output_ex[vals])\n",
        "  integration_biex = gamma_b*np.trapz(master_output_biex[vals], t_list)\n",
        "  max_biex = max(master_output_biex[vals])\n",
        "\n",
        "  int_freq_ex.append(integration_ex)\n",
        "  int_freq_biex.append(integration_biex)\n",
        "  max_val_ex.append(max_ex)\n",
        "  max_val_biex.append(max_biex)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OKlfTWlbFkAF"
      },
      "outputs": [],
      "source": [
        "#plot\n",
        "g = n-1\n",
        "print(n)\n",
        "fig, ax = plt.subplots(figsize=(8,5))\n",
        "ax.plot(omega_max, master_output_ex[g], 'r')\n",
        "ax.plot(omega_max, master_output_biex[g], 'b')\n",
        "#ax.plot(omega_max, final_pulse[g], 'g')\n",
        "ax.set_xlabel('frequency')\n",
        "ax.set_ylabel('Occupation probability')\n",
        "ax.set_title('Vacuum Rabi oscillations')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "AOSfsfyNubUv",
        "outputId": "322cf9cd-ff1f-459d-f5bb-38c1393312b2"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1.5997036641130988\n"
          ]
        }
      ],
      "source": [
        "val = gamma_b*np.trapz(master_output_biex[n-1],t_list)\n",
        "print(val)"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "CDidPnSsbTlE"
      },
      "source": [
        "##References\n",
        "[1] T. Huber et al., \"Coherence and Degree of Time-Bin Entanglement from Quantum Dots,\" *Physical Review B*, vol. 93, *(20)*, 2016.\n",
        "\n",
        "[2] K. Fischer et al., \"Dynamical Modelling of Pulsed Two-Photon Interference,\" *New Journal of Physics*, vol. 18, *(11)*, 2017."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}