{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Optics 2 – Exercise 4 – Dark states\n",
    "\n",
    "## Stimulated Raman adiabatic passage\n",
    "\n",
    "<img src=\"3level.png\">\n",
    "The total Hamiltonian is\n",
    "$$H = -\\hbar\\bigl(\\omega_{01}|g_1\\rangle\\langle g_1|+\\omega_{02}|g_2\\rangle\\langle g_2|\\bigr)+\\frac{\\hbar\\Omega_1}{2}\\bigl(\\sigma_1e^{i\\omega_1t}+\\sigma_1^\\dagger e^{-i\\omega_1t}\\bigr)+\\frac{\\hbar\\Omega_2}{2}\\bigl(\\sigma_2e^{i\\omega_2t}+\\sigma_2^\\dagger e^{-i\\omega_2t}\\bigr).$$\n",
    "We apply two pulses with frequencies of $\\omega_1$ and $\\omega_2$ with a Gaussian shape\n",
    "$$\\Omega_{1,2}(t)=\\Omega_{1,2}\\exp\\biggl(-\\biggl(\\frac{t-t_{1,2}}{T}\\biggr)^2\\biggr).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qutip import *\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" States \"\"\"\n",
    "e = basis(3,0)\n",
    "g1 = basis(3,1)\n",
    "g2 = basis(3,2)\n",
    "n_e = e*e.dag()\n",
    "n_g1 = g1*g1.dag()\n",
    "n_g2 = g2*g2.dag()\n",
    "\n",
    "\"\"\" Energy levels \"\"\"\n",
    "ωe = 0* 2*np.pi\n",
    "ωg1 = -1*2*np.pi\n",
    "ωg2 = -0.8*2*np.pi\n",
    "\n",
    "\"\"\" Laser frequencies and detuning \"\"\"\n",
    "Δ = 0.1* 2*np.pi\n",
    "ω1 = ωe-ωg1-Δ\n",
    "ω2 = ωe-ωg2-Δ\n",
    "\n",
    "\"\"\" Rabi frequencies \"\"\"\n",
    "Ω1 = 1e6* 2*np.pi\n",
    "Ω2 = 1e6* 2*np.pi\n",
    "θ = np.arctan(Ω2/Ω1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" Pulse parameters \"\"\"\n",
    "\"\"\" Tweak these parameters to achieve the adiabaticity condition \"\"\"\n",
    "Δt = 0\n",
    "T = 1e-6\n",
    "\n",
    "def Gauss_1(t):\n",
    "    return np.exp(-(t+Δt)**2/T**2)\n",
    "def Gauss_2(t):    \n",
    "    return np.exp(-(t-Δt)**2/T**2)\n",
    "def Omega_1(t,args):\n",
    "    return Gauss_1(t)*Ω1/2.0*np.exp(-1j*ω1*t)\n",
    "def Omega_1_dag(t,args):\n",
    "    return Gauss_1(t)*Ω1/2.0*np.exp(1j*ω1*t)\n",
    "def Omega_2(t,args):\n",
    "    return Gauss_2(t)*Ω2/2.0*np.exp(-1j*ω2*t)\n",
    "def Omega_2_dag(t,args):\n",
    "    return Gauss_2(t)*Ω2/2.0*np.exp(1j*ω2*t)\n",
    "\n",
    "H0 = ωg1*n_g1+ωg2*n_g2+ωe*n_e\n",
    "H=[H0,[e*g1.dag(),Omega_1],[g1*e.dag(),Omega_1_dag],[e*g2.dag(),Omega_2],[g2*e.dag(),Omega_2_dag]]\n",
    "psi0 = g1\n",
    "\n",
    "t=np.linspace(-5*T,5*T,1000)\n",
    "pulse1 = Ω1/(2*np.pi)*Gauss_1(t)\n",
    "pulse2 = Ω2/(2*np.pi)*Gauss_2(t)\n",
    "result = mesolve(H,psi0,t,[],[n_g1, n_g2, n_e])\n",
    "\n",
    "fig,ax = plt.subplots(2,1,figsize=(10,10),sharex=True,gridspec_kw=dict(hspace=0.1))\n",
    "ax[0].plot(t,pulse1,label='ω1',lw=3.0)\n",
    "ax[0].plot(t,pulse2,label='ω2',lw=3.0)\n",
    "ax[1].plot(t,result.expect[0],label='N_g1',lw=3.0)\n",
    "ax[1].plot(t,result.expect[1],label='N_g2',lw=3.0)\n",
    "ax[1].plot(t,result.expect[2],label='N_e',lw=3.0)\n",
    "ax[0].legend(fontsize=14)\n",
    "ax[1].legend(fontsize=14)\n",
    "ax[0].set_ylabel('Pulse ntensity',fontsize=16)\n",
    "ax[1].set_ylabel('Population',fontsize=16)\n",
    "ax[1].set_xlabel('Time (s)',fontsize=16);"
   ]
  },
  {
   "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}