diff --git a/Exercise_4/3level.png b/Exercise_4/3level.png
new file mode 100644
index 0000000..afdca46
Binary files /dev/null and b/Exercise_4/3level.png differ
diff --git a/Exercise_4/exercise4_dark_states.pdf b/Exercise_4/exercise4_dark_states.pdf
new file mode 100644
index 0000000..65e2bd4
Binary files /dev/null and b/Exercise_4/exercise4_dark_states.pdf differ
diff --git a/Exercise_4/stirap.ipynb b/Exercise_4/stirap.ipynb
new file mode 100644
index 0000000..3754ef5
--- /dev/null
+++ b/Exercise_4/stirap.ipynb
@@ -0,0 +1,133 @@
+{
+ "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",
+    "\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": 1,
+   "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": 2,
+   "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-Δ"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\" Pulse parameters \"\"\"\n",
+    "\"\"\" Tweak these parameters to achieve the adiabaticity condition \"\"\"\n",
+    "Ω1 = 1e6* 2*np.pi   # Rabi frequency between |g1> and |e>\n",
+    "Ω2 = 1e6* 2*np.pi   # Rabi frequency between |g2> and |e>\n",
+    "Δt = 0              # Laser pulse delay\n",
+    "T = 1e-6             # Pulse width\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 (ipykernel)",
+   "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.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}