diff --git a/docs/ampform/LambdaKpi0.ipynb b/docs/ampform/LambdaKpi0.ipynb
new file mode 100644
index 0000000..68a05c2
--- /dev/null
+++ b/docs/ampform/LambdaKpi0.ipynb
@@ -0,0 +1,612 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# pγ → ΛK⁺π⁰ amplitude model"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "PWA study on $p \\gamma \\to \\Lambda K^+ \\pi^0$.\n",
+    "We formulate the helicity amplitude model symbolically using `AmpForm` here."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "mystnb": {
+     "code_prompt_show": "Hide warnings"
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from __future__ import annotations\n",
+    "\n",
+    "import logging\n",
+    "import os\n",
+    "import warnings\n",
+    "\n",
+    "os.environ[\"TF_CPP_MIN_LOG_LEVEL\"] = \"3\"\n",
+    "logging.disable(logging.WARNING)\n",
+    "warnings.filterwarnings(\"ignore\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "mystnb": {
+     "code_prompt_show": "Import Python Libraries"
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from collections import defaultdict\n",
+    "\n",
+    "import ampform\n",
+    "import graphviz\n",
+    "import jax\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import qrules\n",
+    "from ampform.dynamics.builder import RelativisticBreitWignerBuilder\n",
+    "from ampform.io import aslatex, improve_latex_rendering\n",
+    "from IPython.display import Math\n",
+    "from qrules.particle import Particle, Spin, create_particle, load_pdg\n",
+    "from tensorwaves.data import (\n",
+    "    SympyDataTransformer,\n",
+    "    TFPhaseSpaceGenerator,\n",
+    "    TFUniformRealNumberGenerator,\n",
+    ")\n",
+    "from tensorwaves.function.sympy import create_parametrized_function\n",
+    "\n",
+    "improve_latex_rendering()\n",
+    "particle_db = load_pdg()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Decay definition"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Particle definitions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "Lambda"
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "particle_db[\"Lambda\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "K+"
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "particle_db[\"K+\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "pi0"
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "particle_db[\"pi0\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "gamma"
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "particle_db[\"gamma\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "p"
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "particle_db[\"p\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "| Particle         | Lambda            | K+               | pi0              | gamma            | p                |\n",
+    "|------------------|-------------------|------------------|------------------|------------------|------------------|\n",
+    "| **Name**         | Lambda            | K+               | pi0              | gamma            | p                |\n",
+    "| **PID**          | 3122              | 321              | 111              | 22               | 2212             |\n",
+    "| **Latex**        | \\(\\Lambda\\)       | \\(K^{+}\\)        | \\(\\pi^{0}\\)      | \\(\\gamma\\)       | \\(p\\)            |\n",
+    "| **Spin**         | 0.5               | 0.0              | 0.0              | 1.0              | 0.5              |\n",
+    "| **Mass**         | 1.115683          | 0.49367700000000003 | 0.1349768       | 0.0              | 0.93827208816    |\n",
+    "| **Width**        | 2.501e-15         | 5.317e-17        | 7.81e-09         |                  |                  |\n",
+    "| **Charge**       |                   | 1                |                  |                  | 1                |\n",
+    "| **Isospin**      | Spin(0, 0)        | Spin(1/2, +1/2)  | Spin(1, 0)       |                  | Spin(1/2, +1/2)  |\n",
+    "| **Strangeness**  | -1                | 1                |                  |                  |                  |\n",
+    "| **Baryon Number**| 1                 |                  |                  |                  | 1                |\n",
+    "| **Parity**       | +1                | -1               | -1               | -1               | +1               |\n",
+    "| **C Parity**     |                   |                  | +1               | -1               |                  |\n",
+    "| **G Parity**     |                   |                  | -1               |                  |                  |\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Initial state definition"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Mass for $p \\gamma$ system"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "E_lab_gamma = 8.5\n",
+    "m_proton = 0.938\n",
+    "m_0 = np.sqrt(2 * E_lab_gamma * m_proton + m_proton**2)\n",
+    "m_eta = 0.548\n",
+    "m_pi = 0.135\n",
+    "m_0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Add custom particle $p \\gamma$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pgamma1 = Particle(\n",
+    "    name=\"pgamma1\",\n",
+    "    latex=r\"p\\gamma (s1/2)\",\n",
+    "    spin=0.5,\n",
+    "    mass=m_0,\n",
+    "    charge=1,\n",
+    "    isospin=Spin(1 / 2, +1 / 2),\n",
+    "    baryon_number=1,\n",
+    "    parity=-1,\n",
+    "    pid=99990,\n",
+    ")\n",
+    "pgamma2 = create_particle(\n",
+    "    template_particle=pgamma1,\n",
+    "    name=\"pgamma2\",\n",
+    "    latex=R\"p\\gamma (s3/2)\",\n",
+    "    spin=1.5,\n",
+    "    pid=pgamma1.pid + 1,\n",
+    ")\n",
+    "particle_db.update([pgamma1, pgamma2])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Generate transitions"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For simplicity, we use the initial state $p \\gamma$ (with spin-$\\frac{1}{2}$), \n",
+    "and set the allowed interaction type to be strong only,\n",
+    "the formalism is selected to be helicity formalism instead of canonical."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    ":::{seealso}\n",
+    "[Helicity versus canonical](https://ampform.readthedocs.io/stable/usage/helicity/formalism.html)\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "reaction = qrules.generate_transitions(\n",
+    "    initial_state=(\"pgamma1\"),\n",
+    "    final_state=[\"Lambda\", \"K+\", \"pi0\"],\n",
+    "    allowed_interaction_types=[\"strong\"],\n",
+    "    formalism=\"helicity\",\n",
+    "    particle_db=particle_db,\n",
+    "    max_angular_momentum=4,\n",
+    "    max_spin_magnitude=4,\n",
+    "    mass_conservation_factor=0,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "dot = qrules.io.asdot(reaction, collapse_graphs=True)\n",
+    "graphviz.Source(dot)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Formulate amplitude model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_builder = ampform.get_builder(reaction)\n",
+    "model_builder.config.scalar_initial_state_mass = True\n",
+    "model_builder.config.stable_final_state_ids = 0, 1, 2\n",
+    "bw_builder = RelativisticBreitWignerBuilder(\n",
+    "    energy_dependent_width=False,\n",
+    "    form_factor=False,\n",
+    ")\n",
+    "for name in reaction.get_intermediate_particles().names:\n",
+    "    model_builder.dynamics.assign(name, bw_builder)\n",
+    "model = model_builder.formulate()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model.intensity"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The first term in the amplitude model:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input",
+     "scroll-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "(symbol, expr), *_ = model.amplitudes.items()\n",
+    "Math(aslatex({symbol: expr}, terms_per_line=1))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "Model parameters"
+    },
+    "tags": [
+     "scroll-output",
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "Math(aslatex(model.parameter_defaults))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "mystnb": {
+     "code_prompt_show": "Kinematic variable definitions"
+    },
+    "tags": [
+     "hide-input",
+     "scroll-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "Math(aslatex(model.kinematic_variables))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Visualization"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "unfolded_expression = model.expression.doit()\n",
+    "intensity_func = create_parametrized_function(\n",
+    "    expression=unfolded_expression,\n",
+    "    parameters=model.parameter_defaults,\n",
+    "    backend=\"jax\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "phsp_event = 500_000\n",
+    "rng = TFUniformRealNumberGenerator(seed=0)\n",
+    "phsp_generator = TFPhaseSpaceGenerator(\n",
+    "    initial_state_mass=reaction.initial_state[-1].mass,\n",
+    "    final_state_masses={i: p.mass for i, p in reaction.final_state.items()},\n",
+    ")\n",
+    "phsp_momenta = phsp_generator.generate(phsp_event, rng)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "helicity_transformer = SympyDataTransformer.from_sympy(\n",
+    "    model.kinematic_variables,\n",
+    "    backend=\"jax\",\n",
+    ")\n",
+    "phsp = helicity_transformer(phsp_momenta)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%config InlineBackend.figure_formats = ['png']\n",
+    "fig, ax = plt.subplots(dpi=200)\n",
+    "hist = ax.hist2d(\n",
+    "    phsp[\"m_01\"].real ** 2,\n",
+    "    phsp[\"m_12\"].real ** 2,\n",
+    "    bins=200,\n",
+    "    cmin=1e-6,\n",
+    "    density=True,\n",
+    "    cmap=\"jet\",\n",
+    "    vmax=0.15,\n",
+    "    weights=intensity_func(phsp),\n",
+    ")\n",
+    "ax.set_title(\"Model-weighted Phase space Dalitz Plot\")\n",
+    "ax.set_xlabel(R\"$m^2(\\Lambda K^+)\\;\\left[\\mathrm{GeV}^2\\right]$\")\n",
+    "ax.set_ylabel(R\"$m^2(K^+ \\pi^0)\\;\\left[\\mathrm{GeV}^2\\right]$\")\n",
+    "cbar = fig.colorbar(hist[3], ax=ax)\n",
+    "fig.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-output",
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "resonances = defaultdict(set)\n",
+    "for transition in reaction.transitions:\n",
+    "    topology = transition.topology\n",
+    "    top_decay_products = topology.get_edge_ids_outgoing_from_node(0)\n",
+    "    (resonance_id, resonance), *_ = transition.intermediate_states.items()\n",
+    "    recoil_id, *_ = top_decay_products - {resonance_id}\n",
+    "    resonances[recoil_id].add(resonance.particle)\n",
+    "resonances = {k: sorted(v, key=lambda p: p.mass) for k, v in resonances.items()}\n",
+    "{k: [p.name for p in v] for k, v in resonances.items()}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input",
+     "full-width"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "fig, axes = plt.subplots(figsize=(14, 4), ncols=3, sharey=True)\n",
+    "ax1, ax2, ax3 = axes\n",
+    "\n",
+    "max_value = 0\n",
+    "color_id = 0\n",
+    "intensities = intensity_func(phsp)\n",
+    "for recoil_id, ax in enumerate(axes):\n",
+    "    decay_products = sorted({0, 1, 2} - {recoil_id})\n",
+    "    key = f\"m_{''.join(str(i) for i in decay_products)}\"\n",
+    "    bin_values, bin_edges = jax.numpy.histogram(\n",
+    "        phsp[key].real,\n",
+    "        bins=100,\n",
+    "        density=True,\n",
+    "        weights=intensities,\n",
+    "    )\n",
+    "    max_value = max(max_value, bin_values.max())\n",
+    "    ax.fill_between(\n",
+    "        bin_edges[:-1],\n",
+    "        bin_values,\n",
+    "        alpha=0.5,\n",
+    "        step=\"pre\",\n",
+    "    )\n",
+    "    ax.set_ylim(0, None)\n",
+    "    for resonance in resonances[recoil_id]:\n",
+    "        ax.axvline(\n",
+    "            resonance.mass,\n",
+    "            c=f\"C{color_id}\",\n",
+    "            ls=\"dotted\",\n",
+    "            label=resonance.name,\n",
+    "        )\n",
+    "        color_id += 1\n",
+    "\n",
+    "        product_latex = \" \".join(\n",
+    "            [reaction.final_state[i].latex for i in decay_products],\n",
+    "        )\n",
+    "        ax.set_xlabel(f\"$m({product_latex})$ [GeV]\")\n",
+    "\n",
+    "        if recoil_id == 0:\n",
+    "            ax.legend(fontsize=\"small\", loc=\"upper right\")\n",
+    "        elif recoil_id == 1:\n",
+    "            ax.legend(fontsize=\"small\", loc=\"upper center\", bbox_to_anchor=(0.6, 1.0))\n",
+    "        else:\n",
+    "            ax.legend(fontsize=\"small\", loc=\"upper center\")\n",
+    "for ax in axes:\n",
+    "    ax.set_ylim(0, max_value * 1.1)\n",
+    "\n",
+    "fig.tight_layout()\n",
+    "plt.show(fig)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/docs/conf.py b/docs/conf.py
index 71a2ab7..f884a50 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -4,11 +4,12 @@
     "hypothesis": True,
 }
 exclude_patterns = [
-    "**.ipynb_checkpoints",
+    "_build",
     ".DS_Store",
     ".pixi",
+    ".virtual_documents",
+    "**.ipynb_checkpoints",
     "Thumbs.db",
-    "_build",
 ]
 extensions = [
     "myst_nb",
diff --git a/docs/index.md b/docs/index.md
index c3dad8b..ed23ae0 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -20,4 +20,5 @@ symbolics
 pgamma-state
 branching-fraction
 ampform/etapi0p
+ampform/LambdaKpi0
 ```
diff --git a/pixi.toml b/pixi.toml
index 28d6fdd..bfcd8f5 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -75,3 +75,6 @@ sphinx-autobuild \
     docs/ docs/_build/html
 """
 env = {EXECUTE_NB = "yes"}
+
+[tasks.lab]
+cmd = "jupyter lab docs"