docs/usage/spectral-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports for plotting\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "sns.set_style('darkgrid')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectral Analysis\n",
    "\n",
    "We can study the response of a quantised system to an input field using spectral analysis. For linear systems, analytic response functions are known. We can perform spectral analysis by transforming into the frequency domain.\n",
    "\n",
    "Here we'll look at how to use the `spectral` module to obtain the absorption and dispersion profiles by transforming into the frequency domain using a discrete Fourier transform."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll define a typical two-level problem. We'll need to use a high resolution in the time domain if we want to consider high frequencies above and below resonance. So here, note that we use a lot of time steps. We wouldn't usually want this! It will take a long time to save and result in a large savefile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mb_solve_json = \"\"\"\n",
    "{\n",
    "  \"atom\": {\n",
    "    \"decays\": [\n",
    "      {\n",
    "        \"channels\": [[0, 1]],\n",
    "        \"rate\": 1.0\n",
    "      }\n",
    "    ],\n",
    "    \"fields\": [\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 1]],\n",
    "        \"rabi_freq\": 1.0e-3,\n",
    "        \"rabi_freq_t_args\": {\n",
    "          \"ampl\": 1.0,\n",
    "          \"centre\": 0.0,\n",
    "          \"fwhm\": 1.0\n",
    "        },\n",
    "        \"rabi_freq_t_func\": \"gaussian\"\n",
    "      }\n",
    "    ],\n",
    "    \"num_states\": 2\n",
    "  },\n",
    "  \"t_min\": -2.0,\n",
    "  \"t_max\": 10.0,\n",
    "  \"t_steps\": 1000,\n",
    "  \"z_min\": -0.2,\n",
    "  \"z_max\": 1.2,\n",
    "  \"z_steps\": 50,\n",
    "  \"interaction_strengths\": [\n",
    "    1.0\n",
    "  ],\n",
    "  \"savefile\": \"spectral-analysis\"\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from maxwellbloch import mb_solve\n",
    "mbs = mb_solve.MBSolve().from_json_str(mb_solve_json)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Omegas_zt, states_zt = mbs.mbsolve(recalc=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results in the Time Domain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll look at the results in the time domain as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(1, figsize=(16, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "cmap_range = np.linspace(0.0, 1.0e-3, 11)\n",
    "cf = ax.contourf(mbs.tlist, mbs.zlist, \n",
    "                 np.abs(mbs.Omegas_zt[0]/(2*np.pi)), \n",
    "                 cmap_range, cmap=plt.cm.Blues)\n",
    "ax.set_title('Rabi Frequency ($\\Gamma / 2\\pi $)')\n",
    "ax.set_xlabel('Time ($1/\\Gamma$)')\n",
    "ax.set_ylabel('Distance ($L$)')\n",
    "for y in [0.0, 1.0]:\n",
    "    ax.axhline(y, c='grey', lw=1.0, ls='dotted')\n",
    "plt.colorbar(cf);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, it's the typical absorption picture [we've seen before](two-level.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results in the Frequency Domain\n",
    "\n",
    "The envelope of a short Gaussian pulse on resonance is equivalent to a superposition of waves across a wide range of frequencies. We can analyse the spectral response of the medium to the field by transforming the simulated Rabi frequency in the time domain, $\\Omega(z, t)$, to the frequency domain, $\\Omega(z, \\omega)$.\n",
    "\n",
    "We determine the **absoprtion** profile by measuring the intensity of the pulse that arrives at the back of the medium relative to the intensity put in at the front. This we can relate to the imaginary part of the optical susceptibility $\\chi_I$ via\n",
    "$$\n",
    "\\frac{k}{2} \\chi_I(\\omega) z = -\\log \\frac{| \\Omega (z, \\omega) |}{| \\Omega(z = 0, \\omega) |}\n",
    "$$\n",
    "where $k$ is the wavenumber.\n",
    "\n",
    "We determine the **dispersion** profile by measuring the phase shift of the pulse that arrives at the back of the medium relative to the pulse that arrived at the front. This we can relate to the real part of optical susceptibility  $\\chi_R$ via\n",
    "$$\n",
    "\\frac{k}{2} \\chi_R(\\omega) z = \\phi(z, \\omega) - \\phi(z=0, \\omega)\n",
    "$$\n",
    "where $\\phi$ is the complex phase defined via\n",
    "$$\n",
    "\\Omega(z, \\omega) = |\\Omega|\\mathrm{e}^{\\mathrm{i}\\phi}.\n",
    "$$\n",
    "\n",
    "See section 2.5 of [my thesis][thesis] for an explanation of how the susceptibility is derived.\n",
    "\n",
    "Methods to provide the dispersion and absorption of a given field are provided in the `spectral` module.\n",
    "\n",
    "[thesis]: http://etheses.dur.ac.uk/11599/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from maxwellbloch import spectral, utility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pal = sns.color_palette('deep')\n",
    "\n",
    "fig = plt.figure(2, figsize=(16, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "freq_list = spectral.freq_list(mbs)\n",
    "ax.plot(freq_list, spectral.absorption(mbs, field_idx=0), \n",
    "        label='Absorption', c=pal[0], lw=3.0)\n",
    "ax.plot(freq_list, spectral.dispersion(mbs, field_idx=0), \n",
    "        label='Dispersion', c=pal[1], lw=3.0)\n",
    "rabi_freq_abs_0 = np.abs(spectral.rabi_freq(mbs, 0))[0]\n",
    "ax.plot(freq_list, \n",
    "        rabi_freq_abs_0/np.max(rabi_freq_abs_0),\n",
    "        label='Pulse, Normalised', ls='dashed', c=pal[7])\n",
    "\n",
    "ax.set_xlim(-3.0, 3.0)\n",
    "ax.set_ylim(-1.0, 1.0)\n",
    "ax.set_xlabel('Frequency ($\\Gamma$)')\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The input pulse profile (normalised) is shown in grey. In blue we see the familiar Lorentzian absorption profile of homogeneous broadening due to the spontaneous decay. The full-width at half maximum (FWHM) of the peak is the decay rate, $1 \\Gamma$. The associated dispersion profile is shown in orange, describing the phase of the input field relative to the frequency of the oscillator as it passes over resonance from lagging to leading."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing with Analytic Results\n",
    "\n",
    "For a weak input field like this one, the field is in the regime of linear optics and it is possible to use the weak probe approximation to derive an analytic lineshape\n",
    "$$\n",
    "\\frac{k}{2} \\chi({\\omega}) = -Ng \\frac{1}{\\mathrm{i}\\frac{\\Gamma}{2} + \\omega}.\n",
    "$$\n",
    "For details on how this is derived, see equation (2.61) of my thesis.\n",
    "\n",
    "These analytic absorption and dispersion profiles for weak fields are also included in the `spectral` module, so you can check if you are really in the linear regime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interaction_strength = mbs.interaction_strengths[0]\n",
    "decay_rate = mbs.atom.decays[0]['rate']\n",
    "absorption_linear_known = spectral.absorption_two_linear_known(freq_list, \n",
    "    interaction_strength, decay_rate)\n",
    "dispersion_linear_known = spectral.dispersion_two_linear_known(freq_list, \n",
    "    interaction_strength, decay_rate)\n",
    "\n",
    "fig = plt.figure(4, figsize=(16, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "ax.plot(freq_list, spectral.absorption(mbs, 0, -1), \n",
    "        label='Absorption', lw=5.0, c=pal[0])\n",
    "ax.plot(freq_list, spectral.dispersion(mbs, 0, -1), \n",
    "        label='Dispersion', lw=5.0, c=pal[1])\n",
    "\n",
    "ax.plot(freq_list, absorption_linear_known, ls='dashed', c='white', lw=2.0)\n",
    "ax.plot(freq_list, dispersion_linear_known, ls='dashed', c='white', lw=2.0)\n",
    "\n",
    "# Widths\n",
    "hm, r1, r2 = utility.half_max_roots(freq_list, spectral.absorption(mbs, field_idx=0))\n",
    "plt.hlines(y=hm, xmin=r1, xmax=r2, linestyle='dotted', color=pal[0])\n",
    "plt.annotate('FWHM: ' + '%0.2f'%(r2 - r1), xy=(r2, hm), color=pal[0],\n",
    "              xycoords='data', xytext=(5, 5), textcoords='offset points');\n",
    "\n",
    "ax.set_xlim(-3.0, 3.0)\n",
    "ax.set_ylim(-1.0, 1.0)\n",
    "ax.set_xlabel('Frequency ($\\Gamma$)')\n",
    "\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot residuals\n",
    "fig = plt.figure(figsize=(16, 2))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(freq_list, spectral.absorption(mbs, 0, -1) - absorption_linear_known, \n",
    "        label='Absorption', lw=2.0, c=pal[0])\n",
    "ax.set_xlim(-3.0, 3.0)\n",
    "ax.set_ylim(-3e-2, 3e-2)\n",
    "ax.set_xlabel('Frequency ($\\Gamma$)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The known analytic lineshapes are dashed white over the simulated lineshapes. For absorption, we see a maximum absolute residual of 3% and a simualated full width at half maximum (FWHM) of $1 \\Gamma$, which is the decay rate. This tells us we are in the linear regime with negligible additional broadening. At higher field intensities we will see power broadening. Inclusion of effects like inhomogeneous broadening due to the thermal motion and collisions will also lead to broadened lineshapes."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}