Discrete symmetries: representations, characters, and the construction of invariants¶
SUN_Groups_Tutorial does this for continuous groups — generators, the Lie
algebra, irreps of any dimension. This is the discrete counterpart: the finite
flavour symmetries (\(\mathbb{Z}_N\), \(S_3\), \(A_4\)) that model builders impose to forbid
terms and force relations between the ones that survive.
The pedagogy is compute step by step, then develop judgment. Nothing below is typed in from the answer. The character table is derived from generator matrices; the number of invariants is predicted by orthogonality and then checked against explicit construction and against the library’s own enumerator.
Each section runs four moves — set up (look at the raw object), collect (let the structure emerge), recognise (name what appeared and why it had to), check (from group theory, not by re-running the same algebra) — and a blockquote asks you to commit to an expectation before the cell that settles it.
The one question this notebook is really about. You impose a flavour symmetry on a model. How many free parameters does the potential have? Which terms are forbidden? You can answer both before writing a single term, and that is what the character table is for.
import itertools
import sympy as sp
from sympy.printing.pretty.stringpict import prettyForm
sp.init_printing(use_latex="mathjax") # text/latex only -- no PNG mimetype
from feynlag import (S3, ZN, Scalar, SU2, U1, ExternalParameter, dag,
check_discrete_invariance)
from feynlag.groups.discrete import DiscreteSymmetry
from feynlag.suggest import reynolds_project, suggest_potential
def ok(msg):
print("\u2713", msg)
def trace(msg, obj=None):
# One line of commentary, optionally followed by a sympy object. Every
# helper below takes `verbose=False` and routes its narration through
# here, so the tracing style is uniform and switched off by default --
# the cells that assert stay quiet, and the "peek" cells do the talking.
print(msg)
if obj is not None:
display(obj)
class MapsTo(sp.Basic):
# lhs |-> rhs, for display only: no equation semantics, so no stray "=".
def __new__(cls, lhs, rhs):
return super().__new__(cls, sp.sympify(lhs), sp.sympify(rhs))
lhs = property(lambda self: self.args[0])
rhs = property(lambda self: self.args[1])
def _latex(self, printer):
return rf"{printer._print(self.lhs)} \mapsto {printer._print(self.rhs)}"
def _sympystr(self, printer):
return f"{printer._print(self.lhs)} |-> {printer._print(self.rhs)}"
def _pretty(self, printer):
left, right = printer._print(self.lhs), printer._print(self.rhs)
return prettyForm(*left.right(" \u21a6 ", right))
1. The generators are the group¶
A finite group is handed to you as a couple of matrices and a promise: multiply them in every possible order and you get the whole group, nothing more and nothing less. For \(S_3\) (the permutations of three objects) the library ships two generators in the real orthogonal basis — a \(2\pi/3\) rotation \(a\) and a reflection \(b\):
Write them down rather than importing them — you should see the matrices you are reasoning about — and then check the library agrees.
Before running the next cells. \(S_3\) has 6 elements. You have 2 matrices. Ask what relations they must satisfy for the closure to stop at 6 rather than running forever: a rotation by \(2\pi/3\) has some order, a reflection has some order, and their product has a third. Guess all three. Then ask the harder question — if you only multiply generators, do you ever need inverses?
# ---- MOVE 1: set up. The three irreps of S3, as matrices. ----------------
c, s = sp.Rational(-1, 2), sp.sqrt(3) / 2
a2 = sp.Matrix([[c, -s], [s, c]]) # rho_2(a): rotation by 2 pi / 3
b2 = sp.Matrix([[1, 0], [0, -1]]) # rho_2(b): reflection
# Keys are ASCII stand-ins for the printed labels: "1p" is 1', "1pp" is 1''
# (a prime is not a legal identifier character, and these are dict keys).
S3_GENS = {
"1": [sp.Matrix([[1]]), sp.Matrix([[1]])], # trivial: everything -> +1
"1p": [sp.Matrix([[1]]), sp.Matrix([[-1]])], # sign rep: odd permutations -> -1
"2": [a2, b2], # the doublet
}
display(sp.Eq(sp.Symbol("rho_2(a)"), a2, evaluate=False))
display(sp.Eq(sp.Symbol("rho_2(b)"), b2, evaluate=False))
# ---- MOVE 2: the defining relations. --------------------------------------
ID2 = sp.eye(2) # NB: `I2` would collide with section 4.1's invariant I_2
rel = {"a^3 = e": sp.simplify(a2**3), "b^2 = e": sp.simplify(b2**2),
"(ab)^2 = e": sp.simplify((a2 * b2)**2)}
for name, M in rel.items():
print(f" {name:<12} ->", "identity" if M == ID2 else "NOT the identity")
assert M == ID2
ok("a^3 = b^2 = (ab)^2 = e -- the presentation of S_3")
# because a^3 = e, a^{-1} = a^2 is itself a product of generators; likewise
# b^{-1} = b. In a FINITE group every inverse is a positive power, so closing
# under multiplication alone is enough -- no inverses needed.
ok("every inverse is a positive power, so multiplication alone closes the group")
a^3 = e -> identity
b^2 = e -> identity
(ab)^2 = e -> identity
✓ a^3 = b^2 = (ab)^2 = e -- the presentation of S_3
✓ every inverse is a positive power, so multiplication alone closes the group
Now close them. Multiply until nothing new appears — that is the group.
# ---- MOVE 3: close the generators into the whole group. -------------------
def close_group(gens, verbose=False):
# All words in `gens`. A word is a tuple of generator indices, so every
# irrep can later be evaluated on the SAME element. (feynlag does this
# internally too, in suggest._group_elements.)
n = gens[0].shape[0]
key = lambda M: tuple(sp.nsimplify(sp.expand(x)) for x in M)
name = lambda w: "".join("abcd"[i] for i in w) or "e"
seen = {key(sp.eye(n)): ()}
frontier = [()]
while frontier:
if verbose:
trace(f"--- frontier: words of length {len(frontier[0]) + 1} ---")
nxt = []
for word in frontier:
for i in range(len(gens)):
cand = word + (i,)
M = sp.eye(n)
for j in cand:
M = sp.expand(M * gens[j])
if key(M) not in seen:
seen[key(M)] = cand
nxt.append(cand)
if verbose:
trace(f" {name(cand):<4} = {cand!s:<8} NEW", M)
elif verbose:
trace(f" {name(cand):<4} = {cand!s:<8} already seen, "
f"as {name(seen[key(M)])}")
frontier = nxt
return list(seen.values())
WORDS = close_group(S3_GENS["2"])
print("words found:", WORDS)
print("|S_3| =", len(WORDS))
assert len(WORDS) == 6
ok("closing two matrices gives exactly 6 elements -- the generators ARE the group")
words found: [(), (0,), (1,), (0, 0), (0, 1), (1, 0)]
|S_3| = 6
✓ closing two matrices gives exactly 6 elements -- the generators ARE the group
# ---- peek: the same closure, narrating itself. ---------------------------
# The call above stays quiet; this one shows the loop. Watch the LAST
# frontier in particular -- that is where termination comes from.
_ = close_group(S3_GENS["2"], verbose=True)
--- frontier: words of length 1 ---
a = (0,) NEW
b = (1,) NEW
--- frontier: words of length 2 ---
aa = (0, 0) NEW
ab = (0, 1) NEW
ba = (1, 0) NEW
bb = (1, 1) already seen, as e
--- frontier: words of length 3 ---
aaa = (0, 0, 0) already seen, as e
aab = (0, 0, 1) already seen, as ba
aba = (0, 1, 0) already seen, as b
abb = (0, 1, 1) already seen, as a
baa = (1, 0, 0) already seen, as ab
bab = (1, 0, 1) already seen, as aa
Six elements from two matrices, and the closure terminated on its own. The trace shows why it terminated: every length-3 word reproduces an element already found, so the frontier empties and the loop stops without anyone telling it the order of the group. That is the whole content of “a group is generated by”: you never store the group, you store the generators and regenerate it on demand.
Recall:
a = \((123)\) (rotation)
b = \((23)\) (reflection)
Word |
Meaning |
Permutation |
|---|---|---|
|
Identity |
\(()\) |
|
\(a\) |
\((123)\) |
|
\(a^2\) |
\((132)\) |
|
\(b\) |
\((23)\) |
|
\(ab\) |
\((12)\) |
|
\(ba\) |
\((13)\) |
Mind the order. rho multiplies a word left to right, so (0,1) is the matrix
\(ab\) — and a matrix acts on a column vector rightmost factor first. So \(ab\) means
“\(b\), then \(a\)”: \(b\) sends \(1\to1\) and \(a\) sends \(1\to2\), giving \(ab=(12)\). Reading it
the other way round gets you \((13)\), which is \(ba\).
The judgment to keep: checking invariance under the generators is checking
invariance under the whole group. That is why feynlag’s
check_discrete_invariance only ever loops over generators — and why it is cheap.
2. A representation is a homomorphism¶
A representation \(\rho\) assigns a matrix to each group element such that multiplication is preserved:
That single equation is the whole definition. Everything else — characters, invariants, selection rules — is a consequence.
Before running the next cell. We have three candidate representations of the same group, of dimensions 1, 1, 2. Ask whether the one-dimensional ones can really be called representations of \(S_3\) — they map six elements onto just \(\pm1\), so they are badly non-injective. Does the definition above require injectivity? And if not, what is the trivial rep telling you?
# ---- MOVE 1 + 2: does rho(g) rho(h) = rho(gh) hold, in every irrep? -------
def rho(label, word, gens=S3_GENS):
M = sp.eye(gens[label][0].shape[0])
for i in word:
M = M * gens[label][i]
return sp.expand(M)
def compose(w1, w2, gens=S3_GENS, verbose=False):
# The word for g h, found by matching rho_2(g) rho_2(h) in the group.
target = sp.expand(rho("2", w1, gens) * rho("2", w2, gens))
if verbose:
trace(f"rho({w1}) =", rho("2", w1, gens))
trace(f"rho({w2}) =", rho("2", w2, gens))
trace(f"their product, rho({w1}) rho({w2}) =", target)
for w in WORDS:
if sp.simplify(rho("2", w, gens) - target) == sp.zeros(2, 2):
if verbose:
trace(f"found in WORDS: that matrix is the element {w}, "
f"so gh = {''.join('abcd'[k] for k in w) or 'e'}")
return w
raise AssertionError("product left the group -- not closed!")
for label in S3_GENS:
for w1, w2 in itertools.product(WORDS, repeat=2):
lhs = sp.expand(rho(label, w1) * rho(label, w2))
rhs = rho(label, compose(w1, w2))
assert sp.simplify(lhs - rhs) == sp.zeros(*lhs.shape), (label, w1, w2)
ok("rho(g) rho(h) = rho(gh) for all 36 pairs, in all three irreps")
✓ rho(g) rho(h) = rho(gh) for all 36 pairs, in all three irreps
Nothing required \(\rho\) to be injective. The trivial rep sends every element to \(1\) and is a perfectly good representation — it is the statement “this object does not transform at all”, i.e. it is an invariant. Hunting for invariants will turn out to be hunting for copies of the trivial rep, which is what §4 does.
The sign rep \(1'\) is the other non-injective one: it remembers only whether a permutation is even or odd. It is not the trivial rep, so a \(1'\) object is not invariant — it flips sign under any reflection. Forgetting that is a standard way to write down a “invariant” that isn’t.
One pair, worked slowly. The assertion above swept 36 pairs past you at once. Take a
single pair, because the convention hiding inside it is the one that trips people later:
\(w_1=(0,)\) — the generator \(a\), a rotation by \(120^\circ\) — and \(w_2=(1,)\), the reflection
\(b\). Rather than reading the arithmetic off a page, hand it to compose in verbose mode:
it multiplies \(\rho(w_1)\rho(w_2)\) and then searches WORDS for the element carrying that
matrix.
# ---- peek: one pair of the 36, showing its work --------------------------
compose((0,), (1,), verbose=True)
rho((0,)) =
rho((1,)) =
their product, rho((0,)) rho((1,)) =
found in WORDS: that matrix is the element (0, 1), so gh = ab
Compare. The left-hand side is the product matrix; the right-hand side is \(\rho\big((0,1)\big)\), the same matrix reached by a different route. They agree, so \(\rho(g)\rho(h)=\rho(gh)\) holds on this pair — and the loop above simply does this 36 times, in all three irreps.
The same statement in permutations. With \(a=(123)\) and \(b=(23)\),
since the rightmost factor acts first: \(b\) fixes \(1\) and then \(a\) sends \(1\to2\); \(b\) sends \(2\to3\) and then \(a\) sends \(3\to1\); \(3\to2\to3\). So the cell has confirmed that the matrix of \((12)\) really is the product of the matrices of \((123)\) and \((23)\).
Order matters here and is easy to get backwards — \(ba=(13)\), not \((12)\). Both rows are in the table above.
2.1 The same matrices, inside feynlag¶
Before going further, tie the matrices above to what the library actually does. Assigning
a multiplet to an irrep and asking for generator_maps() returns the substitution
\(\phi_i \to \sum_k M[i,k]\,\phi_k\) that invariance checking applies.
SU2L = SU2("SU2L", coupling=ExternalParameter("gw", 0.6535, positive=True))
U1Y = U1("U1Y", coupling=ExternalParameter("g1", 0.3580, positive=True))
higgs = lambda name: Scalar(name, {SU2L: 2, U1Y: sp.Rational(1, 2)})
H1, H2, HS = higgs("H1"), higgs("H2"), higgs("HS")
s3 = S3()
s3.assign("2", H1, H2) # (H1, H2) is an S3 doublet
s3.assign("1", HS) # HS is an S3 singlet
maps = s3.generator_maps()
comp = [list(H1.components), list(H2.components)]
for gi, M in enumerate(S3_GENS["2"]):
for i in range(2):
for j in range(len(comp[0])):
want = sum(M[i, k] * comp[k][j] for k in range(2))
assert sp.simplify(maps[gi][comp[i][j]] - want) == 0
ok("feynlag's generator_maps() is exactly phi_i -> sum_k M[i,k] phi_k, "
"with the matrices written above")
print("\n e.g. under generator a:")
display(sp.Eq(sp.Symbol(r"H_{1,1} \mapsto"),
sp.simplify(maps[0][comp[0][0]]), evaluate=False))
✓ feynlag's generator_maps() is exactly phi_i -> sum_k M[i,k] phi_k, with the matrices written above
e.g. under generator a:
Worth looking at the two objects directly before trusting the assertion above. comp is
the component grid — row \(i\) is the \(SU(2)\times U(1)\) components of the \(i\)-th member of
the \(S_3\) doublet — and maps is one substitution dictionary per generator, keyed by those
same component symbols.
comp, maps
And the matrices those dictionaries encode, in generator order: index \(0\) is \(\rho_2(a)\), index \(1\) is \(\rho_2(b)\) — the same two written by hand in §1.
for gi, M in enumerate(S3_GENS["2"]):
print("generator", gi) # the INDEX is a label, not mathematics
display(M)
generator 0
generator 1
3. Conjugacy classes and characters¶
The character of a representation is just the trace, \(\chi_r(g)=\operatorname{Tr}\rho_r(g)\). It throws away almost all the information in the matrices — and that is precisely why it is useful:
it is basis-independent (trace is cyclic, so \(\operatorname{Tr}(U\rho U^{-1})=\operatorname{Tr}\rho\));
it is constant on conjugacy classes, for the same reason;
and the characters of the irreps are orthonormal, which turns “how many invariants are there?” into an inner product.
Before running the next cells. Derive the conjugacy classes of \(S_3\) by hand first — conjugation is relabelling, so classes should group permutations by shape. How many classes do you get? Then note that \(\chi_r(e)=\dim r\) always, and recall \(\sum_r (\dim r)^2 = |G|\). With \(|S_3|=6\) and three irreps, is the set \(\{1,1,2\}\) forced?
# ---- MOVE 1: conjugacy classes, by conjugating every element. -------------
def conjugacy_classes(words, gens, faithful, verbose=False):
# Orbits of g -> h g h^-1. `faithful` names an irrep that separates the
# elements, so a matrix can be turned back into its word.
d = gens[faithful][0].shape[0]
def word_of(M):
for w in words:
if sp.simplify(rho(faithful, w, gens) - M) == sp.zeros(d, d):
return w
raise AssertionError("not a group element")
out, unassigned = [], set(words)
while unassigned:
g = min(unassigned, key=lambda w: (len(w), w))
if verbose:
trace(f"--- seed g = {g} : conjugate it by every h ---")
orbit = set()
for h in words:
Mh = rho(faithful, h, gens)
img = word_of(sp.expand(Mh * rho(faithful, g, gens) * Mh.inv()))
if verbose:
tag = "new to this class" if img not in orbit else ""
trace(f" h = {h!s:<8} h g h^-1 = {img!s:<8} {tag}")
orbit.add(img)
out.append(sorted(orbit, key=lambda w: (len(w), w)))
if verbose:
trace(f" => class of size {len(orbit)} : {out[-1]}\n")
unassigned -= orbit
return out
classes = conjugacy_classes(WORDS, S3_GENS, "2")
print("conjugacy classes (as words in the generators):")
for cl in classes:
print(f" size {len(cl)} :", cl)
print("\n#classes =", len(classes))
conjugacy classes (as words in the generators):
size 1 : [()]
size 2 : [(0,), (0, 0)]
size 3 : [(1,), (0, 1), (1, 0)]
#classes = 3
The classes as an object. The cell above printed them one line at a time; the value
itself is still worth a look, because two later cells lean on its exact shape. It is a list of lists of words — the
same tuples of generator indices §1 closed the group into — and each inner list arrives
sorted by (len(w), w). That sort is not cosmetic: the character table below takes
cl[0] as the representative of each class, so a deterministic “first” is what makes the
table’s column order reproducible rather than incidental.
classes
Sizes \(1+2+3=6\), and the shapes are the ones conjugation should give: the identity alone, the two 3-cycles together, the three transpositions together. The next cell watches a single orbit close, which is where that grouping comes from.
# ---- peek: watch one orbit close under conjugation -----------------------
# Conjugation is relabelling, so the orbit stops as soon as every h maps the
# seed back into elements already collected.
_ = conjugacy_classes(WORDS, S3_GENS, "2", verbose=True)
--- seed g = () : conjugate it by every h ---
h = () h g h^-1 = () new to this class
h = (0,) h g h^-1 = ()
h = (1,) h g h^-1 = ()
h = (0, 0) h g h^-1 = ()
h = (0, 1) h g h^-1 = ()
h = (1, 0) h g h^-1 = ()
=> class of size 1 : [()]
--- seed g = (0,) : conjugate it by every h ---
h = () h g h^-1 = (0,) new to this class
h = (0,) h g h^-1 = (0,)
h = (1,) h g h^-1 = (0, 0) new to this class
h = (0, 0) h g h^-1 = (0,)
h = (0, 1) h g h^-1 = (0, 0)
h = (1, 0) h g h^-1 = (0, 0)
=> class of size 2 : [(0,), (0, 0)]
--- seed g = (1,) : conjugate it by every h ---
h = () h g h^-1 = (1,) new to this class
h = (0,) h g h^-1 = (1, 0) new to this class
h = (1,) h g h^-1 = (1,)
h = (0, 0) h g h^-1 = (0, 1) new to this class
h = (0, 1) h g h^-1 = (1, 0)
h = (1, 0) h g h^-1 = (0, 1)
=> class of size 3 : [(1,), (0, 1), (1, 0)]
# ---- MOVE 2: collect. The character table emerges. ----------------------
chi = {lab: [sp.nsimplify(sp.expand(sp.trace(rho(lab, w)))) for w in WORDS]
for lab in S3_GENS}
reps = [cl[0] for cl in classes] # one element per class
names = ["e", "a (3-cycle)", "b (transposition)"]
hdr = f"{'':<5}" + "".join(f"{n:<22}" for n in names)
print(hdr); print("-" * len(hdr))
for lab in ("1", "1p", "2"):
row = f"{lab:<5}"
for r in reps:
row += f"{chi[lab][WORDS.index(r)]!s:<22}"
print(row)
# constant on classes -- verify, do not assume
for lab in S3_GENS:
for cl in classes:
vals = {chi[lab][WORDS.index(w)] for w in cl}
assert len(vals) == 1, (lab, cl, vals)
ok("every character is constant on each conjugacy class")
e a (3-cycle) b (transposition)
-----------------------------------------------------------------------
1 1 1 1
1p 1 1 -1
2 2 -1 0
✓ every character is constant on each conjugacy class
The objects behind the table. Three of them, and the indexing is the part worth
fixing in your head. WORDS is the element ordering §1’s closure happened to produce, and
everything in this section is indexed by it — so chi[lab] holds one entry per element,
six of them, not one per class. reps is the other half of the bookkeeping: the first word
of each conjugacy class, taken from the classes above.
reps, WORDS, chi
([(), (0,), (1,)],
[(), (0,), (1,), (0, 0), (0, 1), (1, 0)],
{'1': [1, 1, 1, 1, 1, 1],
'1p': [1, 1, -1, 1, -1, -1],
'2': [2, -1, 0, -1, 0, 0]})
So the printed row is a lookup, not a computation: chi[lab][WORDS.index(r)] picks the
one entry per class that the loop just proved is shared by all of them. Printing three
columns out of six is honest only because of the assertion sitting in the same cell.
Read it: \(\mathbf{1'}\) agrees with \(\mathbf{1}\) on the 3-cycles and differs on the transpositions — it is the even/odd sign. And \(\chi_2(b)=0\) says the doublet’s reflection is traceless, which for a \(2\times2\) orthogonal matrix means eigenvalues \(\pm1\): a genuine reflection, not a rotation.
Now the three checks that make the table trustworthy.
# ---- MOVE 3 + 4: check. Three independent constraints. -----------------
G = len(WORDS)
# (a) orthonormality: <chi_r, chi_s> = delta_rs
inner = lambda p, q: sp.simplify(
sum(sp.conjugate(chi[p][i]) * chi[q][i] for i in range(G)) / G)
for p in S3_GENS:
for q in S3_GENS:
assert inner(p, q) == (1 if p == q else 0), (p, q, inner(p, q))
ok("<chi_r, chi_s> = delta_rs -- the irreps are orthonormal")
# (b) sum of squared dimensions is the group order
dims = {lab: S3_GENS[lab][0].shape[0] for lab in S3_GENS}
assert sum(d**2 for d in dims.values()) == G
ok(f"sum_r (dim r)^2 = {G} = |S_3| -- so 1+1+2 is the ONLY option, and the "
"list of irreps is complete")
# (c) #irreps = #conjugacy classes
assert len(S3_GENS) == len(classes) == 3
ok("#irreps = #conjugacy classes = 3")
✓ <chi_r, chi_s> = delta_rs -- the irreps are orthonormal
✓ sum_r (dim r)^2 = 6 = |S_3| -- so 1+1+2 is the ONLY option, and the list of irreps is complete
✓ #irreps = #conjugacy classes = 3
Check (b) is the one worth internalising. \(\sum_r(\dim r)^2=|G|\) with \(|G|=6\) and three classes forces \(\{1,1,2\}\): there is no fourth irrep of \(S_3\) waiting to be found, and no 3-dimensional one. When someone writes an \(S_3\) triplet, they mean a reducible \(1\oplus 2\) or \(1'\oplus 2\) — worth knowing before you try to assign three generations to “the triplet of \(S_3\).”
4. Counting invariants before building them¶
Here is the payoff. An invariant is a copy of the trivial representation. The number of times \(\mathbf{1}\) appears in a representation with character \(\chi\) is the inner product with the trivial character, which is just an average:
And characters multiply under tensor products, \(\chi_{V\otimes W}=\chi_V\chi_W\). So for \(n\) doublets,
One line, no Clebsch–Gordan tables, no construction.
Before running the next cell. Predict \(n=2\): how many ways are there to contract two doublets into a scalar? Then \(n=3\) — is there a cubic invariant at all? (Think about what \(\mathbf 2\otimes\mathbf 2\) decomposes into, and whether any piece can pair with a third doublet.) Then \(n=4\): guess before you look, and notice you probably guessed the number of pairings \(\{12|34\},\{13|24\},\{14|23\}\).
# ---- MOVE 1 + 2: the averaging formula, for n distinct doublets. ----------
def n_invariants(chars, verbose=False, classes=None):
# Copies of the trivial rep in the tensor product with these characters.
total = sp.simplify(sum(sp.prod([c[i] for c in chars]) for i in range(G)) / G)
if verbose:
# Characters are constant on classes (section 3), so the |G|-term
# average collapses to one row per class, weighted by the class size.
rows = ([(cl[0], len(cl)) for cl in classes] if classes
else [(w, 1) for w in WORDS])
trace(f" {'class rep':<10} {'size':>5} prod_i chi_i(g)")
pieces = []
for rep_, size in rows:
v = sp.simplify(sp.prod([c[WORDS.index(rep_)] for c in chars]))
trace(f" {rep_!s:<10} {size:5d} {v}")
pieces.append(f"{size}*{v}")
trace(f" (1/{G})[{' + '.join(pieces)}] = {total}")
return total
for n in range(2, 7):
print(f" invariants in 2^(x{n}) ({n} DISTINCT doublets) : "
f"{n_invariants([chi['2']] * n)}")
invariants in 2^(x2) (2 DISTINCT doublets) : 1
invariants in 2^(x3) (3 DISTINCT doublets) : 1
invariants in 2^(x4) (4 DISTINCT doublets) : 3
invariants in 2^(x5) (5 DISTINCT doublets) : 5
invariants in 2^(x6) (6 DISTINCT doublets) : 11
Every element, then every class. The n = 4 row above is a sum of \(|G|=6\) terms, one
per element. The first cell below prints all six; the second groups them by conjugacy
class. Same number both times — the characters are constant on classes (§3), so each
class contributes its size times one shared product.
# ---- peek: the same quartic count, all six elements ------------------------
# With no `classes`, the trace walks every element of WORDS -- the raw |G|-term
# average, exactly as the formula is written. The next cell collapses it.
n_invariants([chi["2"]] * 4, verbose=True)
class rep size prod_i chi_i(g)
() 1 16
(0,) 1 1
(1,) 1 0
(0, 0) 1 1
(0, 1) 1 0
(1, 0) 1 0
(1/6)[1*16 + 1*1 + 1*0 + 1*1 + 1*0 + 1*0] = 3
# ---- peek: the quartic count, one class at a time ------------------------
# The n = 4 row above is just this average. Three quartics, and you can see
# which classes contributed.
n_invariants([chi["2"]] * 4, verbose=True, classes=classes)
class rep size prod_i chi_i(g)
() 1 16
(0,) 2 1
(1,) 3 0
(1/6)[1*16 + 2*1 + 3*0] = 3
So: one quadratic, one cubic, three quartic.
The cubic is worth a pause. A single \(S_3\) doublet admits a cubic invariant — that is the
\(\mathbb{Z}_3\)-covariant \(\mathrm{Re}\big[(r_1+ir_2)^3\big]=r_1^3-3r_1r_2^2\) that drives the
whole 3HDM-\(S_3\) scalar potential (see THDM_S3_Tutorial). It exists because \(\chi_2\) is
real and \(\chi_2(a)=-1\), so the average survives.
But there is a trap in the quartic count, and it is the single most common invariant-counting mistake in flavour model building.
4.1 Distinct multiplets versus one multiplet¶
\(n_{\text{inv}}(\mathbf 2^{\otimes4})=3\) counts invariants built from four independent doublets \(x,y,z,w\). If instead you have one doublet \(x\) and want quartic terms \(\sim x^4\), the four slots are no longer independent — the product is symmetrised, and you must count in \(\mathrm{Sym}^4(\mathbf 2)\), not \(\mathbf 2^{\otimes4}\).
The tool is the Molien series, which generates the symmetric-power counts at once:
Before running the next cell. You already know a single \(S_3\) doublet has one quadratic and one cubic invariant. Give them names, since the rest of this section uses them:
\[I_2 \equiv x_1^2+x_2^2, \qquad I_3 \equiv x_1^3-3x_1x_2^2\]— the subscript is the degree, and \(I_3\) is §4’s cubic \(\mathrm{Re}\big[(x_1+ix_2)^3\big]\) with \(r\to x\). (The other cubic combination, \(x_2^3-3x_2x_1^2\), is not invariant: it is even under \(a\) but flips sign under the reflection \(b\), making it a \(\mathbf 1'\) covariant — the same antisymmetry §5.1 returns to.)
If those two are the only generators of the invariant ring, what is the quartic count — and at which degree does something genuinely new first appear? Predict the answer before looking, then check the series against the guess “everything is a polynomial in \(I_2\) and \(I_3\)”.
# ---- MOVE 2: the Molien series for ONE doublet. --------------------------
t = sp.Symbol("t")
molien = sp.simplify(
sum(1 / sp.det(sp.eye(2) - t * rho("2", w)) for w in WORDS) / G)
display(sp.Eq(sp.Symbol("M(t)"), sp.factor(molien), evaluate=False))
# sp.factor groups it as (t-1)^2 (t+1) (t^2+t+1); the form quoted below is
# (1-t^2)(1-t^3). Same series -- so PROVE it rather than asserting by eye.
claimed = 1 / ((1 - t**2) * (1 - t**3))
assert sp.simplify(molien - claimed) == 0
print(" the same series, regrouped as one quadratic x one cubic generator:")
display(sp.Eq(sp.Symbol("M(t)"), claimed, evaluate=False))
ser = sp.expand(sp.series(molien, t, 0, 9).removeO())
counts = [(n, ser.coeff(t, n)) for n in range(9)]
print("\n degree :", "".join(f"{n:4d}" for n, _ in counts))
print(" #inv :", "".join(f"{c!s:>4}" for _, c in counts))
the same series, regrouped as one quadratic x one cubic generator:
degree : 0 1 2 3 4 5 6 7 8
#inv : 1 0 1 1 1 1 2 1 2
# ---- peek: the Molien average, class by class ----------------------------
# Same collapse as above: 1/det(1 - t rho(g)) is a class function too.
acc = 0
for cl in classes:
term = sp.simplify(1 / sp.det(sp.eye(2) - t * rho("2", cl[0])))
acc += len(cl) * term
trace(f" class {cl[0]!s:<8} (size {len(cl)}) contributes", term)
trace(f" the weighted sum over {G} elements, divided by |G| = {G} :",
sp.factor(sp.simplify(acc / G)))
assert sp.simplify(acc / G - molien) == 0
ok("the class-by-class average rebuilds M(t) -- nothing was skipped")
class () (size 1) contributes
class (0,) (size 2) contributes
class (1,) (size 3) contributes
the weighted sum over 6 elements, divided by |G| = 6 :
✓ the class-by-class average rebuilds M(t) -- nothing was skipped
— the invariant ring of one \(S_3\) doublet is freely generated by one quadratic and one cubic. So the degree counts are the number of ways to write \(n=2i+3j\): one invariant at degrees 2, 3, 4, 5, and the first two at degree 6 (\(I_2^3\) and \(I_3^2\)).
The quartic count is therefore 1, not 3. Same group, same irrep — the answer depends entirely on whether the four legs are distinct fields or one field repeated.
# ---- MOVE 3 + 4: recognise, then check against the free-ring claim. ------
free_series = sp.series(1 / ((1 - t**2) * (1 - t**3)), t, 0, 9).removeO()
assert sp.simplify(sp.expand(free_series - ser)) == 0
ok("Molien = 1/((1-t^2)(1-t^3)) exactly -- one quadratic and one cubic "
"generator, freely")
quartic_one = ser.coeff(t, 4)
quartic_four = n_invariants([chi["2"]] * 4)
print("\n quartics from ONE doublet :", quartic_one)
print(" quartics from FOUR doublets :", quartic_four)
assert quartic_one == 1 and quartic_four == 3
ok("1 versus 3 -- distinctness of the legs changes the answer, and neither "
"number is a typo")
✓ Molien = 1/((1-t^2)(1-t^3)) exactly -- one quadratic and one cubic generator, freely
quartics from ONE doublet : 1
quartics from FOUR doublets : 3
✓ 1 versus 3 -- distinctness of the legs changes the answer, and neither number is a typo
5. Clebsch–Gordan: actually building the invariant¶
Counting says how many. To write them down you need the decomposition. For \(S_3\) in the real orthogonal basis,
and feynlag ships the coefficients as S3.doublet_product:
Before running the next cell. Check the dimensions add up: \(2\times2=1+1+2\). ✓. Now the important one — set \(y=x\) (one doublet, not two) and evaluate the \(\mathbf 1'\) piece by hand. What do you get, and what does that imply for how many invariants a single doublet can form?
# ---- MOVE 1: the three channels, and how they transform. -----------------
x = sp.symbols("x1 x2"); y = sp.symbols("y1 y2")
cg = S3.doublet_product
P = cg(x, y)
for ch in ("1", "1p", "2"):
lab = {"1": "1", "1p": "1'", "2": "2"}[ch] # 1p is ASCII for 1'
display(sp.Eq(sp.Symbol(rf"\{{x \otimes y\}}_{{{lab}}}"),
sp.Matrix(P[ch]) if ch == "2" else P[ch], evaluate=False))
What doublet_product hands back. The display above is formatted; the return value is
a plain dictionary keyed by the irrep labels "1", "1p" and "2". The asymmetry in it
matters for everything that follows.
P
{'1': x1*y1 + x2*y2, '1p': x1*y2 - x2*y1, '2': (x1*y1 - x2*y2, -x1*y2 - x2*y1)}
The two singlet channels are bare expressions, while "2" is a tuple of two — a doublet,
in exactly the form cg accepts as an argument. That is why the cell above wraps only the
"2" entry in sp.Matrix to display it, and why §5.2 can write
cg(cg(x, x)["2"], cg(x, x)["2"])["1"]: the doublet channel feeds straight back in as if
it were a fresh pair of fields.
# ---- MOVE 2 + 4: each channel must transform as its label claims. --------
def act(vec, M):
# The doublet (v1, v2) under a generator matrix M.
return tuple(sum(M[i, k] * vec[k] for k in range(2)) for i in range(2))
for gi, M in enumerate(S3_GENS["2"]):
xp, yp = act(x, M), act(y, M)
Pp = cg(xp, yp)
# the singlet is untouched
assert sp.simplify(sp.expand(Pp["1"] - P["1"])) == 0
# the 1' picks up det(M) = +1 for the rotation, -1 for the reflection
sign = sp.simplify(M.det())
assert sp.simplify(sp.expand(Pp["1p"] - sign * P["1p"])) == 0
# the doublet channel transforms as a doublet
want = act(P["2"], M)
assert all(sp.simplify(sp.expand(Pp["2"][i] - want[i])) == 0 for i in range(2))
ok("the three channels transform as 1, 1' (with sign = det) and 2 -- verified, "
"not quoted")
✓ the three channels transform as 1, 1' (with sign = det) and 2 -- verified, not quoted
5.1 The \(\mathbf 1'\) trap¶
Now the thing that catches people.
# ---- the antisymmetric singlet vanishes on a repeated doublet ------------
same = cg(x, x)
display(sp.Eq(sp.Symbol(r"\{x \otimes x\}_{1}"),
sp.expand(same["1"]), evaluate=False))
display(sp.Eq(sp.Symbol(r"\{x \otimes x\}_{1'}"),
sp.expand(same["1p"]), evaluate=False))
print(" ^ identically zero")
assert sp.expand(same["1p"]) == 0
ok("1' is ANTISYMMETRIC, so it vanishes when both legs are the same doublet")
^ identically zero
✓ 1' is ANTISYMMETRIC, so it vanishes when both legs are the same doublet
\(x_1x_2-x_2x_1=0\). Obvious once written — and invisible if you copy a Clebsch–Gordan table into a model with one doublet and start counting terms. It is the same fact §4.1 measured from the other side: with one doublet, channels that need two distinct multiplets simply are not available.
Model-building consequence. The number of terms in your potential depends on how many distinct multiplets you have, not just on the group. Two \(S_3\) doublets \((H_1,H_2)\) and \((H_3,H_4)\) give you \(\mathbf 1'\) couplings that a single doublet cannot form.
Now construct the quartics and confirm both counts from §4.
# ---- MOVE 3: construct, and check the rank against the prediction. -------
z = sp.symbols("z1 z2"); w = sp.symbols("w1 w2")
D = {"x": x, "y": y, "z": z, "w": w}
allv = [*x, *y, *z, *w]
def rank_of(exprs, gens, verbose=False, labels=None):
monos = sorted({m for e in exprs for m in sp.Poly(e, *gens).monoms()})
Mx = sp.Matrix([[sp.Poly(e, *gens).coeff_monomial(m) for m in monos]
for e in exprs])
if verbose:
# One row per candidate, one column per monomial: the rank of that
# matrix IS the number of independent invariants.
trace(f" coefficient matrix: {Mx.rows} candidates x "
f"{Mx.cols} monomials", Mx)
names = labels or [f"row {k}" for k in range(Mx.rows)]
keep, r = [], 0 # a maximal independent subset, in order
for k in range(Mx.rows):
if Mx[keep + [k], :].rank() > r:
keep.append(k); r += 1
trace(" independent : " + ", ".join(names[k] for k in keep))
trace(" dependent : " + ", ".join(n for k, n in enumerate(names)
if k not in keep))
return Mx.rank()
cands, labels = [], []
for (p, q), (r_, u) in itertools.combinations(
list(itertools.combinations("xyzw", 2)), 2):
if set(p + q) | set(r_ + u) != set("xyzw"):
continue
for ch in ("1", "1p"): # pair each channel with itself
cands.append(sp.expand(cg(D[p], D[q])[ch] * cg(D[r_], D[u])[ch]))
labels.append(f"({p}{q})_{ch} ({r_}{u})_{ch}")
d1, d2 = cg(D[p], D[q])["2"], cg(D[r_], D[u])["2"]
cands.append(sp.expand(cg(d1, d2)["1"])) # doublet x doublet -> singlet
labels.append(f"({p}{q})_2 . ({r_}{u})_2")
print(f"built {len(cands)} candidate quartics from four distinct doublets:")
for l in labels:
print(" ", l)
r4 = rank_of(cands, allv)
print("\n independent among them:", r4, f" (characters predicted {quartic_four})")
assert r4 == quartic_four == 3
ok("9 candidates, rank 3 -- six relations, and the count matches the character "
"prediction exactly")
built 9 candidate quartics from four distinct doublets:
(xy)_1 (zw)_1
(xy)_1p (zw)_1p
(xy)_2 . (zw)_2
(xz)_1 (yw)_1
(xz)_1p (yw)_1p
(xz)_2 . (yw)_2
(xw)_1 (yz)_1
(xw)_1p (yz)_1p
(xw)_2 . (yz)_2
independent among them: 3 (characters predicted 3)
✓ 9 candidates, rank 3 -- six relations, and the count matches the character prediction exactly
# ---- peek: WHICH three, and which six are the relations ------------------
_ = rank_of(cands, allv, verbose=True, labels=labels)
coefficient matrix: 9 candidates x 8 monomials
independent : (xy)_1 (zw)_1, (xy)_1p (zw)_1p, (xy)_2 . (zw)_2
dependent : (xz)_1 (yw)_1, (xz)_1p (yw)_1p, (xz)_2 . (yw)_2, (xw)_1 (yz)_1, (xw)_1p (yz)_1p, (xw)_2 . (yz)_2
# ---- and the one-doublet case, against the Molien prediction -------------
one = [sp.expand(cg(x, x)["1"]**2),
sp.expand(cg(cg(x, x)["2"], cg(x, x)["2"])["1"])]
display(sp.Eq(sp.Symbol(r"(x\cdot x)^2"), one[0], evaluate=False))
display(sp.Eq(sp.Symbol(r"\{xx\}_2\cdot\{xx\}_2"), one[1], evaluate=False))
r1 = rank_of(one, list(x))
print("\n independent:", r1, f" (Molien predicted {quartic_one})")
assert r1 == quartic_one == 1
ok("two ways to write it, one invariant -- the doublet channel gives nothing new")
independent: 1 (Molien predicted 1)
✓ two ways to write it, one invariant -- the doublet channel gives nothing new
Both predictions land. This is the loop worth internalising:
count with characters → construct with Clebsch–Gordan → check the rank.
If the rank comes out below the character count you have made an algebra slip; if it comes out above, your “invariants” are not all invariant. Either way the disagreement is informative, which is what makes the check worth doing.
6. The Reynolds projector: invariants without tables¶
Clebsch–Gordan coefficients are elegant but they must be derived per group and they
depend on the basis. There is a blunter tool that always works, needs no tables, and is
what feynlag actually uses — group averaging:
\(P\) is a projector onto the invariant subspace: averaging something already invariant returns it unchanged, and averaging anything else either produces its invariant part or kills it outright.
Before running the next cell. Apply \(P\) to a monomial with a definite \(\mathbb{Z}_N\) charge. The average is \(\frac1N\sum_k \omega^{kq}m\) — a geometric series. For which charges \(q\) is it \(m\), and for which is it \(0\)? Convince yourself this reproduces the familiar selection rule without anyone having stated the rule.
# ---- Z_N: the Reynolds average IS the charge selection rule -------------
Z3 = ZN("Z3", 3)
A, B = Scalar("A", real=True), Scalar("B", real=True)
Z3.assign(1, A) # charge 1
Z3.assign(2, B) # charge 2
sA, sB = A.components[0], B.components[0]
for expr, why in ((sA**3, "charge 3 = 0 mod 3"), (sA * sB, "charge 1+2 = 0 mod 3"),
(sA**2, "charge 2 != 0 mod 3"), (sA * sB**2, "charge 5 = 2 != 0")):
# expand_complex: sympy leaves the omega average as
# (1 - (-1)^(1/3) + (-1)^(2/3))/3, which IS zero but does not display so
P_ = sp.simplify(sp.expand_complex(reynolds_project(expr, [Z3])))
print(f" P({expr!s:<8}) = {P_!s:<8} {why}")
ok("the Reynolds average reproduces 'total charge = 0 mod N' with no rule "
"written down")
P(A**3 ) = A**3 charge 3 = 0 mod 3
P(A*B ) = A*B charge 1+2 = 0 mod 3
P(A**2 ) = 0 charge 2 != 0 mod 3
P(A*B**2 ) = 0 charge 5 = 2 != 0
✓ the Reynolds average reproduces 'total charge = 0 mod N' with no rule written down
# ---- S3: averaging produces genuine invariant COMBINATIONS --------------
X1, X2 = Scalar("X1", real=True), Scalar("X2", real=True)
s3b = S3(); s3b.assign("2", X1, X2)
r1_, r2_ = X1.components[0], X2.components[0]
for expr in (r1_**2, r1_ * r2_, r1_**3, r1_**2 * r2_, r1_**2 + r2_**2):
display(sp.Eq(sp.Function("P")(expr),
sp.simplify(sp.expand(reynolds_project(expr, [s3b]))),
evaluate=False))
The average, done by hand. The cell above asked reynolds_project for \(P(r_1^2)\). Here
is the formula \(P(m)=\frac1{|G|}\sum_g g\cdot m\) with nothing hidden: for each of the six
elements of WORDS, take its doublet matrix, substitute \(r_i\mapsto\sum_k\rho_2(g)_{ik}\,r_k\)
into \(m=r_1^2\), and add the result to a running sum. The substitution must be
simultaneous — sequential subs would apply the \(r_2\) rule to the \(r_2\) that the \(r_1\)
rule just introduced — which is why the cell uses xreplace.
# ---- peek: P(r1^2), one group element at a time ----------------------------
running = sp.S.Zero
for wd in WORDS: # `w` is taken: section 5's fourth doublet
M = rho("2", wd) # the doublet matrix for g
sub = {r1_: M[0, 0] * r1_ + M[0, 1] * r2_,
r2_: M[1, 0] * r1_ + M[1, 1] * r2_}
term = sp.expand((r1_**2).xreplace(sub)) # simultaneous substitution
running += term
trace(f"--- g = {wd} ---", sp.Eq(sp.Symbol(r"\rho(g)"), M, evaluate=False))
display(MapsTo(r1_, sub[r1_]), MapsTo(r2_, sub[r2_]))
display(sp.Eq(sp.Symbol(r"g \cdot m"), term, evaluate=False))
display(sp.Eq(sp.Symbol(r"\mathrm{running\ sum}"), running, evaluate=False))
avg = sp.expand(running / G)
display(sp.Eq(sp.Function("P")(r1_**2), avg, evaluate=False))
assert sp.expand(avg - reynolds_project(r1_**2, [s3b])) == 0
ok("the average by hand matches reynolds_project: P(X1^2) = (X1^2 + X2^2)/2")
--- g = () ---
--- g = (0,) ---
--- g = (1,) ---
--- g = (0, 0) ---
--- g = (0, 1) ---
--- g = (1, 0) ---
✓ the average by hand matches reynolds_project: P(X1^2) = (X1^2 + X2^2)/2
Watch the running sum rather than the last line. The cross term \(\sqrt3\,X_1X_2\) appears, flips sign and cancels between elements, and the \(X_2^2\) piece is fed in by the elements that rotate \(X_1\) into \(X_2\). No single \(g\cdot X_1^2\) is proportional to \(X_1^2+X_2^2\) — only the sum over the whole group is.
Note the difference in character between the two groups. For \(\mathbb{Z}_N\) every monomial is an eigenvector of the projector — it survives whole or dies. For \(S_3\) the projector mixes monomials: \(P(r_1^2)=\tfrac12(r_1^2+r_2^2)\) is a combination that no single monomial realises. That is exactly why an abelian symmetry can be enforced by bookkeeping charges, while a non-abelian one cannot.
Check it against §5’s Clebsch–Gordan answer: the projections must lie in the span the CG construction produced.
# ---- MOVE 4: Reynolds and Clebsch-Gordan must agree ---------------------
I2_cg = cg((r1_, r2_), (r1_, r2_))["1"] # r1^2 + r2^2
I3_cg = sp.expand(cg((r1_, r2_), cg((r1_, r2_), (r1_, r2_))["2"])["1"])
P2 = sp.expand(reynolds_project(r1_**2, [s3b]))
P3 = sp.expand(reynolds_project(r1_**3, [s3b]))
for P_, I_, name in ((P2, I2_cg, "quadratic"), (P3, I3_cg, "cubic")):
ratio = sp.simplify(sp.cancel(P_ / I_))
print(f" {name:<9} : P = ({ratio}) x the CG invariant")
assert ratio.is_number and ratio != 0
ok("the Reynolds projections are proportional to the CG invariants -- two "
"independent routes to the same subspace")
# and the projector really is a projector: P(P(m)) = P(m)
assert sp.simplify(sp.expand(reynolds_project(P2, [s3b]) - P2)) == 0
ok("P^2 = P")
quadratic : P = (1/2) x the CG invariant
cubic : P = (1/4) x the CG invariant
✓ the Reynolds projections are proportional to the CG invariants -- two independent routes to the same subspace
✓ P^2 = P
7. Basis dependence — and why we counted with characters¶
The \(S_3\) doublet has (at least) two bases in common use: the real orthogonal one this library commits to, and a complex one in which the \(2\pi/3\) rotation is diagonal, \(\rho(a)=\mathrm{diag}(\omega,\omega^2)\) with \(\omega=e^{2\pi i/3}\). Papers use both, often without saying which — and every Clebsch–Gordan coefficient differs between them.
Before running the next cell. A change of basis is \(\rho\to U\rho U^{-1}\). Which of the following survive it: the individual matrix entries, the trace, the determinant, the eigenvalues, the number of invariants? Decide before looking, then note which of those the counting in §4 relied on.
# ---- the same rep, a different basis ------------------------------------
U = sp.Matrix([[1, sp.I], [1, -sp.I]]) / sp.sqrt(2)
a2c = sp.simplify(U * a2 * U.inv())
b2c = sp.simplify(U * b2 * U.inv())
print("complex basis:")
display(sp.Eq(sp.Symbol(r"\rho'_2(a)"), a2c, evaluate=False))
display(sp.Eq(sp.Symbol(r"\rho'_2(b)"), b2c, evaluate=False))
COMPLEX = {"1": S3_GENS["1"], "1p": S3_GENS["1p"], "2": [a2c, b2c]}
chi_c = {lab: [sp.nsimplify(sp.expand(sp.trace(rho(lab, w, COMPLEX))))
for w in WORDS] for lab in COMPLEX}
print("\n chi_2 (real basis) :", chi["2"])
print(" chi_2 (complex basis) :", chi_c["2"])
assert all(sp.simplify(chi["2"][i] - chi_c["2"][i]) == 0 for i in range(G))
ok("the matrices changed completely; the character did not")
complex basis:
chi_2 (real basis) : [2, -1, 0, -1, 0, 0]
chi_2 (complex basis) : [2, -1, 0, -1, 0, 0]
✓ the matrices changed completely; the character did not
# ---- so every count from section 4 is unchanged -------------------------
def n_inv_basis(chars, n):
return sp.simplify(sum(chars[i]**n for i in range(G)) / G)
for n in (2, 3, 4):
a_, b_ = n_inv_basis(chi["2"], n), n_inv_basis(chi_c["2"], n)
print(f" 2^(x{n}): real basis {a_} complex basis {b_}")
assert a_ == b_
ok("identical invariant counts in both bases -- the PHYSICS cannot depend on "
"the basis, and characters are the basis-independent way to see it")
2^(x2): real basis 1 complex basis 1
2^(x3): real basis 1 complex basis 1
2^(x4): real basis 3 complex basis 3
✓ identical invariant counts in both bases -- the PHYSICS cannot depend on the basis, and characters are the basis-independent way to see it
The judgment: counts are basis-independent, coefficients are not. When you take a Clebsch–Gordan table from a paper you are importing its basis convention along with it, and mixing two papers’ tables silently gives wrong relative factors. The count is the safe thing to compare across sources; the coefficients are not.
feynlag states its convention explicitly in S3.doublet_product’s docstring — the real
orthogonal basis, with the \(\mathbf 2\) channel derived from \(D=\overline{w\,v}\), where
\(w=x_1+ix_2\) and \(v=y_1+iy_2\) package the two doublets as complex numbers. (These \(w,v\)
are local to that formula — not §5’s fourth doublet \(w\), and not the \(\omega\) of §10.)
8. Complex fields and the conjugate leg¶
Real scalars are the easy case. A realistic model has complex fields, and then a term like \(\bar\psi\,\psi\) has one leg transforming with \(M\) and the other with something else. Getting that “something else” wrong is a silent error — the term looks invariant and is not.
Requiring \(\sum_i \bar\psi'_i\psi'_i = \sum_i\bar\psi_i\psi_i\) under \(\psi_i\to\sum_k M[i,k]\psi_k\) and \(\bar\psi_i\to\sum_k X[i,k]\bar\psi_k\) forces \(X^{\mathsf T}M=\mathbb 1\), i.e.
Before running the next cell. For which kind of \(M\) does \(X\) coincide with \(M\) itself? Check your answer against the two groups in play: \(S_3\)’s irreps are real orthogonal, \(\mathbb{Z}_N\)’s are complex phases. Predict which one has \(X\neq M\) — and what would break if you used \(M\) on both legs there.
# ---- X = (M^-1)^T, computed for both groups ------------------------------
print("S_3 doublet:")
for gi, M in enumerate(S3_GENS["2"]):
X = sp.simplify(M.inv().T)
print(f" generator {gi} : X == M ?", sp.simplify(X - M) == sp.zeros(2, 2),
" (M orthogonal?", sp.simplify(M.T * M) == sp.eye(2), ")")
assert sp.simplify(X - M) == sp.zeros(2, 2)
ok("for S_3, X = M exactly -- because the irreps are REAL ORTHOGONAL")
print("\nZ_3 charge-1 irrep:")
M = sp.Matrix([[sp.exp(2 * sp.pi * sp.I / 3)]])
X = sp.simplify(M.inv().T)
display(sp.Eq(sp.Symbol("M"), M[0], evaluate=False))
display(sp.Eq(sp.Symbol("X"), X[0], evaluate=False))
print(" X == conj(M)?", sp.simplify(X[0] - sp.conjugate(M[0])) == 0)
assert sp.simplify(X[0] - M[0]) != 0
ok("for Z_3, X = conj(M) != M -- using M on both legs would give charge 2, "
"not 0, and the 'invariant' would not be one")
S_3 doublet:
generator 0 : X == M ? True (M orthogonal? True )
generator 1 : X == M ? True (M orthogonal? True )
✓ for S_3, X = M exactly -- because the irreps are REAL ORTHOGONAL
Z_3 charge-1 irrep:
X == conj(M)? True
✓ for Z_3, X = conj(M) != M -- using M on both legs would give charge 2, not 0, and the 'invariant' would not be one
This is not a hypothetical: it is exactly what DiscreteSymmetry.fermion_generator_data
computes, via M.inv().T rather than by assuming. The library derives it because the two
cases genuinely differ, and the \(S_3\) coincidence would otherwise hide the bug until
someone used \(\mathbb{Z}_N\).
Let the invariance checker demonstrate the failure directly.
# ---- a wrong assignment must FAIL the check -----------------------------
C_, D_ = Scalar("C", real=True), Scalar("D", real=True)
z3b = ZN("Z3b", 3); z3b.assign(1, C_); z3b.assign(1, D_) # both charge +1
sC, sD = C_.components[0], D_.components[0]
good, bad = sC * sD**2, sC * sD # charges 3 = 0 mod 3, and 2 != 0
for expr in (good, bad):
inv, residual = check_discrete_invariance(expr, z3b)
display(expr)
print(" invariant?", inv)
if not inv:
print(" residual:", end=" ")
display(sp.simplify(sp.expand_complex(residual[0][1])))
assert check_discrete_invariance(good, z3b)[0]
assert not check_discrete_invariance(bad, z3b)[0]
ok("the checker accepts total charge 0 mod 3 and rejects the rest, with the "
"offending residual shown")
invariant? True
invariant? False
residual:
✓ the checker accepts total charge 0 mod 3 and rejects the rest, with the offending residual shown
9. Into a model: the \(S_3\) three-Higgs-doublet potential¶
Everything so far was group theory. Now the model-building question: assign three Higgs doublets to \(S_3\) as \(\mathbf 2\oplus\mathbf 1\) — the doublet \((H_1,H_2)\) and the singlet \(H_S\) — and ask how many terms the potential has.
The fields are complex \(SU(2)\) doublets, so the potential is built from gauge-invariant bilinears \(H_a^\dagger H_b\), and both conjugate and non-conjugate legs appear. For \(S_3\) that costs us nothing: §8 showed \(X=M\) and §3 showed the characters are real, so \(\bar\chi=\chi\) and the counting is unchanged. (For \(\mathbb{Z}_N\) it would not be.)
Before running the next cell. The library’s enumerator will report a number. Before you see it, ask what it should be: how many quadratic terms can you form (\(H^\dagger H\) is degree 2 in fields), and roughly how many quartic? Then, more usefully, ask what the symmetry has forbidden — a general 3HDM potential has far more parameters than this.
# ---- MOVE 1: enumerate the S3-invariant potential ------------------------
terms = suggest_potential([H1, H2, HS], [SU2L, U1Y], discrete_groups=[s3])
print(f"suggest_potential found {len(terms)} independent operators")
fieldsyms = [c for f in (H1, H2, HS) for c in list(f.components)]
allsyms = fieldsyms + [sp.conjugate(c) for c in fieldsyms]
deg = lambda e: sp.Poly(sp.expand(e), *allsyms).total_degree()
quad = [t for t in terms if deg(t.expr) == 2]
quart = [t for t in terms if deg(t.expr) == 4]
print(" quadratic (mass) terms :", len(quad))
print(" quartic (coupling) :", len(quart))
assert (len(terms), len(quad), len(quart)) == (10, 2, 8)
ok("2 mass + 8 quartic = 10 operators -- the S3 3HDM basis")
suggest_potential found 10 independent operators
quadratic (mass) terms : 2
quartic (coupling) : 8
✓ 2 mass + 8 quartic = 10 operators -- the S3 3HDM basis
Two mass terms and eight quartics. Compare that with an unconstrained three-Higgs-doublet potential, where every \(H_a^\dagger H_b\) pairing is allowed independently — that is where the symmetry earns its keep.
# ---- MOVE 2 + 3: what did the symmetry forbid? --------------------------
free = suggest_potential([H1, H2, HS], [SU2L, U1Y]) # no discrete group
print(f"without S3 : {len(free)} operators")
print(f"with S3 : {len(terms)} operators")
print(f"forbidden : {len(free) - len(terms)}")
ok(f"imposing S3 removes {len(free) - len(terms)} of the {len(free)} parameters "
"-- that is the entire point of a flavour symmetry")
without S3 : 33 operators
with S3 : 10 operators
forbidden : 23
✓ imposing S3 removes 23 of the 33 parameters -- that is the entire point of a flavour symmetry
# ---- MOVE 4: every surviving term really is invariant --------------------
for t in terms:
inv, residual = check_discrete_invariance(sp.expand(t.expr), s3)
assert inv, (t.expr, residual)
ok(f"all {len(terms)} enumerated operators pass check_discrete_invariance "
"independently of how they were produced")
# and a deliberately S3-breaking term must be rejected
bad_term = sp.expand((dag(H1) * H1.mat)[0]**2) # H1 alone, not the doublet
assert not check_discrete_invariance(bad_term, s3)[0]
ok("(H1^dag H1)^2 alone is correctly REJECTED -- the doublet must appear "
"symmetrically, so the check has teeth")
✓ all 10 enumerated operators pass check_discrete_invariance independently of how they were produced
✓ (H1^dag H1)^2 alone is correctly REJECTED -- the doublet must appear symmetrically, so the check has teeth
The 10-operator basis is what THDM_S3_Tutorial and 3HDM-\(S_3\) model building in general
build on: 2 mass parameters \(\mu^2\) and 8 quartics \(\lambda_{1..8}\). You now know where the
number 8 comes from, and — more importantly — that you could have predicted it.
9.1 What suggest_potential actually does¶
That was a black box, and it should not stay one — because it is not a new algorithm. It is the tools of §5, §6 and §8, run in sequence:
# |
stage |
the section it comes from |
|---|---|---|
1 |
build the gauge blocks \(H_a^\dagger H_b\) |
the gauge side — §9’s opening paragraph |
2 |
keep only the charge-neutral monomials |
§6: the \(\mathbb{Z}_N\) average is this rule |
3 |
|
§6 — literally the function demoed there |
4 |
hermitian completion, \(\mathcal O\to\mathcal O+\mathcal O^\dagger\) |
§8, the conjugate leg |
5 |
drop the linearly dependent ones |
§5: construct, then check the rank |
6 |
assert every survivor really is invariant |
§8, and MOVE 4 above |
Nothing on that list is new, so rather than take the “10” on trust, rebuild it. The next two
cells run stages 1–5 with nothing but dag, reynolds_project and §5’s own rank_of, and
finish by comparing against suggest_potential’s own answer.
Before running the next cells. Commit to three numbers first. (i) Three doublets, and a block is \(H_a^\dagger H_b\) — how many blocks? (ii) A potential term is one block (dimension 2) or two (dimension 4), and \(H_a^\dagger H_b\,H_c^\dagger H_d\) is the same operator whichever order you write the factors in — how many monomials? (iii) The one worth pausing on: is the unconstrained potential’s parameter count that monomial number, or something smaller? Look at \(H_1^\dagger H_2\) and ask whether it can appear in a real Lagrangian on its own.
# ---- MOVE 1: stages 1 and 2. Build the blocks, then apply the charge rule.
Hs = [H1, H2, HS]
blocks, dropped = {}, []
for A, B in itertools.product(Hs, repeat=2):
# the TWO gauge-invariant ways to pair two SU(2) doublets, with the
# hypercharge each one carries
pairings = [
(f"({A.name}+{B.name})", # A^dag B
sp.expand((dag(A) * B.mat)[0]),
B.charge(U1Y) - A.charge(U1Y)),
(f"({A.name}.{B.name})", # A^T (i sigma_2) B
sp.expand(A.components[0] * B.components[1]
- A.components[1] * B.components[0]),
A.charge(U1Y) + B.charge(U1Y)),
]
for label, expr, Y in pairings:
if Y == 0:
blocks[label] = expr
else:
dropped.append((label, Y))
print("blocks built :", len(blocks) + len(dropped))
print(" survive the charge rule :", len(blocks), " ", sorted(blocks))
print(" killed by hypercharge :", len(dropped),
f" e.g. {dropped[0][0]} carries Y = {dropped[0][1]}")
# a term is one block (dim 2) or two (dim 4); order does not make a new operator
mono = []
for n in (1, 2):
for combo in itertools.combinations_with_replacement(sorted(blocks), n):
mono.append((sp.expand(sp.Mul(*[blocks[k] for k in combo])),
" ".join(combo)))
print("\ncharge-neutral monomials of dimension 2 and 4 :", len(mono))
assert (len(blocks), len(dropped), len(mono)) == (9, 9, 54)
ok("half the blocks die on hypercharge before S_3 is ever consulted -- 9 of 18 "
"survive, and they generate 54 monomials")
blocks built : 18
survive the charge rule : 9 ['(H1+H1)', '(H1+H2)', '(H1+HS)', '(H2+H1)', '(H2+H2)', '(H2+HS)', '(HS+H1)', '(HS+H2)', '(HS+HS)']
killed by hypercharge : 9 e.g. (H1.H1) carries Y = 1
charge-neutral monomials of dimension 2 and 4 : 54
✓ half the blocks die on hypercharge before S_3 is ever consulted -- 9 of 18 survive, and they generate 54 monomials
Open the funnel’s input. blocks is a dictionary from a label to the expanded
component expression it stands for, and mono is the list of candidate operators, each
paired with a label built by joining its blocks’ names. Three of the fifty-four are enough
to see the shape.
blocks, mono[:3]
({'(H1+H1)': H1_1*conjugate(H1_1) + H1_2*conjugate(H1_2),
'(H1+H2)': H2_1*conjugate(H1_1) + H2_2*conjugate(H1_2),
'(H1+HS)': HS_1*conjugate(H1_1) + HS_2*conjugate(H1_2),
'(H2+H1)': H1_1*conjugate(H2_1) + H1_2*conjugate(H2_2),
'(H2+H2)': H2_1*conjugate(H2_1) + H2_2*conjugate(H2_2),
'(H2+HS)': HS_1*conjugate(H2_1) + HS_2*conjugate(H2_2),
'(HS+H1)': H1_1*conjugate(HS_1) + H1_2*conjugate(HS_2),
'(HS+H2)': H2_1*conjugate(HS_1) + H2_2*conjugate(HS_2),
'(HS+HS)': HS_1*conjugate(HS_1) + HS_2*conjugate(HS_2)},
[(H1_1*conjugate(H1_1) + H1_2*conjugate(H1_2), '(H1+H1)'),
(H2_1*conjugate(H1_1) + H2_2*conjugate(H1_2), '(H1+H2)'),
(HS_1*conjugate(H1_1) + HS_2*conjugate(H1_2), '(H1+HS)')])
Those labels are not decoration. They are what lets the next cells report which monomials the projection killed and which got paired with a conjugate — a rank alone would say \(54\to33\to10\) and leave you no way to check the claim against the operators you would have written by hand.
# ---- MOVE 2 + 3 + 4: project, complete, take the rank -- then check. -----
def independent(exprs):
# Section 5's rank_of is the right tool, but it cannot be pointed at these
# directly: they carry conjugate(...) nodes, and sympy's Poly will not
# accept those as generators. Map each conjugate to a Dummy -- SHARED
# across every candidate, so the same conjugate becomes the same Dummy
# everywhere and the coefficient rows stay comparable. (This is what
# suggest._normalize_atoms does; vacuum/masses.py uses the same trick to
# build charged mass matrices.) Section 5 needed none of it: real doublets.
reg = {}
def norm(e):
e = sp.expand(e)
for node in e.atoms(sp.conjugate):
reg.setdefault(node, sp.Dummy())
return e.xreplace(reg)
rows = [norm(e) for e in exprs]
gens = sorted({s for e in rows for s in e.free_symbols},
key=sp.default_sort_key)
return rank_of(rows, gens) # section 5's helper, unchanged
def funnel(groups, verbose=False, sample=4):
# stage 3: project onto the invariant subspace (a no-op with no group)
projected = [(sp.expand(reynolds_project(e, groups)) if groups else e, l)
for e, l in mono]
kept = [(e, l) for e, l in projected if e != 0]
killed = [l for e, l in projected if e == 0]
# stage 4: hermitian completion. O + O^dag, unless O is already hermitian
completed, hc, paired = [], 0, []
for e, l in kept:
if sp.expand(e - sp.conjugate(e)) == 0:
completed.append(e)
else:
completed.append(sp.expand(e + sp.conjugate(e)))
hc += 1
paired.append(l)
basis = independent(completed) # stage 5
if verbose:
trace(f" stage 3 (project) : {len(kept)} of {len(mono)} monomials "
"survive")
if killed:
trace(" killed, e.g. : " + ", ".join(killed[:sample]))
trace(f" stage 4 (+ h.c.) : {hc} of {len(kept)} are not "
f"self-conjugate and get paired with their conjugate")
if paired:
trace(" paired, e.g. : " + ", ".join(paired[:sample]))
trace(f" stage 5 (rank) : {basis} independent operators")
return len(kept), hc, basis
print(f"{'':<20} {'blocks':>8} {'monomials':>10} {'after P':>9} "
f"{'+ h.c.':>8} {'basis':>7}")
result = {}
for name, groups in (("with S_3", [s3]), ("no discrete group", [])):
after_P, hc, basis = funnel(groups)
result[name] = basis
print(f" {name:<18} {len(blocks):8d} {len(mono):10d} {after_P:9d} "
f"{hc:8d} {basis:7d}")
assert result["with S_3"] == len(terms) == 10
assert result["no discrete group"] == len(free) == 33
ok("the hand-built funnel reproduces suggest_potential exactly -- 10 with S_3 "
"and 33 without -- so it really is sections 5, 6 and 8 composed")
blocks monomials after P + h.c. basis
with S_3 9 54 26 12 10
no discrete group 9 54 54 42 33
✓ the hand-built funnel reproduces suggest_potential exactly -- 10 with S_3 and 33 without -- so it really is sections 5, 6 and 8 composed
# ---- peek: what the funnel throws away, with no discrete group at all ----
# This is the branch the paragraph below is about: hermiticity alone, before
# S_3 is ever consulted, is what takes 54 monomials down to 33 parameters.
_ = funnel([], verbose=True)
stage 3 (project) : 54 of 54 monomials survive
stage 4 (+ h.c.) : 42 of 54 are not self-conjugate and get paired with their conjugate
paired, e.g. : (H1+H2), (H1+HS), (H2+H1), (H2+HS)
stage 5 (rank) : 33 independent operators
Both numbers land, from tools that were on the table by §6. Two things the funnel shows that the bare “10” could not.
Hermiticity is doing more work than the symmetry gets credit for. Look at the + h.c.
column: without a discrete group, \(42\) of the \(54\) monomials are not self-conjugate and get
paired with their conjugates at stage 4. That pairing is what takes \(54\to33\) — and it
happens before any flavour symmetry is imposed. \(H_1^\dagger H_2\) cannot stand alone in a
real Lagrangian; it needs \(H_2^\dagger H_1\) beside it, and the two together are one
parameter, not two. So the reduction is really two steps: hermiticity \(54\to33\), then
\(S_3\) \(33\to10\). The cell in §9 reports only the second.
Stage 2 is not a formality. MOVE 1 built \(18\) blocks, not \(9\) — the other nine are the \(SU(2)\) \(\varepsilon\)-contractions \(A^{\mathsf T}(i\sigma_2)B\), which are every bit as gauge-invariant. Every one carries \(Y_A+Y_B=1\neq0\) for three \(Y=\tfrac12\) doublets, so the charge rule deletes all nine before \(S_3\) is consulted. That is §6’s abelian selection rule removing half the candidates on its own. (In a model with a \(Y=-\tfrac12\) doublet they would survive, and the potential would be larger.)
Stage 6 is also why §9’s MOVE 4 was worth running. The enumerator already applies the
invariance oracle to its own output — so re-checking each term with
check_discrete_invariance is not repetition, it is confirming the answer does not depend
on the route that produced it. If the projector of stage 3 and the checker’s generator
substitution ever disagreed, that cell is where it would surface.
9.2 The same count, for a Yukawa term¶
§9 counted the potential, where every factor was an \(H_a^\dagger H_b\) block built from scalars in the same pair of irreps. A Yukawa term is the more general problem, and the one flavour model building actually runs into:
three different multiplets, each free to sit in \(\mathbf1\) or \(\mathbf2\) independently. Nothing so far has counted that — §4 did \(n\) copies of one irrep, §9 did the potential.
You already have the tool, though. §4’s n_invariants takes a list of characters, one
per factor; it was only ever handed \(n\) copies of the same one. Hand it three different
ones and it answers the Yukawa question unchanged:
Before running the next cell. There are \(2^3=8\) ways to assign \((\bar L,H,e_R)\) to \(\{\mathbf2,\mathbf1\}\). Predict how many carry an invariant. Two hints: \(\mathbf1\otimes\mathbf1\otimes\mathbf1\) obviously does, and any assignment with exactly one \(\mathbf2\) cannot — a lone doublet has no invariant component, which is \(\langle\chi_2,\chi_1\rangle=0\) from §3. How many does that leave, and are they all non-zero?
Nothing about the fermions enters the count. \(\bar L\) contributes a conjugate leg, but §3 showed \(S_3\)’s characters are real, so \(\bar\chi=\chi\) and the arithmetic is the same one §4 ran — exactly the observation §9 opened with.
# ---- the Yukawa problem: THREE multiplets, in MIXED irreps ---------------
# n_invariants already accepts a list of characters -- section 4 simply never
# called it with different ones. No new machinery.
print(f"{'(Lbar, H, eR)':<18} # invariants")
yuk = {}
for combo in itertools.product(("2", "1"), repeat=3):
yuk[combo] = n_invariants([chi[r] for r in combo])
c1, c2, c3 = combo
print(f" ({c1:<2} {c2:<2} {c3:<2}){'':>6} {yuk[combo]}")
total = sum(yuk.values())
print("\n total S3-invariant Yukawa structures :", total)
assert total == 5
assert yuk[("2", "2", "2")] == 1 # the multiplicity 3HDM-S3 models rely on
assert all(v == 0 for k, v in yuk.items() if list(k).count("2") == 1)
# and the five are a model-building CHOICE: allowing 1' as well opens up more
with_1p = sum(n_invariants([chi[r] for r in c])
for c in itertools.product(("1", "1p", "2"), repeat=3))
print(" the same count if 1' is also allowed :", with_1p, "(over 27 assignments)")
ok("five Yukawa structures from 1 and 2 alone -- and the (2,2,2) row is the "
"(1/6)[1*8 + 3*0 + 2*(-1)] = 1 that a 3HDM-S3 lepton sector relies on")
(Lbar, H, eR) # invariants
(2 2 2 ) 1
(2 2 1 ) 1
(2 1 2 ) 1
(2 1 1 ) 0
(1 2 2 ) 1
(1 2 1 ) 0
(1 1 2 ) 0
(1 1 1 ) 1
total S3-invariant Yukawa structures : 5
the same count if 1' is also allowed : 11 (over 27 assignments)
✓ five Yukawa structures from 1 and 2 alone -- and the (2,2,2) row is the (1/6)[1*8 + 3*0 + 2*(-1)] = 1 that a 3HDM-S3 lepton sector relies on
# ---- peek: the (2,2,2) row, the one a 3HDM-S3 lepton sector relies on
n_invariants([chi["2"]] * 3, verbose=True, classes=classes)
class rep size prod_i chi_i(g)
() 1 8
(0,) 2 -1
(1,) 3 0
(1/6)[1*8 + 2*-1 + 3*0] = 1
Five, and each is a term someone has to write down. Before the table, the notation — it packs three different contractions into two symbols, and it is worth unpacking once.
Notation for this table. Everything here acts on the \(S_3\) index only. The \(SU(2)\) contraction and the spinor sandwich are the same in every row, so they are suppressed.
Fields. \(\bar L=(\bar L_1,\bar L_2)\), \(H=(H_1,H_2)\) and \(e_R=(e_{1R},e_{2R})\) are \(S_3\) doublets; \(\bar L_S\), \(H_S\) and \(e_{SR}\) are the singlets (the \(S\) is the one in
HSin the code).One operation, two shorthands — both from §5:
\(\{a\otimes b\}_{\mathbf r}\) is channel \(\mathbf r\) of \(\mathbf2\otimes\mathbf2\), i.e.
S3.doublet_product(a, b)["r"]. The \(\mathbf2\) channel is itself a doublet, \(\{a\otimes b\}_{\mathbf2}=\big(a_1b_1-a_2b_2,\;-(a_1b_2+a_2b_1)\big)\).\(a\cdot b\equiv\{a\otimes b\}_{\mathbf1}=a_1b_1+a_2b_2\). The dot is just the singlet channel, so it never needs a subscript.
The literature often writes \([\,\cdot\,]_{\mathbf r}\) for the same braces, and \((\bar L\cdot H)_{\mathbf1}\) for what is here simply \(\bar L\cdot H\).
\((\bar L, H, e_R)\) |
the invariant |
written out |
|---|---|---|
\((\mathbf2,\mathbf2,\mathbf2)\) |
\(\{\bar L\otimes H\}_{\mathbf2}\cdot e_R\) |
\((\bar L_1H_1-\bar L_2H_2)\,e_{1R}-(\bar L_1H_2+\bar L_2H_1)\,e_{2R}\) |
\((\mathbf2,\mathbf2,\mathbf1)\) |
\((\bar L\cdot H)\,e_{SR}\) |
\((\bar L_1H_1+\bar L_2H_2)\,e_{SR}\) |
\((\mathbf2,\mathbf1,\mathbf2)\) |
\((\bar L\cdot e_R)\,H_S\) |
\(\bar L_1H_Se_{1R}+\bar L_2H_Se_{2R}\) |
\((\mathbf1,\mathbf2,\mathbf2)\) |
\(\bar L_S\,(H\cdot e_R)\) |
\(\bar L_S\,(H_1e_{1R}+H_2e_{2R})\) |
\((\mathbf1,\mathbf1,\mathbf1)\) |
\(\bar L_S H_S e_{SR}\) |
\(\bar L_S H_S e_{SR}\) |
Read the middle column by counting doublets. Two doublets and a singlet: dot the doublets, and the singlet rides along. Three doublets: there is nothing to dot all three with, so fuse two into a new doublet with \(\{\,\otimes\,\}_{\mathbf2}\), then dot that with the third. Only the first row takes two steps — and it raises the obvious question of whether it matters which two you fuse. The next cell answers that, and one more thing.
# ---- peek: the (2,2,2) invariant, taken apart ------------------------------
Lb_, H_, eR_ = sp.symbols("Lb1 Lb2"), sp.symbols("H1 H2"), sp.symbols("e1 e2")
dot = lambda a, b: cg(a, b)["1"]
fuse2 = lambda a, b: cg(a, b)["2"]
# (1) the dot IS the singlet channel
assert sp.expand(dot(Lb_, H_) - (Lb_[0] * H_[0] + Lb_[1] * H_[1])) == 0
# (2) fuse any two first -- the invariant does not care
orders = {"{Lb x H}_2 . eR": dot(fuse2(Lb_, H_), eR_),
"Lb . {H x eR}_2": dot(Lb_, fuse2(H_, eR_)),
"{Lb x eR}_2 . H": dot(fuse2(Lb_, eR_), H_)}
T2 = sp.expand(orders["{Lb x H}_2 . eR"])
for label, expr in orders.items():
print(f" {label:<18}:", sp.expand(expr))
assert sp.expand(expr - T2) == 0
# (3) ...because its coefficient tensor C_ijk is fully symmetric
poly = sp.Poly(T2, *Lb_, *H_, *eR_)
Cijk = {(i + 1, j + 1, k + 1): poly.coeff_monomial(Lb_[i] * H_[j] * eR_[k])
for i, j, k in itertools.product(range(2), repeat=3)}
print("\n nonzero C_ijk :", {ijk: c for ijk, c in Cijk.items() if c != 0})
assert all(Cijk[p] == c for ijk, c in Cijk.items()
for p in itertools.permutations(ijk))
# (4) all three legs on ONE doublet: the cubic invariant section 6 projected
cubic = sp.expand(dot(fuse2((r1_, r2_), (r1_, r2_)), (r1_, r2_)))
display(sp.Eq(sp.Symbol(r"\{X\otimes X\}_{2}\cdot X"), cubic, evaluate=False))
assert sp.expand(cubic - 4 * reynolds_project(r1_**3, [s3b])) == 0
# (5) the dot survives ANY rotation; the cubic only S3's 120-degree ones
def rotated_change(expr, theta):
c, s = sp.cos(theta), sp.sin(theta)
return sp.simplify(sp.expand(
expr.xreplace({r1_: c * r1_ - s * r2_, r2_: s * r1_ + c * r2_}) - expr))
square = sp.expand(dot((r1_, r2_), (r1_, r2_)))
for theta in (2 * sp.pi / 3, sp.pi / 2):
print(f" rotate by {str(theta):<8}: dot changes by {rotated_change(square, theta)},"
f" cubic changes by {rotated_change(cubic, theta)}")
assert rotated_change(square, sp.pi / 2) == 0
assert rotated_change(cubic, 2 * sp.pi / 3) == 0
assert rotated_change(cubic, sp.pi / 2) != 0
ok("the (2,2,2) invariant is one symmetric tensor, fusion order is irrelevant, and "
"it is the cubic Re(z^3) -- invariant under S3's rotations, not under all of them")
{Lb x H}_2 . eR : H1*Lb1*e1 - H1*Lb2*e2 - H2*Lb1*e2 - H2*Lb2*e1
Lb . {H x eR}_2 : H1*Lb1*e1 - H1*Lb2*e2 - H2*Lb1*e2 - H2*Lb2*e1
{Lb x eR}_2 . H : H1*Lb1*e1 - H1*Lb2*e2 - H2*Lb1*e2 - H2*Lb2*e1
nonzero C_ijk : {(1, 1, 1): 1, (1, 2, 2): -1, (2, 1, 2): -1, (2, 2, 1): -1}
rotate by 2*pi/3 : dot changes by 0, cubic changes by 0
rotate by pi/2 : dot changes by 0, cubic changes by -X1**3 + 3*X1**2*X2 + 3*X1*X2**2 - X2**3
✓ the (2,2,2) invariant is one symmetric tensor, fusion order is irrelevant, and it is the cubic Re(z^3) -- invariant under S3's rotations, not under all of them
Three things to read off, each a callback:
The all-doublet row needs §5, but not a choice. The invariant is a single, fully symmetric 3-index tensor, \(C_{111}=1\) and \(C_{122}=C_{212}=C_{221}=-1\). Any pairing of the legs reads off the same tensor, so “fuse two, keep the \(\mathbf2\), dot with the third” is a recipe for finding it, not a decision that changes it. The count above is the other half of the guarantee: multiplicity one means there is nothing else to find.
Only that row is basis-sensitive (§7). Write a doublet as \(z=x_1+ix_2\). The dot is \(\lvert z\rvert^2\), untouched by every rotation, so the four dot rows look the same in any orthogonal doublet basis. The triple structure, with all legs set equal, is \(x_1^3-3x_1x_2^2=\mathrm{Re}\,z^3\) — the very cubic §6 projected — and a rotation by \(\theta\) multiplies \(z^3\) by \(e^{3i\theta}\). It is invariant only for \(S_3\)’s own \(120^\circ\) rotations, so its coefficients depend on how the doublet basis sits relative to them: the fermionic analogue of \(\lambda_4\) in the potential. Counts transfer between papers; these coefficients do not.
\(\mathbf1'\) was a choice. Allowing it too gives 11 invariants over 27 assignments. A \(\mathbf1'\) field would also need the third channel, \(\{a\otimes b\}_{\mathbf1'}=a_1b_2-a_2b_1\), which no dot covers. Assigning every field to \(\mathbf1\) or \(\mathbf2\) is a model-building decision, and the number five is downstream of it.
This is where a real model starts. A 3HDM-\(S_3\) lepton sector writes one Yukawa coupling
per structure — five, not the \(3\times3\times3\) a flavour-blind model would allow. Check
each with check_discrete_invariance, and confirm the count independently with
suggest_yukawa, the fermionic sibling of the pipeline §9.1 walked through. And when
importing the all-doublet structure from a paper, transcribe it through its basis: written
in the other doublet basis, it fails the invariance check here — the second bullet
above, biting in practice.
10. Capstone: \(A_4\) from scratch¶
\(S_3\) is the warm-up. The workhorse of flavour model building is \(A_4\) — the even permutations of four objects, the rotation group of the tetrahedron, order 12 — because it has a 3-dimensional irrep that can hold three generations at once.
feynlag ships only \(\mathbb{Z}_N\) and \(S_3\). But DiscreteSymmetry advertises itself as
extensible, and the contract is small: supply _irrep_generators, one matrix per
generator per irrep. Let us hold it to that.
\(A_4\) is presented by \(S^2=T^3=(ST)^3=e\), with irreps \(\mathbf 1,\mathbf 1',\mathbf 1'', \mathbf 3\) (note \(1+1+1+9=12\) ✓). In the standard basis
and the three singlets are distinguished by \(T\to 1,\omega,\omega^2\).
Before running the next cells. \(A_4\) has four irreps, so four conjugacy classes. Predict \(\chi_3\): \(\chi_3(e)=3\) certainly. \(S\) is a double transposition — what is \(\mathrm{Tr}\,S\)? \(T\) is a 3-cycle, and a cyclic permutation matrix has trace equal to the number of fixed points — how many does it have? Fill in the row before you look.
# ---- MOVE 1: define A4 on top of the shipped base class -----------------
omega = sp.exp(2 * sp.pi * sp.I / 3) # `w` is taken: section 5's fourth doublet
S_ = sp.Matrix([[1, 0, 0], [0, -1, 0], [0, 0, -1]])
T_ = sp.Matrix([[0, 1, 0], [0, 0, 1], [1, 0, 0]])
class A4(DiscreteSymmetry):
# The alternating group A_4: irreps 1, 1', 1'' and 3.
def __init__(self, name="A4"):
super().__init__(name)
self._irrep_generators = {
"1": [sp.Matrix([[1]]), sp.Matrix([[1]])],
"1p": [sp.Matrix([[1]]), sp.Matrix([[omega]])],
"1pp": [sp.Matrix([[1]]), sp.Matrix([[omega**2]])],
"3": [S_, T_],
}
A4_GENS = A4()._irrep_generators
for name, M, n in (("S^2", S_**2, 2), ("T^3", T_**3, 3), ("(ST)^3", (S_ * T_)**3, 3)):
assert sp.simplify(M) == sp.eye(3), name
print(f" {name:<7} = identity ✓")
ok("S^2 = T^3 = (ST)^3 = e -- the presentation of A_4")
S^2 = identity ✓
T^3 = identity ✓
(ST)^3 = identity ✓
✓ S^2 = T^3 = (ST)^3 = e -- the presentation of A_4
# ---- MOVE 2: close it, and derive the character table -------------------
W4 = close_group(A4_GENS["3"])
G4 = len(W4)
print("|A_4| =", G4)
assert G4 == 12
assert sum(A4_GENS[l][0].shape[0]**2 for l in A4_GENS) == G4
ok("order 12, and sum (dim r)^2 = 1+1+1+9 = 12 -- the irrep list is complete")
chi4 = {lab: [sp.nsimplify(sp.expand(sp.trace(rho(lab, wd, A4_GENS))))
for wd in W4] for lab in A4_GENS}
classes4 = conjugacy_classes(W4, A4_GENS, "3")
print("\nconjugacy classes:", [len(cl) for cl in classes4], " (sizes)")
names4 = ["e", "S (double transp.)", "T (3-cycle)", "T^2 (3-cycle)"]
hdr4 = f"{'':<6}" + "".join(f"{n:<24}" for n in names4)
print(); print(hdr4); print("-" * len(hdr4))
for lab in ("1", "1p", "1pp", "3"):
row = f"{lab:<6}"
for cl in classes4:
val = sp.simplify(sp.expand_complex(chi4[lab][W4.index(cl[0])]))
row += f"{val!s:<24}"
print(row)
for lab in A4_GENS: # constant on classes, as before
for cl in classes4:
vals = {sp.simplify(chi4[lab][W4.index(w)] - chi4[lab][W4.index(cl[0])])
for w in cl}
assert vals == {0}, (lab, cl)
ok("A_4's characters are constant on its 4 conjugacy classes")
|A_4| = 12
✓ order 12, and sum (dim r)^2 = 1+1+1+9 = 12 -- the irrep list is complete
conjugacy classes: [1, 3, 4, 4] (sizes)
e S (double transp.) T (3-cycle) T^2 (3-cycle)
------------------------------------------------------------------------------------------------------
1 1 1 1 1
1p 1 1 -1/2 + sqrt(3)*I/2 -1/2 - sqrt(3)*I/2
1pp 1 1 -1/2 - sqrt(3)*I/2 -1/2 + sqrt(3)*I/2
3 3 -1 0 0
✓ A_4's characters are constant on its 4 conjugacy classes
Why the shipped helpers accepted a group they had never seen. Put \(A_4\)’s bookkeeping next to \(S_3\)’s from §3 — the conjugacy classes, and the characters of the three-dimensional irrep.
classes4, chi4["3"]
Same two shapes as before: a list of lists of words, and one character per element,
indexed by the closure order W4 rather than by class. Nothing in close_group, rho
or conjugacy_classes ever asked which group it was holding, so supplying
_irrep_generators really was the whole contract.
# ---- MOVE 4: orthonormality, on a group the library never saw -----------
inner4 = lambda p, q: sp.simplify(
sum(sp.conjugate(chi4[p][i]) * chi4[q][i] for i in range(G4)) / G4)
for p in A4_GENS:
for q in A4_GENS:
assert inner4(p, q) == (1 if p == q else 0), (p, q, inner4(p, q))
ok("<chi_r, chi_s> = delta_rs for all 16 pairs -- A_4's character table is "
"consistent, and nothing about it was typed in")
print()
for n in (2, 3, 4):
n_inv = sp.simplify(sum(chi4["3"][i]**n for i in range(G4)) / G4)
print(f" invariants in 3^(x{n}) (distinct triplets) : {n_inv}")
✓ <chi_r, chi_s> = delta_rs for all 16 pairs -- A_4's character table is consistent, and nothing about it was typed in
invariants in 3^(x2) (distinct triplets) : 1
invariants in 3^(x3) (distinct triplets) : 2
invariants in 3^(x4) (distinct triplets) : 7
\(\chi_3=(3,-1,0,0)\): the double transposition \(S\) has trace \(-1\), and the 3-cycles are traceless because a 3-cycle permutation matrix has no fixed points.
The counts are the reason \(A_4\) is the flavour group of choice: one invariant in \(\mathbf 3\otimes\mathbf 3\) (so a mass term exists), and two in \(\mathbf 3\otimes\mathbf 3\otimes\mathbf 3\) — the symmetric and antisymmetric cubics, which is exactly the freedom a Yukawa sector needs to generate non-trivial mixing.
Finally: do the shipped tools work on this group they have never met?
# ---- the library's machinery, unchanged, on A4 --------------------------
a4 = A4()
F1, F2, F3 = (Scalar(n, real=True) for n in ("F1", "F2", "F3"))
a4.assign("3", F1, F2, F3)
f1, f2, f3 = (f.components[0] for f in (F1, F2, F3))
inv, _ = check_discrete_invariance(f1**2 + f2**2 + f3**2, a4)
notinv, residual = check_discrete_invariance(f1**2, a4)
print(" check_discrete_invariance(F1^2 + F2^2 + F3^2) :", inv)
print(" check_discrete_invariance(F1^2) :", notinv,
" residual:", end=" ")
display(sp.simplify(residual[0][1]))
assert inv and not notinv
display(sp.Eq(sp.Function("P")(f1**2),
sp.simplify(reynolds_project(f1**2, [a4])), evaluate=False))
display(sp.Eq(sp.Function("P")(f1 * f2),
sp.simplify(reynolds_project(f1 * f2, [a4])), evaluate=False))
assert sp.simplify(reynolds_project(f1 * f2, [a4])) == 0
ok("check_discrete_invariance and reynolds_project both work on A_4 with no "
"library change -- the projector finds (F1^2+F2^2+F3^2)/3 and kills F1 F2")
check_discrete_invariance(F1^2 + F2^2 + F3^2) : True
check_discrete_invariance(F1^2) : False residual:
✓ check_discrete_invariance and reynolds_project both work on A_4 with no library change -- the projector finds (F1^2+F2^2+F3^2)/3 and kills F1 F2
# ---- and the count agrees with explicit construction --------------------
quad_inv = sp.expand(reynolds_project(f1**2, [a4]))
n_quad = sp.simplify(sum(chi4["3"][i]**2 for i in range(G4)) / G4)
basis = [sp.expand(reynolds_project(m, [a4]))
for m in (f1**2, f2**2, f3**2, f1*f2, f1*f3, f2*f3)]
nz = [b for b in basis if b != 0]
r = sp.Matrix([[sp.Poly(b, f1, f2, f3).coeff_monomial(m)
for m in sorted({m for e in nz
for m in sp.Poly(e, f1, f2, f3).monoms()})]
for b in nz]).rank()
print(f" projecting all 6 quadratic monomials: {len(nz)} survive, rank {r}")
print(" characters predicted:", n_quad)
assert r == n_quad == 1
ok("one quadratic invariant, predicted and constructed -- the same loop as "
"section 5, on a brand-new group")
projecting all 6 quadratic monomials: 3 survive, rank 1
characters predicted: 1
✓ one quadratic invariant, predicted and constructed -- the same loop as section 5, on a brand-new group
11. Recap¶
The mechanics.
A finite group is its generators; closing them under multiplication regenerates it. Checking invariance under generators is checking it under the group.
A representation is any \(\rho\) with \(\rho(g)\rho(h)=\rho(gh)\) — injectivity not required. The trivial rep is “invariant”; the sign rep \(\mathbf 1'\) is not.
Characters are traces. Basis-independent, constant on classes, orthonormal: \(\sum_r(\dim r)^2=|G|\) and #irreps = #classes tell you when your list is complete.
The tools, in the order you should reach for them.
question |
tool |
|---|---|
how many invariants, \(n\) distinct multiplets? |
\(\frac1{|G|}\sum_g\prod_i\chi_i(g)\) |
how many, from one multiplet? |
the Molien series |
what are they, explicitly? |
Clebsch–Gordan ( |
…without a CG table? |
the Reynolds projector ( |
is this term invariant? |
|
what is the whole basis? |
|
The four traps, each of which produced a wrong answer above before the check caught it:
\(\mathbf 1'\) vanishes on a repeated multiplet. One \(S_3\) doublet forms no antisymmetric singlet with itself.
Distinct legs versus one field. Four distinct doublets give 3 quartic invariants; one doublet gives 1. Use characters for the first, Molien for the second.
Basis conventions. CG coefficients differ between the real and complex \(S_3\) bases; counts do not. Compare counts across papers, never coefficients.
The conjugate leg. \(X=(M^{-1})^{\mathsf T}\), which equals \(M\) only for real orthogonal irreps. True for \(S_3\), false for \(\mathbb{Z}_N\) — so the \(S_3\) coincidence hides the error.
And the model-building punchline. Imposing \(S_3\) on three Higgs doublets cuts the
potential to 2 mass terms and 8 quartics; the forbidden operators are what the symmetry
buys you. You can compute that number from the character table before writing a single
term — and A_4, a group the library does not ship, needed no library change at all.
Where to go next: THDM_S3_Tutorial builds the 3HDM-\(S_3\) potential from these
invariants and breaks the symmetry; ModelBuilding_Tutorial does the same enumeration
for continuous \(U(1)\) charges with anomaly cancellation; SUN_Groups_Tutorial is the
continuous-group counterpart to this notebook.