{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c6dc2355",
   "metadata": {},
   "source": [
    "\n",
    "# EIT Playground: Susceptibility, Group Velocity, and Pulse Propagation\n",
    "\n",
    "This notebook lets you explore **Electromagnetically Induced Transparency (EIT)** in a three–level Λ system and compare it with a simple **gain-doublet** toy model.\n",
    "You can tune the **interaction strength** (control-field Rabi frequency \\(\\Omega_c\\)), visualize the real and imaginary parts of the susceptibility \\(\\chi\\),\n",
    "and compute the **group velocity**. A time‑domain **pulse propagation** demo shows slow‑light delay (and fast‑light advancement in the toy model).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3aa88b2",
   "metadata": {},
   "source": [
    "## Derivation sketch (weak-probe Λ-EIT)\n",
    "\n",
    "We consider a three-level Λ system with probe coupling $|1\\rangle \\leftrightarrow |3\\rangle$ and control coupling $|2\\rangle \\leftrightarrow |3\\rangle$.\n",
    "\n",
    "In the **weak-probe** limit (linear response in the probe), the steady-state optical Bloch equations give a probe polarization  \n",
    "$P = \\epsilon_0 \\chi E_p$ with susceptibility of the form  \n",
    "\n",
    "$$\n",
    "\\chi(\\Delta) \\propto \n",
    "\\frac{\\Delta + i\\gamma_{31}}\n",
    "{(\\Delta + i\\gamma_{31})(\\Delta + i\\gamma_{21}) - \\Omega_c^2/4}\\,.\n",
    "$$\n",
    "\n",
    "Here  \n",
    "- $\\Delta = \\omega - \\omega_{31}$ is the probe detuning (angular),  \n",
    "- $\\gamma_{31}$ is the optical decoherence (HWHM) on the excited state,  \n",
    "- $\\gamma_{21}$ is the ground-state coherence dephasing,  \n",
    "- $\\Omega_c$ is the **control-field Rabi frequency** (our “interaction strength”).\n",
    "\n",
    "---\n",
    "\n",
    "### From susceptibility to refractive index and group velocity\n",
    "For dilute media with small $|\\chi|$,\n",
    "\n",
    "$$\n",
    "n(\\omega) \\approx 1 + \\frac{1}{2}\\chi'(\\omega),\n",
    "\\qquad\n",
    "\\alpha(\\omega) \\propto \\chi''(\\omega),\n",
    "$$\n",
    "\n",
    "where $\\chi'$ and $\\chi''$ are the real (dispersion) and imaginary (absorption/gain) parts.  \n",
    "The **group index** is\n",
    "\n",
    "$$\n",
    "n_g(\\omega) = n(\\omega) + \\omega \\frac{dn}{d\\omega}\n",
    "\\quad \\Rightarrow \\quad\n",
    "v_g = \\frac{c}{n_g}.\n",
    "$$\n",
    "\n",
    "Near two-photon resonance, EIT produces a sharp, normal dispersion slope (large positive $dn/d\\omega$) → **slow light**.  \n",
    "Conversely, engineered anomalous dispersion (e.g., a gain doublet) can give $n_g < 1$ or even negative → **fast light** (no FTL signalling).\n",
    "\n",
    "---\n",
    "\n",
    "### Bandwidth–delay trade-off\n",
    "Increasing the control strength $\\Omega_c$ broadens the EIT window but reduces the slope;  \n",
    "reducing $\\Omega_c$ steepens the slope but narrows the window and increases residual absorption via $\\gamma_{21}$.  \n",
    "Practical slow-light therefore balances $\\Omega_c$, $\\gamma_{21}$, and $\\chi_0$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5104eaf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact, FloatSlider, Dropdown, VBox, HBox, Layout, Checkbox\n",
    "from IPython.display import display, Markdown\n",
    "\n",
    "# Constants\n",
    "c = 299_792_458.0  # m/s\n",
    "pi = np.pi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "eaa9914e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def chi_eit(Delta, chi0, gamma31, gamma21, Oc):\n",
    "    \"\"\"\n",
    "    Weak-probe EIT susceptibility for a 3-level Λ system (angular units).\n",
    "    chi0 is a dimensionless scale factor = N|mu13|^2/(eps0*hbar).\n",
    "    \"\"\"\n",
    "    a = Delta + 1j*gamma31\n",
    "    b = Delta + 1j*gamma21\n",
    "    denom = a*b - (Oc**2)/4.0\n",
    "    return chi0 * (a / denom)\n",
    "\n",
    "\n",
    "def chi_gain_doublet(Delta, chi0, gamma, sep, gain=True):\n",
    "    \"\"\"\n",
    "    Toy 'fast light' model: two Lorentzian features separated by 'sep' (angular).\n",
    "    If gain=True, Im[chi] < 0 between peaks (anomalous dispersion).\n",
    "    \"\"\"\n",
    "    D1 = Delta - sep/2.0\n",
    "    D2 = Delta + sep/2.0\n",
    "    L1 = 1.0/(D1 + 1j*gamma)\n",
    "    L2 = 1.0/(D2 + 1j*gamma)\n",
    "    sgn = -1.0 if gain else 1.0\n",
    "    return chi0 * sgn * (L1 + L2)\n",
    "\n",
    "\n",
    "def n_from_chi(chi_complex):\n",
    "    # For |chi| << 1, n ≈ 1 + chi/2 (complex n includes absorption/gain via Im part)\n",
    "    return 1.0 + 0.5*chi_complex\n",
    "\n",
    "\n",
    "def group_index_from_n(omega, n_vals, Delta):\n",
    "    \"\"\"\n",
    "    n_g = n + omega * dn/domega. Given Δ = ω - ω0, d/dω = d/dΔ.\n",
    "    \"\"\"\n",
    "    dndDelta = np.gradient(np.real(n_vals), Delta)  # dispersion from real part\n",
    "    return np.real(n_vals) + omega * dndDelta\n",
    "\n",
    "\n",
    "def make_axes(title):\n",
    "    fig = plt.figure(figsize=(6,4))\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_title(title)\n",
    "    ax.grid(True, alpha=0.3)\n",
    "    return fig, ax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "116d3d60",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "### Spectral view & group velocity"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "27814aadd96f4596827824a4610a25a9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Dropdown(description='Model:', options=('EIT', 'Gain-doublet'), value='EIT'), FloatSlider(value…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "**EIT parameters**"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4d0d850fab944ccc8015ee30d680fdb4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatSlider(value=6.0, description='γ31 (MHz):', max=20.0, min=1.0, step=0.5), FloatSlider(valu…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "**Gain-doublet parameters**"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e3140ce431ca4ab7afd92495d6fa08b4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Checkbox(value=True, description='Gain (Im χ < 0)'), FloatSlider(value=30.0, description='Separ…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "aa3ce20f265949e8babefeb9ff7d12f5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='Model:', options=('EIT', 'Gain-doublet'), value='EIT'), FloatSlide…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "def demo(model='EIT', \n",
    "         wavelength_nm=795.0, \n",
    "         chi0=1e-2, \n",
    "         gamma31_MHz=6.0, \n",
    "         gamma21_kHz=10.0, \n",
    "         Oc_MHz=5.0,\n",
    "         span_MHz=50.0,\n",
    "         gain_mode=True,\n",
    "         sep_MHz=30.0,\n",
    "         gamma_MHz=5.0):\n",
    "    # Units\n",
    "    MHz = 1e6 * 2*np.pi  # angular\n",
    "    kHz = 1e3 * 2*np.pi\n",
    "\n",
    "    # Frequency grid (angular)\n",
    "    span = span_MHz * MHz\n",
    "    Delta = np.linspace(-span, span, 1401)  # rad/s\n",
    "\n",
    "    # Carrier frequency from wavelength\n",
    "    lam = wavelength_nm * 1e-9\n",
    "    omega0 = 2*np.pi * c / lam  # rad/s\n",
    "\n",
    "    if model == 'EIT':\n",
    "        gamma31 = gamma31_MHz * MHz\n",
    "        gamma21 = gamma21_kHz * kHz\n",
    "        Oc = Oc_MHz * MHz\n",
    "        chi = chi_eit(Delta, chi0=chi0, gamma31=gamma31, gamma21=gamma21, Oc=Oc)\n",
    "    else:\n",
    "        sep = sep_MHz * MHz\n",
    "        gamma = gamma_MHz * MHz\n",
    "        chi = chi_gain_doublet(Delta, chi0=chi0, gamma=gamma, sep=sep, gain=gain_mode)\n",
    "\n",
    "    chi_r = np.real(chi)\n",
    "    chi_i = np.imag(chi)\n",
    "\n",
    "    # Plot χ' and χ''\n",
    "    fig1, ax1 = make_axes(\"Susceptibility χ' (real) vs Δ\")\n",
    "    ax1.plot(Delta/(2*np.pi*1e6), chi_r)\n",
    "    ax1.set_xlabel(\"Detuning Δ (MHz)\")\n",
    "    ax1.set_ylabel(\"χ' (dimensionless)\")\n",
    "    plt.show(fig1)\n",
    "\n",
    "    fig2, ax2 = make_axes(\"Susceptibility χ'' (imag) vs Δ\")\n",
    "    ax2.plot(Delta/(2*np.pi*1e6), chi_i)\n",
    "    ax2.set_xlabel(\"Detuning Δ (MHz)\")\n",
    "    ax2.set_ylabel(\"χ'' (dimensionless)\")\n",
    "    plt.show(fig2)\n",
    "\n",
    "    # Complex refractive index and group velocity\n",
    "    n_vals = n_from_chi(chi)             # complex n\n",
    "    n_real = np.real(n_vals)\n",
    "    ng = group_index_from_n(omega0, n_vals, Delta)\n",
    "    vg = c / ng  # m/s (may be negative)\n",
    "\n",
    "    fig3, ax3 = make_axes(\"Refractive index n (real) vs Δ\")\n",
    "    ax3.plot(Delta/(2*np.pi*1e6), n_real)\n",
    "    ax3.set_xlabel(\"Detuning Δ (MHz)\")\n",
    "    ax3.set_ylabel(\"n (dimensionless)\")\n",
    "    plt.show(fig3)\n",
    "\n",
    "    fig4, ax4 = make_axes(\"Group velocity v_g vs Δ\")\n",
    "    ax4.plot(Delta/(2*np.pi*1e6), vg)\n",
    "    ax4.set_xlabel(\"Detuning Δ (MHz)\")\n",
    "    ax4.set_ylabel(\"v_g (m/s)\")\n",
    "    plt.show(fig4)\n",
    "\n",
    "    # Report values at Δ = 0\n",
    "    idx0 = np.argmin(np.abs(Delta))\n",
    "    display(Markdown(\n",
    "        f\"\"\"**At two-photon resonance (Δ = 0):**  \n",
    "- χ' = {chi_r[idx0]:.3e}, χ'' = {chi_i[idx0]:.3e}  \n",
    "- n ≈ {n_real[idx0]:.9f}  \n",
    "- n_g ≈ {ng[idx0]:.6e}  \n",
    "- v_g ≈ {vg[idx0]:.6e} m/s\"\"\"\n",
    "    ))\n",
    "\n",
    "\n",
    "def ui():\n",
    "    model = Dropdown(options=['EIT', 'Gain-doublet'], value='EIT', description='Model:')\n",
    "    wavelength = FloatSlider(value=795.0, min=400.0, max=1600.0, step=1.0, description='λ (nm):', readout_format='.0f')\n",
    "    chi0 = FloatSlider(value=1e-2, min=1e-3, max=1e0, step=1e-3, readout_format='.1e', description='χ₀:')\n",
    "    span = FloatSlider(value=50.0, min=5.0, max=200.0, step=1.0, description='Span (MHz):')\n",
    "\n",
    "    # EIT parameters\n",
    "    gamma31 = FloatSlider(value=6.0, min=1.0, max=20.0, step=0.5, description='γ31 (MHz):')\n",
    "    gamma21 = FloatSlider(value=10.0, min=0.1, max=200.0, step=0.1, description='γ21 (kHz):')\n",
    "    Oc = FloatSlider(value=5.0, min=0.1, max=50.0, step=0.1, description='Ωc (MHz):')\n",
    "\n",
    "    # Gain-doublet parameters\n",
    "    gain_mode = Checkbox(value=True, description='Gain (Im χ < 0)')\n",
    "    sep = FloatSlider(value=30.0, min=1.0, max=200.0, step=1.0, description='Separation (MHz):')\n",
    "    ggamma = FloatSlider(value=5.0, min=0.5, max=50.0, step=0.5, description='γ (MHz):')\n",
    "\n",
    "    def _render(model, wavelength_nm, chi0_val, span_MHz, \n",
    "                gamma31_MHz, gamma21_kHz, Oc_MHz, \n",
    "                gain_mode_val, sep_MHz, gamma_MHz):\n",
    "        if model == 'EIT':\n",
    "            demo('EIT', wavelength_nm, chi0_val, gamma31_MHz, gamma21_kHz, Oc_MHz, span_MHz,\n",
    "                 gain_mode_val, sep_MHz, gamma_MHz)\n",
    "        else:\n",
    "            demo('Gain-doublet', wavelength_nm, chi0_val, gamma31_MHz, gamma21_kHz, Oc_MHz, span_MHz,\n",
    "                 gain_mode_val, sep_MHz, gamma_MHz)\n",
    "\n",
    "    display(Markdown(\"### Spectral view & group velocity\"))\n",
    "    display(HBox([model, wavelength, chi0, span], layout=Layout(flex_flow='row wrap')))\n",
    "    box_eit = HBox([gamma31, gamma21, Oc], layout=Layout(flex_flow='row wrap'))\n",
    "    box_gain = HBox([gain_mode, sep, ggamma], layout=Layout(flex_flow='row wrap'))\n",
    "    display(Markdown(\"**EIT parameters**\")); display(box_eit)\n",
    "    display(Markdown(\"**Gain-doublet parameters**\")); display(box_gain)\n",
    "\n",
    "    interact(_render,\n",
    "             model=model,\n",
    "             wavelength_nm=wavelength,\n",
    "             chi0_val=chi0,\n",
    "             span_MHz=span,\n",
    "             gamma31_MHz=gamma31,\n",
    "             gamma21_kHz=gamma21,\n",
    "             Oc_MHz=Oc,\n",
    "             gain_mode_val=gain_mode,\n",
    "             sep_MHz=sep,\n",
    "             gamma_MHz=ggamma)\n",
    "\n",
    "ui()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03da6298",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}