docs/report/017.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{autolink-concat}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "physics",
     "sympy"
    ]
   },
   "source": [
    "::::{margin}\n",
    ":::{card} Symbolic phase space boundary for a three-body decay\n",
    "TR-017\n",
    "^^^\n",
    "This reports shows how define the physical phase space region on a Dalitz plot using a Kibble function.\n",
    "+++\n",
    "✅&nbsp;[compwa.github.io#139](https://github.com/ComPWA/compwa.github.io/issues/139)\n",
    ":::\n",
    "::::"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Phase space for a three-body decay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "%pip install -q ampform==0.14.0 sympy==1.10.1 tensorwaves==0.4.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_formats = ['svg']\n",
    "\n",
    "from __future__ import annotations\n",
    "\n",
    "import os\n",
    "import warnings\n",
    "from typing import TYPE_CHECKING\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "from ampform.sympy import (\n",
    "    UnevaluatedExpression,\n",
    "    create_expression,\n",
    "    implement_doit_method,\n",
    "    make_commutative,\n",
    ")\n",
    "from IPython.display import Math\n",
    "from ipywidgets import FloatSlider, VBox, interactive_output\n",
    "from matplotlib.patches import Patch\n",
    "from tensorwaves.function.sympy import create_parametrized_function\n",
    "\n",
    "if TYPE_CHECKING:\n",
    "    from matplotlib.axis import Axis\n",
    "    from matplotlib.contour import QuadContourSet\n",
    "    from matplotlib.lines import Line2D\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "editable": true,
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "if \"COLAB_RELEASE_TAG\" in os.environ:\n",
    "    import subprocess\n",
    "\n",
    "    from google.colab import output\n",
    "\n",
    "    output.enable_custom_widget_manager()\n",
    "    subprocess.run(\"pip install -q ipympl\".split(), check=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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}`Byckling:1971vca`\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{3} = m_{0}^{2} + m_{1}^{2} + m_{2}^{2} + m_{3}^{2} - \\sigma_{1} - \\sigma_{2}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def compute_third_mandelstam(sigma1, sigma2, m0, m1, m2, m3) -> sp.Expr:\n",
    "    return m0**2 + m1**2 + m2**2 + m3**2 - sigma1 - sigma2\n",
    "\n",
    "\n",
    "m0, m1, m2, m3 = sp.symbols(\"m:4\")\n",
    "s1, s2, s3 = sp.symbols(\"sigma1:4\")\n",
    "computed_s3 = compute_third_mandelstam(s1, s2, m0, m1, m2, m3)\n",
    "Math(Rf\"{sp.latex(s3)} = {sp.latex(computed_s3)}\")"
   ]
  },
  {
   "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}`Byckling:1971vca`\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "scroll-input",
     "hide-input",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\phi\\left(\\sigma_{1}, \\sigma_{2}\\right) = \\lambda\\left(\\lambda\\left(\\sigma_{2}, m_{2}^{2}, m_{0}^{2}\\right), \\lambda\\left(\\sigma_{3}, m_{3}^{2}, m_{0}^{2}\\right), \\lambda\\left(\\sigma_{1}, m_{1}^{2}, m_{0}^{2}\\right)\\right)$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@make_commutative\n",
    "@implement_doit_method\n",
    "class Kibble(UnevaluatedExpression):\n",
    "    def __new__(cls, sigma1, sigma2, sigma3, m0, m1, m2, m3, **hints) -> Kibble:\n",
    "        return create_expression(cls, sigma1, sigma2, sigma3, m0, m1, m2, m3, **hints)\n",
    "\n",
    "    def evaluate(self) -> sp.Expr:\n",
    "        sigma1, sigma2, sigma3, m0, m1, m2, m3 = self.args\n",
    "        return Källén(\n",
    "            Källén(sigma2, m2**2, m0**2),\n",
    "            Källén(sigma3, m3**2, m0**2),\n",
    "            Källén(sigma1, m1**2, m0**2),\n",
    "        )\n",
    "\n",
    "    def _latex(self, printer, *args):\n",
    "        sigma1, sigma2, *_ = map(printer._print, self.args)\n",
    "        return Rf\"\\phi\\left({sigma1}, {sigma2}\\right)\"\n",
    "\n",
    "\n",
    "@make_commutative\n",
    "@implement_doit_method\n",
    "class Källén(UnevaluatedExpression):\n",
    "    def __new__(cls, x, y, z, **hints) -> Källén:\n",
    "        return create_expression(cls, x, y, z, **hints)\n",
    "\n",
    "    def evaluate(self) -> sp.Expr:\n",
    "        x, y, z = self.args\n",
    "        return x**2 + y**2 + z**2 - 2 * x * y - 2 * y * z - 2 * z * x\n",
    "\n",
    "    def _latex(self, printer, *args):\n",
    "        x, y, z = map(printer._print, self.args)\n",
    "        return Rf\"\\lambda\\left({x}, {y}, {z}\\right)\"\n",
    "\n",
    "\n",
    "kibble = Kibble(s1, s2, s3, m0, m1, m2, m3)\n",
    "Math(Rf\"{sp.latex(kibble)} = {sp.latex(kibble.doit(deep=False))}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and $\\lambda$ is the **Källén function**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\lambda\\left(x, y, z\\right) = x^{2} - 2 x y - 2 x z + y^{2} - 2 y z + z^{2}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, y, z = sp.symbols(\"x:z\")\n",
    "expr = Källén(x, y, z)\n",
    "Math(f\"{sp.latex(expr)} = {sp.latex(expr.doit())}\")"
   ]
  },
  {
   "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",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\begin{cases} 1 & \\text{for}\\: \\phi\\left(\\sigma_{1}, \\sigma_{2}\\right) \\leq 0 \\\\\\text{NaN} & \\text{otherwise} \\end{cases}$"
      ],
      "text/plain": [
       "Piecewise((1, Kibble(sigma1, sigma2, m0**2 + m1**2 + m2**2 + m3**2 - sigma1 - sigma2, m0, m1, m2, m3) <= 0), (nan, True))"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def is_within_phasespace(\n",
    "    sigma1, sigma2, m0, m1, m2, m3, outside_value=sp.nan\n",
    ") -> sp.Piecewise:\n",
    "    sigma3 = compute_third_mandelstam(sigma1, sigma2, m0, m1, m2, m3)\n",
    "    kibble = Kibble(sigma1, sigma2, sigma3, m0, m1, m2, m3)\n",
    "    return sp.Piecewise(\n",
    "        (1, sp.LessThan(kibble, 0)),\n",
    "        (outside_value, True),\n",
    "    )\n",
    "\n",
    "\n",
    "is_within_phasespace(s1, s2, m0, m1, m2, m3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phsp_expr = is_within_phasespace(s1, s2, m0, m1, m2, m3, outside_value=0)\n",
    "phsp_func = create_parametrized_function(\n",
    "    phsp_expr.doit(),\n",
    "    parameters={m0: 2.2, m1: 0.2, m2: 0.4, m3: 0.4},\n",
    "    backend=\"numpy\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "scroll-input",
     "hide-input",
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "sliders = {\n",
    "    \"m0\": FloatSlider(description=\"m0\", max=3, value=2.1, step=0.01),\n",
    "    \"m1\": FloatSlider(description=\"m1\", max=2, value=0.2, step=0.01),\n",
    "    \"m2\": FloatSlider(description=\"m2\", max=2, value=0.4, step=0.01),\n",
    "    \"m3\": FloatSlider(description=\"m3\", max=2, value=0.4, step=0.01),\n",
    "}\n",
    "\n",
    "resolution = 300\n",
    "X, Y = np.meshgrid(\n",
    "    np.linspace(0, 4, num=resolution),\n",
    "    np.linspace(0, 4, num=resolution),\n",
    ")\n",
    "data = {\"sigma1\": X, \"sigma2\": Y}\n",
    "\n",
    "sidebar_ratio = 0.15\n",
    "fig, ((ax1, _), (ax, ax2)) = plt.subplots(\n",
    "    figsize=(7, 7),\n",
    "    ncols=2,\n",
    "    nrows=2,\n",
    "    gridspec_kw={\n",
    "        \"height_ratios\": [sidebar_ratio, 1],\n",
    "        \"width_ratios\": [1, sidebar_ratio],\n",
    "    },\n",
    ")\n",
    "_.remove()\n",
    "ax.set_xlim(0, 4)\n",
    "ax.set_ylim(0, 4)\n",
    "ax.set_xlabel(R\"$\\sigma_1$\")\n",
    "ax.set_ylabel(R\"$\\sigma_2$\")\n",
    "ax.set_xticks(range(5))\n",
    "ax.set_yticks(range(5))\n",
    "ax1.set_xlim(0, 4)\n",
    "ax2.set_ylim(0, 4)\n",
    "for a in [ax1, ax2]:\n",
    "    a.set_xticks([])\n",
    "    a.set_yticks([])\n",
    "    a.axis(\"off\")\n",
    "fig.tight_layout()\n",
    "fig.subplots_adjust(wspace=0, hspace=0)\n",
    "\n",
    "fig.canvas.toolbar_visible = False\n",
    "fig.canvas.header_visible = False\n",
    "fig.canvas.footer_visible = False\n",
    "\n",
    "MESH: QuadContourSet | None = None\n",
    "PROJECTIONS: tuple[Line2D, Line2D] = None\n",
    "BOUNDARIES: list[Line2D] | None = None\n",
    "\n",
    "\n",
    "def plot(**parameters):\n",
    "    draw_boundaries(\n",
    "        parameters[\"m0\"],\n",
    "        parameters[\"m1\"],\n",
    "        parameters[\"m2\"],\n",
    "        parameters[\"m3\"],\n",
    "    )\n",
    "    global MESH, PROJECTIONS\n",
    "    if MESH is not None:\n",
    "        for coll in MESH.collections:\n",
    "            ax.collections.remove(coll)\n",
    "    phsp_func.update_parameters(parameters)\n",
    "    Z = phsp_func(data)\n",
    "    MESH = ax.contour(X, Y, Z, colors=\"black\")\n",
    "    contour = MESH.collections[0]\n",
    "    contour.set_facecolor(\"lightgray\")\n",
    "    x = X[0]\n",
    "    y = Y[:, 0]\n",
    "    Zx = np.nansum(Z, axis=0)\n",
    "    Zy = np.nansum(Z, axis=1)\n",
    "    if PROJECTIONS is None:\n",
    "        PROJECTIONS = (\n",
    "            ax1.plot(x, Zx, c=\"black\", lw=2)[0],\n",
    "            ax2.plot(Zy, y, c=\"black\", lw=2)[0],\n",
    "        )\n",
    "    else:\n",
    "        PROJECTIONS[0].set_data(x, Zx)\n",
    "        PROJECTIONS[1].set_data(Zy, y)\n",
    "    ax1.relim()\n",
    "    ax2.relim()\n",
    "    ax1.autoscale_view(scalex=False)\n",
    "    ax2.autoscale_view(scaley=False)\n",
    "    create_legend(ax)\n",
    "    fig.canvas.draw()\n",
    "\n",
    "\n",
    "def draw_boundaries(m0, m1, m2, m3) -> None:\n",
    "    global BOUNDARIES\n",
    "    s1_min = (m2 + m3) ** 2\n",
    "    s1_max = (m0 - m1) ** 2\n",
    "    s2_min = (m1 + m3) ** 2\n",
    "    s2_max = (m0 - m2) ** 2\n",
    "    if BOUNDARIES is None:\n",
    "        BOUNDARIES = [\n",
    "            ax.axvline(s1_min, c=\"red\", ls=\"dotted\", label=\"$(m_2+m_3)^2$\"),\n",
    "            ax.axhline(s2_min, c=\"blue\", ls=\"dotted\", label=\"$(m_1+m_3)^2$\"),\n",
    "            ax.axvline(s1_max, c=\"red\", ls=\"dashed\", label=\"$(m_0-m_1)^2$\"),\n",
    "            ax.axhline(s2_max, c=\"blue\", ls=\"dashed\", label=\"$(m_0-m_2)^2$\"),\n",
    "        ]\n",
    "    else:\n",
    "        BOUNDARIES[0].set_data(get_line_data(s1_min))\n",
    "        BOUNDARIES[1].set_data(get_line_data(s2_min, horizontal=True))\n",
    "        BOUNDARIES[2].set_data(get_line_data(s1_max))\n",
    "        BOUNDARIES[3].set_data(get_line_data(s2_max, horizontal=True))\n",
    "\n",
    "\n",
    "def create_legend(ax: Axis):\n",
    "    if ax.get_legend() is not None:\n",
    "        return\n",
    "    label = Rf\"${sp.latex(kibble)}\\leq0$\"\n",
    "    ax.legend(\n",
    "        handles=[\n",
    "            Patch(label=label, ec=\"black\", fc=\"lightgray\"),\n",
    "            *BOUNDARIES,\n",
    "        ],\n",
    "        loc=\"upper right\",\n",
    "        facecolor=\"white\",\n",
    "        framealpha=1,\n",
    "    )\n",
    "\n",
    "\n",
    "def get_line_data(value, horizontal: bool = False) -> np.ndarray:\n",
    "    pair = (value, value)\n",
    "    if horizontal:\n",
    "        return np.array([(0, 1), pair])\n",
    "    return np.array([pair, (0, 1)])\n",
    "\n",
    "\n",
    "output = interactive_output(plot, controls=sliders)\n",
    "VBox([output, *sliders.values()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "{{ run_interactive }}\n",
    "\n",
    "```{image} https://user-images.githubusercontent.com/29308176/165257660-af66c9c1-2da9-475f-b4c4-4a52950acd00.gif\n",
    ":alt: Cell output - interactive Dalitz plot\n",
    ":align: center\n",
    "```\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The phase space boundary can be described analytically in terms of $\\sigma_1$ or $\\sigma_2$, in which case there are two solutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sol1, sol2 = sp.solve(kibble.doit().subs(s3, computed_s3), s2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "remove-input",
     "keep_output",
     "full-width"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\begin{array}{c}\n",
       "  \\frac{- m_{0}^{2} m_{2}^{2} + m_{0}^{2} m_{3}^{2} + m_{0}^{2} \\sigma_{1} + m_{1}^{2} m_{2}^{2} - m_{1}^{2} m_{3}^{2} + m_{1}^{2} \\sigma_{1} + m_{2}^{2} \\sigma_{1} + m_{3}^{2} \\sigma_{1} - \\sigma_{1}^{2} - \\sqrt{\\left(m_{0}^{2} - 2 m_{0} m_{1} + m_{1}^{2} - \\sigma_{1}\\right) \\left(m_{0}^{2} + 2 m_{0} m_{1} + m_{1}^{2} - \\sigma_{1}\\right) \\left(m_{2}^{2} - 2 m_{2} m_{3} + m_{3}^{2} - \\sigma_{1}\\right) \\left(m_{2}^{2} + 2 m_{2} m_{3} + m_{3}^{2} - \\sigma_{1}\\right)}}{2 \\sigma_{1}} \\\\\n",
       "  \\frac{- m_{0}^{2} m_{2}^{2} + m_{0}^{2} m_{3}^{2} + m_{0}^{2} \\sigma_{1} + m_{1}^{2} m_{2}^{2} - m_{1}^{2} m_{3}^{2} + m_{1}^{2} \\sigma_{1} + m_{2}^{2} \\sigma_{1} + m_{3}^{2} \\sigma_{1} - \\sigma_{1}^{2} + \\sqrt{\\left(m_{0}^{2} - 2 m_{0} m_{1} + m_{1}^{2} - \\sigma_{1}\\right) \\left(m_{0}^{2} + 2 m_{0} m_{1} + m_{1}^{2} - \\sigma_{1}\\right) \\left(m_{2}^{2} - 2 m_{2} m_{3} + m_{3}^{2} - \\sigma_{1}\\right) \\left(m_{2}^{2} + 2 m_{2} m_{3} + m_{3}^{2} - \\sigma_{1}\\right)}}{2 \\sigma_{1}} \\\\\n",
       "\\end{array}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latex = R\"\\begin{array}{c}\" + \"\\n\"\n",
    "for sol in [sol1, sol2]:\n",
    "    latex += Rf\"  {sp.latex(sol)} \\\\\" + \"\\n\"\n",
    "latex += R\"\\end{array}\"\n",
    "Math(latex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The boundary cannot be parametrized analytically in polar coordinates, but there is a numeric solution. The idea is to solve the condition $\\phi(\\sigma_1,\\sigma_2)=0$ after the following substitutions:\n",
    "\n",
    "```{margin}\n",
    "See {cite}`Byckling:1971vca`, pp.&nbsp;109–112\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "T0, T1, T2, T3 = sp.symbols(\"T0:4\")\n",
    "r, theta = sp.symbols(\"r theta\", nonnegative=True)\n",
    "substitutions = {\n",
    "    s1: (m2 + m3) ** 2 + T1,\n",
    "    s2: (m1 + m3) ** 2 + T2,\n",
    "    s3: (m1 + m2) ** 2 + T3,\n",
    "    T1: T0 / 3 - r * sp.cos(theta),\n",
    "    T2: T0 / 3 - r * sp.cos(theta + 2 * sp.pi / 3),\n",
    "    T3: T0 / 3 - r * sp.cos(theta + 4 * sp.pi / 3),\n",
    "    T0: (\n",
    "        m0**2 + m1**2 + m2**2 + m3**2 - (m2 + m3) ** 2 - (m1 + m3) ** 2 - (m1 + m2) ** 2\n",
    "    ),\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "remove-input",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\begin{array}{rcl}\n",
       "  \\sigma_{1} &=& T_{1} + \\left(m_{2} + m_{3}\\right)^{2} \\\\\n",
       "  \\sigma_{2} &=& T_{2} + \\left(m_{1} + m_{3}\\right)^{2} \\\\\n",
       "  \\sigma_{3} &=& T_{3} + \\left(m_{1} + m_{2}\\right)^{2} \\\\\n",
       "  T_{1} &=& \\frac{T_{0}}{3} - r \\cos{\\left(\\theta \\right)} \\\\\n",
       "  T_{2} &=& \\frac{T_{0}}{3} + r \\sin{\\left(\\theta + \\frac{\\pi}{6} \\right)} \\\\\n",
       "  T_{3} &=& \\frac{T_{0}}{3} + r \\cos{\\left(\\theta + \\frac{\\pi}{3} \\right)} \\\\\n",
       "  T_{0} &=& m_{0}^{2} + m_{1}^{2} + m_{2}^{2} + m_{3}^{2} - \\left(m_{1} + m_{2}\\right)^{2} - \\left(m_{1} + m_{3}\\right)^{2} - \\left(m_{2} + m_{3}\\right)^{2} \\\\\n",
       "\\end{array}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latex = R\"\\begin{array}{rcl}\" + \"\\n\"\n",
    "for item in substitutions.items():\n",
    "    lhs, rhs = map(sp.latex, item)\n",
    "    latex += Rf\"  {lhs} &=& {rhs} \\\\\" + \"\\n\"\n",
    "latex += R\"\\end{array}\"\n",
    "Math(latex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For every value of $\\theta \\in [0, 2\\pi)$, the value of $r$ can now be found by solving the condition $\\phi(r, \\theta)=0$. Note that $\\phi(r, \\theta)$ is a **cubic polynomial of $r$**. For instance, if we take $m_0=5, m_1=2, m_{2,3}=1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input",
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "phi_r = (\n",
    "    kibble.doit()\n",
    "    .subs(substitutions)  # substitute sigmas\n",
    "    .subs(substitutions)  # substitute T123\n",
    "    .subs(substitutions)  # substitute T0\n",
    "    .subs({m0: 5, m1: 2, m2: 1, m3: 1})\n",
    "    .simplify()\n",
    "    .collect(r)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "remove-input",
     "keep_output"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\begin{eqnarray}\n",
       "\\phi(r, \\theta) & = & r^{3} \\cdot \\left(56 \\sqrt{3} \\sin{\\left(\\theta \\right)} + 400 \\cos^{3}{\\left(\\theta \\right)} - 356 \\cos{\\left(\\theta \\right)} + 112 \\cos{\\left(\\theta + \\frac{\\pi}{3} \\right)}\\right) \\nonumber\\\\\n",
       "& & + r^{2} \\cdot \\left(2000 \\cos^{2}{\\left(\\theta \\right)} + 2100\\right) \\nonumber\\\\\n",
       "& & - 4800 r \\cos{\\left(\\theta \\right)} \\nonumber\\\\\n",
       "& & + -25200 \n",
       "\\end{eqnarray}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lhs = sp.Symbol(R\"\\phi(r, \\theta)\")\n",
    "latex = sp.multiline_latex(lhs, phi_r, environment=\"eqnarray\")\n",
    "Math(latex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lowest value of $r$ that satisfies $\\phi(r,\\theta)=0$ defines the phase space boundary."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "toc_visible": true
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}