{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "intro",
   "metadata": {},
   "source": [
    "# V-Type Three-Level: Simultons\n",
    "\n",
    "In the **V-type three-level system** a probe field couples the ground state\n",
    "|0⟩ to |1⟩ and a coupling field couples the same ground state to |2⟩.\n",
    "Because both fields share the ground state, they interact through the common\n",
    "atomic coherence.\n",
    "\n",
    "When the combined pulse area satisfies the V-system area theorem,\n",
    "\n",
    "$$\n",
    "\\theta_\\mathrm{total} = \\sqrt{\\theta_\\mathrm{probe}^2 + \\theta_\\mathrm{coupling}^2} = 2\\pi n,\n",
    "$$\n",
    "\n",
    "the two fields lock into shape-preserving **simulton** solutions that propagate\n",
    "without absorption — the vector generalisation of the SIT soliton.\n",
    "\n",
    "This notebook surveys the main regimes:\n",
    "\n",
    "1. **Weak simulton** (0.5π + 0.5π): total area below threshold — both fields absorbed.\n",
    "2. **Simulton propagation** (0.5π + 1.5π): sub-threshold redistribution between fields.\n",
    "3. **√2π simulton** (√2π + √2π): threshold simulton — minimum area for lossless propagation.\n",
    "4. **√8π soliton formation** (√8π + √8π): two-soliton simulton.\n",
    "5. **Soliton collision**: two time-displaced solitons merge into a simulton.\n",
    "6. **Optical surfer**: weak CW probe dragged forward by a 2π coupling soliton.\n",
    "7. **Double optical surfer**: weak CW probe driven by a 4π coupling soliton."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "imports",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "\n",
    "from maxwellbloch import mb_solve\n",
    "\n",
    "sns.set_style(\"darkgrid\")\n",
    "\n",
    "sech_fwhm_conv = 1.0 / 2.6339157938\n",
    "t_width = 1.0 * sech_fwhm_conv  # [τ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sec-0505-intro",
   "metadata": {},
   "source": [
    "## Weak simulton — 0.5π + 0.5π\n",
    "\n",
    "With both probe and coupling carrying area 0.5π each, the combined area\n",
    "$\\theta_\\mathrm{total} = \\sqrt{(0.5\\pi)^2 + (0.5\\pi)^2} \\approx 0.71\\pi < 2\\pi$.\n",
    "This is below the simulton threshold so the medium absorbs both fields —\n",
    "neither propagates as a soliton and the total area decreases monotonically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "setup-0505",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 0.5\n",
    "ampl_0505 = n / t_width / (2 * np.pi)\n",
    "\n",
    "mb_solve_json_0505 = \"\"\"\n",
    "{\n",
    "  \"atom\": {\n",
    "    \"fields\": [\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 1]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"probe\",\n",
    "        \"rabi_freq\": 0.20960035913554168,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      },\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 2]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"coupling\",\n",
    "        \"rabi_freq\": 0.20960035913554168,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      }\n",
    "    ],\n",
    "    \"num_states\": 3\n",
    "  },\n",
    "  \"t_min\": -2.0,\n",
    "  \"t_max\": 10.0,\n",
    "  \"t_steps\": 60,\n",
    "  \"z_min\": -0.2,\n",
    "  \"z_max\": 1.2,\n",
    "  \"z_steps\": 70,\n",
    "  \"z_steps_inner\": 2,\n",
    "  \"interaction_strengths\": [10.0, 10.0],\n",
    "  \"savefile\": \"mbs-vee-sech-0.5pi-0.5pi\"\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mbs_0505 = mb_solve.MBSolve().from_json_str(mb_solve_json_0505)\n",
    "mbs_0505.mbsolve(recalc=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-0505-spacetime",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 12))\n",
    "cmap_range = np.linspace(0.0, 0.8, 11)\n",
    "\n",
    "ax = fig.add_subplot(211)\n",
    "cf = ax.contourf(mbs_0505.tlist, mbs_0505.zlist,\n",
    "                 np.abs(mbs_0505.Omegas_zt[0] / (2 * np.pi)), cmap_range, cmap=plt.cm.Blues)\n",
    "ax.set_title(r\"Rabi Frequency ($\\Gamma / 2\\pi$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Probe\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "ax = fig.add_subplot(212)\n",
    "cf = ax.contourf(mbs_0505.tlist, mbs_0505.zlist,\n",
    "                 np.abs(mbs_0505.Omegas_zt[1] / (2 * np.pi)), cmap_range, cmap=plt.cm.Greens)\n",
    "ax.set_xlabel(r\"Time ($1/\\Gamma$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Coupling\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    for y in [0.0, 1.0]:\n",
    "        ax.axhline(y, c=\"grey\", lw=1.0, ls=\"dotted\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-0505-area",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_area = np.sqrt(mbs_0505.fields_area()[0]**2 + mbs_0505.fields_area()[1]**2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "ax.plot(mbs_0505.zlist, mbs_0505.fields_area()[0] / np.pi, label=\"Probe\", clip_on=False)\n",
    "ax.plot(mbs_0505.zlist, mbs_0505.fields_area()[1] / np.pi, label=\"Coupling\", clip_on=False)\n",
    "ax.plot(mbs_0505.zlist, total_area / np.pi, label=\"Total\", ls=\"dashed\", clip_on=False)\n",
    "ax.legend()\n",
    "ax.set_ylim([0.0, 2.0])\n",
    "ax.set_xlabel(\"Distance ($L$)\")\n",
    "ax.set_ylabel(r\"Pulse Area ($\\pi$)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sec-0515-intro",
   "metadata": {},
   "source": [
    "## Simulton propagation — 0.5π probe, 1.5π coupling\n",
    "\n",
    "With probe at 0.5π and coupling at 1.5π the combined area\n",
    "$\\theta_\\mathrm{total} = \\sqrt{(0.5\\pi)^2 + (1.5\\pi)^2} \\approx 1.58\\pi$.\n",
    "Still sub-threshold, but the fields are coupled: the medium redistributes\n",
    "energy from the stronger coupling into the probe as they propagate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "setup-0515",
   "metadata": {},
   "outputs": [],
   "source": [
    "ampl_0515_probe = 0.5 / t_width / (2 * np.pi)\n",
    "ampl_0515_coupling = 1.5 / t_width / (2 * np.pi)\n",
    "\n",
    "mb_solve_json_0515 = \"\"\"\n",
    "{\n",
    "  \"atom\": {\n",
    "    \"fields\": [\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 1]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"probe\",\n",
    "        \"rabi_freq\": 0.20960035913554168,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      },\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 2]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"coupling\",\n",
    "        \"rabi_freq\": 0.628801077406625,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      }\n",
    "    ],\n",
    "    \"num_states\": 3\n",
    "  },\n",
    "  \"t_min\": -2.0,\n",
    "  \"t_max\": 10.0,\n",
    "  \"t_steps\": 60,\n",
    "  \"z_min\": -0.2,\n",
    "  \"z_max\": 1.2,\n",
    "  \"z_steps\": 70,\n",
    "  \"z_steps_inner\": 2,\n",
    "  \"interaction_strengths\": [10.0, 10.0],\n",
    "  \"savefile\": \"mb-solve-vee-sech-0.5pi-1.5pi\"\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mbs_0515 = mb_solve.MBSolve().from_json_str(mb_solve_json_0515)\n",
    "mbs_0515.mbsolve(recalc=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-0515-spacetime",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 12))\n",
    "cmap_range = np.linspace(0.0, 0.8, 11)\n",
    "\n",
    "ax = fig.add_subplot(211)\n",
    "cf = ax.contourf(mbs_0515.tlist, mbs_0515.zlist,\n",
    "                 np.abs(mbs_0515.Omegas_zt[0] / (2 * np.pi)), cmap_range, cmap=plt.cm.Blues)\n",
    "ax.set_title(r\"Rabi Frequency ($\\Gamma / 2\\pi$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Probe\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "ax = fig.add_subplot(212)\n",
    "cf = ax.contourf(mbs_0515.tlist, mbs_0515.zlist,\n",
    "                 np.abs(mbs_0515.Omegas_zt[1] / (2 * np.pi)), cmap_range, cmap=plt.cm.Greens)\n",
    "ax.set_xlabel(r\"Time ($1/\\Gamma$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Coupling\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    for y in [0.0, 1.0]:\n",
    "        ax.axhline(y, c=\"grey\", lw=1.0, ls=\"dotted\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-0515-area",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_area = np.sqrt(mbs_0515.fields_area()[0]**2 + mbs_0515.fields_area()[1]**2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "ax.plot(mbs_0515.zlist, mbs_0515.fields_area()[0] / np.pi, label=\"Probe\", clip_on=False)\n",
    "ax.plot(mbs_0515.zlist, mbs_0515.fields_area()[1] / np.pi, label=\"Coupling\", clip_on=False)\n",
    "ax.plot(mbs_0515.zlist, total_area / np.pi, label=\"Total\", ls=\"dashed\", clip_on=False)\n",
    "ax.legend()\n",
    "ax.set_ylim([0.0, 2.0])\n",
    "ax.set_xlabel(\"Distance ($L$)\")\n",
    "ax.set_ylabel(r\"Pulse Area ($\\pi$)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sec-s2-intro",
   "metadata": {},
   "source": [
    "## √2π simulton formation\n",
    "\n",
    "At $\\theta_\\mathrm{probe} = \\theta_\\mathrm{coupling} = \\sqrt{2}\\pi$ the combined area\n",
    "is exactly $\\theta_\\mathrm{total} = \\sqrt{2\\pi^2 + 2\\pi^2} = 2\\pi$ — the threshold\n",
    "simulton. Both fields propagate without loss, locked together in shape and velocity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "setup-s2",
   "metadata": {},
   "outputs": [],
   "source": [
    "ampl_s2 = np.sqrt(2) / t_width / (2 * np.pi)\n",
    "\n",
    "mb_solve_json_s2 = \"\"\"\n",
    "{\n",
    "  \"atom\": {\n",
    "    \"fields\": [\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 1]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"probe\",\n",
    "        \"rabi_freq\": 0.592839341136,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      },\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 2]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"coupling\",\n",
    "        \"rabi_freq\": 0.592839341136,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      }\n",
    "    ],\n",
    "    \"num_states\": 3\n",
    "  },\n",
    "  \"t_min\": -2.0,\n",
    "  \"t_max\": 10.0,\n",
    "  \"t_steps\": 60,\n",
    "  \"z_min\": -0.2,\n",
    "  \"z_max\": 1.2,\n",
    "  \"z_steps\": 70,\n",
    "  \"z_steps_inner\": 1,\n",
    "  \"num_density_z_func\": \"square\",\n",
    "  \"num_density_z_args\": {\"on\": 0.0, \"off\": 1.0, \"ampl\": 1.0},\n",
    "  \"interaction_strengths\": [10.0, 10.0],\n",
    "  \"savefile\": \"mbs-vee-sech-sqrt2pi-sqrt2pi\"\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mbs_s2 = mb_solve.MBSolve().from_json_str(mb_solve_json_s2)\n",
    "mbs_s2.mbsolve(recalc=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-s2-spacetime",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 12))\n",
    "cmap_range = np.linspace(0.0, 0.8, 11)\n",
    "\n",
    "ax = fig.add_subplot(211)\n",
    "cf = ax.contourf(mbs_s2.tlist, mbs_s2.zlist,\n",
    "                 np.abs(mbs_s2.Omegas_zt[0] / (2 * np.pi)), cmap_range, cmap=plt.cm.Blues)\n",
    "ax.set_title(r\"Rabi Frequency ($\\Gamma / 2\\pi$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Probe\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "ax = fig.add_subplot(212)\n",
    "cf = ax.contourf(mbs_s2.tlist, mbs_s2.zlist,\n",
    "                 np.abs(mbs_s2.Omegas_zt[1] / (2 * np.pi)), cmap_range, cmap=plt.cm.Greens)\n",
    "ax.set_xlabel(r\"Time ($1/\\Gamma$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Coupling\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    for y in [0.0, 1.0]:\n",
    "        ax.axhline(y, c=\"grey\", lw=1.0, ls=\"dotted\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-s2-area",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_area = np.sqrt(mbs_s2.fields_area()[0]**2 + mbs_s2.fields_area()[1]**2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "ax.plot(mbs_s2.zlist, mbs_s2.fields_area()[0] / np.pi, label=\"Probe\", clip_on=False)\n",
    "ax.plot(mbs_s2.zlist, mbs_s2.fields_area()[1] / np.pi, label=\"Coupling\", clip_on=False)\n",
    "ax.plot(mbs_s2.zlist, total_area / np.pi, label=\"Total\", ls=\"dashed\", clip_on=False)\n",
    "ax.legend()\n",
    "ax.set_ylim([0.0, 2.0])\n",
    "ax.set_xlabel(\"Distance ($L$)\")\n",
    "ax.set_ylabel(r\"Pulse Area ($\\pi$)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sec-s8-intro",
   "metadata": {},
   "source": [
    "## √8π soliton formation\n",
    "\n",
    "At $\\theta_\\mathrm{probe} = \\theta_\\mathrm{coupling} = \\sqrt{8}\\pi = 2\\sqrt{2}\\pi$ the\n",
    "combined area is $4\\pi$, supporting a two-soliton simulton. Each field pulse\n",
    "breaks into two SIT-soliton-like sub-pulses that travel together through the medium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "setup-s8",
   "metadata": {},
   "outputs": [],
   "source": [
    "ampl_s8 = np.sqrt(8) / t_width / (2 * np.pi)\n",
    "\n",
    "mb_solve_json_s8 = \"\"\"\n",
    "{\n",
    "  \"atom\": {\n",
    "    \"decays\": [{\"channels\": [[0, 1], [0, 2]], \"rate\": 0.0}],\n",
    "    \"energies\": [],\n",
    "    \"fields\": [\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 1]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"probe\",\n",
    "        \"rabi_freq\": 1.18567868227,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      },\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 2]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"coupling\",\n",
    "        \"rabi_freq\": 1.18567868227,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      }\n",
    "    ],\n",
    "    \"num_states\": 3\n",
    "  },\n",
    "  \"t_min\": -2.0,\n",
    "  \"t_max\": 10.0,\n",
    "  \"t_steps\": 60,\n",
    "  \"z_min\": -0.2,\n",
    "  \"z_max\": 1.2,\n",
    "  \"z_steps\": 70,\n",
    "  \"z_steps_inner\": 1,\n",
    "  \"interaction_strengths\": [10.0, 10.0],\n",
    "  \"savefile\": \"mbs-vee-sech-sqrt8pi-sqrt8pi\"\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mbs_s8 = mb_solve.MBSolve().from_json_str(mb_solve_json_s8)\n",
    "mbs_s8.mbsolve(recalc=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-s8-spacetime",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 12))\n",
    "cmap_range = np.linspace(0.0, 2.0, 11)\n",
    "\n",
    "ax = fig.add_subplot(211)\n",
    "cf = ax.contourf(mbs_s8.tlist, mbs_s8.zlist,\n",
    "                 np.abs(mbs_s8.Omegas_zt[0] / (2 * np.pi)), cmap_range, cmap=plt.cm.Blues)\n",
    "ax.set_title(r\"Rabi Frequency ($\\Gamma / 2\\pi$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Probe\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "ax = fig.add_subplot(212)\n",
    "cf = ax.contourf(mbs_s8.tlist, mbs_s8.zlist,\n",
    "                 np.abs(mbs_s8.Omegas_zt[1] / (2 * np.pi)), cmap_range, cmap=plt.cm.Greens)\n",
    "ax.set_xlabel(r\"Time ($1/\\Gamma$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Coupling\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    for y in [0.0, 1.0]:\n",
    "        ax.axhline(y, c=\"grey\", lw=1.0, ls=\"dotted\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-s8-area",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_area = np.sqrt(mbs_s8.fields_area()[0]**2 + mbs_s8.fields_area()[1]**2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "ax.plot(mbs_s8.zlist, mbs_s8.fields_area()[0] / np.pi, label=\"Probe\", clip_on=False)\n",
    "ax.plot(mbs_s8.zlist, mbs_s8.fields_area()[1] / np.pi, label=\"Coupling\", clip_on=False)\n",
    "ax.plot(mbs_s8.zlist, total_area / np.pi, label=\"Total\", ls=\"dashed\", clip_on=False)\n",
    "ax.legend()\n",
    "ax.set_ylim([0.0, 4.0])\n",
    "ax.set_xlabel(\"Distance ($L$)\")\n",
    "ax.set_ylabel(r\"Pulse Area ($\\pi$)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sec-coll-intro",
   "metadata": {},
   "source": [
    "## Soliton collision — solitons form simulton\n",
    "\n",
    "A $\\sqrt{8}\\pi$ sech pulse on the probe enters at $t = 0$; a second\n",
    "$\\sqrt{8}\\pi$ sech pulse on the coupling enters at $t = 3$. In the medium\n",
    "the two pulses interact via the atoms: the collision reshapes both into a\n",
    "shared simulton — demonstrating that solitons can lock together even when\n",
    "launched at different times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "setup-coll",
   "metadata": {},
   "outputs": [],
   "source": [
    "ampl_coll = np.sqrt(8) / t_width / (2 * np.pi)\n",
    "\n",
    "mb_solve_json_coll = \"\"\"\n",
    "{{\n",
    "  \"atom\": {{\n",
    "    \"fields\": [\n",
    "      {{\n",
    "        \"coupled_levels\": [[0, 1]],\n",
    "        \"label\": \"probe\",\n",
    "        \"rabi_freq\": 1.0,\n",
    "        \"rabi_freq_t_args\": {{\"ampl\": {ampl_coll}, \"centre\": 0.0, \"width\": {t_width}}},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      }},\n",
    "      {{\n",
    "        \"coupled_levels\": [[0, 2]],\n",
    "        \"label\": \"coupling\",\n",
    "        \"rabi_freq\": 1.0,\n",
    "        \"rabi_freq_t_args\": {{\"ampl\": {ampl_coll}, \"centre\": 3.0, \"width\": {t_width}}},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      }}\n",
    "    ],\n",
    "    \"num_states\": 3\n",
    "  }},\n",
    "  \"t_min\": -2.0,\n",
    "  \"t_max\": 14.0,\n",
    "  \"t_steps\": 200,\n",
    "  \"z_min\": -0.5,\n",
    "  \"z_max\": 1.5,\n",
    "  \"z_steps\": 500,\n",
    "  \"z_steps_inner\": 1,\n",
    "  \"interaction_strengths\": [50.0, 10.0],\n",
    "  \"savefile\": \"mbs-vee-sech-sqrt8pi-sqrt8pi-collision\"\n",
    "}}\n",
    "\"\"\".format(ampl_coll=ampl_coll, t_width=t_width)\n",
    "\n",
    "mbs_coll = mb_solve.MBSolve().from_json_str(mb_solve_json_coll)\n",
    "mbs_coll.mbsolve(recalc=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-coll-spacetime",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 12))\n",
    "\n",
    "ax = fig.add_subplot(211)\n",
    "cf = ax.contourf(mbs_coll.tlist, mbs_coll.zlist,\n",
    "                 np.abs(mbs_coll.Omegas_zt[0] / (2 * np.pi)), cmap=plt.cm.Blues)\n",
    "ax.set_title(r\"Rabi Frequency ($\\Gamma / 2\\pi$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Probe\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "ax = fig.add_subplot(212)\n",
    "cf = ax.contourf(mbs_coll.tlist, mbs_coll.zlist,\n",
    "                 np.abs(mbs_coll.Omegas_zt[1] / (2 * np.pi)), cmap=plt.cm.Greens)\n",
    "ax.set_xlabel(r\"Time ($1/\\Gamma$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Coupling\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    for y in [0.0, 1.0]:\n",
    "        ax.axhline(y, c=\"grey\", lw=1.0, ls=\"dotted\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-coll-area",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_area = np.sqrt(mbs_coll.fields_area()[0]**2 + mbs_coll.fields_area()[1]**2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "ax.plot(mbs_coll.zlist, mbs_coll.fields_area()[0] / np.pi, label=\"Probe\", clip_on=False)\n",
    "ax.plot(mbs_coll.zlist, mbs_coll.fields_area()[1] / np.pi, label=\"Coupling\", clip_on=False)\n",
    "ax.plot(mbs_coll.zlist, total_area / np.pi, label=\"Total\", ls=\"dashed\", clip_on=False)\n",
    "ax.legend()\n",
    "ax.set_ylim([0.0, 4.0])\n",
    "ax.set_xlabel(\"Distance ($L$)\")\n",
    "ax.set_ylabel(r\"Pulse Area ($\\pi$)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sec-surf2-intro",
   "metadata": {},
   "source": [
    "## Optical surfer — 2π coupling soliton on CW probe\n",
    "\n",
    "A weak CW probe ($\\Omega_\\mathrm{probe} = 10^{-3}\\,\\Gamma$) is ramped on at $t = -1$.\n",
    "A 2π sech soliton is launched on the coupling field at $t = 0$. The coupling\n",
    "soliton propagates as a SIT soliton; where it overlaps the probe, it opens a\n",
    "local transparency window that drags the probe forward — an **optical surfer**\n",
    "riding on the CW background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "setup-surf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "ampl_surf2 = 2.0 / t_width / (2 * np.pi)\n",
    "\n",
    "mb_solve_json_surf2 = \"\"\"\n",
    "{\n",
    "  \"atom\": {\n",
    "    \"fields\": [\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 1]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"probe\",\n",
    "        \"rabi_freq\": 1.0e-3,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"on\": -1.0, \"fwhm\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"ramp_on\"\n",
    "      },\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 2]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"coupling\",\n",
    "        \"rabi_freq\": 0.8384014365421667,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      }\n",
    "    ],\n",
    "    \"num_states\": 3\n",
    "  },\n",
    "  \"t_min\": -2.0,\n",
    "  \"t_max\": 10.0,\n",
    "  \"t_steps\": 120,\n",
    "  \"z_min\": -0.2,\n",
    "  \"z_max\": 1.2,\n",
    "  \"z_steps\": 140,\n",
    "  \"z_steps_inner\": 2,\n",
    "  \"interaction_strengths\": [10.0, 10.0],\n",
    "  \"savefile\": \"mbs-vee-weak-cw-sech-2pi\"\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mbs_surf2 = mb_solve.MBSolve().from_json_str(mb_solve_json_surf2)\n",
    "mbs_surf2.mbsolve(recalc=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-surf2-spacetime",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 12))\n",
    "\n",
    "ax = fig.add_subplot(211)\n",
    "cf = ax.contourf(mbs_surf2.tlist, mbs_surf2.zlist,\n",
    "                 np.abs(mbs_surf2.Omegas_zt[0] / (2 * np.pi)),\n",
    "                 np.linspace(0.0, 2e-3, 11), cmap=plt.cm.Blues)\n",
    "ax.set_title(r\"Rabi Frequency ($\\Gamma / 2\\pi$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Probe\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "ax = fig.add_subplot(212)\n",
    "cf = ax.contourf(mbs_surf2.tlist, mbs_surf2.zlist,\n",
    "                 np.abs(mbs_surf2.Omegas_zt[1] / (2 * np.pi)),\n",
    "                 np.linspace(0.0, 1.0, 11), cmap=plt.cm.Greens)\n",
    "ax.set_xlabel(r\"Time ($1/\\Gamma$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Coupling\", va=\"top\", ha=\"left\", transform=ax.transAxes, color=\"grey\", fontsize=16)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    for y in [0.0, 1.0]:\n",
    "        ax.axhline(y, c=\"grey\", lw=1.0, ls=\"dotted\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-surf2-area",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_area = np.sqrt(mbs_surf2.fields_area()[0]**2 + mbs_surf2.fields_area()[1]**2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "ax.plot(mbs_surf2.zlist, mbs_surf2.fields_area()[0] / np.pi, label=\"Probe\", clip_on=False)\n",
    "ax.plot(mbs_surf2.zlist, mbs_surf2.fields_area()[1] / np.pi, label=\"Coupling\", clip_on=False)\n",
    "ax.plot(mbs_surf2.zlist, total_area / np.pi, label=\"Total\", ls=\"dashed\", clip_on=False)\n",
    "ax.legend()\n",
    "ax.set_ylim([0.0, 2.0])\n",
    "ax.set_xlabel(\"Distance ($L$)\")\n",
    "ax.set_ylabel(r\"Pulse Area ($\\pi$)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sec-surf4-intro",
   "metadata": {},
   "source": [
    "## Double optical surfer — 4π coupling soliton on CW probe\n",
    "\n",
    "As for the 2π optical surfer, but with a 4π coupling soliton\n",
    "($\\theta_\\mathrm{coupling} = 4\\pi = 2 \\times 2\\pi$).\n",
    "The higher area produces **two** surfer solitons: the coupling pulse\n",
    "generates a pair of probe brightening regions that both propagate\n",
    "through the medium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "setup-surf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ampl_surf4 = 4.0 / t_width / (2 * np.pi)\n",
    "\n",
    "mb_solve_json_surf4 = \"\"\"\n",
    "{\n",
    "  \"atom\": {\n",
    "    \"fields\": [\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 1]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"probe\",\n",
    "        \"rabi_freq\": 1.0e-3,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"on\": -1.0, \"fwhm\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"ramp_on\"\n",
    "      },\n",
    "      {\n",
    "        \"coupled_levels\": [[0, 2]],\n",
    "        \"detuning\": 0.0,\n",
    "        \"detuning_positive\": true,\n",
    "        \"label\": \"coupling\",\n",
    "        \"rabi_freq\": 1.6768028730843334,\n",
    "        \"rabi_freq_t_args\": {\"ampl\": 1.0, \"centre\": 0.0, \"width\": 0.3796628587572578},\n",
    "        \"rabi_freq_t_func\": \"sech\"\n",
    "      }\n",
    "    ],\n",
    "    \"num_states\": 3\n",
    "  },\n",
    "  \"t_min\": -2.0,\n",
    "  \"t_max\": 10.0,\n",
    "  \"t_steps\": 120,\n",
    "  \"z_min\": -0.2,\n",
    "  \"z_max\": 1.2,\n",
    "  \"z_steps\": 140,\n",
    "  \"z_steps_inner\": 2,\n",
    "  \"interaction_strengths\": [10.0, 10.0],\n",
    "  \"savefile\": \"mbs-vee-weak-cw-sech-4pi\"\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "mbs_surf4 = mb_solve.MBSolve().from_json_str(mb_solve_json_surf4)\n",
    "mbs_surf4.mbsolve(recalc=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-surf4-spacetime",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 12))\n",
    "\n",
    "ax = fig.add_subplot(211)\n",
    "cf = ax.contourf(mbs_surf4.tlist, mbs_surf4.zlist,\n",
    "                 np.abs(mbs_surf4.Omegas_zt[0] / (2 * np.pi)),\n",
    "                 np.linspace(0.0, 2.5e-3, 11), cmap=plt.cm.Blues)\n",
    "ax.set_title(r\"Rabi Frequency ($\\Gamma / 2\\pi$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Probe\", va=\"top\", ha=\"left\", transform=ax.transAxes,\n",
    "        color=\"k\", fontsize=16, alpha=0.5)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "ax = fig.add_subplot(212)\n",
    "cf = ax.contourf(mbs_surf4.tlist, mbs_surf4.zlist,\n",
    "                 np.abs(mbs_surf4.Omegas_zt[1] / (2 * np.pi)),\n",
    "                 np.linspace(0.0, 2.5, 11), cmap=plt.cm.Greens)\n",
    "ax.set_xlabel(r\"Time ($1/\\Gamma$)\")\n",
    "ax.set_ylabel(\"Distance ($L$)\")\n",
    "ax.text(0.02, 0.95, \"Coupling\", va=\"top\", ha=\"left\", transform=ax.transAxes,\n",
    "        color=\"k\", fontsize=15, alpha=0.5)\n",
    "plt.colorbar(cf)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    for y in [0.0, 1.0]:\n",
    "        ax.axhline(y, c=\"grey\", lw=1.0, ls=\"dotted\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-surf4-area",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_area = np.sqrt(mbs_surf4.fields_area()[0]**2 + mbs_surf4.fields_area()[1]**2)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "ax.plot(mbs_surf4.zlist, mbs_surf4.fields_area()[0] / np.pi, label=\"Probe\", clip_on=False)\n",
    "ax.plot(mbs_surf4.zlist, mbs_surf4.fields_area()[1] / np.pi, label=\"Coupling\", clip_on=False)\n",
    "ax.plot(mbs_surf4.zlist, total_area / np.pi, label=\"Total\", ls=\"dashed\", clip_on=False)\n",
    "ax.legend()\n",
    "ax.set_ylim([0.0, 4.0])\n",
    "ax.set_xlabel(\"Distance ($L$)\")\n",
    "ax.set_ylabel(r\"Pulse Area ($\\pi$)\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}