docs/usage/two-level.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports for plotting\n",
    "\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": [
    "# Getting Started with the Two-Level Atom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first model we'll consider is a medium of two-level atoms, coupled by a weak pulse.\n",
    "\n",
    "```\n",
    "____ |1>\n",
    " |\n",
    " | Ω\n",
    " |\n",
    "____ |0>\n",
    "```\n",
    "\n",
    "We'll need to define a field, and the two-level atom for the field to address."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the Field and Atom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from maxwellbloch import field, ob_atom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll define the field. The quantity we care about is the input Rabi frequency as a function of time, $\\Omega(t)$. Let's say our input pulse has a max $\\Omega_0 = 2\\pi \\cdot 0.001 \\mathrm{~ MHz}$, and is a Gaussian with a full-width at half-maximum (FWHM) of $1 \\mathrm{~ \\mu s}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "the_field = field.Field(\n",
    "    coupled_levels=[[0, 1]],\n",
    "    rabi_freq=1.0e-3,               # [2π MHz]\n",
    "    rabi_freq_t_func='gaussian',\n",
    "    rabi_freq_t_args={\n",
    "        'ampl': 1.0,                \n",
    "        'centre': 0.0,              # [μs]\n",
    "        'fwhm': 1.0},               # [μs]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first property `coupled_levels` is a list of pairs to be coupled by the field. In this case we only have two levels, indexed `0` (the ground state) and `1` (the excited state). The properties `rabi_freq`, `rabi_freq_t_func` and `rabi_freq_t_args` together define the field profile entering the medium. What does the profile look like in this case?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from maxwellbloch import t_funcs\n",
    "\n",
    "tlist = np.linspace(-2, 10, 100)\n",
    "plt.plot(tlist, t_funcs.gaussian(1)(tlist, args={ 'ampl_1': 1.0, 'centre_1': 0.0, 'fwhm_1': 1.0}));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's a Gaussian-shaped pulse, centred at $t=0$. Many more pre-defined profiles are defined in the `t_funcs` module, see [REF t_funcs doc] for details. \n",
    "\n",
    "So that's the field. Now we need to define an atom to address. Here we just need to state that `num_states=2`. The atom also takes a list of fields, in this case we've only got one, `the_field`, so we'll add that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "the_atom = ob_atom.OBAtom(\n",
    "    num_states=2,\n",
    "    )\n",
    "the_atom.fields.append(the_field)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the Maxwell-Bloch Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have everything required to solve the problem. We just need to define the space $z$ and time $t$ dimesions we want to solve the problem over. Let's say that our medium is $1 \\mathrm{~ cm}$ long and has constant number density, for example a vapour cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from maxwellbloch import mb_solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mbs = mb_solve.MBSolve(\n",
    "    t_min=-2.0,                     # [μs]\n",
    "    t_max=10.0,                     # [μs]\n",
    "    t_steps=100, \n",
    "    z_min=-0.5,                     # [cm]\n",
    "    z_max=1.5,                      # [cm]\n",
    "    z_steps=10,\n",
    "    interaction_strengths=[0.1]     # [2π MHz /cm] \n",
    ")\n",
    "mbs.atom = the_atom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we've set `z_min=-0.5` and `z_max=1.5` so that we can see how the pulse travels before and after it is in the medium, i.e. as it travels in vacuum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interaction Strength\n",
    "\n",
    "The interaction strength ($Ng$) for each field describes the strength of the coupling of the field to the medium. You can think of it as the amount of absorption per unit length of medium, so in this system it has units of $\\mathrm{2\\pi~MHz~/cm}$.\n",
    "\n",
    "Now we can solve the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Omegas_zt, states_zt = mbs.mbsolve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solver indicates progress as it goes, in 10% increments by default."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting and Analysing Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solver returns two objects. `Omegas_zt` contains the complex value of each field, at each $z$-point and each $t$-point. In this case we have only one field, $50+1$ $z$-points (for 50 `z_steps`) and $100+1$ $t$-points (for 100 `t_steps`). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Omegas_zt.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`states_zt` contains the density matrix representation of the states of the atoms at each $z$-point and each $t$-point. The density matrix $\\rho$ is of size $n \\times n$ where $n$ is the number of states in the atom, in this case $n=2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states_zt.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both the Rabi frequency $\\Omega$ and the density matrix $\\rho$ are complex-valued, so the elements of both `Omegas_zt` and `states_zt`\n",
    "are of `np.complex` datatype."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 5))\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 ($\\mathrm{MHz} / 2\\pi $), Speed-of-Light Reference Frame')\n",
    "ax.set_xlabel('Time ($\\mathrm{\\mu s}$)')\n",
    "ax.set_ylabel('Distance ($\\mathrm{cm}$)')\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": [
    "The plot above shows the absolute value of `Omegas_zt[0]`, the simulated result of propagation of the weak pulse through the two-level medium. Time $t$ is shown on the $x$-axis and the propagation distance $z$ is shown on the $y$-axis such that the field enters at the bottom of the plot. The\n",
    "horizontal slice at $z = 0$ thus represents the input field pulse. \n",
    "\n",
    "This results is presented in the speed-of-light reference frame. That means that a pulse travelling at the speed of light $c$ moves vertically in a plot like the one above. See [REF reference frame] for how we shift into the lab frame."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis\n",
    "\n",
    "What happens in this simulation? We can see that the pulse is attenuated by the medium, in that the maximum of the peak leaving the medium at $z=1~\\mathrm{cm}$ is lower than the maximum peak entering at $z=0~\\mathrm{cm}$. \n",
    "\n",
    "We also see that the pulse is slightly fast, such that the peak arrives at the back of the medium, $z=1~\\mathrm{cm}$, at\n",
    "a time $t < 0~\\mathrm{\\mu s}$ in this speed-of-light reference frame. This is typical of propagation through a medium with a normal dispersion profile. [REF spectral analysis]. There's a little bit of ringing after (to the right) of the pulse, seen in the lighter region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the System with JSON\n",
    "\n",
    "We went through defining `Field`, `Atom` and `MBSolve` objects separately so I could introduce these concepts, but in practice it is rare that we will need to interact with these directly. The whole problem can be defined via a JSON file or string. Here is the same system above defined in JSON."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mbs_json = \"\"\"\n",
    "{\n",
    "  \"atom\": {\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\": 100,\n",
    "  \"z_min\": -0.2,\n",
    "  \"z_max\": 1.2,\n",
    "  \"z_steps\": 10,\n",
    "  \"interaction_strengths\": [\n",
    "    0.1\n",
    "  ]\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can load these using the `from_json(fname)` method for a JSON file or the `from_json_str(json_str)` method for a string like the one above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mbs = mb_solve.MBSolve().from_json_str(mbs_json)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can solve the problem exactly as we did above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Omegas_zt, states_zt = mbs.mbsolve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From now on we'll nearly always define MaxwellBloch problems like this with JSON. And just to check we get the same result as above, we'll plot the simulated field again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 5))\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 ($\\mathrm{MHz} / 2\\pi $), Speed-of-Light Reference Frame')\n",
    "ax.set_xlabel('Time ($\\mathrm{\\mu s}$)')\n",
    "ax.set_ylabel('Distance ($\\mathrm{cm}$)')\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": [
    "### Adding More Atoms\n",
    "\n",
    "What happens if we want to increase the number density $N$ of the medium? We'll define the same problem but increase the `interaction_strength` by a factor of 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mbs_json = \"\"\"\n",
    "{\n",
    "  \"atom\": {\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\": 100,\n",
    "  \"z_min\": -0.5,\n",
    "  \"z_max\": 1.5,\n",
    "  \"z_steps\": 20,\n",
    "  \"interaction_strengths\": [\n",
    "    1.0\n",
    "  ]\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mbs = mb_solve.MBSolve().from_json_str(mbs_json)\n",
    "Omegas_zt, states_zt = mbs.mbsolve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And again we'll plot the absolute value of the Rabi frequency, $| \\Omega(z,t) |$, over the medium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 5))\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 ($\\mathrm{MHz} / 2\\pi $), Speed-of-Light Reference Frame')\n",
    "ax.set_xlabel('Time ($\\mathrm{\\mu s}$)')\n",
    "ax.set_ylabel('Distance ($\\mathrm{cm}$)')\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": [
    "## Spontaneous Decay\n",
    " \n",
    "To include spontaneous decay in the Lindblad master equation we have to add collapse operators. We add a list of `decays` into the `atom`, each with a list of channels describing the decays (ordered `lower` then `upper`) and a decay rate in the same units as the rabi frequency (in this case $2\\pi \\mathrm{~ MHz}$). Lets say in our system the decay rate is $2\\pi \\cdot 1 \\mathrm{~ MHz}$.\n",
    "\n",
    "Otherwise we'll leave the system the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mbs_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\": 100,\n",
    "  \"z_min\": -0.5,\n",
    "  \"z_max\": 1.5,\n",
    "  \"z_steps\": 20,\n",
    "  \"interaction_strengths\": [\n",
    "    1.0\n",
    "  ]\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mbs = mb_solve.MBSolve().from_json_str(mbs_json)\n",
    "Omegas_zt, states_zt = mbs.mbsolve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 5))\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 ($\\mathrm{MHz} / 2\\pi $), Speed-of-Light Reference Frame')\n",
    "ax.set_xlabel('Time ($\\mathrm{\\mu s}$)')\n",
    "ax.set_ylabel('Distance ($\\mathrm{cm}$)')\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": [
    "We see that the ringing in the system has been damped. In other solvers, decay is neglected and the Schrödinger equation is solved rather than the Lindblad master equation. In real atomic systems, unless the times we're looking at are very short compared to the decay lifetime, inhomogeneous broadening by spontaneous decay has considerable effect and is important to include in a model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using a Natural Unit System"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to use any reciprocal angular momentum and time units (e.g. $2\\pi \\mathrm{~ MHz}$ and $\\mathrm{~ \\mu s}$ in the above example). The length unit is also arbitrary (e.g. $\\mathrm{cm}$) as long as the interaction strength $Ng$ is given in the same units (e.g. $\\mathrm{2\\pi~MHz~/cm}$).\n",
    "\n",
    "For a two-level system we have a single natural linewidth, and so it is convenient to introduce a natural unit system, with frequencies in units of the natural linewidth $\\Gamma$, times in units of the reciprocal spontaneous decay lifetime $\\tau = 1/\\Gamma$ and distances in units of the length of the medium $L$.\n",
    "\n",
    "By introducing this natural unit system we are able to reduce the number of parameters involved in the mathematical problem. For example, it becomes clear that increasing the length of the medium ten times is equivalent to raising the number density by the same scale, or by choosing a system with a suitably higher dipole moment.\n",
    "\n",
    "For most examples in this documentation from now on, we will use this natural unit system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shifting to the Fixed Frame of Reference\n",
    "\n",
    "We solve the problem in a frame of reference that moves with the speed of light across the medium. This means we can keep time in the system separate from distance, which makes the coupled equations easier to solve. We can also see all the important details of the solution in this frame of reference, as in the plots above. But if we want to look at how a field actually propagates in time, we need to shift back to a fixed (or laboratory) frame of reference. To do that we need to connect the time and space dimensions via the speed of light in the system. The utilities for this shift are in the `fixed` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from maxwellbloch import fixed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the speed of light based on the units in the system. For example, in the units above, $c \\approx 3 \\times 10^4 ~ \\mathrm{cm / \\mu s }$. If I put that speed in to the simulation results above, you will see no difference in the output. It will be too fast, and appear that the pulse progresses through the medium instantanously. Instead we can put in a slow speed-of-light to get an idea of what the shift does, say $10^{-5} c$. Something like looking at the pulse in slow motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "speed_of_light = 0.3 # [cm /μs] THIS IS 10^5 SLOWER THAN C!\n",
    "\n",
    "tlist_fixed_frame = fixed.t_list(mbs, speed_of_light)\n",
    "field_fixed_frame = fixed.rabi_freq_abs(mbs, 0, speed_of_light)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 5))\n",
    "cmap_range = np.linspace(0.0, 1.0e-3, 11)\n",
    "cf = ax.contourf(tlist_fixed_frame, mbs.zlist, \n",
    "                 np.abs(field_fixed_frame/(2*np.pi)), \n",
    "                 cmap_range, cmap=plt.cm.Blues)\n",
    "ax.set_title('Rabi Frequency ($\\mathrm{MHz} / 2\\pi $), Fixed Reference Frame')\n",
    "ax.set_xlabel('Fixed Time ($\\mathrm{\\mu s}$)')\n",
    "ax.set_ylabel('Distance ($\\mathrm{cm}$)')\n",
    "for x in [0.0, 1/speed_of_light]:\n",
    "    ax.axvline(x, c='red', lw=1.0, ls='dashed')\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": [
    "So we see that in this simulation with a slow speed-of-light, the pulse would arrive at the back of the medium around $3.33 \\mathrm{~ \\mu s}$ after it hit the front of the medium. The shift to the fixed reference frame has the effect of skewing the results rightward as we move through the medium, representing progression through time as well as space."
   ]
  }
 ],
 "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
}