diff --git a/docs/_extend_docstrings.py b/docs/_extend_docstrings.py
index f01d1203d..6a7bd4327 100644
--- a/docs/_extend_docstrings.py
+++ b/docs/_extend_docstrings.py
@@ -295,16 +295,16 @@ def extend_EuclideanNorm() -> None:
     _append_code_rendering(expr)
 
 
-def extend_formulate_form_factor() -> None:
-    from ampform.dynamics import formulate_form_factor
+def extend_FormFactor() -> None:
+    from ampform.dynamics.form_factor import FormFactor
 
     s, m_a, m_b, L, d = sp.symbols("s m_a m_b L d")
-    form_factor = formulate_form_factor(s, m_a, m_b, angular_momentum=L, meson_radius=d)
+    form_factor = FormFactor(s, m_a, m_b, angular_momentum=L, meson_radius=d)
     _append_to_docstring(
-        formulate_form_factor,
+        FormFactor,
         f"""
     .. math:: {sp.latex(form_factor)}
-        :label: formulate_form_factor
+        :label: FormFactor
     """,
     )
 
diff --git a/docs/usage/dynamics.ipynb b/docs/usage/dynamics.ipynb
index 629aaa20d..a19e15b87 100644
--- a/docs/usage/dynamics.ipynb
+++ b/docs/usage/dynamics.ipynb
@@ -91,15 +91,6 @@
     "```"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "%matplotlib widget"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -116,8 +107,6 @@
    },
    "outputs": [],
    "source": [
-    "%config InlineBackend.figure_formats = ['svg']\n",
-    "\n",
     "import logging\n",
     "import warnings\n",
     "\n",
@@ -159,11 +148,11 @@
    },
    "outputs": [],
    "source": [
-    "from ampform.dynamics import BlattWeisskopfSquared\n",
+    "from ampform.dynamics.form_factor import BlattWeisskopfSquared\n",
     "\n",
     "L = sp.Symbol(\"L\", integer=True, nonnegative=True)\n",
     "z = sp.Symbol(\"z\", nonnegative=True, real=True)\n",
-    "ff2 = BlattWeisskopfSquared(z, L)"
+    "bl2 = BlattWeisskopfSquared(z, L)"
    ]
   },
   {
@@ -187,54 +176,90 @@
     "from ampform.io import aslatex\n",
     "\n",
     "ell = sp.Symbol(R\"\\ell\", integer=True, nonnegative=True)\n",
-    "exprs = [\n",
-    "    BlattWeisskopfSquared(z, L),\n",
-    "    SphericalHankel1(ell, z),\n",
-    "]\n",
+    "exprs = [bl2, SphericalHankel1(ell, z)]\n",
     "Math(aslatex({e: e.doit(deep=False) for e in exprs}))"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "editable": true,
+    "jupyter": {
+     "source_hidden": true
+    },
+    "slideshow": {
+     "slide_type": ""
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "%matplotlib inline\n",
+    "\n",
+    "bl2_func = sp.lambdify((z, L), bl2.doit())\n",
+    "x_values = np.linspace(0.0, 5.0, num=500)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.set_xlabel(f\"${sp.latex(z)}$\")\n",
+    "ax.set_ylabel(f\"${sp.latex(bl2)}$\")\n",
+    "ax.set_ylim(0, 10)\n",
+    "\n",
+    "for i in range(5):\n",
+    "    y_values = bl2_func(x_values, L=i)\n",
+    "    ax.plot(x_values, y_values, color=f\"C{i}\", label=f\"$L={i}$\")\n",
+    "ax.legend()\n",
+    "fig.show()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The Blatt-Weisskopf form factor is used to 'dampen' the breakup-momentum at threshold and when going to infinity. A usual choice for $z$ is therefore $z=q^2d^2$ with $q^2$ the {class}`.BreakupMomentumSquared` and $d$ the impact parameter (also called meson radius):"
+    "The Blatt-Weisskopf form factor is used to 'dampen' the breakup-momentum at threshold and when going to infinity. A usual choice for $z$ is therefore $z=q^2d^2$ with $q^2$ the {class}`.BreakupMomentumSquared` and $d$ the impact parameter (also called meson radius). The {class}`.FormFactor` expression class can be used for this:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
-    "tags": []
+    "jupyter": {
+     "source_hidden": true
+    }
    },
    "outputs": [],
    "source": [
-    "from ampform.dynamics import BreakupMomentumSquared\n",
+    "from ampform.dynamics.form_factor import FormFactor\n",
     "\n",
-    "m, m_a, m_b, d = sp.symbols(\"m, m_a, m_b, d\", nonnegative=True)\n",
-    "s = m**2\n",
-    "q_squared = BreakupMomentumSquared(s, m_a, m_b)\n",
-    "ff2 = BlattWeisskopfSquared(q_squared * d**2, angular_momentum=L)"
+    "s, m1, m2, d = sp.symbols(\"s m1 m2 d\", nonnegative=True)\n",
+    "ff2 = FormFactor(s, m1, m2, angular_momentum=L, meson_radius=d)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
+    "editable": true,
     "jupyter": {
      "source_hidden": true
     },
+    "slideshow": {
+     "slide_type": ""
+    },
     "tags": [
-     "hide-cell"
+     "hide-input"
     ]
    },
    "outputs": [],
    "source": [
-    "np_blatt_weisskopf, sliders = symplot.prepare_sliders(\n",
-    "    plot_symbol=m,\n",
-    "    expression=ff2.doit(),\n",
-    ")\n",
-    "np_breakup_momentum = sp.lambdify((m, L, d, m_a, m_b), q_squared.doit())"
+    "from ampform.dynamics.form_factor import BreakupMomentumSquared\n",
+    "\n",
+    "q2 = BreakupMomentumSquared(s, m1, m2)\n",
+    "exprs = [ff2, q2]\n",
+    "Math(aslatex({e: e.doit(deep=False) for e in exprs}))"
    ]
   },
   {
@@ -242,24 +267,12 @@
    "execution_count": null,
    "metadata": {
     "tags": [
-     "remove-cell"
+     "remove-input"
     ]
    },
    "outputs": [],
    "source": [
-    "plot_domain = np.linspace(0.01, 4, 500)\n",
-    "sliders.set_ranges(\n",
-    "    L=(0, 10),\n",
-    "    m_a=(0, 2, 200),\n",
-    "    m_b=(0, 2, 200),\n",
-    "    d=(0, 5),\n",
-    ")\n",
-    "sliders.set_values(\n",
-    "    L=1,\n",
-    "    m_a=0.3,\n",
-    "    m_b=0.2,\n",
-    "    d=3,\n",
-    ")"
+    "%matplotlib widget"
    ]
   },
   {
@@ -270,38 +283,64 @@
      "source_hidden": true
     },
     "tags": [
-     "remove-cell"
+     "hide-input",
+     "remove-output",
+     "scroll-input"
     ]
    },
    "outputs": [],
    "source": [
+    "import ipywidgets as w\n",
+    "\n",
+    "ff2_func = sp.lambdify((s, m1, m2, L, d), ff2.doit())\n",
+    "q2_func = sp.lambdify((s, m1, m2, L, d), q2.doit())\n",
+    "\n",
+    "x = np.linspace(0.01, 4, 500)\n",
+    "sliders = dict(\n",
+    "    m1=w.FloatSlider(description=\"$m_1$\", min=0, max=2, value=0.3),\n",
+    "    m2=w.FloatSlider(description=\"$m_2$\", min=0, max=2, value=0.2),\n",
+    "    L=w.IntSlider(description=\"$L$\", min=0, max=10, value=1),\n",
+    "    d=w.FloatSlider(description=\"$d$\", min=0.1, max=5, value=1),\n",
+    ")\n",
+    "\n",
     "fig, ax = plt.subplots(figsize=(8, 5), tight_layout=True)\n",
     "ax.set_xlabel(\"$m$\")\n",
     "ax.axhline(0, c=\"black\", linewidth=0.5)\n",
-    "controls = iplt.plot(\n",
-    "    plot_domain,\n",
-    "    np_blatt_weisskopf,\n",
-    "    **sliders,\n",
-    "    ylim=\"auto\",\n",
-    "    label=R\"$B_L^2\\left(q(m)\\right)$\",\n",
-    "    linestyle=\"dotted\",\n",
-    "    ax=ax,\n",
-    ")\n",
-    "iplt.plot(\n",
-    "    plot_domain,\n",
-    "    np_breakup_momentum,\n",
-    "    controls=controls,\n",
-    "    ylim=\"auto\",\n",
-    "    label=\"$q^2(m)$\",\n",
-    "    linestyle=\"dashed\",\n",
-    "    ax=ax,\n",
-    ")\n",
-    "iplt.title(\n",
-    "    \"Effect of Blatt-Weisskopf factor\\n$L={L}, m_a={m_a:.1f}, m_b={m_b:.1f}$\",\n",
-    "    controls=controls,\n",
-    ")\n",
-    "plt.legend()\n",
-    "plt.show()"
+    "LINES = None\n",
+    "\n",
+    "\n",
+    "def plot(m1, m2, L, d):\n",
+    "    global LINES\n",
+    "    s = x**2\n",
+    "    y_ff2 = ff2_func(s, m1, m2, L, d)\n",
+    "    y_q2 = q2_func(s, m1, m2, L, d)\n",
+    "    m_thr = m1 + m2\n",
+    "    left = x < m_thr\n",
+    "    right = x > m_thr\n",
+    "    if LINES is None:\n",
+    "        LINES = (\n",
+    "            ax.axvline(m_thr, c=\"black\", ls=\"dotted\", label=\"$m_1+m_2$\"),\n",
+    "            ax.plot(x[left], y_ff2[left], c=\"C0\", ls=\"dashed\")[0],\n",
+    "            ax.plot(x[left], y_q2[left], c=\"C1\", ls=\"dashed\")[0],\n",
+    "            ax.plot(x[right], y_ff2[right], c=\"C0\", label=f\"${sp.latex(ff2)}$\")[0],\n",
+    "            ax.plot(x[right], y_q2[right], c=\"C1\", label=f\"${sp.latex(q2)}$\")[0],\n",
+    "        )\n",
+    "    else:\n",
+    "        LINES[0].set_xdata(m_thr)\n",
+    "        LINES[1].set_data(x[left], y_ff2[left])\n",
+    "        LINES[2].set_data(x[left], y_q2[left])\n",
+    "        LINES[3].set_data(x[right], y_ff2[right])\n",
+    "        LINES[4].set_data(x[right], y_q2[right])\n",
+    "    y_min = np.nanmin(y_q2)\n",
+    "    y_max = np.nanmax(y_q2[right])\n",
+    "    ax.set_ylim(y_min, y_max)\n",
+    "    fig.canvas.draw()\n",
+    "\n",
+    "\n",
+    "UI = w.VBox(list(sliders.values()))\n",
+    "OUTPUT = w.interactive_output(plot, controls=sliders)\n",
+    "ax.legend(loc=\"upper right\")\n",
+    "display(UI, OUTPUT)"
    ]
   },
   {
@@ -321,7 +360,7 @@
     "    from IPython.display import SVG\n",
     "\n",
     "    output_file = \"blatt-weisskopf.svg\"\n",
-    "    plt.savefig(output_file)\n",
+    "    fig.savefig(output_file)\n",
     "    display(SVG(output_file))"
    ]
   },
@@ -381,32 +420,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The relativistic Breit-Wigner can be adapted slightly, so that its amplitude goes to zero at threshold ($m_0 = m_a + m_b$) and that it becomes normalizable. This is done with {ref}`form factors <usage/dynamics:Form factor>` and can be obtained with the function {func}`.relativistic_breit_wigner_with_ff`:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "jupyter": {
-     "source_hidden": true
-    },
-    "tags": [
-     "hide-cell"
-    ]
-   },
-   "outputs": [],
-   "source": [
-    "from ampform.dynamics import (\n",
-    "    BreakupMomentumSquared,\n",
-    "    PhaseSpaceFactor,  # noqa: F401\n",
-    ")\n",
-    "from ampform.sympy.math import ComplexSqrt\n",
-    "\n",
-    "\n",
-    "def breakup_momentum(s: sp.Symbol, m_a: sp.Symbol, m_b: sp.Symbol) -> sp.Expr:\n",
-    "    q_squared = BreakupMomentumSquared(s, m_a, m_b)\n",
-    "    return ComplexSqrt(q_squared)"
+    "The relativistic Breit-Wigner can be adapted slightly, so that its amplitude goes to zero at threshold ($m_0 = m1 + m2$) and that it becomes normalizable. This is done with {ref}`form factors <usage/dynamics:Form factor>` and can be obtained with the function {func}`.relativistic_breit_wigner_with_ff`:"
    ]
   },
   {
@@ -423,8 +437,8 @@
     "    s=s,\n",
     "    mass0=m0,\n",
     "    gamma0=w0,\n",
-    "    m_a=m_a,\n",
-    "    m_b=m_b,\n",
+    "    m_a=m1,\n",
+    "    m_b=m2,\n",
     "    angular_momentum=L,\n",
     "    meson_radius=1,\n",
     "    phsp_factor=PhaseSpaceFactorSWave,\n",
@@ -443,6 +457,9 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
     "tags": []
    },
    "outputs": [],
@@ -454,8 +471,8 @@
     "    s=s,\n",
     "    mass0=m0,\n",
     "    gamma0=w0,\n",
-    "    m_a=m_a,\n",
-    "    m_b=m_b,\n",
+    "    m_a=m1,\n",
+    "    m_b=m2,\n",
     "    angular_momentum=L,\n",
     "    meson_radius=1,\n",
     "    phsp_factor=PhaseSpaceFactorSWave,\n",
@@ -504,21 +521,21 @@
     "\n",
     "# Two types of relativistic Breit-Wigners\n",
     "rel_bw_with_ff = relativistic_breit_wigner_with_ff(\n",
-    "    s=s,\n",
+    "    s=m**2,\n",
     "    mass0=m0,\n",
     "    gamma0=w0,\n",
-    "    m_a=m_a,\n",
-    "    m_b=m_b,\n",
+    "    m_a=m1,\n",
+    "    m_b=m2,\n",
     "    angular_momentum=L,\n",
     "    meson_radius=d,\n",
     "    phsp_factor=PhaseSpaceFactorComplex,\n",
     ")\n",
     "rel_bw_with_ff_ac = relativistic_breit_wigner_with_ff(\n",
-    "    s=s,\n",
+    "    s=m**2,\n",
     "    mass0=m0,\n",
     "    gamma0=w0,\n",
-    "    m_a=m_a,\n",
-    "    m_b=m_b,\n",
+    "    m_a=m1,\n",
+    "    m_b=m2,\n",
     "    angular_momentum=L,\n",
     "    meson_radius=d,\n",
     "    phsp_factor=PhaseSpaceFactorSWave,\n",
@@ -530,30 +547,30 @@
     "    expression=rel_bw_with_ff.doit(),\n",
     ")\n",
     "np_rel_bw_with_ff_ac = sp.lambdify(\n",
-    "    args=(m, w0, L, d, m0, m_a, m_b),\n",
+    "    args=(m, w0, L, d, m0, m1, m2),\n",
     "    expr=rel_bw_with_ff_ac.doit(),\n",
     ")\n",
     "np_rel_bw = sp.lambdify(\n",
-    "    args=(m, w0, L, d, m0, m_a, m_b),\n",
+    "    args=(m, w0, L, d, m0, m1, m2),\n",
     "    expr=rel_bw.doit(),\n",
     ")\n",
     "\n",
     "# Set sliders\n",
-    "plot_domain = np.linspace(start=0, stop=4, num=500)\n",
+    "plot_domain = np.linspace(0, 4, num=500)\n",
     "sliders.set_ranges(\n",
     "    m0=(0, 4, 200),\n",
     "    Gamma0=(0, 1, 100),\n",
     "    L=(0, 10),\n",
-    "    m_a=(0, 2, 200),\n",
-    "    m_b=(0, 2, 200),\n",
+    "    m1=(0, 2, 200),\n",
+    "    m2=(0, 2, 200),\n",
     "    d=(0, 5),\n",
     ")\n",
     "sliders.set_values(\n",
-    "    m0=1.6,\n",
+    "    m0=1.5,\n",
     "    Gamma0=0.6,\n",
     "    L=0,\n",
-    "    m_a=0.6,\n",
-    "    m_b=0.6,\n",
+    "    m1=0.6,\n",
+    "    m2=0.6,\n",
     "    d=1,\n",
     ")\n",
     "\n",
@@ -567,7 +584,7 @@
     "for ax in axes:\n",
     "    ax.axhline(0, c=\"gray\", linewidth=0.5)\n",
     "\n",
-    "rho_c = PhaseSpaceFactorComplex(m**2, m_a, m_b)\n",
+    "rho_c = PhaseSpaceFactorComplex(m**2, m1, m2)\n",
     "ax_ff.set_title(f\"'Complex' phase space factor: ${sp.latex(rho_c.evaluate())}$\")\n",
     "ax_ac.set_title(\"$S$-wave Chew-Mandelstam as phase space factor\")\n",
     "\n",
@@ -581,7 +598,7 @@
     "    ax=ax_ff,\n",
     ")\n",
     "iplt.plot(\n",
-    "    plot_domain.real,\n",
+    "    plot_domain,\n",
     "    lambda *args, **kwargs: np_rel_bw_with_ff(*args, **kwargs).imag,\n",
     "    label=\"imaginary\",\n",
     "    controls=controls,\n",
@@ -589,7 +606,7 @@
     "    ax=ax_ff,\n",
     ")\n",
     "iplt.plot(\n",
-    "    plot_domain.real,\n",
+    "    plot_domain,\n",
     "    lambda *args, **kwargs: np.abs(np_rel_bw_with_ff(*args, **kwargs)) ** 2,\n",
     "    label=\"absolute\",\n",
     "    controls=controls,\n",
@@ -601,7 +618,7 @@
     "\n",
     "\n",
     "iplt.plot(\n",
-    "    plot_domain.real,\n",
+    "    plot_domain,\n",
     "    lambda *args, **kwargs: np_rel_bw_with_ff_ac(*args, **kwargs).real,\n",
     "    label=\"real\",\n",
     "    controls=controls,\n",
@@ -609,7 +626,7 @@
     "    ax=ax_ac,\n",
     ")\n",
     "iplt.plot(\n",
-    "    plot_domain.real,\n",
+    "    plot_domain,\n",
     "    lambda *args, **kwargs: np_rel_bw_with_ff_ac(*args, **kwargs).imag,\n",
     "    label=\"imaginary\",\n",
     "    controls=controls,\n",
@@ -617,7 +634,7 @@
     "    ax=ax_ac,\n",
     ")\n",
     "iplt.plot(\n",
-    "    plot_domain.real,\n",
+    "    plot_domain,\n",
     "    lambda *args, **kwargs: np.abs(np_rel_bw_with_ff_ac(*args, **kwargs)) ** 2,\n",
     "    label=\"absolute\",\n",
     "    controls=controls,\n",
@@ -636,12 +653,12 @@
     "        alpha=0.3,\n",
     "    )\n",
     "    iplt.axvline(\n",
-    "        lambda m_a, m_b, **kwargs: m_a + m_b,\n",
+    "        lambda m1, m2, **kwargs: m1 + m2,\n",
     "        controls=controls,\n",
     "        ax=ax,\n",
     "        c=\"black\",\n",
     "        alpha=0.3,\n",
-    "        label=f\"${sp.latex(m_a)} + {sp.latex(m_b)}$\",\n",
+    "        label=f\"${sp.latex(m1)} + {sp.latex(m2)}$\",\n",
     "    )\n",
     "ax_ac.legend(loc=\"upper right\")\n",
     "fig.tight_layout()\n",
@@ -666,7 +683,7 @@
     "    from IPython.display import SVG\n",
     "\n",
     "    output_file = \"relativistic-breit-wigner-with-form-factor.svg\"\n",
-    "    plt.savefig(output_file)\n",
+    "    fig.savefig(output_file)\n",
     "    display(SVG(output_file))"
    ]
   }
@@ -690,7 +707,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.8"
+   "version": "3.11.9"
   }
  },
  "nbformat": 4,
diff --git a/docs/usage/dynamics/k-matrix.ipynb b/docs/usage/dynamics/k-matrix.ipynb
index 0492224d0..b401194e2 100644
--- a/docs/usage/dynamics/k-matrix.ipynb
+++ b/docs/usage/dynamics/k-matrix.ipynb
@@ -128,19 +128,12 @@
     "from mpl_interactions.controller import Controls\n",
     "\n",
     "import symplot\n",
-    "from ampform.dynamics import (\n",
-    "    BlattWeisskopfSquared,\n",
-    "    PhaseSpaceFactor,\n",
-    "    kmatrix,\n",
-    "    relativistic_breit_wigner,\n",
-    ")\n",
+    "from ampform.dynamics import PhaseSpaceFactor, kmatrix, relativistic_breit_wigner\n",
     "\n",
     "logging.basicConfig()\n",
     "logging.getLogger().setLevel(logging.ERROR)\n",
     "\n",
-    "warnings.filterwarnings(\"ignore\")\n",
-    "\n",
-    "BlattWeisskopfSquared.max_angular_momentum = 3"
+    "warnings.filterwarnings(\"ignore\")"
    ]
   },
   {
@@ -536,7 +529,7 @@
     "\\end{eqnarray}\n",
     "$$ (P-vector parametrization)\n",
     "\n",
-    "with $B_{R,i}(q(s))$ the **centrifugal damping factor** (see {class}`.BlattWeisskopfSquared`) for channel $i$ and $\\beta_R^0$ some (generally complex) constants that describe the production information of the decaying state $R$. Usually, these constants are rescaled just like the residue functions in {eq}`residue-function`:\n",
+    "with $B_{R,i}(q(s))$ the **centrifugal damping factor** (see {class}`.FormFactor` and {class}`.BlattWeisskopfSquared`) for channel $i$ and $\\beta_R^0$ some (generally complex) constants that describe the production information of the decaying state $R$. Usually, these constants are rescaled just like the residue functions in {eq}`residue-function`:\n",
     "\n",
     "[^damping-factor-P-parametrization]: Just as with [^phase-space-factor-normalization], we have smuggled a bit in the last equation in order to be able to reproduce Equation (50.23) in {pdg-review}`2021; Resonances; p.9` in the case $n=1,n_R=1$, on which {func}`.relativistic_breit_wigner_with_ff` is based.\n",
     "\n",
@@ -1513,11 +1506,8 @@
     "    # Set slider values and ranges\n",
     "    m0_values = np.linspace(x_min, x_max, num=n_poles + 2)\n",
     "    m0_values = m0_values[1:-1]\n",
-    "    max_angular_momentum = BlattWeisskopfSquared.max_angular_momentum\n",
-    "    if max_angular_momentum is None:\n",
-    "        max_angular_momentum = 8\n",
     "    if \"L\" in sliders:\n",
-    "        sliders.set_ranges(L=(0, max_angular_momentum))\n",
+    "        sliders.set_ranges(L=(0, 10))\n",
     "    sliders.set_ranges(i=(0, n_channels - 1))\n",
     "    for R in range(1, n_poles + 1):\n",
     "        for i in range(n_channels):\n",
@@ -1843,7 +1833,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.8"
+   "version": "3.11.9"
   }
  },
  "nbformat": 4,
diff --git a/pyproject.toml b/pyproject.toml
index e6116b9ef..247afa9a9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -345,6 +345,7 @@ split-on-trailing-comma = false
     "E402",
     "E703",
     "F404",
+    "N803",
     "N806",
     "N816",
     "PLC2401",
diff --git a/src/ampform/dynamics/__init__.py b/src/ampform/dynamics/__init__.py
index fcc3d40ea..671c91bcb 100644
--- a/src/ampform/dynamics/__init__.py
+++ b/src/ampform/dynamics/__init__.py
@@ -6,13 +6,17 @@
 from __future__ import annotations
 
 from typing import TYPE_CHECKING, Any
+from warnings import warn
 
 import sympy as sp
 
 # pyright: reportUnusedImport=false
-from ampform.dynamics.form_factor import BlattWeisskopfSquared
+from ampform.dynamics.form_factor import (
+    BlattWeisskopfSquared,  # noqa: F401
+    FormFactor,
+)
 from ampform.dynamics.phasespace import (
-    BreakupMomentumSquared,
+    BreakupMomentumSquared,  # noqa: F401
     EqualMassPhaseSpaceFactor,  # noqa: F401
     PhaseSpaceFactor,
     PhaseSpaceFactorAbs,  # noqa: F401
@@ -36,10 +40,10 @@ class EnergyDependentWidth(sp.Expr):
     :cite:`ParticleDataGroup:2020ssz`, equation (6). Default value for
     :code:`phsp_factor` is `.PhaseSpaceFactor`.
 
-    Note that the `.BlattWeisskopfSquared` of AmpForm is normalized in the sense that
-    equal powers of :math:`z` appear in the nominator and the denominator, while the
-    definition in the PDG (as well as some other sources), always have :math:`1` in the
-    nominator of the Blatt-Weisskopf. In that case, one needs an additional factor
+    Note that the `.FormFactor` of AmpForm is normalized in the sense that equal powers
+    of :math:`z` appear in the nominator and the denominator, while the definition in
+    the PDG (as well as some other sources), always have :math:`1` in the nominator of
+    the Blatt-Weisskopf. In that case, one needs an additional factor
     :math:`\left(q/q_0\right)^{2L}` in the definition for :math:`\Gamma(m)`.
     """
 
@@ -56,25 +60,17 @@ class EnergyDependentWidth(sp.Expr):
     name: str | None = argument(default=None, sympify=False)
 
     def evaluate(self) -> sp.Expr:
-        s, mass0, gamma0, m_a, m_b, angular_momentum, meson_radius = self.args
-        q_squared = BreakupMomentumSquared(s, m_a, m_b)
-        q0_squared = BreakupMomentumSquared(mass0**2, m_a, m_b)  # type: ignore[operator]
-        form_factor_sq = BlattWeisskopfSquared(
-            q_squared * meson_radius**2,  # type: ignore[operator]
-            angular_momentum,
-        )
-        form_factor0_sq = BlattWeisskopfSquared(
-            q0_squared * meson_radius**2,  # type: ignore[operator]
-            angular_momentum,
-        )
-        rho = self.phsp_factor(s, m_a, m_b)
-        rho0 = self.phsp_factor(mass0**2, m_a, m_b)  # type: ignore[operator]
-        return gamma0 * (form_factor_sq / form_factor0_sq) * (rho / rho0)
+        s, m0, width0, m1, m2, angular_momentum, meson_radius = self.args
+        ff2 = FormFactor(s, m1, m2, angular_momentum, meson_radius)
+        ff2_0 = FormFactor(m0**2, m1, m2, angular_momentum, meson_radius)  # type: ignore[operator]
+        rho = self.phsp_factor(s, m1, m2)
+        rho0 = self.phsp_factor(m0**2, m1, m2)  # type: ignore[operator]
+        return width0 * (ff2 / ff2_0) * (rho / rho0)
 
     def _latex_repr_(self, printer: LatexPrinter, *args) -> str:
         s = printer._print(self.args[0])
-        gamma0 = self.args[2]
-        subscript = _indices_to_subscript(determine_indices(gamma0))
+        width0 = self.args[2]
+        subscript = _indices_to_subscript(determine_indices(width0))
         name = Rf"\Gamma{subscript}" if self.name is None else self.name
         return Rf"{name}\left({s}\right)"
 
@@ -98,12 +94,12 @@ def relativistic_breit_wigner_with_ff(  # noqa: PLR0917
     meson_radius,
     phsp_factor: PhaseSpaceFactorProtocol = PhaseSpaceFactor,
 ) -> sp.Expr:
-    """Relativistic Breit-Wigner with `.BlattWeisskopfSquared` factor.
+    """Relativistic Breit-Wigner with `.FormFactor`.
 
     See :ref:`usage/dynamics:_With_ form factor` and :pdg-review:`2021; Resonances;
     p.9`.
     """
-    form_factor = formulate_form_factor(s, m_a, m_b, angular_momentum, meson_radius)
+    form_factor = FormFactor(s, m_a, m_b, angular_momentum, meson_radius)
     energy_dependent_width = EnergyDependentWidth(
         s, mass0, gamma0, m_a, m_b, angular_momentum, meson_radius, phsp_factor
     )
@@ -115,10 +111,11 @@ def relativistic_breit_wigner_with_ff(  # noqa: PLR0917
 def formulate_form_factor(s, m_a, m_b, angular_momentum, meson_radius) -> sp.Expr:
     """Formulate a Blatt-Weisskopf form factor.
 
-    Returns the production process factor :math:`n_a` from Equation (50.26) in
-    :pdg-review:`2021; Resonances; p.9`, which features the
-    `~sympy.functions.elementary.miscellaneous.sqrt` of a `.BlattWeisskopfSquared`.
+    .. deprecated:: 0.16
     """
-    q_squared = BreakupMomentumSquared(s, m_a, m_b)
-    ff_squared = BlattWeisskopfSquared(q_squared * meson_radius**2, angular_momentum)
-    return sp.sqrt(ff_squared)
+    warn(
+        message="Use the FormFactor expression class instead.",
+        category=DeprecationWarning,
+        stacklevel=1,
+    )
+    return FormFactor(s, m_a, m_b, angular_momentum, meson_radius)
diff --git a/src/ampform/dynamics/builder.py b/src/ampform/dynamics/builder.py
index 775823952..7e717e661 100644
--- a/src/ampform/dynamics/builder.py
+++ b/src/ampform/dynamics/builder.py
@@ -9,11 +9,7 @@
 from attrs import field, frozen
 from attrs.validators import instance_of
 
-from ampform.dynamics import (
-    EnergyDependentWidth,
-    formulate_form_factor,
-    relativistic_breit_wigner,
-)
+from ampform.dynamics import EnergyDependentWidth, FormFactor, relativistic_breit_wigner
 from ampform.dynamics.phasespace import (
     EqualMassPhaseSpaceFactor,
     PhaseSpaceFactor,
@@ -83,17 +79,17 @@ def create_non_dynamic_with_ff(
 ) -> BuilderReturnType:
     """Generate (only) a Blatt-Weisskopf form factor for a two-body decay.
 
-    See also :func:`.formulate_form_factor`.
+    See also :class:`.FormFactor`.
     """
     if variable_pool.angular_momentum is None:
         msg = "Angular momentum is not defined but is required in the form factor!"
         raise ValueError(msg)
     res_identifier = resonance.latex if resonance.latex else resonance.name
     meson_radius = sp.Symbol(f"d_{{{res_identifier}}}", positive=True)
-    form_factor = formulate_form_factor(
+    form_factor = FormFactor(
         s=variable_pool.incoming_state_mass**2,
-        m_a=variable_pool.outgoing_state_mass1,
-        m_b=variable_pool.outgoing_state_mass2,
+        m1=variable_pool.outgoing_state_mass1,
+        m2=variable_pool.outgoing_state_mass2,
         angular_momentum=variable_pool.angular_momentum,
         meson_radius=meson_radius,
     )
@@ -111,8 +107,8 @@ class RelativisticBreitWignerBuilder:
 
     Args:
         form_factor: Formulate a relativistic Breit-Wigner function multiplied
-            by a Blatt-Weisskopf form factor (`.BlattWeisskopfSquared`), like in
-            Equation (50.26) on :pdg-review:`2021; Resonances; p.9`.
+            by a Blatt-Weisskopf form factor (`.FormFactor`), like in Equation (50.26)
+            on :pdg-review:`2021; Resonances; p.9`.
         energy_dependent_width: Use an `.EnergyDependentWidth` in the
             denominator of the Breit-Wigner.
         phsp_factor: A class that complies with the
@@ -214,10 +210,10 @@ def __create_form_factor(
 
         inv_mass = variable_pool.incoming_state_mass
         _, __, meson_radius = self.__create_symbols(resonance)
-        form_factor = formulate_form_factor(
+        form_factor = FormFactor(
             s=inv_mass**2,
-            m_a=variable_pool.outgoing_state_mass1,
-            m_b=variable_pool.outgoing_state_mass2,
+            m1=variable_pool.outgoing_state_mass1,
+            m2=variable_pool.outgoing_state_mass2,
             angular_momentum=variable_pool.angular_momentum,
             meson_radius=meson_radius,
         )
diff --git a/src/ampform/dynamics/form_factor.py b/src/ampform/dynamics/form_factor.py
index e54b0b6ce..3b5b5edcf 100644
--- a/src/ampform/dynamics/form_factor.py
+++ b/src/ampform/dynamics/form_factor.py
@@ -7,9 +7,34 @@
 
 import sympy as sp
 
+from ampform.dynamics.phasespace import BreakupMomentumSquared
 from ampform.sympy import unevaluated
 
 
+@unevaluated
+class FormFactor(sp.Expr):
+    """Formulate a Blatt-Weisskopf form factor.
+
+    Returns the production process factor :math:`n_a` from Equation (50.26) in
+    :pdg-review:`2021; Resonances; p.9`, which features the
+    `~sympy.functions.elementary.miscellaneous.sqrt` of a `.BlattWeisskopfSquared`.
+    """
+
+    s: Any
+    m1: Any
+    m2: Any
+    angular_momentum: Any
+    meson_radius: Any = 1
+
+    _latex_repr_ = R"\mathcal{{F}}_{{{angular_momentum}}}\left({s}, {m1}, {m2}\right)"
+
+    def evaluate(self):
+        s, m1, m2, angular_momentum, meson_radius = self.args
+        q2 = BreakupMomentumSquared(s, m1, m2)
+        ff_squared = BlattWeisskopfSquared(q2 * meson_radius**2, angular_momentum)
+        return sp.sqrt(ff_squared)
+
+
 @unevaluated
 class BlattWeisskopfSquared(sp.Expr):
     r"""Normalized Blatt-Weisskopf function :math:`B_L^2(z)`, with :math:`B_L^2(1)=1`.
diff --git a/src/ampform/dynamics/kmatrix.py b/src/ampform/dynamics/kmatrix.py
index defa8071a..da2e3bf63 100644
--- a/src/ampform/dynamics/kmatrix.py
+++ b/src/ampform/dynamics/kmatrix.py
@@ -19,8 +19,8 @@
     EnergyDependentWidth,
     PhaseSpaceFactor,
     PhaseSpaceFactorProtocol,
-    formulate_form_factor,
 )
+from ampform.dynamics.form_factor import FormFactor
 from ampform.sympy import create_symbol_matrix
 
 
@@ -384,9 +384,7 @@ def parametrization(  # noqa: PLR0917
         gamma = residue_constant[pole_id, i]
         mass0 = pole_position[pole_id]
         width = pole_width[pole_id, i]
-        form_factor = formulate_form_factor(
-            s, m_a[i], m_b[i], angular_momentum, meson_radius
-        )
+        form_factor = FormFactor(s, m_a[i], m_b[i], angular_momentum, meson_radius)
         return sp.Sum(
             beta * gamma * mass0 * width * form_factor / (mass0**2 - s),
             (pole_id, 1, n_poles),
diff --git a/tests/dynamics/test_builder.py b/tests/dynamics/test_builder.py
index a06338df8..a12a4f732 100644
--- a/tests/dynamics/test_builder.py
+++ b/tests/dynamics/test_builder.py
@@ -2,11 +2,12 @@
 import sympy as sp
 from qrules.particle import Particle
 
-from ampform.dynamics import EnergyDependentWidth, formulate_form_factor
+from ampform.dynamics import EnergyDependentWidth
 from ampform.dynamics.builder import (
     RelativisticBreitWignerBuilder,
     TwoBodyKinematicVariableSet,
 )
+from ampform.dynamics.form_factor import FormFactor
 
 
 class TestRelativisticBreitWignerBuilder:
@@ -53,9 +54,7 @@ def test_simple_breit_wigner(
         m2 = variable_set.outgoing_state_mass2
         L = variable_set.angular_momentum  # noqa: N806
         d = sp.Symbol(R"d_{N}", positive=True)
-        form_factor = formulate_form_factor(
-            s, m1, m2, angular_momentum=L, meson_radius=d
-        )
+        form_factor = FormFactor(s, m1, m2, angular_momentum=L, meson_radius=d)
         assert bw_with_ff / bw == form_factor
         assert set(parameters) == {m0, w0, d}
         assert parameters[m0] == particle.mass
@@ -89,9 +88,7 @@ def test_breit_wigner_with_energy_dependent_width(
         builder.form_factor = True
         bw_with_ff, parameters = builder(particle, variable_set)
         ang_mom = variable_set.angular_momentum
-        form_factor = formulate_form_factor(
-            s, m1, m2, angular_momentum=ang_mom, meson_radius=d
-        )
+        form_factor = FormFactor(s, m1, m2, angular_momentum=ang_mom, meson_radius=d)
         assert bw_with_ff / bw == form_factor
         assert set(parameters) == {m0, w0, d}
         assert parameters[m0] == particle.mass
diff --git a/tests/sympy/test_cache.py b/tests/sympy/test_cache.py
index 41fd50906..a99ba42c0 100644
--- a/tests/sympy/test_cache.py
+++ b/tests/sympy/test_cache.py
@@ -111,17 +111,17 @@ def test_get_readable_hash_large(amplitude_model: tuple[str, HelicityModel]):
         # https://github.com/ComPWA/ampform/actions/runs/3277058875/jobs/5393849802
         # https://github.com/ComPWA/ampform/actions/runs/3277143883/jobs/5394043014
         expected_hash = {
-            "canonical-helicity": "pythonhashseed-0-5383158178341674913",
-            "helicity": "pythonhashseed-0-7110179181475816071",
+            "canonical-helicity": "pythonhashseed-0-4409019767276782833",
+            "helicity": "pythonhashseed-0+8495836064961054249",
         }[formalism]
     elif sys.version_info >= (3, 11):
         expected_hash = {
-            "canonical-helicity": "pythonhashseed-0-9019524491135114657",
-            "helicity": "pythonhashseed-0+220334524130979941",
+            "canonical-helicity": "pythonhashseed-0-8140852268928771574",
+            "helicity": "pythonhashseed-0-991855900379383849",
         }[formalism]
     else:
         expected_hash = {
-            "canonical-helicity": "pythonhashseed-0-8505502895987205495",
-            "helicity": "pythonhashseed-0-1430245260241162669",
+            "canonical-helicity": "pythonhashseed-0+3166036244969111461",
+            "helicity": "pythonhashseed-0+4247688887304834148",
         }[formalism]
     assert get_readable_hash(model.expression) == expected_hash