{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Velocity Classes for Modelling Doppler Broadening in Thermal Systems\n", "\n", "So far the models we have considered are appropriate for ultracold systems with close-to-stationary atoms. For thermal atoms we might want to consider the averaging effect of motion in the line of propagation.\n", "\n", "An atom moving with a velocity component $v$ in the $z$-direction will interact with a Doppler-shifted field frequency $\\omega - kv$. This shift is effected over a 1D Maxwell-Boltzmann probability distribution function of velocity\n", "\n", "$$\n", "f(v) = \\frac{1}{u\\sqrt{\\pi}} \\mathrm{e}^{-(kv/u)^2}\n", "$$\n", "\n", "where the thermal width $u = kv_w$, $k$ is the wavenumber of the monochromatic field and $v_w = 2k_b T/ M$ is the most probable speed of the Maxwell-Boltzmann distribution for a temperature $T$ and atomic mass $m$.\n", "\n", "In MaxwellBloch we model this by solving the system over a range of `velocity_classes`, each detuning the system from resonance. The results are convoluted over the Maxwell-Boltzmann distribution. \n", "The thermal width is provided in `thermal_width` in the same $2\\pi~\\Gamma$ units as the `decay`s and `rabi_freq`s.\n", "We provide velocity classes from `thermal_delta_min` to `thermal_delta_max`, again in units of $2\\pi~\\Gamma$. The number of classes we choose to solve is given in `thermal_delta_steps`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Imports for plotting\n", "import matplotlib.pyplot as plt\n", "\n", "%matplotlib inline\n", "import numpy as np\n", "import seaborn as sns\n", "\n", "sns.set_style(\"darkgrid\")" ] }, { "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", " \"velocity_classes\": {\n", " \"thermal_delta_min\": -0.1,\n", " \"thermal_delta_max\": 0.1,\n", " \"thermal_delta_steps\": 4,\n", " \"thermal_width\": 0.05\n", " },\n", " \"savefile\": \"velocity-classes\"\n", "}\n", "\"\"\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from maxwellbloch import mb_solve\n", "\n", "mbs = mb_solve.MBSolve().from_json_str(mb_solve_json)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check the set of velocity classes we've defined:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mbs.thermal_delta_list / (2 * np.pi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The weights of the Maxwell-Boltzmann distribution at these deltas is given by:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mbs.thermal_weights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And so we can plot the numerical approximation to the Gaussian Maxwell-Boltzmann distribution:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "maxboltz = mb_solve.maxwell_boltzmann(\n", " mbs.thermal_delta_list, 2 * np.pi * mbs.velocity_classes[\"thermal_width\"]\n", ")\n", "plt.plot(mbs.thermal_delta_list, maxboltz, marker=\"o\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is useful to look at what the numerical integration looks like for these velocity classes. If it is close to 1, the thermal distribution should be well covered." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.trapezoid(mbs.thermal_weights, mbs.thermal_delta_list)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can solve as before. Now at each $z$-step, the system will be solved `thermal_delta_steps` times, once for each velocity class, and so the time taken to solve scales linearly." ] }, { "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": "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(\n", " mbs.tlist,\n", " mbs.zlist,\n", " np.abs(mbs.Omegas_zt[0] / (2 * np.pi)),\n", " cmap_range,\n", " cmap=plt.cm.Blues,\n", ")\n", "ax.set_title(r\"Rabi Frequency ($\\Gamma / 2\\pi $)\")\n", "ax.set_xlabel(r\"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": [ "## Results in the Frequency Domain\n", "\n", "For a two-level system in the linear regime, this convolution of a Lorentzian function with a Gaussian can be determined analytically and is known as a Voigt proļ¬le. It is provided in MaxwellBloch as `spectral.voigt_two_linear_known(freq_list, decay_rate, thermal_width)` so we can compare the simulation solution with the known analytic result." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from maxwellbloch import spectral, utility" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interaction_strength = mbs.interaction_strengths[0]\n", "decay_rate = mbs.atom.decays[0][\"rate\"]\n", "freq_list = spectral.freq_list(mbs)\n", "absorption_linear_known = spectral.absorption_two_linear_known(\n", " freq_list, interaction_strength, decay_rate\n", ")\n", "dispersion_linear_known = spectral.dispersion_two_linear_known(\n", " freq_list, interaction_strength, decay_rate\n", ")\n", "\n", "fig = plt.figure(4, figsize=(16, 6))\n", "ax = fig.add_subplot(111)\n", "pal = sns.color_palette(\"deep\")\n", "\n", "ax.plot(\n", " freq_list, spectral.absorption(mbs, 0, -1), label=\"Absorption\", lw=5.0, c=pal[0]\n", ")\n", "ax.plot(\n", " freq_list, spectral.dispersion(mbs, 0, -1), label=\"Dispersion\", lw=5.0, c=pal[1]\n", ")\n", "\n", "ax.plot(\n", " freq_list,\n", " absorption_linear_known,\n", " ls=\"dotted\",\n", " c=pal[0],\n", " lw=2.0,\n", " label=\"Absorption, No Thermal\",\n", ")\n", "ax.plot(\n", " freq_list,\n", " dispersion_linear_known,\n", " ls=\"dotted\",\n", " c=pal[1],\n", " lw=2.0,\n", " label=\"Dispersion, No Thermal\",\n", ")\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(\n", " \"FWHM: \" + \"%0.2f\" % (r2 - r1),\n", " xy=(r2, hm),\n", " color=pal[0],\n", " xycoords=\"data\",\n", " xytext=(5, 5),\n", " textcoords=\"offset points\",\n", ")\n", "\n", "voigt = spectral.voigt_two_linear_known(freq_list, 1.0, 0.05).imag\n", "ax.plot(\n", " freq_list,\n", " voigt,\n", " c=\"white\",\n", " ls=\"dashed\",\n", " lw=2.0,\n", " label=\"Known Absorption, Voigt Profile\",\n", ")\n", "\n", "ax.set_xlim(-3.0, 3.0)\n", "ax.set_ylim(-1.0, 1.0)\n", "ax.set_xlabel(r\"Frequency ($\\Gamma$)\")\n", "\n", "ax.legend();" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Plot residuals\n", "fig = plt.figure(figsize=(16, 2))\n", "ax = fig.add_subplot(111)\n", "ax.plot(\n", " freq_list,\n", " spectral.absorption(mbs, 0, -1) - voigt,\n", " label=\"Absorption\",\n", " lw=2.0,\n", " c=pal[0],\n", ")\n", "ax.set_xlim(-3.0, 3.0)\n", "ax.set_ylim(-3e-2, 3e-2)\n", "ax.set_xlabel(r\"Frequency ($\\Gamma$)\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding Inner Steps\n", "\n", "If the thermal width is much larger than the decay rate, you may wish to add a narrower band around the resonance frequency to sample the Lorentzian sufficiently while covering the Maxwell-Boltzmann distribution efficiently. This can be done by providing an inner range defined by `thermal_delta_inner_min`, `thermal_delta_inner_min` and `thermal_delta_inner_steps`. Any duplicated velocity classes (in the constructed `thermal_delta_list`) will only be counted once. Below is an example." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vc = {\n", " \"thermal_delta_min\": -0.1,\n", " \"thermal_delta_max\": 0.1,\n", " \"thermal_delta_steps\": 4,\n", " \"thermal_delta_inner_min\": -0.05,\n", " \"thermal_delta_inner_max\": 0.05,\n", " \"thermal_delta_inner_steps\": 10,\n", " \"thermal_width\": 0.05,\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the thermal delta range is $[-0.1, -0.05, 0.05, 1.0]$ and the inner range is $[-0.05, -0.04, \\dots, 0.04, 0.05]$. These are combined to form:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mbs.build_velocity_classes(velocity_classes=vc)\n", "print(mbs.thermal_delta_list / (2 * np.pi))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "maxboltz = mb_solve.maxwell_boltzmann(\n", " mbs.thermal_delta_list, 2 * np.pi * mbs.velocity_classes[\"thermal_width\"]\n", ")\n", "plt.plot(mbs.thermal_delta_list, maxboltz, marker=\"o\");" ] } ], "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 }