Extending the SM by a gauged U(1)\(_X\): a Z′ tutorial in feynlag

A Z′ boson is what you get when the SM gauge group grows an extra abelian factor U(1)\(_X\). The SM fermion content (plus three right-handed neutrinos) admits exactly two independent family-universal anomaly-free abelian charges — hypercharge \(Y\) and \(B\!-\!L\) — so the most general consistent assignment is a linear combination

\[X = a\,Y + b\,(B\!-\!L)\]

with \(a\), \(b\) free parameters. This notebook keeps them (and the singlet’s charge \(q_S\)) symbolic all the way through the pipeline — declaration, invariance checking, gauge currents, vertex extraction.

Because the Higgs carries hypercharge, it carries \(X\)-charge \(a/2\), and the \(Z\) and \(Z'\) mix at tree level — the physically interesting (and technically instructive) regime. What’s new relative to SM_Feynman_Rules_Tutorial.ipynb (read that first — ExternalParameter, Dmu, Scalar, Model are not re-explained here):

• two independent U(1) factors in one model (groups are distinct objects; a field’s reps dict just gains a second U(1) entry),

• an SM-singlet complex scalar \(S\) whose VEV \(v_X\) breaks U(1)\(_X\), with a Higgs-portal coupling \(\lambda_{HS}|H|^2|S|^2\) mixing \(h\) and \(s\),

• three right-handed neutrinos (required for the \(B\!-\!L\) anomaly to cancel) with a Dirac Yukawa via \(\tilde H\),

• a 3×3 neutral gauge mass matrix diagonalized by two chained 2×2 rotations — Weinberg first, then \(Z\)–\(Z'\).

The headline formulas we’ll derive and verify:

\[\tan 2\theta' = \frac{2\,a\,g_Z g_X v^2}{a^2 g_X^2 v^2 + 4 g_X^2 q_S^2 v_X^2 - g_Z^2 v^2}, \qquad \tan 2\theta_s = \frac{\lambda_{HS}\, v\, v_X}{\lambda v^2 - \lambda_S v_X^2}, \qquad m_Z^2\, m_{Z'}^2 = \tfrac{1}{4} g_Z^2 g_X^2 q_S^2 v^2 v_X^2 .\]

\(a \to 0\) is the pure gauged-\(B\!-\!L\) limit: no mixing, \(m_{Z'} = g_X q_S v_X\).

import sympy as sp

from feynlag import (
    Bilinear, DiracGamma, Dmu, ExternalParameter, InternalParameter,
    Lagrangian, Model, Rotation, SU2, SU3, Scalar, U1, WeylFermion,
    conjugate_pair, dag, diagonalize_orthogonal_2x2, diracPL, diracPR,
    extract_fermion_vertices, fermion_feynman_rule, fermion_gauge_current,
    fermion_mass_matrix, rotation_2x2,
)

1. Symmetries and parameters

Two things distinguish this from the plain-SM setup:

• A second U(1). U1X = U1("U1X", coupling=gX) — groups compare by identity, so U1Y and U1X coexist as separate keys in every field’s reps dict, each carrying its own coupling and its own gauge boson.

• Symbolic charges. aX, bX, qS are ordinary ExternalParameters (dimensionless, so the default unit_dim=0 is right). A U(1) “generator” in feynlag is just the 1×1 matrix [[charge]], and nothing downstream assumes the charge is a number.

The charge map qcharge(Y, B−L) = a·Y + b·(B−L) below is the entire flavor structure of the extension — every field’s \(X\)-charge follows from its known \(Y\) and \(B\!-\!L\).

gw = ExternalParameter("g_w", 0.6535, positive=True)
g1 = ExternalParameter("g_1", 0.3580, positive=True)
gX = ExternalParameter("g_X", 0.5, positive=True)
gs = ExternalParameter("g_s", 1.22, positive=True)
SU2L, U1Y = SU2("SU2L", coupling=gw), U1("U1Y", coupling=g1)
U1X = U1("U1X", coupling=gX)
SU3c = SU3("SU3c", coupling=gs)

aX = ExternalParameter("a_X", 0.4)      # X = aX*Y + bX*(B-L), kept symbolic
bX = ExternalParameter("b_X", -0.7)
qS = ExternalParameter("q_S", 1.0)      # singlet X-charge, also symbolic

v = ExternalParameter("v", 246.0, positive=True, unit_dim=1)
vX = ExternalParameter("v_X", 2000.0, positive=True, unit_dim=1)
lam = ExternalParameter("lambda", 0.129)
lamS = ExternalParameter("lambda_S", 0.2)
lamHS = ExternalParameter("lambda_HS", 0.05)
mu2 = InternalParameter("mu_2", unit_dim=2)
muS2 = InternalParameter("mu_{S2}", unit_dim=2)


def qcharge(Y, BmL):
    """The anomaly-free charge map X = a*Y + b*(B-L)."""
    return aX.s * Y + bX.s * BmL

2. Fields

Two scalars now:

• the SM Higgs doublet H, which picks up \(X\)-charge \(a/2\) (it has \(Y=1/2\), \(B\!-\!L=0\)) — this is what turns on the tree-level \(Z\)–\(Z'\) mixing;

• an SM-singlet complex scalar S with the free charge \(q_S\) and VEV \(v_X\).

S is a one-component complex scalar; expand_vev on a complex component automatically creates the real fluctuation symbols and registers \(S \to (v_X + S0_r + i\,S0_i)/\sqrt{2}\) — same \(1/\sqrt2\) convention as the Higgs. (\(S0_r\) is the CP-even “dark Higgs” \(s\); \(S0_i\) will turn out to be the Goldstone the \(Z'\) eats.)

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

S = Scalar("S", reps={U1X: qS.s}, component_names=["S0"])
S.expand_vev({S.components[0]: vX})

S.shift_map     # the automatic complex-VEV expansion

{S0: sqrt(2)*(I*S0_i + S0_r + v_X)/2}

The fermions are the SM’s three generations plus three right-handed neutrinos (nuR, with \(Y=0\), \(B\!-\!L=-1\), hence \(X\)-charge \(-b\) — they carry only the new gauge charge). Without them the \(B\!-\!L\) part of \(X\) would be anomalous; with them, any \((a, b)\) is anomaly-free. Every SM field’s reps dict simply gains its U1X: qcharge(...) entry — nothing else about the declarations changes.

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

QL = WeylFermion("QL", reps={SU2L: 2, U1Y: sp.Rational(1, 6), SU3c: 3,
                             U1X: qcharge(sp.Rational(1, 6), sp.Rational(1, 3))},
                 chirality="L", nflavors=3,
                 component_names=["uL_1", "uL_2", "uL_3",
                                  "dL_1", "dL_2", "dL_3"])
uR = WeylFermion("uR", reps={U1Y: sp.Rational(2, 3), SU3c: 3,
                             U1X: qcharge(sp.Rational(2, 3), sp.Rational(1, 3))},
                 chirality="R", nflavors=3,
                 component_names=["uR_1", "uR_2", "uR_3"])
dR = WeylFermion("dR", reps={U1Y: -sp.Rational(1, 3), SU3c: 3,
                             U1X: qcharge(-sp.Rational(1, 3), sp.Rational(1, 3))},
                 chirality="R", nflavors=3,
                 component_names=["dR_1", "dR_2", "dR_3"])

W, B = SU2L.bosons("W"), U1Y.bosons("B")
X = U1X.bosons("X")
G = SU3c.bosons("G")

Gp, H0 = H.components
nuL, eL = Ll.components
nuLbar, eLbar = Ll.bar_components
eRc, eRbar = eR.components[0], eR.bar_components[0]
nuRc, nuRbar = nuR.components[0], nuR.bar_components[0]

3. The Lagrangian

Scalar potential with a Higgs portal

\[V = -\mu^2|H|^2 + \lambda|H|^4 - \mu_S^2|S|^2 + \lambda_S|S|^4 + \lambda_{HS}|H|^2|S|^2\]

The portal term is the only renormalizable, gauge-invariant coupling between the two scalar sectors — it will mix \(h\) and \(s\) after both VEVs turn on. The kinetic terms are built with Dmu, which gauges every group in a field’s reps automatically: Dmu(H) includes \(W\), \(B\) and \(X\); Dmu(S) includes only \(X\).

HdH = (dag(H) * H.mat)[0]
SdS = (dag(S) * S.mat)[0]
V = (-mu2.s * HdH + lam.s * HdH**2
     - muS2.s * SdS + lamS.s * SdS**2
     + lamHS.s * HdH * SdS)
DH, DS = Dmu(H), Dmu(S)

Yukawas: charged leptons, and neutrinos via \(\tilde H\)

The charged-lepton Yukawa is the standard \(-Y_e^{ij}\,\bar L_i H\, e_{Rj} + \text{h.c.}\) The neutrino Dirac Yukawa \(-Y_\nu^{ij}\,\bar L_i \tilde H\, \nu_{Rj} + \text{h.c.}\) needs the conjugate doublet \(\tilde H = i\sigma_2 H^* = (H^{0*}, -G^{+*})\), written inline component-by-component — the same pattern the up-quark Yukawa uses in the SM tutorial. Check the \(X\)-charge balance: \(-q_L + q_{\tilde H} + q_{\nu_R} = (a/2 + b) + (-a/2) + (-b) = 0\) — invariant by construction, for any \(a, b\). (check_invariance() will verify this symbolically in §4.)

i, j = sp.symbols("i j", integer=True)
Ye, Yv = sp.IndexedBase("Ye"), sp.IndexedBase("Yv")

LYuk = -(Ye[i, j] * Gp * Bilinear(nuLbar[i], diracPR, eRc[j])
         + Ye[i, j] * H0 * Bilinear(eLbar[i], diracPR, eRc[j]))
LYuk += -(sp.conjugate(Ye[i, j]) * sp.conjugate(Gp)
          * Bilinear(eRbar[j], diracPL, nuL[i])
          + sp.conjugate(Ye[i, j]) * sp.conjugate(H0)
          * Bilinear(eRbar[j], diracPL, eL[i]))

LYukN = -(Yv[i, j] * sp.conjugate(H0)
          * Bilinear(nuLbar[i], diracPR, nuRc[j])
          + Yv[i, j] * (-sp.conjugate(Gp))
          * Bilinear(eLbar[i], diracPR, nuRc[j]))
LYukN += -(sp.conjugate(Yv[i, j]) * H0
           * Bilinear(nuRbar[j], diracPL, nuL[i])
           + sp.conjugate(Yv[i, j]) * (-Gp)
           * Bilinear(nuRbar[j], diracPL, eL[i]))

The quark Yukawas are verbatim from the SM tutorial (color summed in Python, flavor-generic, no CKM), and the gauge currents come from fermion_gauge_current — which, like Dmu, is generic over the reps dict, so every fermion’s \(X\)-boson current (with its symbolic charge) appears automatically.

Yu, Yd = sp.IndexedBase("Yu"), sp.IndexedBase("Yd")
qL_u, qL_d = QL.components[:3], QL.components[3:]
qLbar_u, qLbar_d = QL.bar_components[:3], QL.bar_components[3:]
uRc, uRbar = uR.components, uR.bar_components
dRc, dRbar = dR.components, dR.bar_components


def qcolor(bar_list, gamma, field_list, idx1, idx2):
    return sum(Bilinear(bar_list[c][idx1], gamma, field_list[c][idx2])
               for c in range(3))


LYuk_d = -(Yd[i, j] * Gp * qcolor(qLbar_u, diracPR, dRc, i, j)
           + Yd[i, j] * H0 * qcolor(qLbar_d, diracPR, dRc, i, j))
LYuk_d += -(sp.conjugate(Yd[i, j]) * sp.conjugate(Gp)
            * qcolor(dRbar, diracPL, qL_u, j, i)
            + sp.conjugate(Yd[i, j]) * sp.conjugate(H0)
            * qcolor(dRbar, diracPL, qL_d, j, i))

LYuk_u = -(Yu[i, j] * sp.conjugate(H0) * qcolor(qLbar_u, diracPR, uRc, i, j)
           + Yu[i, j] * (-sp.conjugate(Gp)) * qcolor(qLbar_d, diracPR, uRc, i, j))
LYuk_u += -(sp.conjugate(Yu[i, j]) * H0 * qcolor(uRbar, diracPL, qL_u, j, i)
            + sp.conjugate(Yu[i, j]) * (-Gp) * qcolor(uRbar, diracPL, qL_d, j, i))

current = (fermion_gauge_current(Ll, i) + fermion_gauge_current(eR, i)
           + fermion_gauge_current(nuR, i))
current_quarks = (fermion_gauge_current(QL, i)
                  + fermion_gauge_current(uR, i)
                  + fermion_gauge_current(dR, i))

L = Lagrangian()
L.add((dag(DH) * DH)[0] + (dag(DS) * DS)[0], sector="kinetic")
L.add(-V, sector="potential")
L.add(LYuk + LYukN, sector="yukawa")
L.add(current, sector="yukawa")
L.add(LYuk_d + LYuk_u, sector="yukawa")
L.add(current_quarks, sector="yukawa");

4. Assemble the Model and check invariance

This is the first model in feynlag with symbol-valued gauge charges going through the invariance checker: the infinitesimal-transformation residuals now carry factors of \(a\), \(b\), \(q_S\), and invariance means they cancel identically in those symbols — a strictly stronger statement than for numeric charges (it certifies the whole \((a,b)\) family at once, and would catch any term that’s only invariant for a special charge choice).

model = Model("SM-U1X", gauge_groups=[SU2L, U1Y, U1X, SU3c],
              fields=[H, S, Ll, eR, nuR, QL, uR, dR, W, B, X, G],
              parameters=[gw, g1, gX, gs, aX, bX, qS, v, vX,
                          lam, lamS, lamHS, mu2, muS2], lagrangian=L)

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

InvarianceReport(33 checks, OK)

5. EWSB: two tadpoles

Two VEVs, two tadpole conditions, two internal parameters (\(\mu^2\), \(\mu_S^2\)) to solve for. The portal coupling links them: each \(\mu^2\) picks up a \(\lambda_{HS}\)-term proportional to the other sector’s VEV.

model.solve_tadpoles([mu2, muS2])

{mu_2: (2*lambda*v**2 + lambda_HS*v_X**2)/2,
 mu_{S2}: (lambda_HS*v**2 + 2*lambda_S*v_X**2)/2}

6. The scalar sector: portal mixing and Goldstone counting

The CP-even block over \((h_0, s_0) = (H0_r, S0_r)\):

h0, G0 = sp.Symbol("H0_r", real=True), sp.Symbol("H0_i", real=True)
s0, aX0 = sp.Symbol("S0_r", real=True), sp.Symbol("S0_i", real=True)

Meven = model.mass_matrix([h0, s0])
Meven

\[\begin{split}\displaystyle \left[\begin{matrix}2 \lambda v^{2} & \lambda_{HS} v v_{X}\\\lambda_{HS} v v_{X} & 2 \lambda_{S} v_{X}^{2}\end{matrix}\right]\end{split}\]

Goldstone counting is a sharp consistency check: two neutral gauge bosons are about to become massive (\(Z\), \(Z'\)), so exactly two neutral pseudoscalars must be exactly massless — and indeed the whole CP-odd \((G^0, \mathrm{Im}\,S)\) block vanishes identically after the tadpole substitution (as does the charged-Goldstone mass, eaten by \(W^\pm\)):

Modd = model.mass_matrix([G0, aX0])
Mch = model.mass_matrix([Gp], charged=True)
assert Modd == sp.zeros(2, 2) and sp.simplify(Mch[0, 0]) == 0
Modd

\[\begin{split}\displaystyle \left[\begin{matrix}0 & 0\\0 & 0\end{matrix}\right]\end{split}\]

diagonalize_orthogonal_2x2 with a free angle symbol thetas gives the mass-basis rotation \((h_1, h_2)\) with the defining relation attached (\(h_1 \to h\), the SM Higgs, as \(\theta_s \to 0\)):

h1, h2 = sp.symbols("h1 h2", real=True)
thetas = sp.Symbol("theta_s", real=True)
rot_s = diagonalize_orthogonal_2x2(Meven, [h0, s0], [h1, h2], angle=thetas)
model.rotate(rot_s)
rot_s.angle_relation

\[\displaystyle \tan{\left(2 \theta_{s} \right)} = \frac{\lambda_{HS} v v_{X}}{\lambda v^{2} - \lambda_{S} v_{X}^{2}}\]

7. The neutral gauge sector: one 3×3, two chained 2×2 rotations

With the Higgs charged under all three neutral generators, the \((W_3, B, X)\) mass matrix is a genuine 3×3 — the \(a\)-proportional entries in the third row/column are the Higgs’s doing, and \(\det = 0\) promises exactly one massless photon:

W1, W2, W3 = W.components
B0 = B.components[0]
X0 = X.components[0]

M3 = model.gauge_mass_matrix([W3, B0, X0])
assert sp.simplify(M3.det()) == 0
M3

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{g_{w}^{2} v^{2}}{4} & - \frac{g_{1} g_{w} v^{2}}{4} & - \frac{a_{X} g_{X} g_{w} v^{2}}{4}\\- \frac{g_{1} g_{w} v^{2}}{4} & \frac{g_{1}^{2} v^{2}}{4} & \frac{a_{X} g_{1} g_{X} v^{2}}{4}\\- \frac{a_{X} g_{X} g_{w} v^{2}}{4} & \frac{a_{X} g_{1} g_{X} v^{2}}{4} & \frac{g_{X}^{2} \left(a_{X}^{2} v^{2} + 4 q_{S}^{2} v_{X}^{2}\right)}{4}\end{matrix}\right]\end{split}\]

There is no symbolic 3×3 eigen-diagonalizer worth using (nested radicals — see vacuum/diagonalize.py’s docstring), and none is needed: the photon direction is fixed by symmetry alone. Rotation 1 is the standard Weinberg rotation, untouched by the extension — and it decouples the photon exactly, because the unbroken electric charge doesn’t care about \(X\). What survives is an intermediate \(Z_0\) that still mixes with \(X\):

Z0, A, Z, Zp = sp.symbols("Z0 A Z Zp", real=True)
thetaW = sp.atan(g1.s / gw.s)
model.rotate(Rotation([W3, B0], [Z0, A], rotation_2x2(-thetaW)))

cw, sw = sp.cos(thetaW), sp.sin(thetaW)
R1e = sp.Matrix([[cw, -sw, 0], [sw, cw, 0], [0, 0, 1]])   # (W3,B,X) -> (Z0,A,X)
M3r = sp.simplify(R1e * M3 * R1e.T)
assert M3r[1, :] == sp.zeros(1, 3)   # photon row exactly zero
M3r

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{v^{2} \left(g_{1}^{2} + g_{w}^{2}\right)}{4} & 0 & - \frac{a_{X} g_{X} v^{2} \sqrt{g_{1}^{2} + g_{w}^{2}}}{4}\\0 & 0 & 0\\- \frac{a_{X} g_{X} v^{2} \sqrt{g_{1}^{2} + g_{w}^{2}}}{4} & 0 & g_{X}^{2} \left(\frac{a_{X}^{2} v^{2}}{4} + q_{S}^{2} v_{X}^{2}\right)\end{matrix}\right]\end{split}\]

Rotation 2 diagonalizes the surviving \((Z_0, X)\) block. Note what makes this chaining legal: physical_lagrangian applies registered rotations sequentially, in registration order — each one is a single xreplace pass. The Weinberg rotation leaves the plain symbol Z0 in the Lagrangian, and the second rotation’s substitution consumes exactly that symbol. (Order matters: register Weinberg first. If the rotations were applied as one simultaneous substitution instead, Z0 — absent from the original Lagrangian — would never match anything.)

With a free angle symbol thetap, every downstream coupling comes out in terms of \(\theta'\), and the defining \(\tan 2\theta'\) relation is attached to the rotation:

blockZX = M3r.extract([0, 2], [0, 2])
thetap = sp.Symbol("theta_p", real=True)
rot_zzp = diagonalize_orthogonal_2x2(blockZX, [Z0, X0], [Z, Zp], angle=thetap)
model.rotate(rot_zzp)
rot_zzp.angle_relation

\[\displaystyle \tan{\left(2 \theta_{p} \right)} = \frac{2 a_{X} g_{X} v^{2} \sqrt{g_{1}^{2} + g_{w}^{2}}}{a_{X}^{2} g_{X}^{2} v^{2} - g_{1}^{2} v^{2} + 4 g_{X}^{2} q_{S}^{2} v_{X}^{2} - g_{w}^{2} v^{2}}\]

8. Dual verification, masses, and the W± rotation

Per CONVENTIONS.md, every derived rotation gets verified twice: exactly (substitute the closed-form angle and check the off-diagonals vanish) and numerically at the benchmark point. The benchmark lands where it should — \(m_Z \approx 91.6\) GeV, and \(m_{Z'} \approx g_X q_S v_X = 1000\) GeV up to a small mixing correction:

check_rot = Rotation([Z0, X0], [Z, Zp], rotation_2x2(rot_zzp.angle_solution))
ok, residuals = check_rot.check(blockZX)
assert ok, residuals

bench = {gw.s: 0.6535, g1.s: 0.3580, gX.s: 0.5, aX.s: 0.4,
         qS.s: 1.0, v.s: 246.0, vX.s: 2000.0}
mZ2, mZp2 = [complex(m.subs(bench)).real
             for m in check_rot.masses_squared(blockZX.subs(bench),
                                               simplifier=sp.nsimplify)]
print(f"m_Z  = {mZ2**0.5:.2f} GeV")
print(f"m_Z' = {mZp2**0.5:.2f} GeV")

m_Z  = 91.62 GeV
m_Z' = 1000.31 GeV

Wp, Wm = sp.symbols("Wp Wm")
Umix = sp.Matrix([[1, -sp.I], [1, sp.I]]) / sp.sqrt(2)
model.rotate(Rotation([W1, W2], [Wp, Wm], Umix, kind="unitary"))

\[\begin{split}\displaystyle \left[\begin{matrix}Wp\\Wm\end{matrix}\right] = \left[\begin{matrix}\frac{\sqrt{2} W_{1}}{2} - \frac{\sqrt{2} i W_{2}}{2}\\\frac{\sqrt{2} W_{1}}{2} + \frac{\sqrt{2} i W_{2}}{2}\end{matrix}\right]\end{split}\]

9. Extracting the physical couplings

The boson list now holds eleven physical fields — two CP-even scalars, two pseudoscalars, the charged Goldstones, and five gauge bosons. First a bosonic highlight from the kinetic sector: the \(h_1 Z Z'\) vertex, which exists only because of the mixings (it vanishes as \(\theta' \to 0\) with \(\theta_s \to 0\)… every factor of it is BSM):

Gm, cmap = conjugate_pair(Gp, "Gm")
boson_fields = [h1, h2, G0, aX0, Gp, Gm, Z, A, Zp, Wp, Wm]

rules = model.feynman_rules(boson_fields, sector="kinetic",
                            conjugate_map=cmap, simplifier=sp.simplify)
key = tuple(sorted((h1, Z, Zp), key=lambda f: f.sort_key()))
rules[key]

\[\displaystyle \frac{i \left(- a_{X} g_{X} v \left(g_{1}^{2} + g_{w}^{2}\right)^{3} \cos{\left(2 \theta_{p} \right)} \cos{\left(\theta_{s} \right)} - \frac{g_{1}^{2} g_{w}^{2} v \left(g_{1}^{2} + g_{w}^{2}\right)^{\frac{3}{2}} \left(\sin{\left(2 \theta_{p} - \theta_{s} \right)} + \sin{\left(2 \theta_{p} + \theta_{s} \right)}\right)}{2} + g_{X}^{2} \left(g_{1}^{2} + g_{w}^{2}\right)^{\frac{5}{2}} \left(a_{X}^{2} v \cos{\left(\theta_{s} \right)} + 4 q_{S}^{2} v_{X} \sin{\left(\theta_{s} \right)}\right) \sin{\left(\theta_{p} \right)} \cos{\left(\theta_{p} \right)} - \frac{v \left(g_{1}^{2} + g_{w}^{2}\right)^{\frac{3}{2}} \left(g_{1}^{4} + g_{w}^{4}\right) \left(\sin{\left(2 \theta_{p} - \theta_{s} \right)} + \sin{\left(2 \theta_{p} + \theta_{s} \right)}\right)}{4}\right)}{2 \left(g_{1}^{2} + g_{w}^{2}\right)^{\frac{5}{2}}}\]

Z′ couplings — the headline

physical_lagrangian(sector="yukawa") applies the vacuum shift, both tadpole solutions, and all five registered rotations; extract_fermion_vertices groups the result by (bar, gamma, field). Every \(Z' f f\) coupling should come out as

\[g_{Z'ff} = -\sin\theta'\; g_Z\,(T^3 - Q\,s_W^2)\; +\; \cos\theta'\; g_X\, q_f\]

— the weak-basis \(Z_0\) and \(X\) couplings reshuffled by the mixing angle.

LYuk_phys = model.physical_lagrangian(sector="yukawa")
table = extract_fermion_vertices(LYuk_phys, boson_fields)

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)

sp.simplify(coeff(eRbar[i], gammaR, eRc[i], Zp))   # Z' e_R e_R

\[\displaystyle - a_{X} g_{X} \cos{\left(\theta_{p} \right)} - b_{X} g_{X} \cos{\left(\theta_{p} \right)} - \frac{g_{1}^{2} \sin{\left(\theta_{p} \right)}}{\sqrt{g_{1}^{2} + g_{w}^{2}}}\]

sp.simplify(coeff(nuRbar[i], gammaR, nuRc[i], Zp))  # Z' nu_R nu_R: pure gX*(-b)*cos

\[\displaystyle - b_{X} g_{X} \cos{\left(\theta_{p} \right)}\]

\(\nu_R\) has \(T^3 = Q = Y = 0\): its whole gauge life is the \(X\) charge \(-b\), so its \(Z'\) coupling is pure \(\cos\theta'\, g_X\,(-b)\) — a Z′ portal to otherwise-sterile neutrinos.

Modified Z couplings

The same reshuffle gives every SM fermion a \(\sin\theta' \cdot g_X q_f\) correction to its \(Z\) coupling (\(\theta' \to 0\) recovers the SM — the practical handle electroweak precision data uses to bound \(\theta'\)):

zeR = sp.simplify(coeff(eRbar[i], gammaR, eRc[i], Z))
zeR

\[\displaystyle - a_{X} g_{X} \sin{\left(\theta_{p} \right)} - b_{X} g_{X} \sin{\left(\theta_{p} \right)} + \frac{g_{1}^{2} \cos{\left(\theta_{p} \right)}}{\sqrt{g_{1}^{2} + g_{w}^{2}}}\]

zeR.subs(thetap, 0)    # SM limit: gZ * sw^2

\[\displaystyle \frac{g_{1}^{2}}{\sqrt{g_{1}^{2} + g_{w}^{2}}}\]

Photon couplings — a sanity check

The second rotation touched \(Z_0\) and \(X\) but never \(A\), so photon couplings are exactly the SM electric charges, for any \(a, b, \theta'\) — and the neutrinos (both chiralities) stay dark:

sp.simplify(coeff(eLbar[i], gammaL, eL[i], A))     # -e, exactly

\[\displaystyle - \frac{g_{1} g_{w}}{\sqrt{g_{1}^{2} + g_{w}^{2}}}\]

coeff(nuRbar[i], gammaR, nuRc[i], A)               # 0 identically

\[\displaystyle 0\]

Higgs-portal couplings

The scalar rotation splits every SM Yukawa coupling between the two mass eigenstates: \(h_1 \bar f f \propto \cos\theta_s\) (the SM value, suppressed) and \(h_2 \bar f f \propto \sin\theta_s\) (how the mostly-singlet scalar talks to SM fermions at all):

sp.simplify(coeff(eLbar[i], diracPR, eRc[j], h1))

\[\displaystyle - \frac{\sqrt{2} \cos{\left(\theta_{s} \right)} {Ye}_{i,j}}{2}\]

sp.simplify(coeff(eLbar[i], diracPR, eRc[j], h2))

\[\displaystyle \frac{\sqrt{2} \sin{\left(\theta_{s} \right)} {Ye}_{i,j}}{2}\]

Neutrino Dirac masses

The \(\tilde H\) Yukawa delivers a flavor-generic Dirac mass matrix \(M_\nu = Y_\nu\, v/\sqrt2\) (a Majorana \(\nu_R\) mass from \(S\,\bar\nu_R^c \nu_R\) — i.e. a type-I seesaw — would need charge-conjugation machinery feynlag doesn’t have yet):

M_nu = fermion_mass_matrix(LYukN, nuLbar, nuRc, model.vacuum, 3, (i, j),
                           gamma=diracPR)
M_nu

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{\sqrt{2} v {Yv}_{0,0}}{2} & \frac{\sqrt{2} v {Yv}_{0,1}}{2} & \frac{\sqrt{2} v {Yv}_{0,2}}{2}\\\frac{\sqrt{2} v {Yv}_{1,0}}{2} & \frac{\sqrt{2} v {Yv}_{1,1}}{2} & \frac{\sqrt{2} v {Yv}_{1,2}}{2}\\\frac{\sqrt{2} v {Yv}_{2,0}}{2} & \frac{\sqrt{2} v {Yv}_{2,1}}{2} & \frac{\sqrt{2} v {Yv}_{2,2}}{2}\end{matrix}\right]\end{split}\]

10. Recap

From one charge map \(X = aY + b(B\!-\!L)\) and one singlet scalar, feynlag mechanically derived: symbolic-charge gauge invariance of the whole Lagrangian, two coupled tadpoles, the portal-mixed scalar spectrum with its two exactly-massless Goldstones, the 3×3 neutral gauge mass matrix diagonalized by chained Weinberg and \(Z\)–\(Z'\) rotations with a closed-form \(\tan 2\theta'\), and the full physical coupling table — \(Z'ff\) for every species, the \(\theta'\)-shifted \(Z\) couplings, untouched photon couplings, and \(\cos\theta_s/\sin\theta_s\)-split Higgs couplings.

• Exact pinned formulas for everything here (plus the \(a \to 0\) B−L limit, the W-sector non-renormalization, and the mass invariants): tests/test_u1x.py.

• The full flat script: examples/sm_u1x.py.

• The plain-SM baseline: SM_Feynman_Rules_Tutorial.ipynb.

• Not attempted here: kinetic mixing \(\varepsilon B_{\mu\nu}X^{\mu\nu}\) (no field-strength primitive in the library yet) and UFO export (the UFO particle table needs numeric charges — substitute values for \(a\), \(b\), \(q_S\) first).