diff --git a/.cspell.json b/.cspell.json
index bb5795c08..1f4730b40 100644
--- a/.cspell.json
+++ b/.cspell.json
@@ -68,6 +68,7 @@
         "bdist",
         "bgcolor",
         "boldsymbol",
+        "byckling",
         "cahn",
         "cano",
         "celltoolbar",
@@ -254,6 +255,7 @@
         "isospin",
         "itertools",
         "jupyter",
+        "Källén",
         "lambdification",
         "lambdified",
         "lambdify",
diff --git a/docs/_extend_docstrings.py b/docs/_extend_docstrings.py
index f8fdd9727..17cf7f6f9 100644
--- a/docs/_extend_docstrings.py
+++ b/docs/_extend_docstrings.py
@@ -178,6 +178,24 @@ def extend_ComplexSqrt() -> None:
     )
 
 
+def extend_compute_third_mandelstam() -> None:
+    from ampform.kinematics.phasespace import compute_third_mandelstam
+
+    m0, m1, m2, m3 = sp.symbols("m:4")
+    s1, s2 = sp.symbols("sigma1 sigma2")
+    expr = compute_third_mandelstam(s1, s2, m0, m1, m2, m3)
+    _append_to_docstring(
+        compute_third_mandelstam,
+        Rf"""
+
+    .. math:: \sigma_3 = {sp.latex(expr)}
+        :label: compute_third_mandelstam
+
+    Note that this expression is symmetric in :math:`\sigma_{{1,2,3}}`.
+    """,
+    )
+
+
 def extend_EqualMassPhaseSpaceFactor() -> None:
     from ampform.dynamics.phasespace import (
         EqualMassPhaseSpaceFactor,
@@ -276,6 +294,35 @@ def extend_formulate_form_factor() -> None:
     )
 
 
+def extend_Kallen() -> None:
+    from ampform.kinematics.phasespace import Kallen
+
+    x, y, z = sp.symbols("x:z")
+    expr = Kallen(x, y, z)
+    _append_latex_doit_definition(expr)
+    _append_to_docstring(
+        Kallen,
+        """
+    .. seealso:: `.BreakupMomentumSquared`
+    """,
+    )
+
+
+def extend_Kibble() -> None:
+    from ampform.kinematics.phasespace import Kibble
+
+    m0, m1, m2, m3 = sp.symbols("m:4")
+    s1, s2, s3 = sp.symbols("sigma1:4")
+    expr = Kibble(s1, s2, s3, m0, m1, m2, m3)
+    _append_latex_doit_definition(expr)
+    _append_to_docstring(
+        Kibble,
+        R"""
+    with :math:`\lambda` defined by :eq:`Kallen`.
+    """,
+    )
+
+
 def extend_InvariantMass() -> None:
     from ampform.kinematics import InvariantMass
 
@@ -284,6 +331,24 @@ def extend_InvariantMass() -> None:
     _append_latex_doit_definition(expr)
 
 
+def extend_is_within_phasespace() -> None:
+    from ampform.kinematics.phasespace import is_within_phasespace
+
+    m0, m1, m2, m3 = sp.symbols("m:4")
+    s1, s2 = sp.symbols("sigma1 sigma2")
+    expr = is_within_phasespace(s1, s2, m0, m1, m2, m3)
+    _append_to_docstring(
+        is_within_phasespace,
+        Rf"""
+
+    .. math:: {sp.latex(expr)}
+        :label: is_within_phasespace
+
+    with :math:`\phi` defined by :eq:`Kibble`.
+    """,
+    )
+
+
 def extend_PhaseSpaceFactor() -> None:
     from ampform.dynamics.phasespace import PhaseSpaceFactor
 
diff --git a/docs/bibliography.bib b/docs/bibliography.bib
index 1237d9809..8cd2d2ed6 100755
--- a/docs/bibliography.bib
+++ b/docs/bibliography.bib
@@ -38,6 +38,16 @@ @article{aubertDalitzPlotAnalysis2005
   url = {https://link.aps.org/doi/10.1103/PhysRevD.72.052008}
 }
 
+@book{bycklingParticleKinematics1973,
+  title = {Particle {{Kinematics}}},
+  author = {Byckling, Eero and Kajantie, Keijo},
+  year = {1973},
+  publisher = {{Wiley}},
+  address = {{London, New York}},
+  isbn = {978-0-471-12885-4},
+  lccn = {QC794.6.K5 B95}
+}
+
 @article{cahnMystery9801986,
   title = {Mystery of the δ(980)},
   author = {Cahn, R.N. and Landshoff, P.V.},
diff --git a/docs/usage.ipynb b/docs/usage.ipynb
index 16170c20d..943255660 100644
--- a/docs/usage.ipynb
+++ b/docs/usage.ipynb
@@ -335,6 +335,7 @@
     "usage/dynamics\n",
     "usage/helicity/formalism\n",
     "usage/helicity/spin-alignment\n",
+    "usage/kinematics\n",
     "```"
    ]
   }
diff --git a/docs/usage/kinematics.ipynb b/docs/usage/kinematics.ipynb
new file mode 100644
index 000000000..b6852aecd
--- /dev/null
+++ b/docs/usage/kinematics.ipynb
@@ -0,0 +1,248 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "hideCode": true,
+    "hideOutput": true,
+    "hidePrompt": true,
+    "jupyter": {
+     "source_hidden": true
+    },
+    "slideshow": {
+     "slide_type": "skip"
+    },
+    "tags": [
+     "remove-cell",
+     "skip-execution"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# WARNING: advised to install a specific version, e.g. ampform==0.1.2\n",
+    "%pip install -q ampform[doc,viz] IPython"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "hideCode": true,
+    "hideOutput": true,
+    "hidePrompt": true,
+    "jupyter": {
+     "source_hidden": true
+    },
+    "slideshow": {
+     "slide_type": "skip"
+    },
+    "tags": [
+     "remove-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "import os\n",
+    "\n",
+    "from IPython.display import display  # noqa: F401\n",
+    "\n",
+    "STATIC_WEB_PAGE = {\"EXECUTE_NB\", \"READTHEDOCS\"}.intersection(os.environ)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "```{autolink-concat}\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "# Kinematics"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "import sympy as sp\n",
+    "from IPython.display import Math\n",
+    "\n",
+    "from ampform.io import aslatex"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    ":::{margin}\n",
+    "This notebook originates from {doc}`compwa-org:report/017`.\n",
+    ":::\n",
+    "\n",
+    "Kinematics for a three-body decay $0 \\to 123$ can be fully described by two **Mandelstam variables** $\\sigma_1, \\sigma_2$, because the third variable $\\sigma_3$ can be expressed in terms $\\sigma_1, \\sigma_2$, the mass $m_0$ of the initial state, and the masses $m_1, m_2, m_3$ of the final state. As can be seen, the roles of $\\sigma_1, \\sigma_2, \\sigma_3$ are interchangeable.\n",
+    "\n",
+    "```{margin}\n",
+    "See Eq. (1.2) in {cite}`bycklingParticleKinematics1973`\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from ampform.kinematics.phasespace import compute_third_mandelstam\n",
+    "\n",
+    "m0, m1, m2, m3 = sp.symbols(\"m:4\")\n",
+    "s1, s2, s3 = sp.symbols(\"sigma1:4\")\n",
+    "s3_expr = compute_third_mandelstam(s1, s2, m0, m1, m2, m3)\n",
+    "\n",
+    "latex = aslatex({s3: s3_expr})\n",
+    "Math(latex)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The phase space is defined by the closed area that satisfies the condition $\\phi(\\sigma_1,\\sigma_2) \\leq 0$, where $\\phi$ is a **Kibble function**:\n",
+    "\n",
+    "\n",
+    "```{margin}\n",
+    "See §V.2 in {cite}`bycklingParticleKinematics1973`\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from ampform.kinematics.phasespace import Kibble\n",
+    "\n",
+    "kibble = Kibble(s1, s2, s3, m0, m1, m2, m3)\n",
+    "\n",
+    "latex = aslatex({kibble: kibble.evaluate()})\n",
+    "Math(latex)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "and $\\lambda$ is the **Källén function**:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from ampform.kinematics.phasespace import Kallen\n",
+    "\n",
+    "x, y, z = sp.symbols(\"x:z\")\n",
+    "kallen = Kallen(x, y, z)\n",
+    "\n",
+    "latex = aslatex({kallen: kallen.evaluate()})\n",
+    "Math(latex)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Any distribution over the phase space can now be defined using a two-dimensional grid over a Mandelstam pair $\\sigma_1,\\sigma_2$ of choice, with the condition $\\phi(\\sigma_1,\\sigma_2)<0$ selecting the values that are physically allowed."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from ampform.kinematics.phasespace import is_within_phasespace\n",
+    "\n",
+    "is_within_phasespace(s1, s2, m0, m1, m2, m3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "See {doc}`compwa-org:report/017` for an interactive visualization of the phase space region and an analytic expression for the phase space boundary."
+   ]
+  }
+ ],
+ "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.8.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/src/ampform/kinematics/phasespace.py b/src/ampform/kinematics/phasespace.py
new file mode 100644
index 000000000..708de9374
--- /dev/null
+++ b/src/ampform/kinematics/phasespace.py
@@ -0,0 +1,75 @@
+# pylint: disable=arguments-differ, abstract-method, protected-access,
+# pylint: disable=too-many-arguments, unbalanced-tuple-unpacking
+"""Functions for determining phase space boundaries.
+
+.. seealso:: :doc:`/usage/kinematics`
+"""
+from __future__ import annotations
+
+import sympy as sp
+
+from ampform.sympy import (
+    UnevaluatedExpression,
+    create_expression,
+    implement_doit_method,
+    make_commutative,
+)
+
+
+@make_commutative
+@implement_doit_method
+class Kibble(UnevaluatedExpression):
+    """Kibble function for determining the phase space region."""
+
+    def __new__(
+        cls, sigma1, sigma2, sigma3, m0, m1, m2, m3, **hints
+    ) -> Kibble:
+        return create_expression(
+            cls, sigma1, sigma2, sigma3, m0, m1, m2, m3, **hints
+        )
+
+    def evaluate(self) -> Kallen:
+        sigma1, sigma2, sigma3, m0, m1, m2, m3 = self.args
+        return Kallen(
+            Kallen(sigma2, m2**2, m0**2),  # type: ignore[operator]
+            Kallen(sigma3, m3**2, m0**2),  # type: ignore[operator]
+            Kallen(sigma1, m1**2, m0**2),  # type: ignore[operator]
+        )
+
+    def _latex(self, printer, *args):
+        sigma1, sigma2, *_ = map(printer._print, self.args)
+        return Rf"\phi\left({sigma1}, {sigma2}\right)"
+
+
+@make_commutative
+@implement_doit_method
+class Kallen(UnevaluatedExpression):
+    """Källén function, used for computing break-up momenta."""
+
+    def __new__(cls, x, y, z, **hints) -> Kallen:
+        return create_expression(cls, x, y, z, **hints)
+
+    def evaluate(self) -> sp.Expr:
+        x, y, z = self.args
+        return x**2 + y**2 + z**2 - 2 * x * y - 2 * y * z - 2 * z * x  # type: ignore[operator]
+
+    def _latex(self, printer, *args):
+        x, y, z = map(printer._print, self.args)
+        return Rf"\lambda\left({x}, {y}, {z}\right)"
+
+
+def is_within_phasespace(
+    sigma1, sigma2, m0, m1, m2, m3, outside_value=sp.nan
+) -> sp.Piecewise:
+    """Determine whether a set of masses lie within phase space."""
+    sigma3 = compute_third_mandelstam(sigma1, sigma2, m0, m1, m2, m3)
+    kibble = Kibble(sigma1, sigma2, sigma3, m0, m1, m2, m3)
+    return sp.Piecewise(
+        (1, sp.LessThan(kibble, 0)),
+        (outside_value, True),
+    )
+
+
+def compute_third_mandelstam(sigma1, sigma2, m0, m1, m2, m3) -> sp.Add:
+    """Compute the third Mandelstam variable in a three-body decay."""
+    return m0**2 + m1**2 + m2**2 + m3**2 - sigma1 - sigma2
diff --git a/tests/test_dynamics.py b/tests/dynamics/test_dynamics.py
similarity index 100%
rename from tests/test_dynamics.py
rename to tests/dynamics/test_dynamics.py
diff --git a/tests/test_angular_distributions.py b/tests/helicity/test_angular_distributions.py
similarity index 100%
rename from tests/test_angular_distributions.py
rename to tests/helicity/test_angular_distributions.py
diff --git a/tests/test_parity_prefactor.py b/tests/helicity/test_parity_prefactor.py
similarity index 100%
rename from tests/test_parity_prefactor.py
rename to tests/helicity/test_parity_prefactor.py
diff --git a/tests/kinematics/test_phasespace.py b/tests/kinematics/test_phasespace.py
new file mode 100644
index 000000000..0c8d4b6ed
--- /dev/null
+++ b/tests/kinematics/test_phasespace.py
@@ -0,0 +1,21 @@
+import pytest
+
+from ampform.kinematics.phasespace import is_within_phasespace
+
+
+@pytest.mark.parametrize(
+    ("s1", "s2", "expected"),
+    [
+        (0.0, 3.0, 0),
+        (1.0, 1.0, 1),
+        (2.0, 2.0, 1),
+    ],
+)
+def test_is_within_phasespace(s1, s2, expected):
+    # See widget https://compwa-org.rtfd.io/report/017.html
+    m0 = 2.1
+    m1 = 0.2
+    m2 = 0.4
+    m3 = 0.4
+    computed = is_within_phasespace(s1, s2, m0, m1, m2, m3, outside_value=0)
+    assert computed.doit() == expected