Vector-like fermions in feynlag: SM + a vector-like lepton doublet

A vector-like fermion has both chiralities in the same gauge representation, so — unlike a chiral SM fermion — it can carry a bare, gauge-invariant Dirac mass all by itself. Mixed with an ordinary SM fermion through the Higgs, it produces a textbook flavor-mixing problem: a 2×2 mass matrix diagonalized by two independent angles (θ_L, θ_R, one per chirality), and — because the heavy and light states are no longer pure gauge eigenstates — flavor-changing neutral currents (FCNCs) that don’t exist in the SM.

This notebook walks through examples/sm_vll.py stage by stage, reusing the same declare → write → check → break-symmetry → diagonalize → extract pipeline as SM_Feynman_Rules_Tutorial.ipynb (read that one first if ExternalParameter, Dmu, Scalar, or Model are unfamiliar — this notebook won’t re-explain them). What’s new here:

  • a vector-like doublet Ψ = (N, E), built as two WeylFermions with identical reps (not a DiracFermion — see §2),

  • a bare mass term -M Ψ̄_L Ψ_R that needs no scalar VEV to exist,

  • biunitary (2×2 SVD) diagonalization with two physically distinct angles, via the analytic diagonalize_svd_2x2 helper,

  • rotating fermion fields into their mass basis and re-extracting Feynman rules — which exercises expand_bilinear, a piece of machinery this model is the first to need.

The headline physics we’ll derive and verify symbolically:

  1. the Z-FCNC lives only in the right-handed current (the left-handed Z coupling is exactly protected),

  2. a clean Higgs-coupling sum rule for the deviation from -m_e/v,

  3. a new right-handed W current with no SM analogue.

import sympy as sp
from feynlag import (
    Bilinear, DiracGamma, Dmu, ExternalParameter, InternalParameter,
    Lagrangian, Model, Rotation, SU2, Scalar, U1, WeylFermion, dag,
    diagonalize_svd_2x2, diracPL, diracPR, expand_bilinear,
    extract_fermion_vertices, fermion_gauge_current, fermion_mass_matrix,
    rotation_2x2, weinberg_rotation, charged_current_rotation,
)

1. Symmetries and parameters

Same SU2L/U1Y electroweak setup as the plain-SM tutorial, plus three new couplings for the VLL sector:

  • ye — the ordinary (tiny) SM electron Yukawa,

  • lamE — the mixing Yukawa between the heavy doublet and e_R,

  • MPsi — the bare vector-like mass. It needs unit_dim=1 (mass dimension 1) so Model.check_invariance()’s power-counting check can see that MPsi · (dimension-3 Bilinear) is a dimension-4 operator — exactly like v needs unit_dim=1 for the Yukawa terms below.

gw = ExternalParameter("gw", 0.6535, positive=True)
g1 = ExternalParameter("g1", 0.3580, positive=True)
SU2L, U1Y = SU2("SU2L", coupling=gw), U1("U1Y", coupling=g1)

v = ExternalParameter("v", 246.0, positive=True, unit_dim=1)
lam = ExternalParameter("lam", 0.129, tex="\\lambda")
mu2 = InternalParameter("mu2", unit_dim=2, tex="\\mu^2")

ye = ExternalParameter("ye", 0.01, positive=True, tex="y_e")
lamE = ExternalParameter("lamE", 0.4, positive=True, tex="\\lambda_E")
MPsi = ExternalParameter("MPsi", 1000.0, positive=True, unit_dim=1, tex="M_\\Psi")
MPsi.unit_dim

1

2. Fields

The Higgs and the SM lepton generation are declared exactly as in the SM tutorial (one flavor: nflavors=1, so the mixing physics isn’t obscured by 3-generation flavor indices).

H = Scalar("H", reps={SU2L: 2, U1Y: sp.Rational(1, 2)},
           component_names=["Gp", "H0"])
H.expand_vev({H.components[1]: v})

Ll = WeylFermion("Ll", reps={SU2L: 2, U1Y: -sp.Rational(1, 2)},
                 chirality="L", nflavors=1, component_names=["nuL", "eL"])
eRf = WeylFermion("eRf", reps={U1Y: -1}, chirality="R", nflavors=1,
                  component_names=["eR"])
Ll.components, eRf.components

([nuL, eL], [eR])

Now the vector-like doublet. PsiL and PsiR are given the same reps (SU2L: 2, U1Y: -1/2) as Ll — both chiralities in one representation is precisely what makes a bare mass gauge invariant (a chiral SM fermion can’t have one: L_L and e_R transform differently, so nothing gauge-invariant connects them without a scalar).

Why two WeylFermions and not a single DiracFermion? DiracFermion’s chirality=None default routes through fermion_gauge_current with the identity projector, producing a bare DiracGamma(mu)*diracI structure that dirac_conjugate doesn’t know how to Dirac-conjugate — it raises NotImplementedError, silently breaking the hermiticity check. That path is untested (see fields.py’s DiracFermion docstring); modeling the VLL as two WeylFermions with identical reps sidesteps it entirely and is a well-exercised pattern.

PsiL = WeylFermion("PsiL", reps={SU2L: 2, U1Y: -sp.Rational(1, 2)},
                   chirality="L", nflavors=1, component_names=["NL", "EL"])
PsiR = WeylFermion("PsiR", reps={SU2L: 2, U1Y: -sp.Rational(1, 2)},
                   chirality="R", nflavors=1, component_names=["NR", "ER"])
W, B = SU2L.bosons("W"), U1Y.bosons("B")
PsiL.components, PsiR.components

([NL, EL], [NR, ER])

Gp, H0 = H.components
nuL, eL = Ll.components
nuLbar, eLbar = Ll.bar_components
eR, eRbar = eRf.components[0], eRf.bar_components[0]
NL, EL = PsiL.components
NLbar, ELbar = PsiL.bar_components
NR, ER = PsiR.components
NRbar, ERbar = PsiR.bar_components

i = sp.Symbol("fl_i", integer=True)

3. The Lagrangian

Three Yukawa-sector pieces, built up one at a time so the + h.c. pattern is visible at each step (Bilinear(bar, gamma, field) is the opaque fermion-sandwich atom; diracPR/diracPL are the chiral projectors — see the SM tutorial for the full explanation of this convention).

SM Yukawa: -y_e L̄_L H e_R + h.c.

Identical in form to the plain-SM lepton Yukawa.

LYuk = -(ye.s * Gp * Bilinear(nuLbar[i], diracPR, eR[i])
         + ye.s * H0 * Bilinear(eLbar[i], diracPR, eR[i]))
LYuk += -(ye.s * sp.conjugate(Gp) * Bilinear(eRbar[i], diracPL, nuL[i])
          + ye.s * sp.conjugate(H0) * Bilinear(eRbar[i], diracPL, eL[i]))
LYuk
\[\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)\]

Mixing Yukawa: -λ_E Ψ̄_L H e_R + h.c.

Ψ_L carries exactly L_L’s reps, so the component pattern is identical to the SM Yukawa above — just with NL/EL in place of nuL/eL. This single term is what ties the heavy doublet to the light SM lepton.

LYuk += -(lamE.s * Gp * Bilinear(NLbar[i], diracPR, eR[i])
          + lamE.s * H0 * Bilinear(ELbar[i], diracPR, eR[i]))
LYuk += -(lamE.s * sp.conjugate(Gp) * Bilinear(eRbar[i], diracPL, NL[i])
          + lamE.s * sp.conjugate(H0) * Bilinear(eRbar[i], diracPL, EL[i]))
LYuk
\[\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)\]

Bare vector-like mass: -M_Ψ Ψ̄_L Ψ_R + h.c.

No Higgs field appears — this term needs no VEV to exist, which is the whole point of “vector-like.” Both N and E get this same mass M_Ψ.

LYuk += -MPsi.s * (Bilinear(NLbar[i], diracPR, NR[i])
                   + Bilinear(ELbar[i], diracPR, ER[i]))
LYuk += -MPsi.s * (Bilinear(NRbar[i], diracPL, NL[i])
                   + Bilinear(ERbar[i], diracPL, EL[i]))
LYuk
\[\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)\]

Finally the gauge currents (fermion_gauge_current is fully generic — one call per fermion, regardless of how many gauge groups it’s charged under):

current = (fermion_gauge_current(Ll, i) + fermion_gauge_current(eRf, i)
           + fermion_gauge_current(PsiL, i)
           + fermion_gauge_current(PsiR, i))
current
\[\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}\]

4. Assemble the Model and check invariance

The kinetic/potential Higgs sector is the standard SM piece. Assembling everything into a Model and calling check_invariance() is the correctness gate: this is the first model in feynlag where a bare fermion mass and a mixing Yukawa between two different Fermion objects get checked for gauge invariance, hermiticity (via Bilinear._eval_conjugate and the bar_partner registry — which has to correctly pair Ll/eRf legs with PsiL/PsiR legs), and mass dimension all at once.

HdH = (dag(H) * H.mat)[0]
V = -mu2.s * HdH + lam.s * HdH**2
DH = Dmu(H)

L = Lagrangian()
L.add((dag(DH) * DH)[0], sector="kinetic")
L.add(-V, sector="potential")
L.add(LYuk, sector="yukawa")
L.add(current, sector="yukawa")

model = Model("SM-VLL", gauge_groups=[SU2L, U1Y],
              fields=[H, Ll, eRf, PsiL, PsiR, W, B],
              parameters=[gw, g1, v, lam, mu2, ye, lamE, MPsi],
              lagrangian=L)

report = model.check_invariance()
report.raise_on_failure()
report

InvarianceReport(15 checks, OK)

5. EWSB and gauge rotations

Identical to the SM tutorial: solve the tadpole for mu2, then rotate (W3, B) -> (Z, A) (Weinberg angle) and (W1, W2) -> (W+, W-).

model.solve_tadpoles([mu2])

{mu2: lam*v**2}

# physical basis: the standard Weinberg + W± rotations, from feynlag.models
Z, A = weinberg_rotation(model, SU2L, U1Y)
Wp, Wm = charged_current_rotation(model, SU2L, wp="W^+", wm="W^-")

6. The charged-lepton mass matrix

fermion_mass_matrix(L, bar_base, field_base, ...) extracts the coefficient of one (bar, field) Bilinear pair from the vacuum-shifted Lagrangian — it makes no assumption that bar_base and field_base belong to the same Fermion object, so calling it 4× over (ē_L, Ē_L) × (e_R, E_R) and assembling a 2×2 by hand gives the full mixing mass matrix:

j = sp.Symbol("fl_j", integer=True)
bars, rights = (eLbar, ELbar), (eR, ER)
M2 = sp.Matrix(2, 2, lambda a, b: fermion_mass_matrix(
    LYuk, bars[a], rights[b], model.vacuum, 1, (i, j), gamma=diracPR)[0, 0])
M2
\[\begin{split}\displaystyle \left[\begin{matrix}\frac{\sqrt{2} v ye}{2} & 0\\\frac{\sqrt{2} lamE v}{2} & MPsi\end{matrix}\right]\end{split}\]

The lower-triangular shape (zero in the top-right) says something physical: there is no bare L̄_L Ψ_R mass in our Lagrangian (it’s removable by a field redefinition, so we never wrote one) — the only way e_R talks to the doublet sector is through λ_E.

The neutral sector has no such mixing at all: N is a pure Dirac state of mass exactly M_Ψ, and the SM neutrino stays exactly massless (no Bilinear anywhere pairs nuL with a right-handed state).

m_N = fermion_mass_matrix(LYuk, NLbar, NR, model.vacuum, 1, (i, j),
                          gamma=diracPR)[0, 0]
m_nu = fermion_mass_matrix(LYuk, nuLbar, NR, model.vacuum, 1, (i, j),
                           gamma=diracPR)[0, 0]
m_N, m_nu

(MPsi, 0)

7. Biunitary diagonalization

A non-symmetric mass matrix needs two rotations — one for the bar/left-handed fields, one for the field/right-handed ones — the biunitary (SVD) construction: θ_L diagonalizes M·Mᵀ, θ_R diagonalizes Mᵀ·M. diagonalize_svd’s generic route (Matrix.diagonalize(normalize=True)) returns unusable nested-radical expressions for symbolic entries, so diagonalize_svd_2x2 instead calls the same analytic tan 2θ machinery used for scalar mixing (solve_mixing_angle_2x2) on M·Mᵀ and Mᵀ·M separately.

Given the matrix’s shape above (M[0,1]=0), expect θ_L to be small — it can only arise from the product y_e λ_E, both small numbers — while θ_R picks up the direct λ_E v / M_Ψ mixing and is the physically large angle. That’s the defining signature of a vector-like doublet (for a vector-like singlet it’s the other way around).

thL, thR = sp.symbols("thL thR", real=True)
e1L, e2L = sp.IndexedBase("e1L"), sp.IndexedBase("e2L")
e1R, e2R = sp.IndexedBase("e1R"), sp.IndexedBase("e2R")
e1Lbar, e2Lbar = sp.IndexedBase("e1Lbar"), sp.IndexedBase("e2Lbar")
e1Rbar, e2Rbar = sp.IndexedBase("e1Rbar"), sp.IndexedBase("e2Rbar")

rotL, rotR = diagonalize_svd_2x2(M2, [eL[i], EL[i]], [eR[i], ER[i]],
                                 [e1L[i], e2L[i]], [e1R[i], e2R[i]],
                                 angle_left=thL, angle_right=thR)
rotL.angle_relation, rotR.angle_relation

(Eq(tan(2*thL), -2*lamE*v**2*ye/(2*MPsi**2 + lamE**2*v**2 - v**2*ye**2)),
 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)))

Dual verification (symbolic derivation + numeric check, per CONVENTIONS.md): plug in a benchmark point and confirm the rotation really diagonalizes M, with the light state landing in slot 1.

bench = {ye.s: 0.01, lamE.s: 0.4, v.s: 246.0, MPsi.s: 1000.0}
thL_num = float(rotL.angle_solution.subs(bench))
thR_num = float(rotR.angle_solution.subs(bench))
D = rotation_2x2(thL_num) * M2.subs(bench) * rotation_2x2(thR_num).T

assert abs(D[0, 1]) < 1e-9 and abs(D[1, 0]) < 1e-9
assert abs(D[0, 0]) < abs(D[1, 1])   # field 1 = light state
thL_num, thR_num, float(D[0, 0]), float(D[1, 1])

(-0.00012044923364816906,
 -0.06946755651195814,
 1.735287235094696,
 1002.417724593002)

8. Registering the fermion rotations

Rewriting the Lagrangian in the mass basis uses the same Rotation machinery as the boson case above — Rotation([eL[i], EL[i]], [e1L[i], e2L[i]], rotation_2x2(θ)) — but with Indexed old-fields instead of plain Symbols. Two things are specific to fermions:

  1. Register two rotations per chirality, not one: the field-side rotation (eL, EL -> e1L, e2L) and an independent bar-side rotation (eLbar, ELbar -> e1Lbar, e2Lbar) with the same matrix — a Bilinear has a bar leg and a field leg, and both need rewriting.

  2. Everything needs a shared flavor-index symbol. Every term above was written with the same i, so a single xreplace reaches every occurrence at once.

Here’s the subtlety this newly exercises. After xreplace, a Bilinear slot holds cosθ·e1L[i] + sinθ·e2L[i] — a SymPy Add — trapped inside an opaque, custom Function that has no linearity rule built in. Watch what a raw Bilinear looks like right after the substitution, before anything distributes it:

demo_bar, demo_field = sp.IndexedBase("demo_bar"), sp.IndexedBase("demo_field")
demo_new1, demo_new2 = sp.IndexedBase("demo_n1"), sp.IndexedBase("demo_n2")
th_demo = sp.Symbol("th_{demo}", real=True)

demo_term = Bilinear(demo_bar[i], diracPR, demo_field[i])
# what Rotation.substitution() would build for field[i] -> R^-1 * new[i];
# written out directly here since a real Rotation needs 2 old fields:
demo_sub = {demo_field[i]: sp.cos(th_demo) * demo_new1[i]
                          + sp.sin(th_demo) * demo_new2[i]}
demo_rotated = demo_term.xreplace(demo_sub)
demo_rotated, demo_rotated.field

(Bilinear(demo_bar[fl_i], PR, sin(th_{demo})*demo_n2[fl_i] + cos(th_{demo})*demo_n1[fl_i]),
 sin(th_{demo})*demo_n2[fl_i] + cos(th_{demo})*demo_n1[fl_i])

The field slot is now an unresolved Add — extract_fermion_vertices would group by this entire composite key rather than splitting it into two clean (coefficient, Indexed) contributions. expand_bilinear (added to vertices/bilinear.py alongside this example) distributes Bilinear over Add-valued legs, exploiting that it’s linear in each of its bar/field slots:

expand_bilinear(demo_rotated)
\[\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}}\]

Two clean terms, each with a plain Indexed leg and the right cosθ/sinθ coefficient. extract_fermion_vertices and fermion_mass_matrix both call expand_bilinear automatically now, so the real model below doesn’t need to do this by hand — but this is exactly what would silently go wrong without it.

Now register the four real rotations:

model.rotate(rotL)
model.rotate(Rotation([eLbar[i], ELbar[i]], [e1Lbar[i], e2Lbar[i]],
                      rotation_2x2(thL)))
model.rotate(rotR)
model.rotate(Rotation([eRbar[i], ERbar[i]], [e1Rbar[i], e2Rbar[i]],
                      rotation_2x2(thR)))
\[\begin{split}\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]\end{split}\]

9. Extracting the physical couplings

physical_lagrangian applies the vacuum shift, tadpole solution, and every registered rotation (bosonic and fermionic) in one pass; extract_fermion_vertices then groups by (bar, gamma, field) and peels off the boson legs — exactly as in the SM tutorial.

h = sp.Symbol("H0_r", real=True)
LYuk_phys = model.physical_lagrangian(sector="yukawa")
table = extract_fermion_vertices(LYuk_phys, [Z, A, Wp, Wm, h])

mu = sp.Symbol("mu", integer=True)
gammaL = DiracGamma(mu) * diracPL
gammaR = DiracGamma(mu) * diracPR


def coeff(bar, gamma, field, boson):
    entry = table.get((bar, gamma, field))
    if entry is None:
        return sp.S.Zero
    return entry.get(1, {}).get((boson,), sp.S.Zero)

Z couplings — the headline result

e_L and E_L are both T³ = -1/2 doublet members, so rotating between them by θ_L doesn’t change the diagonal Z coupling and produces exactly zero FCNC — the left-handed current is protected by the doublet structure, independent of any mixing angle.

e_R (T³ = 0, an SM singlet) and E_R (T³ = -1/2, from the doublet) are different under SU(2)_L, so rotating them by θ_R both shifts the light state’s diagonal coupling and generates a genuine FCNC.

sp.simplify(coeff(e1Lbar[i], gammaL, e1L[i], Z))   # LH diagonal: unchanged
\[\displaystyle \frac{g_{1}^{2} - gw^{2}}{2 \sqrt{g_{1}^{2} + gw^{2}}}\]
sp.simplify(coeff(e1Lbar[i], gammaL, e2L[i], Z))   # LH FCNC: exactly zero
\[\displaystyle 0\]
sp.simplify(coeff(e1Rbar[i], gammaR, e1R[i], Z))   # RH diagonal: shifted
\[\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}}}\]
sp.simplify(coeff(e1Rbar[i], gammaR, e2R[i], Z))   # RH FCNC: -(gZ/2) sinθR cosθR
\[\displaystyle - \frac{\sqrt{g_{1}^{2} + gw^{2}} \sin{\left(2 thR \right)}}{4}\]

Photon couplings — a sanity check

U(1)_EM is unbroken, so the rotation must not touch the electric charge: diagonal coupling -e in both mass eigenstates, and zero photon FCNC (a mismatch here would mean a bug in the rotation, not new physics).

sp.simplify(coeff(e1Lbar[i], gammaL, e1L[i], A))
\[\displaystyle - \frac{g_{1} gw}{\sqrt{g_{1}^{2} + gw^{2}}}\]
sp.simplify(coeff(e1Lbar[i], gammaL, e2L[i], A))
\[\displaystyle 0\]

W couplings — a genuinely new right-handed current

The left-handed current splits between the light and heavy states by the usual cosθ_L/sinθ_L pattern. But because Ψ_R is also an SU(2)_L doublet (unlike the SM’s e_R singlet), there’s a right-handed charged current W⁺ N̄_R e_R ∝ sinθ_R that has no SM analogue at all.

sp.simplify(coeff(nuLbar[i], gammaL, e1L[i], Wp))
\[\displaystyle \frac{\sqrt{2} gw \cos{\left(thL \right)}}{2}\]
sp.simplify(coeff(nuLbar[i], gammaL, e2L[i], Wp))
\[\displaystyle - \frac{\sqrt{2} gw \sin{\left(thL \right)}}{2}\]
sp.simplify(coeff(NRbar[i], gammaR, e1R[i], Wp))   # the new RH current
\[\displaystyle \frac{\sqrt{2} gw \sin{\left(thR \right)}}{2}\]

Higgs couplings — a sum rule

The light lepton’s Higgs coupling deviates from the SM value -m_e/v by a clean, closed-form amount: sinθ_L sinθ_R M_Ψ/v. Rather than just printing the coupling, build the SM mass formula for this mixed system and verify the deviation symbolically — the residue below should simplify to exactly zero.

cL, sL = sp.cos(thL), sp.sin(thL)
cR, sR = sp.cos(thR), sp.sin(thR)
m_e_expr = (cL * ye.s + sL * lamE.s) * v.s / sp.sqrt(2) * cR + sL * sR * MPsi.s

h_e1e1 = coeff(e1Lbar[i], diracPR, e1R[i], h)
residue = sp.simplify(h_e1e1 + (m_e_expr - sL * sR * MPsi.s) / v.s)
residue
\[\displaystyle 0\]

10. Recap

Starting from one bare mass term and one mixing Yukawa, feynlag mechanically derived: a non-trivial 2×2 mixing matrix, two independent diagonalization angles with the doublet’s characteristic θ_L-small/θ_R-large hierarchy, an exactly-zero left-handed Z-FCNC alongside a nonzero right-handed one, a new right-handed W current, and a closed-form Higgs-coupling sum rule — all derived, not assumed, and all cross-checked symbolically above.