Implement second order concatenation

doc/source/examples/advanced_concatenation.ipynb

@@ -239,7 +239,9 @@
     "    for i in range(3):\n",
     "        for j in range(3):\n",
     "            t = ax[i, j].get_title()[7:]\n",
-    "            ax[i, j].set_title('$F^{(' + pulses[i] + pulses[j] + t)"
+    "            ax[i, j].set_title('$F^{(' + pulses[i] + pulses[j] + t)\n",
+    "    fig.suptitle(f'{gate_type} gates')\n",
+    "    fig.tight_layout()"
    ]
   },
   {
@@ -254,12 +256,81 @@
     "\n",
     "While the imaginary part cancels out when calculating fidelities, $\\mathcal{I}\\propto\\sum_{ij} \\int\\mathrm{d}\\omega S(\\omega)F^{(ij)}(\\omega)$, the real part does not and the offdiagonals therefore lead to corrections in the total fidelity of a composite pulse, $\\mathcal{I}_\\text{tot}\\neq\\sum_g\\mathcal{I}^{(g)}$ with $\\mathcal{I}^{(g)}$ the infidelities of the individual pulses. These corrections can thus in principle also be negative, leading to improved fidelities for composite pulses."
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Second-order filter functions\n",
+    "We can also compute the second-order filter function as the concatenation of individual pulses. Although the calculation is not completely independent of the internal dynamics of the individual pulses as in the first-order case due to the nested time integral, a significant chunk can still be obtained from previously computed quantities. \n",
+    "\n",
+    "To do this, we need to make sure that some intermediate results are cached when computing the first- and second-order filter functions of the individual pulses, setting the `cache_intermediates` flag. We then just set the `calc_second_order_FF` argument to `True`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Clean up from before to have a clean slate for demonstration purposes\n",
+    "from itertools import chain\n",
+    "\n",
+    "for key, pulse in chain(X2.items(), Y2.items()):\n",
+    "    pulse.cleanup('frequency dependent')\n",
+    "    pulse.cache_filter_function(omega[key], cache_intermediates=True,\n",
+    "                               order=1)\n",
+    "    pulse.cache_filter_function(omega[key], cache_intermediates=True,\n",
+    "                               order=2)\n",
+    "\n",
+    "H = {key: ff.concatenate((Y2, X2, X2), calc_second_order_FF=True, which='generalized')\n",
+    "     for (key, X2), (key, Y2) in zip(X2.items(), Y2.items())}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Visualize the FFs:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "for key, pulse in H.items():\n",
+    "    fig, axes = plt.subplots(3, 3, sharex=True, sharey=False,\n",
+    "                            figsize=(9, 6))\n",
+    "    for i in range(1, 4):\n",
+    "        for j in range(1, 4):\n",
+    "            FF1 = pulse.get_filter_function(omega[key], order=1,\n",
+    "                                            which='generalized')\n",
+    "            FF2 = pulse.get_filter_function(omega[key], order=2,\n",
+    "                                            which='generalized')\n",
+    "            axes[i-1, j-1].semilogx(omega[key], FF1[0,0,i,j].real,\n",
+    "                                    color='tab:blue', ls='-')\n",
+    "            axes[i-1, j-1].semilogx(omega[key], FF1[0,0,i,j].imag,\n",
+    "                                    color='tab:blue', ls='--')\n",
+    "            axes[i-1, j-1].semilogx(omega[key], FF2[0,0,i,j].real,\n",
+    "                                    color='tab:orange', ls='-')\n",
+    "            axes[i-1, j-1].semilogx(omega[key], FF2[0,0,i,j].imag,\n",
+    "                                    color='tab:orange', ls='--')\n",
+    "            axes[i-1, j-1].margins(x=0)\n",
+    "    fig.suptitle(f'{key} gates')\n",
+    "    fig.supxlabel(r'$\\omega$')\n",
+    "    fig.supylabel(r'$\\mathcal{F}(\\omega)$')\n",
+    "    fig.tight_layout()"
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (Spyder)",
-   "language": "python3",
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
    "name": "python3"
   },
   "language_info": {
@@ -272,7 +343,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.13.5"
   },
   "widgets": {
    "application/vnd.jupyter.widget-state+json": {