Nonlinear propagation tutorial: Added soliton propagation example

docs/source/tutorials/nonlinear_propagation_split_step.ipynb

@@ -156,12 +156,13 @@
     "    \"\"\"Sellmeier equation for fused silica.\"\"\"\n",
     "    x = wavelength * 1e6  #  convert wavelength to µm\n",
     "\n",
-    "    return np.sqrt(\n",
+    "    n = np.sqrt(\n",
     "        1\n",
     "        + 0.6961663 / (1 - (0.0684043 / x) ** 2)\n",
     "        + 0.4079426 / (1 - (0.1162414 / x) ** 2)\n",
     "        + 0.8974794 / (1 - (9.896161 / x) ** 2)\n",
     "    )\n",
+    "    return n\n",
     "\n",
     "\n",
     "# intensity dependent refractive index\n",
@@ -351,6 +352,379 @@
     "# show laser after propagation\n",
     "laser.show(envelope_type=\"intensity\")"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fceec76b",
+   "metadata": {},
+   "source": [
+    "## Soliton propagation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0ebfbff1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "from scipy.constants import c\n",
+    "\n",
+    "from lasy.laser import Laser\n",
+    "from lasy.profiles import CombinedLongitudinalTransverseProfile\n",
+    "from lasy.profiles.longitudinal import Sech2Profile\n",
+    "from lasy.profiles.transverse import PlaneWaveProfile\n",
+    "from lasy.utils.laser_utils import get_duration, get_laser_power"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0b01f2b6",
+   "metadata": {},
+   "source": [
+    "To verify the validity of the split-step approach, we can simulate the propagation of a soliton pulse, which is a special case where a $\\mathrm{sech}^2$ shaped pulse keeps its temporal shape due to a balance of group velocity dispersion and self-phase modulation. This is typically found in fiber lasers, where the pulses are spatially confined and pulse durations are moderate, such that higher order dispersion plays a subordinate role.\n",
+    "In simple fused silica based fibers, this works only at relatively long wavelengths such as the 1.55µm commonly found in telecom fiber systems. This is because only at these wavelengths does fused silica introduce negative dispersion, that is required to balance the positive dispersion introduced by self-phase modulation.\n",
+    "\n",
+    "We'll start by defining a laser pulse at this wavelength. To minimize the influence of higher order dispersion, we use a moderately long pulse duration of 1ps."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "525a97c3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Physical Parameters\n",
+    "wavelength = 1550e-9  # wavelength in meters\n",
+    "tau = 1000e-15 / (\n",
+    "    2 * np.arccosh(np.sqrt(2))\n",
+    ")  # FWHM pulse duration (1ps) converted to the sech pulse duration\n",
+    "t_peak = 0.0  # peak position in time\n",
+    "pol = (1, 0)  # polarisation vector of the laser"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8a73f6a3",
+   "metadata": {},
+   "source": [
+    "To mimic the spatially confined propagation, we run our simulation in 1D, which can be done by combining `PlaneWaveProfile` in the transverse plane and a `Sech2Profile` in time."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cd743cdb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# longitudinal profile\n",
+    "long_prof = Sech2Profile(wavelength, tau, t_peak)\n",
+    "\n",
+    "# transverse profile\n",
+    "tran_prof = PlaneWaveProfile()\n",
+    "\n",
+    "# combined profile\n",
+    "profile = CombinedLongitudinalTransverseProfile(\n",
+    "    wavelength, pol, long_prof, tran_prof, peak_power=1.0\n",
+    ")\n",
+    "\n",
+    "# Computational Grid\n",
+    "dim = \"xyt\"\n",
+    "lo = (None, None, -10 * tau)\n",
+    "hi = (None, None, 10 * tau)\n",
+    "npoints = (1, 1, 1000)\n",
+    "\n",
+    "laser = Laser(dim, lo, hi, npoints, profile)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "714f9d50",
+   "metadata": {},
+   "source": [
+    "As a next step, we can define the parmeters of the fiber. In this case we take a 1km long fiber, with the correct intensity dependent refractive index $n_2$ at the 1.55µm wavelength."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7824edf1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n2 = 2.41e-20  # intensity dependent refractive index\n",
+    "d_FS = 1000  # meters, total propagation distance in fused silica"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "28df40dc",
+   "metadata": {},
+   "source": [
+    "Now that we have set the pulse and material properties, we can calculate the peak power of the pulse at which the soliton condition is met.\n",
+    "For this we need to calculate the dispersion parameter $\\frac{\\partial^2 k}{\\partial \\omega} = \\beta_2$ of the fiber from its refractive index."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6116fcf9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.interpolate import interp1d\n",
+    "\n",
+    "\n",
+    "def get_beta2(n, omega, omega0):\n",
+    "    \"\"\"Calculate the second order dispersion coefficient.\"\"\"\n",
+    "    # ensure that omega is sorted\n",
+    "    order = np.argsort(omega)\n",
+    "    omega = omega[order]\n",
+    "    n = n[order]\n",
+    "\n",
+    "    # calculate the wavevector k\n",
+    "    k = n * omega / c\n",
+    "\n",
+    "    # second derivative of k with respect to omega\n",
+    "    beta2 = np.gradient(np.gradient(k, omega), omega)\n",
+    "\n",
+    "    # evaluate at omega0\n",
+    "    beta2_0 = interp1d(omega, beta2, bounds_error=True)(omega0)\n",
+    "    return beta2_0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "20bddde2",
+   "metadata": {},
+   "source": [
+    "With this function we can calculate the dispersion parameter beta2 at the central frequency of the laser"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "763479fd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create frequency axis\n",
+    "_, omega = laser.grid.get_spectral_field()\n",
+    "omega += laser.profile.omega0\n",
+    "\n",
+    "# get refractive index\n",
+    "refractive_index = n_fusedsilica(2 * np.pi * c / omega)\n",
+    "\n",
+    "# calculate the dispersion parameter\n",
+    "beta2 = get_beta2(refractive_index, omega, laser.profile.omega0)\n",
+    "print(f\"Dispersion parameter: {beta2 * 1e27:.3f} fs^2/mm\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "289d8f75",
+   "metadata": {},
+   "source": [
+    "With this we can then calculate the soliton power at which SPM and dispersion are perfectly balanced, and normalize the laser power to this value."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9cef2882",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "soliton_power = abs(beta2) * c / (tau**2 * n2 * laser.profile.omega0)\n",
+    "laser.normalize(soliton_power, kind=\"peak_power\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "856eed03",
+   "metadata": {},
+   "source": [
+    "Check that the temporal shape is as expected."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e8884997",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "power = get_laser_power(laser.dim, laser.grid)\n",
+    "\n",
+    "plt.figure()\n",
+    "plt.plot(laser.grid.axes[-1] * 1e15, power / 1e12)\n",
+    "plt.xlim(laser.grid.axes[-1][0] * 1e15, laser.grid.axes[-1][-1] * 1e15)\n",
+    "plt.ylim(0, None)\n",
+    "plt.xlabel(\"Time (fs)\")\n",
+    "plt.ylabel(\"Instantaneous Power (TW)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d995b47d",
+   "metadata": {},
+   "source": [
+    "Define our propagator and nonlinear step:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8a68d5a0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "linear_propagator = AngularSpectrumPropagator(\n",
+    "    omega0=profile.omega0, n=n_fusedsilica, dim=laser.dim\n",
+    ")\n",
+    "nonlinear_phase = NonlinearKerrStep(n2=n2, k0=laser.profile.omega0 / c)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b8067956",
+   "metadata": {},
+   "source": [
+    "Run the actual split-step propagation:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ce8f875a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dz = 100e-3  # meters, step size\n",
+    "\n",
+    "Nsteps = int(d_FS / dz)  # number of propagation steps\n",
+    "\n",
+    "# add propagators to the laser\n",
+    "laser.add_propagator(linear_propagator)\n",
+    "\n",
+    "peak_power = []\n",
+    "duration = []\n",
+    "shape = []\n",
+    "err_rel = []\n",
+    "temporal_field_initial = laser.grid.get_temporal_field()\n",
+    "\n",
+    "# propagate using second order split step method\n",
+    "for i in range(Nsteps):\n",
+    "    laser.propagate(dz / 2)\n",
+    "    laser.grid = nonlinear_phase.apply(grid_in=laser.grid, distance=dz)\n",
+    "    laser.propagate(dz / 2)\n",
+    "\n",
+    "    duration.append(get_duration(laser.grid, laser.dim, level=0.5))\n",
+    "    peak_power.append(np.max(get_laser_power(laser.dim, laser.grid)))\n",
+    "    shape.append(get_laser_power(laser.dim, laser.grid))\n",
+    "\n",
+    "    err_abs = np.abs(\n",
+    "        abs(laser.grid.get_temporal_field()) - abs(temporal_field_initial)\n",
+    "    ).sum()\n",
+    "    err_norm = max(\n",
+    "        np.abs(laser.grid.get_temporal_field()).sum(),\n",
+    "        (np.abs(temporal_field_initial)).sum(),\n",
+    "    )\n",
+    "\n",
+    "    if err_norm > 0:\n",
+    "        err_rel.append(err_abs / err_norm)\n",
+    "    else:\n",
+    "        err_rel.append(0.0)\n",
+    "\n",
+    "    print(f\"Progress: {100 * (i + 1) / Nsteps:.1f}%\", end=\"\\r\")\n",
+    "\n",
+    "duration = np.array(duration)\n",
+    "peak_power = np.array(peak_power)\n",
+    "\n",
+    "shape = np.array(shape)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3506536b",
+   "metadata": {},
+   "source": [
+    "If we now plot the result, we can see that the pulse shape remains constant during the propagation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9c36484a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.imshow(\n",
+    "    shape.T[::10] * 1e-12,\n",
+    "    aspect=\"auto\",\n",
+    "    extent=[0, d_FS, lo[-1] * 1e15, hi[-1] * 1e15],\n",
+    "    cmap=\"Reds\",\n",
+    ")\n",
+    "plt.xlabel(\"Propagation distance (m)\")\n",
+    "plt.ylabel(\"Time (fs)\")\n",
+    "plt.colorbar(label=\"Peak intensity (TW/m$^2$)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "285585a5",
+   "metadata": {},
+   "source": [
+    "We can also look at the relative changes of the pulse duration and peak power during the kilometer of propagation, and see that they are constant to the sub-permille-level, as one would expect from the propagation of a perfect soliton pulse."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "928a3411",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z = np.linspace(0, d_FS, Nsteps)\n",
+    "fig, ax = plt.subplots(2, 1)\n",
+    "\n",
+    "ax[0].plot(z, 1e15 * duration, label=\"Pulse duration\", linestyle=\"-\", color=\"C0\")\n",
+    "ax[0].set_ylabel(r\"Duration [fs]\")\n",
+    "ax[0].legend(loc=\"upper left\")\n",
+    "ax[0].set_ylim(duration.min() * 0.91 * 1e15, duration.max() * 1.11 * 1e15)\n",
+    "ax01 = ax[0].twinx()\n",
+    "\n",
+    "ax01.plot(z, peak_power * 1e-12, label=\"Peak power\", color=\"C1\", linestyle=\"-\")\n",
+    "ax01.set_ylabel(\"P [TW/m$^2$]\")\n",
+    "ax01.legend(loc=\"upper right\")\n",
+    "ax01.set_ylim(peak_power.min() * 0.9 * 1e-12, peak_power.max() * 1.1 * 1e-12)\n",
+    "ax01.set_xlim(z.min(), z.max())\n",
+    "\n",
+    "ax[1].plot(z, duration / duration[0], label=\"Pulse duration\")\n",
+    "ax[1].set_ylabel(r\"$\\tau/\\tau_0$\")\n",
+    "ax[1].legend(loc=\"upper left\")\n",
+    "ax[1].set_ylim(1 - 2e-3, 1 + 2e-3)\n",
+    "ax[1].set_xlabel(\"Propagation distance (m)\")\n",
+    "ax11 = ax[1].twinx()\n",
+    "\n",
+    "ax11.plot(z, peak_power / peak_power[0], label=\"Peak power\", color=\"C1\")\n",
+    "ax11.set_ylabel(\"P/P$_0$\")\n",
+    "ax11.legend(loc=\"upper right\")\n",
+    "ax11.set_ylim(1 - 2e-3, 1 + 2e-3)\n",
+    "ax11.set_xlim(z.min(), z.max())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a4e30cfd",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {