{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Vector-like fermions in `feynlag`: SM + a vector-like lepton doublet\n",
    "\n",
    "A **vector-like fermion** has both chiralities in the *same* gauge\n",
    "representation, so — unlike a chiral SM fermion — it can carry a bare,\n",
    "gauge-invariant Dirac mass all by itself. Mixed with an ordinary SM fermion\n",
    "through the Higgs, it produces a textbook flavor-mixing problem: a 2×2 mass\n",
    "matrix diagonalized by *two independent* angles (θ_L, θ_R, one per\n",
    "chirality), and — because the heavy and light states are no longer pure\n",
    "gauge eigenstates — flavor-changing neutral currents (FCNCs) that don't\n",
    "exist in the SM.\n",
    "\n",
    "This notebook walks through `examples/sm_vll.py` stage by stage, reusing\n",
    "the same declare → write → check → break-symmetry → diagonalize → extract\n",
    "pipeline as **[`SM_Feynman_Rules_Tutorial.ipynb`](SM_Feynman_Rules_Tutorial.ipynb)**\n",
    "(read that one first if `ExternalParameter`, `Dmu`, `Scalar`, or `Model`\n",
    "are unfamiliar — this notebook won't re-explain them). What's new here:\n",
    "\n",
    "- a **vector-like doublet** Ψ = (N, E), built as two `WeylFermion`s with\n",
    "  identical reps (not a `DiracFermion` — see §2),\n",
    "- a **bare mass term** `-M Ψ̄_L Ψ_R` that needs no scalar VEV to exist,\n",
    "- **biunitary (2×2 SVD) diagonalization** with two physically distinct\n",
    "  angles, via the analytic `diagonalize_svd_2x2` helper,\n",
    "- **rotating fermion fields** into their mass basis and re-extracting\n",
    "  Feynman rules — which exercises `expand_bilinear`, a piece of machinery\n",
    "  this model is the first to need.\n",
    "\n",
    "The headline physics we'll derive and verify symbolically:\n",
    "\n",
    "1. the **Z-FCNC lives only in the right-handed current** (the left-handed\n",
    "   Z coupling is exactly protected),\n",
    "2. a clean **Higgs-coupling sum rule** for the deviation from `-m_e/v`,\n",
    "3. a **new right-handed W current** with no SM analogue.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "from feynlag import (\n",
    "    Bilinear, DiracGamma, Dmu, ExternalParameter, InternalParameter,\n",
    "    Lagrangian, Model, Rotation, SU2, Scalar, U1, WeylFermion, dag,\n",
    "    diagonalize_svd_2x2, diracPL, diracPR, expand_bilinear,\n",
    "    extract_fermion_vertices, fermion_gauge_current, fermion_mass_matrix,\n",
    "    rotation_2x2, weinberg_rotation, charged_current_rotation,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## 1. Symmetries and parameters\n",
    "\n",
    "Same `SU2L`/`U1Y` electroweak setup as the plain-SM tutorial, plus three\n",
    "new couplings for the VLL sector:\n",
    "\n",
    "- `ye` — the ordinary (tiny) SM electron Yukawa,\n",
    "- `lamE` — the mixing Yukawa between the heavy doublet and `e_R`,\n",
    "- `MPsi` — the bare vector-like mass. It needs `unit_dim=1` (mass\n",
    "  dimension 1) so `Model.check_invariance()`'s power-counting check can\n",
    "  see that `MPsi · (dimension-3 Bilinear)` is a dimension-4 operator —\n",
    "  exactly like `v` needs `unit_dim=1` for the Yukawa terms below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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",
    "\n",
    "v = ExternalParameter(\"v\", 246.0, positive=True, unit_dim=1)\n",
    "lam = ExternalParameter(\"lam\", 0.129, tex=\"\\\\lambda\")\n",
    "mu2 = InternalParameter(\"mu2\", unit_dim=2, tex=\"\\\\mu^2\")\n",
    "\n",
    "ye = ExternalParameter(\"ye\", 0.01, positive=True, tex=\"y_e\")\n",
    "lamE = ExternalParameter(\"lamE\", 0.4, positive=True, tex=\"\\\\lambda_E\")\n",
    "MPsi = ExternalParameter(\"MPsi\", 1000.0, positive=True, unit_dim=1, tex=\"M_\\\\Psi\")\n",
    "MPsi.unit_dim\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## 2. Fields\n",
    "\n",
    "The Higgs and the SM lepton generation are declared exactly as in the SM\n",
    "tutorial (one flavor: `nflavors=1`, so the mixing physics isn't obscured\n",
    "by 3-generation flavor indices).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([nuL, eL], [eR])"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = Scalar(\"H\", reps={SU2L: 2, U1Y: sp.Rational(1, 2)},\n",
    "           component_names=[\"Gp\", \"H0\"])\n",
    "H.expand_vev({H.components[1]: v})\n",
    "\n",
    "Ll = WeylFermion(\"Ll\", reps={SU2L: 2, U1Y: -sp.Rational(1, 2)},\n",
    "                 chirality=\"L\", nflavors=1, component_names=[\"nuL\", \"eL\"])\n",
    "eRf = WeylFermion(\"eRf\", reps={U1Y: -1}, chirality=\"R\", nflavors=1,\n",
    "                  component_names=[\"eR\"])\n",
    "Ll.components, eRf.components\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "Now the vector-like doublet. `PsiL` and `PsiR` are given the **same**\n",
    "reps `(SU2L: 2, U1Y: -1/2)` as `Ll` — both chiralities in one\n",
    "representation is precisely what makes a bare mass gauge invariant (a\n",
    "chiral SM fermion can't have one: `L_L` and `e_R` transform differently,\n",
    "so nothing gauge-invariant connects them without a scalar).\n",
    "\n",
    "Why two `WeylFermion`s and not a single `DiracFermion`? `DiracFermion`'s\n",
    "`chirality=None` default routes through `fermion_gauge_current` with the\n",
    "identity projector, producing a bare `DiracGamma(mu)*diracI` structure\n",
    "that `dirac_conjugate` doesn't know how to Dirac-conjugate — it raises\n",
    "`NotImplementedError`, silently breaking the hermiticity check. That path\n",
    "is untested (see `fields.py`'s `DiracFermion` docstring); modeling the VLL\n",
    "as two `WeylFermion`s with identical reps sidesteps it entirely and is a\n",
    "well-exercised pattern.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([NL, EL], [NR, ER])"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PsiL = WeylFermion(\"PsiL\", reps={SU2L: 2, U1Y: -sp.Rational(1, 2)},\n",
    "                   chirality=\"L\", nflavors=1, component_names=[\"NL\", \"EL\"])\n",
    "PsiR = WeylFermion(\"PsiR\", reps={SU2L: 2, U1Y: -sp.Rational(1, 2)},\n",
    "                   chirality=\"R\", nflavors=1, component_names=[\"NR\", \"ER\"])\n",
    "W, B = SU2L.bosons(\"W\"), U1Y.bosons(\"B\")\n",
    "PsiL.components, PsiR.components\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gp, H0 = H.components\n",
    "nuL, eL = Ll.components\n",
    "nuLbar, eLbar = Ll.bar_components\n",
    "eR, eRbar = eRf.components[0], eRf.bar_components[0]\n",
    "NL, EL = PsiL.components\n",
    "NLbar, ELbar = PsiL.bar_components\n",
    "NR, ER = PsiR.components\n",
    "NRbar, ERbar = PsiR.bar_components\n",
    "\n",
    "i = sp.Symbol(\"fl_i\", integer=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "## 3. The Lagrangian\n",
    "\n",
    "Three Yukawa-sector pieces, built up one at a time so the `+ h.c.`\n",
    "pattern is visible at each step (`Bilinear(bar, gamma, field)` is the\n",
    "opaque fermion-sandwich atom; `diracPR`/`diracPL` are the chiral\n",
    "projectors — see the SM tutorial for the full explanation of this\n",
    "convention).\n",
    "\n",
    "### SM Yukawa: `-y_e L̄_L H e_R + h.c.`\n",
    "\n",
    "Identical in form to the plain-SM lepton Yukawa.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - (Gp ye \\bar{{\\bar{nuL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}} + H_{0} ye \\bar{{\\bar{eL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}}) - \\left(ye \\overline{Gp} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{nuL}_{fl_{i}} + ye \\overline{H_{0}} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{eL}_{fl_{i}}\\right)$"
      ],
      "text/plain": [
       "-(Gp*ye*Bilinear(nuLbar[fl_i], PR, eR[fl_i]) + H0*ye*Bilinear(eLbar[fl_i], PR, eR[fl_i])) - (ye*conjugate(Gp)*Bilinear(eRbar[fl_i], PL, nuL[fl_i]) + ye*conjugate(H0)*Bilinear(eRbar[fl_i], PL, eL[fl_i]))"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LYuk = -(ye.s * Gp * Bilinear(nuLbar[i], diracPR, eR[i])\n",
    "         + ye.s * H0 * Bilinear(eLbar[i], diracPR, eR[i]))\n",
    "LYuk += -(ye.s * sp.conjugate(Gp) * Bilinear(eRbar[i], diracPL, nuL[i])\n",
    "          + ye.s * sp.conjugate(H0) * Bilinear(eRbar[i], diracPL, eL[i]))\n",
    "LYuk\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "### Mixing Yukawa: `-λ_E Ψ̄_L H e_R + h.c.`\n",
    "\n",
    "`Ψ_L` carries exactly `L_L`'s reps, so the component pattern is identical\n",
    "to the SM Yukawa above — just with `NL`/`EL` in place of `nuL`/`eL`. This\n",
    "single term is what ties the heavy doublet to the light SM lepton.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - (Gp lamE \\bar{{\\bar{NL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}} + H_{0} lamE \\bar{{\\bar{EL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}}) - \\left(Gp ye \\bar{{\\bar{nuL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}} + H_{0} ye \\bar{{\\bar{eL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}}\\right) - \\left(lamE \\overline{Gp} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{NL}_{fl_{i}} + lamE \\overline{H_{0}} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{EL}_{fl_{i}}\\right) - \\left(ye \\overline{Gp} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{nuL}_{fl_{i}} + ye \\overline{H_{0}} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{eL}_{fl_{i}}\\right)$"
      ],
      "text/plain": [
       "-(Gp*lamE*Bilinear(NLbar[fl_i], PR, eR[fl_i]) + H0*lamE*Bilinear(ELbar[fl_i], PR, eR[fl_i])) - (Gp*ye*Bilinear(nuLbar[fl_i], PR, eR[fl_i]) + H0*ye*Bilinear(eLbar[fl_i], PR, eR[fl_i])) - (lamE*conjugate(Gp)*Bilinear(eRbar[fl_i], PL, NL[fl_i]) + lamE*conjugate(H0)*Bilinear(eRbar[fl_i], PL, EL[fl_i])) - (ye*conjugate(Gp)*Bilinear(eRbar[fl_i], PL, nuL[fl_i]) + ye*conjugate(H0)*Bilinear(eRbar[fl_i], PL, eL[fl_i]))"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LYuk += -(lamE.s * Gp * Bilinear(NLbar[i], diracPR, eR[i])\n",
    "          + lamE.s * H0 * Bilinear(ELbar[i], diracPR, eR[i]))\n",
    "LYuk += -(lamE.s * sp.conjugate(Gp) * Bilinear(eRbar[i], diracPL, NL[i])\n",
    "          + lamE.s * sp.conjugate(H0) * Bilinear(eRbar[i], diracPL, EL[i]))\n",
    "LYuk\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "### Bare vector-like mass: `-M_Ψ Ψ̄_L Ψ_R + h.c.`\n",
    "\n",
    "No Higgs field appears — this term needs no VEV to exist, which is the\n",
    "whole point of \"vector-like.\" Both `N` and `E` get this same mass `M_Ψ`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - MPsi \\left(\\bar{{\\bar{EL}}_{fl_{i}}}\\,P_R\\,{ER}_{fl_{i}} + \\bar{{\\bar{NL}}_{fl_{i}}}\\,P_R\\,{NR}_{fl_{i}}\\right) - MPsi \\left(\\bar{{\\bar{ER}}_{fl_{i}}}\\,P_L\\,{EL}_{fl_{i}} + \\bar{{\\bar{NR}}_{fl_{i}}}\\,P_L\\,{NL}_{fl_{i}}\\right) - \\left(Gp lamE \\bar{{\\bar{NL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}} + H_{0} lamE \\bar{{\\bar{EL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}}\\right) - \\left(Gp ye \\bar{{\\bar{nuL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}} + H_{0} ye \\bar{{\\bar{eL}}_{fl_{i}}}\\,P_R\\,{eR}_{fl_{i}}\\right) - \\left(lamE \\overline{Gp} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{NL}_{fl_{i}} + lamE \\overline{H_{0}} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{EL}_{fl_{i}}\\right) - \\left(ye \\overline{Gp} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{nuL}_{fl_{i}} + ye \\overline{H_{0}} \\bar{{\\bar{eR}}_{fl_{i}}}\\,P_L\\,{eL}_{fl_{i}}\\right)$"
      ],
      "text/plain": [
       "-MPsi*(Bilinear(ELbar[fl_i], PR, ER[fl_i]) + Bilinear(NLbar[fl_i], PR, NR[fl_i])) - MPsi*(Bilinear(ERbar[fl_i], PL, EL[fl_i]) + Bilinear(NRbar[fl_i], PL, NL[fl_i])) - (Gp*lamE*Bilinear(NLbar[fl_i], PR, eR[fl_i]) + H0*lamE*Bilinear(ELbar[fl_i], PR, eR[fl_i])) - (Gp*ye*Bilinear(nuLbar[fl_i], PR, eR[fl_i]) + H0*ye*Bilinear(eLbar[fl_i], PR, eR[fl_i])) - (lamE*conjugate(Gp)*Bilinear(eRbar[fl_i], PL, NL[fl_i]) + lamE*conjugate(H0)*Bilinear(eRbar[fl_i], PL, EL[fl_i])) - (ye*conjugate(Gp)*Bilinear(eRbar[fl_i], PL, nuL[fl_i]) + ye*conjugate(H0)*Bilinear(eRbar[fl_i], PL, eL[fl_i]))"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LYuk += -MPsi.s * (Bilinear(NLbar[i], diracPR, NR[i])\n",
    "                   + Bilinear(ELbar[i], diracPR, ER[i]))\n",
    "LYuk += -MPsi.s * (Bilinear(NRbar[i], diracPL, NL[i])\n",
    "                   + Bilinear(ERbar[i], diracPL, EL[i]))\n",
    "LYuk\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "Finally the gauge currents (`fermion_gauge_current` is fully generic —\n",
    "one call per fermion, regardless of how many gauge groups it's charged\n",
    "under):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{B g_{1} \\bar{{\\bar{EL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{EL}_{fl_{i}}}{2} - \\frac{B g_{1} \\bar{{\\bar{ER}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{ER}_{fl_{i}}}{2} - \\frac{B g_{1} \\bar{{\\bar{NL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{NL}_{fl_{i}}}{2} - \\frac{B g_{1} \\bar{{\\bar{NR}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{NR}_{fl_{i}}}{2} - \\frac{B g_{1} \\bar{{\\bar{eL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{eL}_{fl_{i}}}{2} - B g_{1} \\bar{{\\bar{eR}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{eR}_{fl_{i}} - \\frac{B g_{1} \\bar{{\\bar{nuL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{nuL}_{fl_{i}}}{2} + \\frac{W_{1} gw \\bar{{\\bar{EL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{NL}_{fl_{i}}}{2} + \\frac{W_{1} gw \\bar{{\\bar{ER}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{NR}_{fl_{i}}}{2} + \\frac{W_{1} gw \\bar{{\\bar{NL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{EL}_{fl_{i}}}{2} + \\frac{W_{1} gw \\bar{{\\bar{NR}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{ER}_{fl_{i}}}{2} + \\frac{W_{1} gw \\bar{{\\bar{eL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{nuL}_{fl_{i}}}{2} + \\frac{W_{1} gw \\bar{{\\bar{nuL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{eL}_{fl_{i}}}{2} + \\frac{i W_{2} gw \\bar{{\\bar{EL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{NL}_{fl_{i}}}{2} + \\frac{i W_{2} gw \\bar{{\\bar{ER}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{NR}_{fl_{i}}}{2} - \\frac{i W_{2} gw \\bar{{\\bar{NL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{EL}_{fl_{i}}}{2} - \\frac{i W_{2} gw \\bar{{\\bar{NR}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{ER}_{fl_{i}}}{2} + \\frac{i W_{2} gw \\bar{{\\bar{eL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{nuL}_{fl_{i}}}{2} - \\frac{i W_{2} gw \\bar{{\\bar{nuL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{eL}_{fl_{i}}}{2} - \\frac{W_{3} gw \\bar{{\\bar{EL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{EL}_{fl_{i}}}{2} - \\frac{W_{3} gw \\bar{{\\bar{ER}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{ER}_{fl_{i}}}{2} + \\frac{W_{3} gw \\bar{{\\bar{NL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{NL}_{fl_{i}}}{2} + \\frac{W_{3} gw \\bar{{\\bar{NR}}_{fl_{i}}}\\,\\gamma^{\\mu} P_R\\,{NR}_{fl_{i}}}{2} - \\frac{W_{3} gw \\bar{{\\bar{eL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{eL}_{fl_{i}}}{2} + \\frac{W_{3} gw \\bar{{\\bar{nuL}}_{fl_{i}}}\\,\\gamma^{\\mu} P_L\\,{nuL}_{fl_{i}}}{2}$"
      ],
      "text/plain": [
       "-B*g1*Bilinear(ELbar[fl_i], gamma(mu)*PL, EL[fl_i])/2 - B*g1*Bilinear(ERbar[fl_i], gamma(mu)*PR, ER[fl_i])/2 - B*g1*Bilinear(NLbar[fl_i], gamma(mu)*PL, NL[fl_i])/2 - B*g1*Bilinear(NRbar[fl_i], gamma(mu)*PR, NR[fl_i])/2 - B*g1*Bilinear(eLbar[fl_i], gamma(mu)*PL, eL[fl_i])/2 - B*g1*Bilinear(eRbar[fl_i], gamma(mu)*PR, eR[fl_i]) - B*g1*Bilinear(nuLbar[fl_i], gamma(mu)*PL, nuL[fl_i])/2 + W_1*gw*Bilinear(ELbar[fl_i], gamma(mu)*PL, NL[fl_i])/2 + W_1*gw*Bilinear(ERbar[fl_i], gamma(mu)*PR, NR[fl_i])/2 + W_1*gw*Bilinear(NLbar[fl_i], gamma(mu)*PL, EL[fl_i])/2 + W_1*gw*Bilinear(NRbar[fl_i], gamma(mu)*PR, ER[fl_i])/2 + W_1*gw*Bilinear(eLbar[fl_i], gamma(mu)*PL, nuL[fl_i])/2 + W_1*gw*Bilinear(nuLbar[fl_i], gamma(mu)*PL, eL[fl_i])/2 + I*W_2*gw*Bilinear(ELbar[fl_i], gamma(mu)*PL, NL[fl_i])/2 + I*W_2*gw*Bilinear(ERbar[fl_i], gamma(mu)*PR, NR[fl_i])/2 - I*W_2*gw*Bilinear(NLbar[fl_i], gamma(mu)*PL, EL[fl_i])/2 - I*W_2*gw*Bilinear(NRbar[fl_i], gamma(mu)*PR, ER[fl_i])/2 + I*W_2*gw*Bilinear(eLbar[fl_i], gamma(mu)*PL, nuL[fl_i])/2 - I*W_2*gw*Bilinear(nuLbar[fl_i], gamma(mu)*PL, eL[fl_i])/2 - W_3*gw*Bilinear(ELbar[fl_i], gamma(mu)*PL, EL[fl_i])/2 - W_3*gw*Bilinear(ERbar[fl_i], gamma(mu)*PR, ER[fl_i])/2 + W_3*gw*Bilinear(NLbar[fl_i], gamma(mu)*PL, NL[fl_i])/2 + W_3*gw*Bilinear(NRbar[fl_i], gamma(mu)*PR, NR[fl_i])/2 - W_3*gw*Bilinear(eLbar[fl_i], gamma(mu)*PL, eL[fl_i])/2 + W_3*gw*Bilinear(nuLbar[fl_i], gamma(mu)*PL, nuL[fl_i])/2"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "current = (fermion_gauge_current(Ll, i) + fermion_gauge_current(eRf, i)\n",
    "           + fermion_gauge_current(PsiL, i)\n",
    "           + fermion_gauge_current(PsiR, i))\n",
    "current\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {},
   "source": [
    "## 4. Assemble the Model and check invariance\n",
    "\n",
    "The kinetic/potential Higgs sector is the standard SM piece. Assembling\n",
    "everything into a `Model` and calling `check_invariance()` is the\n",
    "correctness gate: this is the first model in `feynlag` where a bare\n",
    "fermion mass *and* a mixing Yukawa **between two different `Fermion`\n",
    "objects** get checked for gauge invariance, hermiticity (via\n",
    "`Bilinear._eval_conjugate` and the `bar_partner` registry — which has to\n",
    "correctly pair `Ll`/`eRf` legs with `PsiL`/`PsiR` legs), and mass\n",
    "dimension all at once.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "InvarianceReport(15 checks, OK)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "HdH = (dag(H) * H.mat)[0]\n",
    "V = -mu2.s * HdH + lam.s * HdH**2\n",
    "DH = Dmu(H)\n",
    "\n",
    "L = Lagrangian()\n",
    "L.add((dag(DH) * DH)[0], sector=\"kinetic\")\n",
    "L.add(-V, sector=\"potential\")\n",
    "L.add(LYuk, sector=\"yukawa\")\n",
    "L.add(current, sector=\"yukawa\")\n",
    "\n",
    "model = Model(\"SM-VLL\", gauge_groups=[SU2L, U1Y],\n",
    "              fields=[H, Ll, eRf, PsiL, PsiR, W, B],\n",
    "              parameters=[gw, g1, v, lam, mu2, ye, lamE, MPsi],\n",
    "              lagrangian=L)\n",
    "\n",
    "report = model.check_invariance()\n",
    "report.raise_on_failure()\n",
    "report\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {},
   "source": [
    "## 5. EWSB and gauge rotations\n",
    "\n",
    "Identical to the SM tutorial: solve the tadpole for `mu2`, then rotate\n",
    "`(W3, B) -> (Z, A)` (Weinberg angle) and `(W1, W2) -> (W+, W-)`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{mu2: lam*v**2}"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.solve_tadpoles([mu2])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "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, wp=\"W^+\", wm=\"W^-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "## 6. The charged-lepton mass matrix\n",
    "\n",
    "`fermion_mass_matrix(L, bar_base, field_base, ...)` extracts the\n",
    "coefficient of one `(bar, field)` `Bilinear` pair from the vacuum-shifted\n",
    "Lagrangian — it makes no assumption that `bar_base` and `field_base`\n",
    "belong to the same `Fermion` object, so calling it 4× over\n",
    "`(ē_L, Ē_L) × (e_R, E_R)` and assembling a 2×2 by hand gives the full\n",
    "mixing mass matrix:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\sqrt{2} v ye}{2} & 0\\\\\\frac{\\sqrt{2} lamE v}{2} & MPsi\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[  sqrt(2)*v*ye/2,    0],\n",
       "[sqrt(2)*lamE*v/2, MPsi]])"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j = sp.Symbol(\"fl_j\", integer=True)\n",
    "bars, rights = (eLbar, ELbar), (eR, ER)\n",
    "M2 = sp.Matrix(2, 2, lambda a, b: fermion_mass_matrix(\n",
    "    LYuk, bars[a], rights[b], model.vacuum, 1, (i, j), gamma=diracPR)[0, 0])\n",
    "M2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "The lower-triangular shape (zero in the top-right) says something\n",
    "physical: there is no bare `L̄_L Ψ_R` mass in our Lagrangian (it's\n",
    "removable by a field redefinition, so we never wrote one) — the only way\n",
    "`e_R` talks to the doublet sector is through `λ_E`.\n",
    "\n",
    "The neutral sector has no such mixing at all: `N` is a pure Dirac state\n",
    "of mass exactly `M_Ψ`, and the SM neutrino stays exactly massless (no\n",
    "`Bilinear` anywhere pairs `nuL` with a right-handed state).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(MPsi, 0)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m_N = fermion_mass_matrix(LYuk, NLbar, NR, model.vacuum, 1, (i, j),\n",
    "                          gamma=diracPR)[0, 0]\n",
    "m_nu = fermion_mass_matrix(LYuk, nuLbar, NR, model.vacuum, 1, (i, j),\n",
    "                           gamma=diracPR)[0, 0]\n",
    "m_N, m_nu\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "## 7. Biunitary diagonalization\n",
    "\n",
    "A non-symmetric mass matrix needs *two* rotations — one for the\n",
    "bar/left-handed fields, one for the field/right-handed ones — the\n",
    "biunitary (SVD) construction: `θ_L` diagonalizes `M·Mᵀ`, `θ_R`\n",
    "diagonalizes `Mᵀ·M`. `diagonalize_svd`'s generic route\n",
    "(`Matrix.diagonalize(normalize=True)`) returns unusable nested-radical\n",
    "expressions for symbolic entries, so `diagonalize_svd_2x2` instead calls\n",
    "the same analytic `tan 2θ` machinery used for scalar mixing\n",
    "(`solve_mixing_angle_2x2`) on `M·Mᵀ` and `Mᵀ·M` separately.\n",
    "\n",
    "Given the matrix's shape above (`M[0,1]=0`), expect `θ_L` to be small —\n",
    "it can only arise from the *product* `y_e λ_E`, both small numbers —\n",
    "while `θ_R` picks up the direct `λ_E v / M_Ψ` mixing and is the\n",
    "physically large angle. That's the defining signature of a vector-like\n",
    "**doublet** (for a vector-like **singlet** it's the other way around).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Eq(tan(2*thL), -2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2)),\n",
       " Eq(2*tan(thR)/(1 - tan(thR)**2), (2*sqrt(2)*MPsi*lamE*v*sin(atan(2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2))/2)**2 - 2*sqrt(2)*MPsi*v*ye*sin(atan(2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2))/2)*cos(atan(2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2))/2))/(-2*MPsi**2*sin(atan(2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2))/2)**2 + lamE**2*v**2*sin(atan(2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2))/2)**2 - 2*lamE*v**2*ye*sin(atan(2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2))/2)*cos(atan(2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2))/2) + v**2*ye**2*cos(atan(2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2))/2)**2)))"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thL, thR = sp.symbols(\"thL thR\", real=True)\n",
    "e1L, e2L = sp.IndexedBase(\"e1L\"), sp.IndexedBase(\"e2L\")\n",
    "e1R, e2R = sp.IndexedBase(\"e1R\"), sp.IndexedBase(\"e2R\")\n",
    "e1Lbar, e2Lbar = sp.IndexedBase(\"e1Lbar\"), sp.IndexedBase(\"e2Lbar\")\n",
    "e1Rbar, e2Rbar = sp.IndexedBase(\"e1Rbar\"), sp.IndexedBase(\"e2Rbar\")\n",
    "\n",
    "rotL, rotR = diagonalize_svd_2x2(M2, [eL[i], EL[i]], [eR[i], ER[i]],\n",
    "                                 [e1L[i], e2L[i]], [e1R[i], e2R[i]],\n",
    "                                 angle_left=thL, angle_right=thR)\n",
    "rotL.angle_relation, rotR.angle_relation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28",
   "metadata": {},
   "source": [
    "Dual verification (symbolic derivation + numeric check, per\n",
    "`CONVENTIONS.md`): plug in a benchmark point and confirm the rotation\n",
    "really diagonalizes `M`, with the light state landing in slot 1.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.00012044923364816906,\n",
       " -0.06946755651195814,\n",
       " 1.735287235094696,\n",
       " 1002.417724593002)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bench = {ye.s: 0.01, lamE.s: 0.4, v.s: 246.0, MPsi.s: 1000.0}\n",
    "thL_num = float(rotL.angle_solution.subs(bench))\n",
    "thR_num = float(rotR.angle_solution.subs(bench))\n",
    "D = rotation_2x2(thL_num) * M2.subs(bench) * rotation_2x2(thR_num).T\n",
    "\n",
    "assert abs(D[0, 1]) < 1e-9 and abs(D[1, 0]) < 1e-9\n",
    "assert abs(D[0, 0]) < abs(D[1, 1])   # field 1 = light state\n",
    "thL_num, thR_num, float(D[0, 0]), float(D[1, 1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30",
   "metadata": {},
   "source": [
    "## 8. Registering the fermion rotations\n",
    "\n",
    "Rewriting the Lagrangian in the mass basis uses the *same* `Rotation`\n",
    "machinery as the boson case above — `Rotation([eL[i], EL[i]], [e1L[i],\n",
    "e2L[i]], rotation_2x2(θ))` — but with `Indexed` old-fields instead of\n",
    "plain `Symbol`s. Two things are specific to fermions:\n",
    "\n",
    "1. **Register two rotations per chirality**, not one: the field-side\n",
    "   rotation (`eL, EL -> e1L, e2L`) and an *independent* bar-side rotation\n",
    "   (`eLbar, ELbar -> e1Lbar, e2Lbar`) with the same matrix — a `Bilinear`\n",
    "   has a bar leg and a field leg, and both need rewriting.\n",
    "2. **Everything needs a shared flavor-index symbol.** Every term above\n",
    "   was written with the same `i`, so a single `xreplace` reaches every\n",
    "   occurrence at once.\n",
    "\n",
    "Here's the subtlety this newly exercises. After `xreplace`, a `Bilinear`\n",
    "slot holds `cosθ·e1L[i] + sinθ·e2L[i]` — a SymPy `Add` — trapped inside an\n",
    "opaque, custom `Function` that has no linearity rule built in. Watch what\n",
    "a raw `Bilinear` looks like right after the substitution, before anything\n",
    "distributes it:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Bilinear(demo_bar[fl_i], PR, sin(th_{demo})*demo_n2[fl_i] + cos(th_{demo})*demo_n1[fl_i]),\n",
       " sin(th_{demo})*demo_n2[fl_i] + cos(th_{demo})*demo_n1[fl_i])"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demo_bar, demo_field = sp.IndexedBase(\"demo_bar\"), sp.IndexedBase(\"demo_field\")\n",
    "demo_new1, demo_new2 = sp.IndexedBase(\"demo_n1\"), sp.IndexedBase(\"demo_n2\")\n",
    "th_demo = sp.Symbol(\"th_{demo}\", real=True)\n",
    "\n",
    "demo_term = Bilinear(demo_bar[i], diracPR, demo_field[i])\n",
    "# what Rotation.substitution() would build for field[i] -> R^-1 * new[i];\n",
    "# written out directly here since a real Rotation needs 2 old fields:\n",
    "demo_sub = {demo_field[i]: sp.cos(th_demo) * demo_new1[i]\n",
    "                          + sp.sin(th_demo) * demo_new2[i]}\n",
    "demo_rotated = demo_term.xreplace(demo_sub)\n",
    "demo_rotated, demo_rotated.field\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32",
   "metadata": {},
   "source": [
    "The `field` slot is now an unresolved `Add` — `extract_fermion_vertices`\n",
    "would group by this *entire* composite key rather than splitting it into\n",
    "two clean `(coefficient, Indexed)` contributions. `expand_bilinear`\n",
    "(added to `vertices/bilinear.py` alongside this example) distributes\n",
    "`Bilinear` over `Add`-valued legs, exploiting that it's linear in each of\n",
    "its bar/field slots:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sin{\\left(th_{demo} \\right)} \\bar{{demo_{bar}}_{fl_{i}}}\\,P_R\\,{demo_{n2}}_{fl_{i}} + \\cos{\\left(th_{demo} \\right)} \\bar{{demo_{bar}}_{fl_{i}}}\\,P_R\\,{demo_{n1}}_{fl_{i}}$"
      ],
      "text/plain": [
       "sin(th_{demo})*Bilinear(demo_bar[fl_i], PR, demo_n2[fl_i]) + cos(th_{demo})*Bilinear(demo_bar[fl_i], PR, demo_n1[fl_i])"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expand_bilinear(demo_rotated)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34",
   "metadata": {},
   "source": [
    "Two clean terms, each with a plain `Indexed` leg and the right\n",
    "`cosθ`/`sinθ` coefficient. `extract_fermion_vertices` and\n",
    "`fermion_mass_matrix` both call `expand_bilinear` automatically now, so\n",
    "the real model below doesn't need to do this by hand — but this is\n",
    "exactly what would silently go wrong without it.\n",
    "\n",
    "Now register the four real rotations:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}{\\bar{e1R}}_{fl_{i}}\\\\{\\bar{e2R}}_{fl_{i}}\\end{matrix}\\right] = \\left[\\begin{matrix}\\sin{\\left(thR \\right)} {\\bar{ER}}_{fl_{i}} + \\cos{\\left(thR \\right)} {\\bar{eR}}_{fl_{i}}\\\\- \\sin{\\left(thR \\right)} {\\bar{eR}}_{fl_{i}} + \\cos{\\left(thR \\right)} {\\bar{ER}}_{fl_{i}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Rotation(eRbar[fl_i]→e1Rbar[fl_i], ERbar[fl_i]→e2Rbar[fl_i], kind='orthogonal')"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.rotate(rotL)\n",
    "model.rotate(Rotation([eLbar[i], ELbar[i]], [e1Lbar[i], e2Lbar[i]],\n",
    "                      rotation_2x2(thL)))\n",
    "model.rotate(rotR)\n",
    "model.rotate(Rotation([eRbar[i], ERbar[i]], [e1Rbar[i], e2Rbar[i]],\n",
    "                      rotation_2x2(thR)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36",
   "metadata": {},
   "source": [
    "## 9. Extracting the physical couplings\n",
    "\n",
    "`physical_lagrangian` applies the vacuum shift, tadpole solution, and\n",
    "every registered rotation (bosonic *and* fermionic) in one pass;\n",
    "`extract_fermion_vertices` then groups by `(bar, gamma, field)` and peels\n",
    "off the boson legs — exactly as in the SM tutorial.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37",
   "metadata": {},
   "outputs": [],
   "source": [
    "h = sp.Symbol(\"H0_r\", real=True)\n",
    "LYuk_phys = model.physical_lagrangian(sector=\"yukawa\")\n",
    "table = extract_fermion_vertices(LYuk_phys, [Z, A, Wp, Wm, h])\n",
    "\n",
    "mu = sp.Symbol(\"mu\", integer=True)\n",
    "gammaL = DiracGamma(mu) * diracPL\n",
    "gammaR = DiracGamma(mu) * diracPR\n",
    "\n",
    "\n",
    "def coeff(bar, gamma, field, boson):\n",
    "    entry = table.get((bar, gamma, field))\n",
    "    if entry is None:\n",
    "        return sp.S.Zero\n",
    "    return entry.get(1, {}).get((boson,), sp.S.Zero)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38",
   "metadata": {},
   "source": [
    "### Z couplings — the headline result\n",
    "\n",
    "`e_L` and `E_L` are **both** `T³ = -1/2` doublet members, so rotating\n",
    "between them by `θ_L` doesn't change the diagonal Z coupling and produces\n",
    "**exactly zero** FCNC — the left-handed current is protected by the\n",
    "doublet structure, independent of any mixing angle.\n",
    "\n",
    "`e_R` (`T³ = 0`, an SM singlet) and `E_R` (`T³ = -1/2`, from the doublet)\n",
    "are *different* under `SU(2)_L`, so rotating them by `θ_R` both shifts\n",
    "the light state's diagonal coupling **and** generates a genuine FCNC.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{g_{1}^{2} - gw^{2}}{2 \\sqrt{g_{1}^{2} + gw^{2}}}$"
      ],
      "text/plain": [
       "(g1**2 - gw**2)/(2*sqrt(g1**2 + gw**2))"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(e1Lbar[i], gammaL, e1L[i], Z))   # LH diagonal: unchanged\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(e1Lbar[i], gammaL, e2L[i], Z))   # LH FCNC: exactly zero\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{- g_{1}^{2} \\sin^{2}{\\left(thR \\right)} + 2 g_{1}^{2} - gw^{2} \\sin^{2}{\\left(thR \\right)}}{2 \\sqrt{g_{1}^{2} + gw^{2}}}$"
      ],
      "text/plain": [
       "(-g1**2*sin(thR)**2 + 2*g1**2 - gw**2*sin(thR)**2)/(2*sqrt(g1**2 + gw**2))"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(e1Rbar[i], gammaR, e1R[i], Z))   # RH diagonal: shifted\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{\\sqrt{g_{1}^{2} + gw^{2}} \\sin{\\left(2 thR \\right)}}{4}$"
      ],
      "text/plain": [
       "-sqrt(g1**2 + gw**2)*sin(2*thR)/4"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(e1Rbar[i], gammaR, e2R[i], Z))   # RH FCNC: -(gZ/2) sinθR cosθR\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43",
   "metadata": {},
   "source": [
    "### Photon couplings — a sanity check\n",
    "\n",
    "`U(1)_EM` is unbroken, so the rotation must not touch the electric\n",
    "charge: diagonal coupling `-e` in both mass eigenstates, and **zero**\n",
    "photon FCNC (a mismatch here would mean a bug in the rotation, not new\n",
    "physics).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{g_{1} gw}{\\sqrt{g_{1}^{2} + gw^{2}}}$"
      ],
      "text/plain": [
       "-g1*gw/sqrt(g1**2 + gw**2)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(e1Lbar[i], gammaL, e1L[i], A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(e1Lbar[i], gammaL, e2L[i], A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46",
   "metadata": {},
   "source": [
    "### W couplings — a genuinely new right-handed current\n",
    "\n",
    "The left-handed current splits between the light and heavy states by the\n",
    "usual `cosθ_L`/`sinθ_L` pattern. But because `Ψ_R` is *also* an `SU(2)_L`\n",
    "doublet (unlike the SM's `e_R` singlet), there's a **right-handed** charged\n",
    "current `W⁺ N̄_R e_R ∝ sinθ_R` that has no SM analogue at all.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\sqrt{2} gw \\cos{\\left(thL \\right)}}{2}$"
      ],
      "text/plain": [
       "sqrt(2)*gw*cos(thL)/2"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(nuLbar[i], gammaL, e1L[i], Wp))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{\\sqrt{2} gw \\sin{\\left(thL \\right)}}{2}$"
      ],
      "text/plain": [
       "-sqrt(2)*gw*sin(thL)/2"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(nuLbar[i], gammaL, e2L[i], Wp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\sqrt{2} gw \\sin{\\left(thR \\right)}}{2}$"
      ],
      "text/plain": [
       "sqrt(2)*gw*sin(thR)/2"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(coeff(NRbar[i], gammaR, e1R[i], Wp))   # the new RH current\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50",
   "metadata": {},
   "source": [
    "### Higgs couplings — a sum rule\n",
    "\n",
    "The light lepton's Higgs coupling deviates from the SM value `-m_e/v`\n",
    "by a clean, closed-form amount: `sinθ_L sinθ_R M_Ψ/v`. Rather than just\n",
    "printing the coupling, build the SM mass formula for this mixed system\n",
    "and verify the deviation symbolically — the residue below should\n",
    "`simplify` to exactly zero.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cL, sL = sp.cos(thL), sp.sin(thL)\n",
    "cR, sR = sp.cos(thR), sp.sin(thR)\n",
    "m_e_expr = (cL * ye.s + sL * lamE.s) * v.s / sp.sqrt(2) * cR + sL * sR * MPsi.s\n",
    "\n",
    "h_e1e1 = coeff(e1Lbar[i], diracPR, e1R[i], h)\n",
    "residue = sp.simplify(h_e1e1 + (m_e_expr - sL * sR * MPsi.s) / v.s)\n",
    "residue\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52",
   "metadata": {},
   "source": [
    "## 10. Recap\n",
    "\n",
    "Starting from one bare mass term and one mixing Yukawa, `feynlag`\n",
    "mechanically derived: a non-trivial 2×2 mixing matrix, two independent\n",
    "diagonalization angles with the doublet's characteristic\n",
    "θ_L-small/θ_R-large hierarchy, an exactly-zero left-handed Z-FCNC\n",
    "alongside a nonzero right-handed one, a new right-handed W current, and a\n",
    "closed-form Higgs-coupling sum rule — all *derived*, not assumed, and all\n",
    "cross-checked symbolically above.\n",
    "\n",
    "- Exact pinned formulas for every coupling here (including the ones we\n",
    "  only spot-checked): [`tests/test_vll.py`](../tests/test_vll.py).\n",
    "- The full flat script this notebook walks through:\n",
    "  [`examples/sm_vll.py`](sm_vll.py).\n",
    "- For the plain-SM baseline this model extends:\n",
    "  [`SM_Feynman_Rules_Tutorial.ipynb`](SM_Feynman_Rules_Tutorial.ipynb) /\n",
    "  [`examples/sm_scalar_gauge.py`](sm_scalar_gauge.py).\n"
   ]
  }
 ],
 "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
}