Model building with feynlag: suggest + anomaly cancellation

The other tutorials (SM_Feynman_Rules_Tutorial, SM_VLL_Tutorial, SM_U1X_Tutorial) take a Lagrangian you have already written by hand and push it through the analysis pipeline. This notebook goes one step earlier — it shows the model-building tools that help you write a consistent Lagrangian in the first place:

• feynlag.suggest — enumerate every gauge- and discrete-invariant operator your field content admits, up to a mass dimension, so you review a basis instead of hand-deriving invariants;

• feynlag.anomalies — compute the triangle-anomaly coefficients symbolically, turning “is this chiral gauge theory even consistent?” into a set of polynomial conditions on the charges.

We build a small dark sector: a gauged \(U(1)_D\) dark force with a complex dark Higgs \(S\) and a dark fermion \(\chi\). Along the way the tools will tell us which charge assignment is consistent, which operators are allowed, and catch a deliberate mistake.

1. Why enumerate?

The FeynRules-style workflow inverts model building: you write down every term you think is allowed, then hope the invariance checker flags whatever you got wrong or forgot. For anything past the SM that is error-prone — a missing \((S^\dagger S)^2\) quartic, a Yukawa whose charges don’t actually close, a bare mass that a symmetry forbids.

feynlag.suggest turns it around: declare the fields and symmetries, and let the library enumerate the invariant operator basis. The engine is feynlag’s own, already-tested invariance machinery used as an oracle — candidates are built from gauge-invariant blocks, filtered for \(U(1)\) neutrality, projected onto the discrete-invariant subspace, reduced to a linearly independent basis, and then every emitted term is re-verified against the full check_* battery.

import sympy as sp
sp.init_printing()

from feynlag import (
    ExternalParameter, InternalParameter, Lagrangian, Model, Scalar, U1,
    WeylFermion, dag, Dmu, diracPR, diracPL, DiracGamma, Bilinear,
    extract_fermion_vertices, fermion_feynman_rule, fermion_gauge_current,
    fermion_mass_matrix, check_anomaly_free, anomaly_coefficients,
    suggest_potential, suggest_yukawa, suggest_kinetic, build_lagrangian,
)

2. Fields with symbolic charges

We gauge a single abelian factor \(U(1)_D\) with coupling \(g_D\). The dark Higgs \(S\) carries charge \(q_S\) and the dark fermion is a pair of Weyl fields \(\chi_L,\chi_R\) carrying charge \(q_\chi\). We keep both charges symbolic — every pipeline stage (generators, invariance checking, gauge currents, mass matrices) handles symbol-valued \(U(1)\) charges; the one exception is UFO export, whose particle table needs numeric charges.

Declaring \(\chi_L\) and \(\chi_R\) with the same charge \(q_\chi\) makes the dark fermion vector-like (a Dirac fermion). Whether that is forced or a choice is exactly what the anomaly analysis in §3 will tell us.

gD = ExternalParameter("gD", 0.5, positive=True)
U1D = U1("U1D", coupling=gD)

q_chi = sp.Symbol("q_chi", real=True)   # dark-fermion charge (symbolic)
q_S   = sp.Symbol("q_S", real=True)     # dark-Higgs charge (symbolic)

chiL = WeylFermion("chiL", reps={U1D: q_chi}, chirality="L",
                   component_names=["chiL"])
chiR = WeylFermion("chiR", reps={U1D: q_chi}, chirality="R",
                   component_names=["chiR"])
S = Scalar("S", reps={U1D: q_S}, component_names=["S0"])

print("chiL charge:", chiL.reps[U1D])
print("chiR charge:", chiR.reps[U1D], " (same -> vector-like Dirac fermion)")
print("S    charge:", S.reps[U1D])

chiL charge: q_chi
chiR charge: q_chi  (same -> vector-like Dirac fermion)
S    charge: q_S

3. Anomaly cancellation constrains the charges

A chiral \(U(1)\) gauge theory is only consistent if its triangle anomalies vanish. With a single abelian factor the relevant coefficients are the cubic \([U(1)_D]^3\) and the mixed gravitational–gauge \(\mathrm{grav}^2\text{-}U(1)_D\). anomaly_coefficients reduces every fermion to a left-handed Weyl field (a right-handed field contributes with its charge sign flipped) and returns the coefficients as polynomials in the charge symbols.

For our vector-like assignment (\(q_{\chi_L}=q_{\chi_R}=q_\chi\)) the right-handed contribution exactly cancels the left-handed one:

model_vec = Model("darkU1_vec", gauge_groups=[U1D], fields=[chiL, chiR])
coeffs_vec = anomaly_coefficients(model_vec)
coeffs_vec

{'[U1D][U1D][U1D]': 0, 'grav^2-U1D': 0}

report_vec = check_anomaly_free(model_vec)
print(report_vec)
print("nonzero coefficients:", report_vec.nonzero or "none — anomaly-free for ANY q_chi")

AnomalyReport(2 coefficients, anomaly-free)
nonzero coefficients: none — anomaly-free for ANY q_chi

Now suppose we had tried a chiral assignment instead — \(\chi_R\) with a different charge \(q_R\). The anomaly coefficients no longer cancel automatically; they become conditions that pin \(q_R\):

q_R = sp.Symbol("q_R", real=True)
chiR_chiral = WeylFermion("chiRx", reps={U1D: q_R}, chirality="R",
                          component_names=["chiRx"])
model_chiral = Model("darkU1_chiral", gauge_groups=[U1D],
                     fields=[chiL, chiR_chiral])
coeffs_chiral = anomaly_coefficients(model_chiral)
coeffs_chiral

{'[U1D][U1D][U1D]': -q_R**3 + q_chi**3, 'grav^2-U1D': -q_R + q_chi}

solution = sp.solve([coeffs_chiral["[U1D][U1D][U1D]"],
                     coeffs_chiral["grav^2-U1D"]], [q_R], dict=True)
print("anomaly-free requires q_R =", solution[0][q_R])

anomaly-free requires q_R = q_chi

So for a single dark fermion pair, anomaly freedom forces \(q_R = q_\chi\) — the fermion must be vector-like. (A genuinely chiral dark \(U(1)\) needs extra fermions arranged so the cubes and the linear charges each sum to zero, exactly as the SM hypercharges do.) We take the minimal consistent option: a vector-like \(\chi\).

4. Enumerate the invariant operators

With a consistent charge assignment in hand, let suggest write the Lagrangian’s operator content.

Scalar potential. Only \(S^\dagger S\) is a neutral quadratic (a bare \(S\), \(S^2\), \(S^3\) all carry net \(U(1)_D\) charge), so the potential is the familiar Mexican-hat pair:

potential = suggest_potential([S], [U1D])
for t in potential:
    print(f"  dim {t.dim}: ", end="")
    display(t)
potential

  dim 2: 

\[\displaystyle S_{0} \overline{S_{0}}\]

  dim 4: 

\[\displaystyle S_{0}^{2} \overline{S_{0}}^{2}\]

[SuggestedTerm(S0 S0† [potential, dim 2]),
 SuggestedTerm(S0 S0 S0† S0† [potential, dim 4])]

Fermion mass / Yukawa. Because \(\chi\) is vector-like, the gauge-invariant operator is a bare Dirac mass \(\bar\chi_L\chi_R\) (dimension 3) — the dark Higgs cannot give \(\chi\) a Yukawa mass here, since \(\bar\chi_L S\,\chi_R\) would need \(q_S = q_{\chi_L}-q_{\chi_R}=0\). suggest_yukawa returns exactly that bare mass, with its \(+\,\text{h.c.}\) completed:

yukawa = suggest_yukawa([chiL, chiR], [S], [U1D])
for t in yukawa:
    print(f"  dim {t.dim}: ", end="")
    display(t)

  dim 3: 

\[\displaystyle \bar{{\bar{chiL}}_{i}}\,P_R\,{chiR}_{j} + \bar{{\bar{chiR}}_{j}}\,P_L\,{chiL}_{i}\]

The tool catches a mistake. Suppose we fat-fingered \(\chi_R\)’s charge to \(q_\chi + 1\). Now no gauge-invariant mass term exists at all — the enumeration comes back empty, telling us the dark fermion would be massless and the model is incomplete:

chiR_bad = WeylFermion("chiRb", reps={U1D: q_chi + 1}, chirality="R",
                       component_names=["chiRb"])
yukawa_bad = suggest_yukawa([chiL, chiR_bad], [S], [U1D])
print("suggested mass terms for the mistuned charge:", yukawa_bad,
      "\n-> empty: no invariant mass, the model as declared is inconsistent")

suggested mass terms for the mistuned charge: [] 
-> empty: no invariant mass, the model as declared is inconsistent

5. Attach couplings with build_lagrangian

build_lagrangian takes the suggested basis, attaches a fresh coupling to each operator (potential/Yukawa terms enter as \(-c\cdot\mathcal{O}\), matching feynlag’s \(\mathcal{L}\supset -V\) convention; generated couplings carry the right mass dimension \(4-\dim\)), and returns a filled Lagrangian. suggest_kinetic adds the canonical kinetic scaffolding — \((D_\mu S)^\dagger(D^\mu S)\) and the \(\chi\) gauge currents.

kinetic = suggest_kinetic([S, chiL, chiR])
L_auto, gen_params = build_lagrangian(potential + yukawa + kinetic)
print("generated couplings:", [p.name for p in gen_params],
      "  (dims:", [p.unit_dim for p in gen_params], ")")
print("kinetic scaffolding terms:", [t.label for t in kinetic])

generated couplings: ['c1', 'c2', 'c3']   (dims: [2, 0, 1] )
kinetic scaffolding terms: ['(DμS)†(DμS)', 'chiL gauge current', 'chiR gauge current']

That L_auto already forms a consistent model. For the physical pipeline below, though, we adopt feynlag’s standard parametrization (as in every examples/ model): the \(S^\dagger S\) mass coefficient is an InternalParameter fixed by the tadpole condition, rather than a free external. The operator content is identical to what suggest enumerated — we are only choosing which parameter the vacuum equation solves for.

6. Run the pipeline

We build the physical model (charges kept symbolic), break \(U(1)_D\) with \(\langle S\rangle = v_D/\sqrt2\), solve the tadpole, and read off the dark-photon mass, the \(\chi\) Dirac mass, and the \(Z_D\chi\chi\) gauge coupling.

vD   = ExternalParameter("vD", 1000.0, positive=True, unit_dim=1)
lamS = ExternalParameter("lamS", 0.2, positive=True)
muS2 = InternalParameter("muS2", unit_dim=2)     # fixed by the tadpole
Mchi = ExternalParameter("Mchi", 300.0, positive=True, unit_dim=1)

Sp = Scalar("S", reps={U1D: q_S}, component_names=["S0"])
Sp.expand_vev({Sp.components[0]: vD})            # S0 -> (vD + S0_r + i S0_i)/sqrt(2)
cL = WeylFermion("chiL", reps={U1D: q_chi}, chirality="L", component_names=["chiL"])
cR = WeylFermion("chiR", reps={U1D: q_chi}, chirality="R", component_names=["chiR"])
Xb = U1D.bosons("X")                              # the dark photon field

SdS = (dag(Sp) * Sp.mat)[0]
V   = -muS2.s * SdS + lamS.s * SdS**2
DS  = Dmu(Sp)

i, j = sp.symbols("i j", integer=True)
clc, crc = cL.components[0], cR.components[0]
clb, crb = cL.bar_components[0], cR.bar_components[0]
# bare Dirac mass  -M (chibar_L chi_R + h.c.)  -- the operator suggest_yukawa found
Lmass = -(Mchi.s * Bilinear(clb[i], diracPR, crc[i])
          + Mchi.s * Bilinear(crb[i], diracPL, clc[i]))
current = fermion_gauge_current(cL, i) + fermion_gauge_current(cR, i)

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

model = Model("darkU1", gauge_groups=[U1D], fields=[Sp, cL, cR, Xb],
              parameters=[gD, vD, lamS, muS2, Mchi], lagrangian=L)

report = model.check_invariance()
print(report)
report.raise_on_failure()
print("tadpole solution:", model.solve_tadpoles([muS2]))

InvarianceReport(11 checks, OK)
tadpole solution: {muS2: lamS*vD**2}

Dark-photon mass. The dark photon eats the CP-odd part of \(S\) and picks up a mass \(m_{Z_D}^2 = g_D^2\,q_S^2\,v_D^2\) — note the explicit \(q_S\) dependence, carried symbolically straight through the mass matrix:

X0 = Xb.components[0]
M_ZD2 = sp.simplify(model.gauge_mass_matrix([X0])[0, 0])
print("m_ZD^2 =", M_ZD2)

m_ZD^2 = gD**2*q_S**2*vD**2

Dark-fermion mass. The bare Dirac term gives \(\chi\) the mass \(M_\chi\) directly:

M_chi = fermion_mass_matrix(Lmass, clb, crc, model.vacuum, 1, (i, j), gamma=diracPR)
print("chi Dirac mass matrix:", M_chi)

chi Dirac mass matrix: Matrix([[Mchi]])

Dark gauge coupling. Extracting the fermion vertices, the \(Z_D\chi\chi\) coupling is the expected vector current \(i\,g_D\,q_\chi\,\gamma^\mu\) (same on both chiralities — \(\chi\) is vector-like), with the symbolic charge intact:

Lg = sp.expand(model.physical_lagrangian(sector="yukawa"))
table = extract_fermion_vertices(Lg, [X0])
mu = sp.Symbol("mu", integer=True)
gL, gR = DiracGamma(mu) * diracPL, DiracGamma(mu) * diracPR

cLcoup = table.get((clb[i], gL, clc[i]), {}).get(1, {}).get((X0,), 0)
cRcoup = table.get((crb[i], gR, crc[i]), {}).get(1, {}).get((X0,), 0)
print("Z_D chiL chiL :", sp.simplify(fermion_feynman_rule(cLcoup, gL, (X0,))))
print("Z_D chiR chiR :", sp.simplify(fermion_feynman_rule(cRcoup, gR, (X0,))))

Z_D chiL chiL : I*gD*q_chi*gamma(mu)*PL
Z_D chiR chiR : I*gD*q_chi*gamma(mu)*PR

7. Recap

Starting from nothing but field declarations with symbolic charges, the model-building tools:

Step	Tool	What it gave us
Consistency	anomaly_coefficients / check_anomaly_free	anomaly-free ⟺ \(\chi\) vector-like (\(q_R=q_\chi\))
Operator basis	suggest_potential, suggest_yukawa	\(\{S^\dagger S,\,(S^\dagger S)^2\}\) and the bare Dirac mass — and an empty result flagging a mistuned charge
Assembly	build_lagrangian, suggest_kinetic	a filled, coupling-tagged Lagrangian
Physics	the standard pipeline	\(m_{Z_D}^2=g_D^2q_S^2v_D^2\), \(m_\chi=M_\chi\), \(Z_D\chi\chi = i g_D q_\chi\gamma^\mu\)

Every suggested operator is re-verified against the full invariance battery (verify=True by default), so the basis you review is guaranteed consistent.

Where to go next: SM_Feynman_Rules_Tutorial walks the full analysis pipeline for the Standard Model; SM_VLL_Tutorial covers vector-like fermions and mass-basis rotations; SM_U1X_Tutorial extends the SM by a \(U(1)_X\) with tree-level \(Z\)–\(Z'\) mixing. Each hand- writes its Lagrangian — the suggest tools here can bootstrap that first step for a new model of your own.