{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# The type-I seesaw in `feynlag`: Majorana neutrino masses, step by step\n",
    "\n",
    "Neutrinos are massless in the Standard Model.  The **type-I seesaw** is the\n",
    "simplest way to give them mass: add right-handed neutrinos $\\nu_R$ (complete\n",
    "gauge singlets) with\n",
    "\n",
    "- a **Dirac Yukawa** $-Y_\\nu \\bar L \\tilde H\\, \\nu_R$ (like every other fermion), and\n",
    "- a **Majorana mass** $-\\tfrac12 M_R\\, \\nu_R^T C \\nu_R$ — *allowed* precisely\n",
    "  because $\\nu_R$ carries no gauge charge, and the origin of the mechanism.\n",
    "\n",
    "After electroweak symmetry breaking the neutral-fermion mass matrix, in the\n",
    "left-handed basis $n = (\\nu_L,\\ \\nu_R^c)$, is\n",
    "\n",
    "$$M_\\nu = \\begin{pmatrix} 0 & m_D \\\\ m_D^T & M_R \\end{pmatrix},\\qquad m_D = Y_\\nu v/\\sqrt2.$$\n",
    "\n",
    "For $M_R \\gg m_D$ this **seesaws**: a naturally tiny light neutrino\n",
    "$m_\\nu \\approx -m_D M_R^{-1} m_D^T$ and a heavy state $\\approx M_R$, mixed by\n",
    "$V \\approx m_D M_R^{-1}$.  This notebook walks the whole construction in feynlag\n",
    "— including extracting the physical heavy-neutrino couplings, which requires the\n",
    "Majorana machinery (`diracC`, `MajoranaBilinear`) and a charge-conjugation-aware\n",
    "mass-basis rotation, because the physical neutrinos are Majorana fields mixing\n",
    "$\\nu_L$ with the charge-conjugate $\\nu_R^c$.\n",
    "\n",
    "We use one generation for transparency; the machinery is generation-general."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## 1. Symmetries, parameters, fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()\n",
    "\n",
    "from feynlag import (\n",
    "    Bilinear, DiracGamma, Dmu, ExternalParameter, InternalParameter, Lagrangian,\n",
    "    MajoranaBilinear, MajoranaRotation, Model, Rotation, SU2, Scalar, U1,\n",
    "    WeylFermion, dag, diagonalize_takagi, diracC, diracPL, diracPR,\n",
    "    extract_fermion_vertices, fermion_gauge_current, fermion_mass_matrix,\n",
    "    majorana_mass_matrix, rotation_2x2, seesaw_light_mass, seesaw_mass_matrix,\n",
    "    weinberg_rotation, charged_current_rotation,\n",
    ")\n",
    "\n",
    "gw = ExternalParameter(\"gw\", 0.6535, positive=True)\n",
    "g1 = ExternalParameter(\"g1\", 0.3580, positive=True)\n",
    "SU2L, U1Y = SU2(\"SU2L\", coupling=gw), U1(\"U1Y\", coupling=g1)\n",
    "v = ExternalParameter(\"v\", 246.0, positive=True, unit_dim=1)\n",
    "lam = ExternalParameter(\"lam\", 0.129)\n",
    "mu2 = InternalParameter(\"mu2\", unit_dim=2)\n",
    "yv = ExternalParameter(\"yv\", 0.01, positive=True)              # Dirac Yukawa\n",
    "MR = ExternalParameter(\"MR\", 1.0e3, positive=True, unit_dim=1)  # seesaw scale"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "The right-handed neutrino is a **total singlet** — `reps={}` — which is exactly\n",
    "what lets it carry a bare Majorana mass (nothing forbids $\\nu_R^T C \\nu_R$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ν_R reps (total singlet): {}\n"
     ]
    }
   ],
   "source": [
    "H = Scalar(\"H\", reps={SU2L: 2, U1Y: sp.Rational(1, 2)}, component_names=[\"Gp\", \"H0\"])\n",
    "H.expand_vev({H.components[1]: v})\n",
    "Ll = WeylFermion(\"Ll\", reps={SU2L: 2, U1Y: -sp.Rational(1, 2)}, chirality=\"L\",\n",
    "                 nflavors=1, component_names=[\"nuL\", \"eL\"])\n",
    "eR = WeylFermion(\"eR\", reps={U1Y: -1}, chirality=\"R\", nflavors=1, component_names=[\"eR\"])\n",
    "nuR = WeylFermion(\"nuR\", reps={}, chirality=\"R\", nflavors=1, component_names=[\"nuR\"])\n",
    "print(\"ν_R reps (total singlet):\", nuR.reps)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "## 2. The neutrino Lagrangian\n",
    "\n",
    "The Dirac Yukawa uses the conjugate doublet $\\tilde H = (H^{0*},\\,-H^{+*})$\n",
    "(the same pattern as the up-type quarks).  The Majorana mass is written with the\n",
    "`MajoranaBilinear` $\\nu_R^T C\\,P_L\\,\\nu_R$ — feynlag's charge-conjugation-aware\n",
    "sandwich (the middle carries `diracC`, the charge-conjugation matrix $C=i\\gamma^2\\gamma^0$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJEAAAAZCAYAAAAxO8yWAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGXUlEQVRoBe2Z4XFUOQyAN5ktgMtVAHRwQAehA3KpAOjgGP7lHwMdJFRwkA6gAyAdkA4ulw5y3+dYHvvtey/Jwu2+ZVYzXtuyZMuSLMtvd66urmZTh6OjoxfI+IrymXJJeUDZp3zI/Xu5/wraU9pbWKEG5itc60eWegrzIxxEB5plp3pG/TImzbg0HrhtvRoN7K5mmeVXwTmMOp+oawfRqYxKNVzQOa8R2/ZqNDB5J1INONBJRx1eZZ86uDPotk7UUcoqujubkBPVisBRjEzfKV5vZ/XYtr0eDaw8J8Lwf7FVryOjiRDJsslxwDF0Qwly4vsZDsQcz1jQvErHtJyCO6BuANyQzNIp9x7lK8XEvr52QU0TRvZ0WzuUja0tErEJn4WfqXWoAvQj0gwZ9CPE97p8ZYJbNOAN53Euja8TxLw7Q1PAp8xem4+6NOC8Xh9T7tPeCEdyD8i6lB3kDdiNxiprBP8jr9fNa9yUeY3RydeXDtWFvnyoSzPYZ863DL6kfko5oegUrmdi3k3WyzzQhMx/F2TbOKbrKfZzxEZAtadl7FD2uPLrLK9cX2VFmKphZFgANq1TaahBYy8wVQj4NfA+dRNJ8rzOvXCVVew3yRwO/3vFM/XmTXvqtUN3U+tyokMFwXgLiTE4T7zFqy69tjJOHvGCkeQ75d119+ZfaDWy0aK5PjOneZFzLsiTx60S3whNzDsUqaqpJtO8kx2GpE5OhGL0SA30hPKcskdRsf9QVI7faYrBaEfYfsjYN/rNE5y++cVz6qHcIDkJNA1Ar6HlbfIh8Bp3zMDNPAMdHSiuroaE+f0afhOoo94ICL/7cdzEemk54Z20HYYUNM8D5gcqwLvxPeWLfceodZB/qTVsfId5TVt6E1TpixOB88own+m9GsCrqAS0ffUIXgE60AXlAPzShoB/AZjPuV23fOFeIBpBwK+TCE3ukOd1TvWg3EMvSnlvAxtlB/VCOZvzo3JDOSr7HFyJOvSNSoLRZ8aYCv1iGzBKafgaxI05geOC10dcV86tM6b1HfzJoJGF3khyPTT6G47/BJlNzAP8XmWUbiIZ/XBY19WxPJS1ToO/1Jln0+zwEbkP5uziK41LiobUiN3TGqcwopBOFk7yJ/RvKDWowC6uO+56Md+MtteeJ9nn5mtKY5RgZlwZ96gLb4zdUHtNz5bgi2mT48PfG12DKGrozOfU2/5teaDdGDtU+zSdme2ySQ0opNOmAq675feQlrlEoqtq6TVqfZXphOK6c4AqoHIXxplXPsE5hkAHi6gyRNOHT3NXa/TRjOHcaxycMbp6TMdb2GdNULdDr+A2wQ616LN51UvOUvWjqdFSdNII1WY9lUalyyCkTgYG16tw8ElB0EXYrlivlQfiokbWbfh7I1RNM9AOGR8z3mtY5jY/87tR0Kap6Eck7uUbWE+0e+1G9RHyMjRpO6APD6R5r1HIq/x0t4h+velGURDEqf+Q6YwEAUaMrrOU0wfvi7xg0FvHdRDz1WPpygFhnpEA/jBgoJatw2nrfKbMxTppnLpxoEww5vhljrpRyd3os6YZabtew8d8U7LDW+Qxv1NnyR/mbgZkXEPdbxzJURg3h9H76nHzkj35BcY1uCc9rreH8jkm0HYuc6ghKM4jAfQqU6/3KrXt/Ca24Yh0bwfwGGGMYr4mvlGbs3kAnNNoMfZfnZFhBl9jWHEjoLxN3jdCW4ZYY+p20B5xG3jo/ctolv47Q3g3rSJTouSAQF+m9xS97oJ+ecJWYzqT35OsNYyn3debyWWKVNR9EcCPheFwkKf1jql0RJ1VI2h8ZTCqvaOYeP9GfUl9J8jzKEc4sgpovn/VE0KvLBpV3Qg60SB9osg/8Kov5W8cPsvwmLrXIcFP1g71/mwjq4dRnzlZ2x+wXaGG+gipIS8oKthQ2jj6EN868cios/sJo3tIdOI34O98CNa5n7618x69bc7nfQRTwimk8lCna2dKsvXJgpxekUITbcCb1zyg/hUcyD2W63ryTqQ1ULxXmpHIb0npK6n4qQGy+cJLORS1kUgRld0rWsXHh1aaGw3aohySjXAiBFbo9K0KwxiRLJMDZPPVYvnVwYNizphgU5xIrz/MJ70vSc/b2Vb/lwbQvYm0LzPzU6/lku9thBMhsHlE89Khv4XVasDXrFeyD5v79dL/Ae+nJWRWiAtGAAAAAElFTkSuQmCC",
      "text/latex": [
       "$\\displaystyle {nuR}_{i}^{T}\\,\\mathcal{C} P_L\\,{nuR}_{j}$"
      ],
      "text/plain": [
       "MajoranaBilinear(nuR[i], C⋅PL, nuR[j])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "i, j = sp.symbols(\"i j\", integer=True)\n",
    "Gp, H0 = H.components\n",
    "nuL, eL = Ll.components\n",
    "nuLbar, eLbar = Ll.bar_components\n",
    "nR, nRbar = nuR.components[0], nuR.bar_components[0]\n",
    "CPL = diracC * diracPL\n",
    "\n",
    "# Dirac Yukawa via H̃:  −Y_ν L̄ H̃ ν_R + h.c.\n",
    "LYukD = -(yv.s * sp.conjugate(H0) * Bilinear(nuLbar[i], diracPR, nR[j])\n",
    "          + yv.s * (-sp.conjugate(Gp)) * Bilinear(eLbar[i], diracPR, nR[j]))\n",
    "LYukD += -(sp.conjugate(yv.s) * H0 * Bilinear(nRbar[j], diracPL, nuL[i])\n",
    "           + sp.conjugate(yv.s) * (-Gp) * Bilinear(nRbar[j], diracPL, eL[i]))\n",
    "\n",
    "# Right-handed Majorana mass:  −½ M_R ν_Rᵀ C ν_R + h.c.\n",
    "op = -sp.Rational(1, 2) * MR.s * MajoranaBilinear(nR[i], CPL, nR[j])\n",
    "LMaj = op + sp.conjugate(op)\n",
    "\n",
    "current = fermion_gauge_current(Ll, i) + fermion_gauge_current(eR, i)\n",
    "display(MajoranaBilinear(nR[i], CPL, nR[j]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## 3. Assemble the model and check invariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "InvarianceReport(16 checks, OK)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tadpole: {mu2: lam*v**2}\n"
     ]
    }
   ],
   "source": [
    "L = Lagrangian()\n",
    "L.add((dag(Dmu(H)) * Dmu(H))[0], sector=\"kinetic\")\n",
    "L.add(-(-mu2.s * (dag(H) * H.mat)[0] + lam.s * (dag(H) * H.mat)[0]**2), sector=\"potential\")\n",
    "L.add(LYukD + current, sector=\"yukawa\")\n",
    "L.add(LMaj, sector=\"other\")      # the Majorana mass — a mass term, not a vertex\n",
    "\n",
    "model = Model(\"SM-seesaw\", gauge_groups=[SU2L, U1Y],\n",
    "              fields=[H, Ll, eR, nuR, SU2L.bosons(\"W\"), U1Y.bosons(\"B\")],\n",
    "              parameters=[gw, g1, v, lam, mu2, yv, MR], lagrangian=L)\n",
    "print(model.check_invariance())\n",
    "model.check_invariance().raise_on_failure()\n",
    "print(\"tadpole:\", model.solve_tadpoles([mu2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "Gauge invariance passes: the Dirac Yukawa is a standard singlet-completing\n",
    "coupling, and the Majorana mass is invariant because $\\nu_R$ is a total singlet.\n",
    "Now rotate the gauge bosons to the mass basis (Weinberg + $W^\\pm$), as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "# physical basis: the standard Weinberg + W± rotations, from feynlag.models\n",
    "Z, A = weinberg_rotation(model, SU2L, U1Y)\n",
    "Wp, Wm = charged_current_rotation(model, SU2L)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "## 4. The seesaw mass matrix\n",
    "\n",
    "The Dirac block comes from the Yukawa (`fermion_mass_matrix`), the Majorana\n",
    "block from `majorana_mass_matrix` (the symmetric $M_R$ — the D.2 Majorana\n",
    "machinery), and `seesaw_mass_matrix` assembles the symmetric $6N\\times6N$ block\n",
    "form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dirac mass m_D = Y_ν v/√2: sqrt(2)*v*yv/2    Majorana mass M_R: MR\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI4AAAA/CAYAAAA2a/DYAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJ1ElEQVR4Ae2c7ZHUOBCGDbUBcEsGexnAEcENGbBsBCwZQPFr+UdBBkAEfGQAFwELGXAZQJHB3fto1UajkT/GY1ueHXWVR7LUklqtV60Pa3Tj4uJiVVXVJz0x/fv8+fM/48DyfjgaUPt/VW3vJGp8ehQE3g/8eH9G7+V1Ig2ogf7rylo8N7p4Joh/pDyPo3ydkamBI8E+RwzldQYNSO/nKuax3DczFLdVEZLpW5xAYb8Iq4ETM5T32TRwqsaIrf1shQ8tqABnqOZGSCfAPFE2r8lK/ltynuEXnVw51SOFux7u3xfjFODkbYr7AoZZm5fyPzZx5AdQTE4XuUC5aYIWd14NCBi1tfElnytsFUjxUv4ThaVWNQFbHm+xOBn0LjAwLIXWBimwNpd4msgDCyDdk/8UPrkADPepdwHfGo8P/yAeS7PBE+dDmjYqFqdNOzvEqSGwFqyYUsRcxjW4RYr3jZ5wPgOQ2EtzKxu5gO2O3FdyH/h3kn/R4+ZETTwKB0jOmjXxhPnI30kFOJ0q2p5BjUMj0fCvfUPVmQQN17j9IR4a+oGeu3XCq/0UwEU4gHIgk/uRd893LHeDR2HIY9YsyRPl47Nrdo6ao+aJkcDW836oRCaCTBJNEfMIMXIpkh9QfJZLg2Fd3DDii6G+4bsPvnKUBusBz135awskv9OJ3DPFARZHegdkblOuiUfx97p4wnxcxh0/WS2OhGXV8EXuUz2YYBT6SX5nejtk34foFxLyierDMFN5lyFsY2PNx1NvdMH855ceeGNdAMbwE9FKPLH1qsFEviLSxDx98nGJUz/ZgKPKMv7fkhv2HnoY7yxF95583bAUWB2o0dqIF4BQb4Y35jI0Ph0p/vTDuwvzaZIgDHgACHqO+frmo+SblA04EoUZflwZJGSyRy9yvZSAPSca/9wD4TjRgFY9rC+NjGsP6erhyjMydzpTOB0vZW1go0w+Y8CDP6XnPvkoaZqO0sGzhKKk1PcZm98QX1ujWSSaoBA13kc9WJp/9IST3bXSxPPHWkDDi/jc/KkhulI8He6bXKc7uXTQDQuu8NZ8mvK38CwWx1fOZGhymf1fF6LXM1m2TjFlvd4qc1Zelcqj82HlUh0UlsGUy+IYKGIzHFbkugxVNCC9fy7riXVhUs3O9G25jVYuVPa2/lzA6SPn7T5MhWddAwJKvHpaZxjpLRdwmNE3kVmjH00MSw5Xw3Ueylqy/MimOnQeGssCHAnGHgUypoYjC5tjPoAMo1IfpY9aYKbMskyOfV0xqSeJepvFmcXkJsovQT00kMXieLk+yGWZGhOTOZaTbRPnOM21fpcusMK2iWidLeshr2wWR8pgifhTrls60vJeQQ/l5ZB0od8a4PsdnyJ42JdhjsgmYTbKBhxfY6wL32VQDMtH9iD+lj+10+mTHKTDDvIqqDmWmiX3nSBsVm/OoapSxRmO2Pou1K4BdHTZzjJvbFbgzFvV/S1NHSze+QVI9SGvHDXLPVTlqPNel+mHp/iQ1+x1KsCZXeXDCxRoTpSa+c3aIa/hOQ5PWYaq4bqbNaUHjTvkRcH+HTfLRmkBzqzNP6wwDxI+XgIcW0kxz+GrexYqwMmi9q0LZc+GTcC1vRuBKNuKtABn6zacP4EA0uuQ15ySTQocVbjzS7F4Or/EzqmQUlY/DUwGHAGC866LvL6jn2oKV5sGJgOOCt3L6zvalFXifmtgEuDI2tR/qJd/cV92f1e/+IZqYBLgSJjwD/V7dX3HUEUeWrrRd45Da+OVubgvu4fWyFPUd1SL44el0Nogc+eXXaVbiY+NreT1Hal4Mlb4qFd3kGehfhrY2uKosTgHwoopRZxSWzvVJ96dru9QeuZIG9d7KBygAbiqiUdR9RUg8BUaTwNbWRw1EA3F1WP8kf69Hs7TOAoar3EbXDw0dvxl91hhG1dziJd/QHLLQjJe4chiZ1SSPEEeYl0OSS46Hif5HPDxI2ubhIpnwUGnROecx34Xp/E8tI/lCx/8dD4j/pveWpYxtrlHbZFxnApEkFmv71CZ7iOe3DOVHVYYULVe76E0ADW82UGv+Ulycb6GzsK5azoSX70bydeD+kKALHmQX+Hc+PFKLhuvtJPdL0i6Su+U810unRLgDqathypf0gu5c1/fEYOAXhUrkLAQKCspKObxVcjr+EZEVixC1wWRf5m0XfVRPHqCQj24AMXRCdEHN3q1gtUlaPkZBBwVSs9HCDt5jwlNDlFeQL7sjnF9h/sjn/IEIKNf3aE85yQDPnqsgRELoLoyrL3XY/wxS/wOH9TUYRjWd6ZBwPGlApRzVQyEz3F9B+VNenWHr9dcDoexAM2lnmTvVzzh8NDYzFM2rIjCYmJIr5R248C/wmgrnp0vQDiikCEkIWa7vkNlobTJr+4YoocR0nxXHuHkNcySIYU5i61im6xImMYBIwzArzwAIXOqnec35DcYOCQWYQXOJBS9Ykp6q8xZWqNETPEkV3dMWYEwb9+IdrbG6Y6wUI/yM2m2Q+pukquwDSsS5WvDVCVeVmEQlzcAGoZ5JtateYinF+0EHAnBXCdc6fQqdAATc6TJr+4YINfQJDSwWQ/rdFgKA5GzQNLvL19AyN9Wpq2iGNLDvOh4Nuy1pe8d1xs4EqTzbE3vUndklCzWm7bKSemWcvaH+Y2zJnIZgqkHDWvE3JGldSWXcIDUZ34DwLjQwcBIesCHpaH9WMwwSuxM2wBnKUrfudILzIDGdUtyNTBDVGjFAQNkFurqLf2L1QrTOi7lCfCgEJxXIQN/bw5MV5IN1IAakca1+Y3lgoVgKHYNLLe2GApn+MGKtM5NFG8AS1kmi2OeMwr1tjijlBZl4hVle0HWG7LewhCJOMUrjRhbBYBCeD1EBQUT3sfa2E4wez4x2a4zKzhH0j3f/1rBaLwpNytwJNAhntVh4urmL0GDsGKMh6iKxlU4VihlRerk4qPTPawDNj01YIgSP2Ak38HAyT1U0cOohBE70JhsFHatSHXiCAgNSP2+6rG9GepJA/KfKTdEyV3pYc+FByKO9DZXuQrVr8LQGUdXmGwDyGcKC/OGx30bUxxbJywskCG2egruT7ktTudZnf5VWTanGsqGkg1BFcdQVA9H8ftGgiBAvL1WSeJD16NRVuCoMrbBZRWicllvYTBBituugdxDVS2dQMTwxDjPZUuFFq6BRQBHoDmRnhirs9/CsPD2Wox4WYcqtOBBw+TPvscAIsLDvQyCCi1IA1mB40HDdyiAYysp5jm9JnwL0uPBiZIVONL24m5hODgEDKzwjYuLi5XSpjaYWN10HWkcWGxJtg8aUPvTsW0kCEU+xeJc6rHP8WHkaN81wkyLf6808EjSHickvvwfuBjfZ1w9RucAAAAASUVORK5CYII=",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & \\frac{\\sqrt{2} v yv}{2}\\\\\\frac{\\sqrt{2} v yv}{2} & MR\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡         √2⋅v⋅yv⎤\n",
       "⎢   0     ───────⎥\n",
       "⎢            2   ⎥\n",
       "⎢                ⎥\n",
       "⎢√2⋅v⋅yv         ⎥\n",
       "⎢───────    MR   ⎥\n",
       "⎣   2            ⎦"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mD = fermion_mass_matrix(LYukD, nuLbar, nR, model.vacuum, 1, (i, j), gamma=diracPR)\n",
    "MRmat = majorana_mass_matrix(LMaj, nR, model.vacuum, 1, (i, j), gamma=CPL)\n",
    "Mnu = seesaw_mass_matrix(mD, MRmat)\n",
    "print(\"Dirac mass m_D = Y_ν v/√2:\", mD[0, 0], \"   Majorana mass M_R:\", MRmat[0, 0])\n",
    "Mnu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "## 5. The seesaw formula and the spectrum\n",
    "\n",
    "Block-diagonalising for $M_R \\gg m_D$ gives the light neutrino mass\n",
    "$m_\\nu \\approx -m_D M_R^{-1} m_D^T = -m_D^2/M_R$ — the **seesaw**: heavy $M_R$ in\n",
    "the denominator makes $m_\\nu$ tiny.  We check it against the *exact* Takagi\n",
    "factorisation of the full matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seesaw formula  m_ν = −m_D²/M_R  = -v**2*yv**2/(2*MR)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact spectrum:  m_light = 3.0258e-03 GeV,  m_heavy = 1000.0 GeV\n",
      "seesaw formula:  m_light = 3.0258e-03 GeV\n"
     ]
    }
   ],
   "source": [
    "print(\"seesaw formula  m_ν = −m_D²/M_R  =\", sp.simplify(seesaw_light_mass(mD, MRmat)[0, 0]))\n",
    "\n",
    "bench = {yv.s: 0.01, v.s: 246.0, MR.s: 1000.0, gw.s: 0.6535, g1.s: 0.3580}\n",
    "Mn = sp.Matrix(2, 2, lambda a, b: sp.nsimplify(Mnu[a, b].subs(bench)))\n",
    "U, D = diagonalize_takagi(Mn)     # M = U D Uᵀ, D ≥ 0 (Takagi, for Majorana)\n",
    "masses = [abs(complex(sp.N(D[k, k]))) for k in range(2)]\n",
    "light = 0 if masses[0] < masses[1] else 1\n",
    "heavy = 1 - light\n",
    "print(f\"exact spectrum:  m_light = {masses[light]:.4e} GeV,  m_heavy = {masses[heavy]:.1f} GeV\")\n",
    "print(f\"seesaw formula:  m_light = {abs(complex(sp.N(seesaw_light_mass(mD, MRmat)[0,0].subs(bench)))):.4e} GeV\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "The exact light mass and the seesaw formula agree. The spectrum is one light\n",
    "state and one heavy state $\\approx M_R$.\n",
    "\n",
    "> **The seesaw tension.** Here $M_R = 1$ TeV gives $m_\\nu \\approx 3$ MeV — far\n",
    "> too heavy.  A *realistic* $m_\\nu \\sim 0.05$ eV needs $M_R \\sim 10^{14}$ GeV (with\n",
    "> $Y_\\nu \\sim 1$), which pushes the heavy neutrino — and its couplings — out of\n",
    "> experimental reach.  We keep $M_R = 1$ TeV here so the mixing below is visible."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "## 6. The light–heavy mixing\n",
    "\n",
    "The physical neutrinos are Majorana mass eigenstates $\\chi_k$ mixing $\\nu_L$ with\n",
    "$\\nu_R^c$.  The **light–heavy mixing** $V$ — the $\\nu_L$ content of the heavy\n",
    "state — is the $(0,\\text{heavy})$ entry of $U^*$, and equals $\\approx m_D/M_R$.\n",
    "It controls every heavy-neutrino coupling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "light–heavy mixing V = |U*[0,heavy]| = 0.00174\n",
      "           m_D/M_R               = 0.00174\n"
     ]
    }
   ],
   "source": [
    "V = abs(complex(sp.N(U.conjugate()[0, heavy])))\n",
    "print(f\"light–heavy mixing V = |U*[0,heavy]| = {V:.5f}\")\n",
    "print(f\"           m_D/M_R               = {float((mD[0,0]/MR.s).subs(bench)):.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## 7. Rotate to the physical Majorana basis\n",
    "\n",
    "`MajoranaRotation` substitutes the weak neutrinos into the physical Majorana\n",
    "fields $\\chi_k$ using the Takagi $U$.  Because $\\chi$ is self-conjugate\n",
    "($\\chi^c=\\chi$), the four Weyl fields map as\n",
    "\n",
    "$$\\nu_L\\to U^*\\,P_L\\chi,\\quad \\bar\\nu_L\\to U\\,\\bar\\chi P_R,\\quad\n",
    "  \\nu_R\\to U\\,P_R\\chi,\\quad \\bar\\nu_R\\to U^*\\,\\bar\\chi P_L,$$\n",
    "\n",
    "which is exactly the charge-conjugation bookkeeping feynlag's `Rotation` cannot\n",
    "express (it is square and conjugation-blind).  We apply it to the current sector\n",
    "and extract the physical vertices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "extracted 22 neutrino vertex groups in the physical basis\n"
     ]
    }
   ],
   "source": [
    "chiL, chiR = sp.IndexedBase(\"chiL\"), sp.IndexedBase(\"chiR\")\n",
    "chiLbar, chiRbar = sp.IndexedBase(\"chiLbar\"), sp.IndexedBase(\"chiRbar\")\n",
    "rot = MajoranaRotation(U, nuL, nR, nuLbar, nRbar, chiL, chiR, chiLbar, chiRbar, n_L=1)\n",
    "Lchi = rot.apply(model.physical_lagrangian(sector=\"yukawa\"), (i, j), 1)\n",
    "tab = extract_fermion_vertices(Lchi, [Wp, Wm, Z])\n",
    "print(\"extracted\", len(tab), \"neutrino vertex groups in the physical basis\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "## 8. The heavy-neutrino couplings\n",
    "\n",
    "The smoking gun: the charged-current **production** coupling\n",
    "$W\\,\\bar\\ell\\,N = (g/\\sqrt2)\\,V$ — the heavy neutrino couples to $W$ only through\n",
    "the mixing $V$.  The light neutrino keeps the full SM coupling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  W⁻ e light ν:  |coupling|/(g/√2) = 1.00000\n",
      "  W⁻ e heavy N:  |coupling|/(g/√2) = 0.00174\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Z ν light ν:   |coupling|/g_Z    = 0.50000\n",
      "  Z ν heavy N:   |coupling|/g_Z    = 0.00087\n"
     ]
    }
   ],
   "source": [
    "mu = sp.Symbol(\"mu\", integer=True)\n",
    "gL = DiracGamma(mu) * diracPL\n",
    "gsq = float((gw.s / sp.sqrt(2)).subs(bench))\n",
    "gz = float(sp.sqrt(gw.s**2 + g1.s**2).subs(bench))\n",
    "\n",
    "def coup(key, boson, norm):\n",
    "    c = tab.get(key, {}).get(1, {}).get((boson,), 0)\n",
    "    return abs(complex(sp.N(sp.simplify(c).subs(bench)))) / norm if c != 0 else 0.0\n",
    "\n",
    "for k, name in ((light, \"light ν\"), (heavy, \"heavy N\")):\n",
    "    cW = coup((eLbar[0], gL, chiL[k]), Wm, gsq)\n",
    "    print(f\"  W⁻ e {name}:  |coupling|/(g/√2) = {cW:.5f}\")\n",
    "for k, name in ((light, \"light ν\"), (heavy, \"heavy N\")):\n",
    "    cZ = coup((chiLbar[light], gL, chiL[k]), Z, gz)\n",
    "    print(f\"  Z ν {name}:   |coupling|/g_Z    = {cZ:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "The heavy-neutrino couplings are suppressed by the mixing $V\\approx m_D/M_R$:\n",
    "$W\\bar eN = (g/\\sqrt2)\\,V$, $Z\\bar\\nu N = (g_Z/2)\\,V$ (up to the neutral-current\n",
    "$\\tfrac12$).  The light neutrino keeps the SM values $g/\\sqrt2$ and $g_Z/2$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "## 9. Decoupling\n",
    "\n",
    "Raising $M_R$ drives the mixing — and the heavy-neutrino production coupling — to\n",
    "zero as $1/M_R$: the seesaw *decouples*.  This is the same computation at three\n",
    "scales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  M_R = 1e+03 GeV:  m_light = 3.026e-03 GeV,  W e N /(g/√2) = 1.739e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  M_R = 1e+04 GeV:  m_light = 3.026e-04 GeV,  W e N /(g/√2) = 1.739e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  M_R = 1e+05 GeV:  m_light = 3.026e-05 GeV,  W e N /(g/√2) = 1.739e-05\n"
     ]
    }
   ],
   "source": [
    "for MRval in (1.0e3, 1.0e4, 1.0e5):\n",
    "    b = dict(bench); b[MR.s] = MRval\n",
    "    Mnk = sp.Matrix(2, 2, lambda a, bb: sp.nsimplify(Mnu[a, bb].subs(b)))\n",
    "    Uk, Dk = diagonalize_takagi(Mnk)\n",
    "    mk = [abs(complex(sp.N(Dk[c, c]))) for c in range(2)]\n",
    "    hv = 0 if mk[0] > mk[1] else 1\n",
    "    rk = MajoranaRotation(Uk, nuL, nR, nuLbar, nRbar, chiL, chiR, chiLbar, chiRbar, n_L=1)\n",
    "    tk = extract_fermion_vertices(rk.apply(model.physical_lagrangian(sector=\"yukawa\"), (i, j), 1), [Wp, Wm, Z])\n",
    "    cW = tk.get((eLbar[0], gL, chiL[hv]), {}).get(1, {}).get((Wm,), 0)\n",
    "    val = abs(complex(sp.N(sp.simplify(cW).subs(b)))) / gsq if cW != 0 else 0.0\n",
    "    print(f\"  M_R = {MRval:.0e} GeV:  m_light = {min(mk):.3e} GeV,  W e N /(g/√2) = {val:.3e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25",
   "metadata": {},
   "source": [
    "## 10. Recap\n",
    "\n",
    "Starting from the SM plus right-handed neutrinos, feynlag delivered the type-I\n",
    "seesaw end to end:\n",
    "\n",
    "| Step | Tool | Result |\n",
    "|------|------|--------|\n",
    "| Majorana mass | `MajoranaBilinear`, `majorana_mass_matrix` | the symmetric $M_R$ block |\n",
    "| Mass matrix | `seesaw_mass_matrix`, `diagonalize_takagi` | light + heavy spectrum |\n",
    "| Seesaw formula | `seesaw_light_mass` | $m_\\nu = -m_D^2/M_R$, matching the exact Takagi |\n",
    "| Physical basis | `MajoranaRotation` | $\\nu_L,\\nu_R \\to$ Majorana $\\chi$ (with $\\nu_R^c$) |\n",
    "| Couplings | `extract_fermion_vertices` | $W\\bar\\ell N=(g/\\sqrt2)V$, decoupling as $M_R\\to\\infty$ |\n",
    "\n",
    "The heavy neutrino's couplings are fixed by a single number, the light–heavy\n",
    "mixing $V \\approx m_D/M_R$ — which is also what makes a high-scale seesaw\n",
    "experimentally invisible.\n",
    "\n",
    "**Where to go next:** `SM_Feynman_Rules_Tutorial` for the full SM pipeline,\n",
    "`SM_U1X_Tutorial` for a $Z'$, and `ModelBuilding_Tutorial` for the\n",
    "`suggest`/`anomalies` model-building tools."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (lagrangian)",
   "language": "python",
   "name": "lagrangian"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}