diff --git a/docs/eta-pi-p/manual.ipynb b/docs/eta-pi-p/manual.ipynb
index cd4b6c0..5e4cce5 100644
--- a/docs/eta-pi-p/manual.ipynb
+++ b/docs/eta-pi-p/manual.ipynb
@@ -11,7 +11,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "This section is a follow-up of previous chapter:[Reaction and Models](reaction-model.md), to formulate the amplitude model for the $\\gamma p \\to \\eta\\pi^0 p$ channel example symbolically. See **[TR‑033](https://compwa.github.io/report/033)** for a purely numerical tutorial.\n",
+    "This section is a follow-up of the previous chapter, [](reaction-model.md), to formulate the amplitude model for the $\\gamma p \\to \\eta\\pi^0 p$ channel example symbolically. See **[TR‑033](https://compwa.github.io/report/033)** for a purely numerical tutorial.\n",
     "\n",
     "The model we want to implement is"
    ]
@@ -103,6 +103,15 @@
     "## Model implementation"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "l_max = 2"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -116,11 +125,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "s12, m_a2, gamma_a2 = sp.symbols(r\"s_{12} m_{a_2} \\Gamma_{a_2}\")\n",
+    "s12, m_a2, gamma_a2, l12 = sp.symbols(r\"s_{12} m_{a_2} \\Gamma_{a_2} l_{12}\")\n",
     "theta1, phi1 = sp.symbols(\"theta_1 phi_1\")\n",
     "a = sp.IndexedBase(\"a\")\n",
     "m = sp.symbols(\"m\", cls=sp.Idx)\n",
-    "A12 = sp.Sum(a[m] * sp.Ynm(2, m, theta1, phi1), (m, -2, 2)) / (\n",
+    "A12 = sp.Sum(a[m] * sp.Ynm(l12, m, theta1, phi1), (m, -l12, l12)) / (\n",
     "    s12 - m_a2**2 + sp.I * m_a2 * gamma_a2\n",
     ")\n",
     "A12"
@@ -131,12 +140,12 @@
    "execution_count": null,
    "metadata": {
     "tags": [
-     "hide-cell"
+     "remove-cell"
     ]
    },
    "outputs": [],
    "source": [
-    "sp.Ynm(2, 1, 1, 1).expand(func=False)"
+    "A12.expand(func=True)"
    ]
   },
   {
@@ -144,18 +153,22 @@
    "execution_count": null,
    "metadata": {
     "tags": [
-     "hide-cell"
+     "remove-cell"
     ]
    },
    "outputs": [],
    "source": [
-    "sp.Ynm(2, 1, 1, 1).expand(func=True)"
+    "A12.expand(func=False)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "tags": [
+     "remove-cell"
+    ]
+   },
    "outputs": [],
    "source": [
     "A12.free_symbols"
@@ -167,11 +180,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "A12_func = sp.lambdify(\n",
-    "    [s12, a[-2], a[-1], a[0], a[1], a[2], m_a2, gamma_a2, theta1, phi1],\n",
-    "    A12.doit().expand(func=True),\n",
-    ")\n",
-    "A12_func"
+    "A12_funcs = [\n",
+    "    sp.lambdify(\n",
+    "        [s12, *(a[j] for j in range(-l_max, l_max + 1)), m_a2, gamma_a2, theta1, phi1],\n",
+    "        expr=A12.subs(l12, i).doit().expand(func=True),\n",
+    "    )\n",
+    "    for i in range(l_max + 1)\n",
+    "]"
    ]
   },
   {
@@ -187,11 +202,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "s23, m_delta, gamma_delta = sp.symbols(r\"s_{23} m_{\\Delta^+} \\Gamma_{\\Delta^+}\")\n",
+    "s23, m_delta, gamma_delta, l23 = sp.symbols(\n",
+    "    r\"s_{23} m_{\\Delta^+} \\Gamma_{\\Delta^+} l_{23}\"\n",
+    ")\n",
     "b = sp.IndexedBase(\"b\")\n",
     "m = sp.symbols(\"m\", cls=sp.Idx)\n",
     "theta2, phi2 = sp.symbols(\"theta_2 phi_2\")\n",
-    "A23 = sp.Sum(b[m] * sp.Ynm(1, m, theta2, phi2), (m, -1, 1)) / (\n",
+    "A23 = sp.Sum(b[m] * sp.Ynm(l23, m, theta2, phi2), (m, -l23, l23)) / (\n",
     "    s23 - m_delta**2 + sp.I * m_delta * gamma_delta\n",
     ")\n",
     "A23"
@@ -203,10 +220,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "A23_func = sp.lambdify(\n",
-    "    [s23, b[-1], b[0], b[1], m_delta, gamma_delta, theta2, phi2],\n",
-    "    A23.doit().expand(func=True),\n",
-    ")"
+    "A23_funcs = [\n",
+    "    sp.lambdify(\n",
+    "        [\n",
+    "            s23,\n",
+    "            *(b[j] for j in range(-l_max, l_max + 1)),\n",
+    "            m_delta,\n",
+    "            gamma_delta,\n",
+    "            theta2,\n",
+    "            phi2,\n",
+    "        ],\n",
+    "        A23.subs(l23, i).doit().expand(func=True),\n",
+    "    )\n",
+    "    for i in range(l_max + 1)\n",
+    "]"
    ]
   },
   {
@@ -222,35 +249,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "s31, m_nstar, gamma_nstar = sp.symbols(r\"s_{31} m_{N^*} \\Gamma_{N^*}\")\n",
     "c = sp.IndexedBase(\"c\")\n",
-    "A31 = c[0] / (s31 - m_nstar**2 + sp.I * m_nstar * gamma_nstar)\n",
-    "A31"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    ":::{tip}\n",
-    "The expression above is originally deduced from $Y_l^m (\\Omega_3)$ when $l=0$ as shown below, so that the numerator is absorbed and becomes $c_0$.\n",
-    ":::"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "tags": [
-     "hide-input"
-    ]
-   },
-   "outputs": [],
-   "source": [
-    "theta3, phi3 = sp.symbols(\"theta_3 phi_3\")\n",
-    "sp.Sum(c[m] * sp.Ynm(0, m, theta3, phi3), (m, 0, 0)) / (\n",
+    "s31, m_nstar, gamma_nstar = sp.symbols(r\"s_{31} m_{N^*} \\Gamma_{N^*}\")\n",
+    "theta3, phi3, l31 = sp.symbols(\"theta_3 phi_3 l_{31}\")\n",
+    "A31 = sp.Sum(c[m] * sp.Ynm(l31, m, theta3, phi3), (m, -l31, l31)) / (\n",
     "    s31 - m_nstar**2 + sp.I * m_nstar * gamma_nstar\n",
-    ")"
+    ")\n",
+    "A31"
    ]
   },
   {
@@ -259,7 +264,20 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "A31_func = sp.lambdify([s31, c[0], m_nstar, gamma_nstar], A31)"
+    "A31_funcs = [\n",
+    "    sp.lambdify(\n",
+    "        [\n",
+    "            s31,\n",
+    "            *(c[j] for j in range(-l_max, l_max + 1)),\n",
+    "            m_nstar,\n",
+    "            gamma_nstar,\n",
+    "            theta3,\n",
+    "            phi3,\n",
+    "        ],\n",
+    "        A31.subs(l31, i).doit().expand(func=True),\n",
+    "    )\n",
+    "    for i in range(l_max + 1)\n",
+    "]"
    ]
   },
   {
@@ -339,6 +357,9 @@
     "jupyter": {
      "source_hidden": true
     },
+    "mystnb": {
+     "code_prompt_show": "Helper functions for formulating symbolic kinematics"
+    },
     "tags": [
      "hide-cell"
     ]
@@ -451,39 +472,45 @@
      "source_hidden": true
     },
     "tags": [
-     "hide-input"
+     "hide-input",
+     "hide-output",
+     "scroll-output"
     ]
    },
    "outputs": [],
    "source": [
     "a_vals = [0, 0.5, 3.5, 4, 2.5]\n",
-    "b_vals = [-1.5, 4, 0.5]\n",
-    "c0_val = 2.5\n",
+    "b_vals = [0, -1.5, 4, 0.5, 0]\n",
+    "c_vals = [0, 0, 2.5, 0, 0]\n",
     "m_a2_val = 1.32\n",
     "m_delta_val = 1.54\n",
     "m_nstar_val = 1.87\n",
     "gamma_a2_val = 0.1\n",
     "gamma_delta_val = 0.1\n",
     "gamma_nstar_val = 0.1\n",
+    "l12_val = 2\n",
+    "l23_val = 1\n",
+    "l31_val = 0\n",
     "\n",
     "parameters_default = {\n",
-    "    a[-2]: a_vals[0],\n",
-    "    a[-1]: a_vals[1],\n",
-    "    a[0]: a_vals[2],\n",
-    "    a[1]: a_vals[3],\n",
-    "    a[2]: a_vals[4],\n",
-    "    b[-1]: b_vals[0],\n",
-    "    b[0]: b_vals[1],\n",
-    "    b[1]: b_vals[2],\n",
-    "    c[0]: c0_val,\n",
     "    m_a2: m_a2_val,\n",
     "    m_delta: m_delta_val,\n",
     "    m_nstar: m_nstar_val,\n",
     "    gamma_a2: gamma_a2_val,\n",
     "    gamma_delta: gamma_delta_val,\n",
     "    gamma_nstar: gamma_nstar_val,\n",
+    "    l12: l12_val,\n",
+    "    l23: l23_val,\n",
+    "    l31: l31_val,\n",
     "}\n",
     "\n",
+    "a_dict = {a[i]: a_vals[i + l_max] for i in range(-l_max, l_max + 1)}\n",
+    "b_dict = {b[i]: b_vals[i + l_max] for i in range(-l_max, l_max + 1)}\n",
+    "c_dict = {c[i]: c_vals[i + l_max] for i in range(-l_max, l_max + 1)}\n",
+    "parameters_default.update(a_dict)\n",
+    "parameters_default.update(b_dict)\n",
+    "parameters_default.update(c_dict)\n",
+    "\n",
     "Latex(aslatex(parameters_default))"
    ]
   },
@@ -519,9 +546,9 @@
     "phi = np.pi / 4\n",
     "theta = np.pi / 4\n",
     "s_values = np.linspace(0, 5, num=500)\n",
-    "A12_values = A12_func(s_values, *a_vals, m_a2_val, gamma_a2_val, theta, phi)\n",
-    "A23_values = A23_func(s_values, *b_vals, m_delta_val, gamma_delta_val, theta, phi)\n",
-    "A31_values = A31_func(s_values, c0_val, m_nstar_val, gamma_nstar_val)"
+    "A12_values = A12_funcs[2](s_values, *a_vals, m_a2_val, gamma_a2_val, theta, phi)\n",
+    "A23_values = A23_funcs[1](s_values, *b_vals, m_delta_val, gamma_delta_val, theta, phi)\n",
+    "A31_values = A31_funcs[0](s_values, *c_vals, m_nstar_val, gamma_nstar_val, theta, phi)"
    ]
   },
   {
@@ -633,6 +660,9 @@
     "jupyter": {
      "source_hidden": true
     },
+    "mystnb": {
+     "code_prompt_show": "Slider definitions"
+    },
     "tags": [
      "hide-input",
      "scroll-input"
@@ -640,10 +670,25 @@
    },
    "outputs": [],
    "source": [
+    "def format_l_as_unicode(name: str) -> str:\n",
+    "    replacements = {\n",
+    "        \"l\": \"ℓ\",\n",
+    "        \"1\": \"₁\",\n",
+    "        \"2\": \"₂\",\n",
+    "        \"3\": \"₃\",\n",
+    "        \"_{\": \"\",\n",
+    "        \"}\": \"\",\n",
+    "    }\n",
+    "    for old, new in replacements.items():\n",
+    "        name = name.replace(old, new)\n",
+    "    return name\n",
+    "\n",
+    "\n",
     "sliders = {}\n",
     "categorized_sliders_m = defaultdict(list)\n",
     "categorized_sliders_gamma = defaultdict(list)\n",
     "categorized_cphi_pair = defaultdict(list)\n",
+    "categorized_sliders_l = []\n",
     "\n",
     "for symbol, value in parameters_default.items():\n",
     "    if symbol.name.startswith(R\"\\Gamma_{\"):\n",
@@ -680,7 +725,7 @@
     "        elif symbol.name.startswith(\"m_{a\"):\n",
     "            categorized_sliders_m[2].append(slider)\n",
     "\n",
-    "    else:\n",
+    "    elif isinstance(symbol, sp.Indexed):\n",
     "        c_latex = sp.latex(symbol)\n",
     "        phi_latex = Rf\"\\phi_{{{c_latex}}}\"\n",
     "        slider_c = w.FloatSlider(\n",
@@ -712,15 +757,32 @@
     "        else:\n",
     "            raise NotImplementedError(symbol.name)\n",
     "\n",
+    "    elif symbol.name.startswith(\"l_{\"):\n",
+    "        slider = w.IntSlider(\n",
+    "            value=1,\n",
+    "            min=0,\n",
+    "            max=2,\n",
+    "            description=format_l_as_unicode(symbol.name),\n",
+    "            continuous_update=False,\n",
+    "        )\n",
+    "        sliders[symbol.name] = slider\n",
+    "        categorized_sliders_l.append(slider)\n",
+    "\n",
+    "    else:\n",
+    "        raise NotImplementedError(symbol.name)\n",
+    "\n",
     "tab_contents = []\n",
-    "resonances_name = [\"N*\", \"Δ⁺\", \"a₂\"]\n",
+    "resonances_name = [\"N* [pη (31)]\", \"Δ⁺ [π⁰p(23)]\", \"a₂ [ηπ⁰(12)]\"]\n",
     "for i in range(len(resonances_name)):\n",
     "    tab_content = w.VBox([\n",
     "        w.HBox(categorized_sliders_m[i] + categorized_sliders_gamma[i]),\n",
     "        w.VBox(categorized_cphi_pair[i]),\n",
     "    ])\n",
     "    tab_contents.append(tab_content)\n",
-    "UI = w.Tab(tab_contents, titles=resonances_name)"
+    "UI = w.VBox([\n",
+    "    w.HBox(categorized_sliders_l),\n",
+    "    w.Tab(tab_contents, titles=resonances_name),\n",
+    "])"
    ]
   },
   {
@@ -729,21 +791,23 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unfolded_expression = intensity_expr.doit().expand(func=True)\n",
-    "intensity_func = create_parametrized_function(\n",
-    "    expression=unfolded_expression,\n",
-    "    parameters=parameters_default,\n",
-    "    backend=\"jax\",\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "intensities = intensity_func(phsp)"
+    "intensity_funcs = np.array([\n",
+    "    [\n",
+    "        [\n",
+    "            create_parametrized_function(\n",
+    "                expression=intensity_expr.xreplace({l12: i, l23: j, l31: k})\n",
+    "                .doit()\n",
+    "                .expand(func=True),\n",
+    "                parameters=parameters_default,\n",
+    "                backend=\"jax\",\n",
+    "            )\n",
+    "            for i in range(l_max + 1)\n",
+    "        ]\n",
+    "        for j in range(l_max + 1)\n",
+    "    ]\n",
+    "    for k in range(l_max + 1)\n",
+    "])\n",
+    "intensity_funcs.shape"
    ]
   },
   {
@@ -753,6 +817,9 @@
     "jupyter": {
      "source_hidden": true
     },
+    "mystnb": {
+     "code_prompt_show": "Function for making merging sliders into complex values"
+    },
     "tags": [
      "hide-input"
     ]
@@ -783,7 +850,8 @@
     },
     "tags": [
      "hide-input",
-     "full-width"
+     "full-width",
+     "scroll-input"
     ]
    },
    "outputs": [],
@@ -800,6 +868,10 @@
     "\n",
     "def update_histogram(**parameters):\n",
     "    global mesh\n",
+    "    l12 = parameters[\"l_{12}\"]\n",
+    "    l23 = parameters[\"l_{23}\"]\n",
+    "    l31 = parameters[\"l_{31}\"]\n",
+    "    intensity_func = intensity_funcs[l12, l23, l31]\n",
     "    parameters = insert_phi(parameters)\n",
     "    intensity_func.update_parameters(parameters)\n",
     "    intensity_weights = intensity_func(phsp)\n",
@@ -810,7 +882,7 @@
     "        weights=intensity_weights,\n",
     "        density=True,\n",
     "    )\n",
-    "    bin_values = jnp.where(bin_values < 1e-6, jnp.nan, bin_values)\n",
+    "    bin_values = jnp.where(bin_values < 1e-10, jnp.nan, bin_values)\n",
     "    x, y = jnp.meshgrid(xedges[:-1], yedges[:-1])\n",
     "    if mesh is None:\n",
     "        mesh = ax_2d.pcolormesh(x, y, bin_values.T, cmap=\"jet\", vmax=0.15)\n",
@@ -848,7 +920,8 @@
     },
     "tags": [
      "hide-input",
-     "full-width"
+     "full-width",
+     "scroll-input"
     ]
    },
    "outputs": [],
@@ -898,6 +971,10 @@
     "\n",
     "\n",
     "def update_plot(**parameters):  # noqa: C901, PLR0912, PLR0914\n",
+    "    l12 = parameters[\"l_{12}\"]\n",
+    "    l23 = parameters[\"l_{23}\"]\n",
+    "    l31 = parameters[\"l_{31}\"]\n",
+    "    intensity_func = intensity_funcs[l12, l23, l31]\n",
     "    parameters = insert_phi(parameters)\n",
     "    intensity_func.update_parameters(parameters)\n",
     "    intensities = intensity_func(phsp)\n",
diff --git a/pyproject.toml b/pyproject.toml
index 00c3b80..db7ad2a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -119,6 +119,7 @@ split-on-trailing-comma = false
     "PLC2701",
     "PLR2004",
     "PLW0603",
+    "RUF001",
     "S101",
     "T201",
     "TCH",