From Symmetry to Theory: A Computational Engine

Introduction

In the dynamical similarity and allometry posts, we started with a group acting on a system and asked what survived the group action. For dynamical similarity, the survivors were quantities shared by whole families of similar trajectories, somewhat like the conserved quantities in Noether’s theorem. For allometry, they were the power laws relating an organism’s traits to its body size. Under the right conditions, we can turn a symmetry into either something the system conserves or a restriction on the forms it can take.

What if we run the calculation in the other direction? Suppose we specify the fields, their symmetries, and how complicated an answer we are willing to consider. Can we list every theory consistent with those choices? Can we then derive the equations, charges, or selection rules without starting over by hand each time?

This post builds an engine for doing that procedure. The larger goal is to treat a family of symmetry-allowed theories as something we can actually search. Given a finite jet representation, a specified linear or affine symmetry action, bounded polynomial and derivative orders, and explicitly declared equivalences, the engine computes the corresponding space of local polynomial expressions. Furthermore, if these derivations can be put under a single algorithm, the resulting catalog could function as a kind of “periodic table” for (some parts of) mathematical physics. Instead of searching the literature for each symmetry calculation, we could query an API. Finally, since we can machine-verify these, it’s a good system to work in conjunction with LLMs.

Setup

Let’s first make the question finite. We need three inputs:

1. a vector space \(V\) containing the fields or variables
2. a group \(G\) acting on \(V\) through a representation \(\rho: G \to \mathrm{GL}(V)\)
3. a truncation rule: a polynomial degree \(d\) and, when derivatives are included, a maximum derivative order \(k\)

The output is the finite-dimensional space of polynomial terms allowed at that truncation. Let’s look at how the computation works.

Invariants at Fixed Degree

Consider a group \(G\) acting on a finite-dimensional real vector space \(V\) via a representation \(\rho \colon G \to \mathrm{GL}(V)\). A polynomial \(f \in \mathbb{R}[V]\) is invariant if

\[ f(\rho(g)x) = f(x) \]

for all \(g \in G\) and \(x \in V\).

The invariants form a graded subalgebra

\[ \mathbb{R}[V]^G \subseteq \mathbb{R}[V] \]

where the degree-\(d\) component \(\mathbb{R}[V]^G_d\) is finite-dimensional. Fixing a basis \(I_1, \dots, I_m\) of \(\mathbb{R}[V]^G_d\), every invariant of degree \(d\) is a linear combination

\[ \sum_{i=1}^m c_i I_i \] \[ c_i \in \mathbb{R} \]

So at a fixed degree, finding the invariants is a finite linear-algebra problem.

Choose coordinates \(x_1,\dots,x_n\) on \(V\). The degree-\(d\) monomials

\[ x^\alpha = x_1^{\alpha_1}\cdots x_n^{\alpha_n} \] \[ |\alpha| = \alpha_1+\cdots+\alpha_n=d \]

form a basis of \(\mathbb{R}[V]_d\). Write them as a column vector

\[ M_d(x) = \begin{pmatrix} m_1(x)\\ \vdots\\ m_N(x) \end{pmatrix} \]

Every degree-\(d\) polynomial can then be written as

\[ f_c(x)=c^\top M_d(x) \]

for a coefficient vector \(c\in\mathbb{R}^N\).

A group element acts on the monomial basis by substitution. Define the matrix \(A_d(g)\) by

\[ M_d(\rho(g)x)=A_d(g)M_d(x) \]

Then

\[ f_c(\rho(g)x) = c^\top A_d(g)M_d(x) \]

The invariance condition \(f_c(\rho(g)x)=f_c(x)\) is therefore the linear equation

\[ \big(A_d(g)^\top-I\big)c=0 \]

Each symmetry element gives a block of equations for the coefficients. Stack the blocks and we get one system

\[ B_d c=0 \]

whose null space is the degree-\(d\) invariant space

\[ \mathbb{R}[V]^G_d=\ker B_d \]

This degree-by-degree construction assumes that the action is linear. Translations are affine: \(x\mapsto x+a\) mixes homogeneous degrees. The engine therefore handles translations separately through the infinitesimal derivations \(\partial/\partial x_i\). A direct finite affine action would instead require monomials through degree \(d\), rather than only those of degree \(d\), or an equivalent homogeneous-coordinate construction.

Finite Groups: Projection onto the Invariants

Here is the smallest useful example. Consider the group \(G = \{e, i\}\) of order two¹, where \(i\) is spatial inversion

\[ i \cdot (v_x, v_y, v_z) = (-v_x, -v_y, -v_z) \]

on \(\mathbb{R}^3\). Invariance means \(f(-\mathbf{v}) = f(\mathbf{v})\). Replacing \(\mathbf{v}\) by \(-\mathbf{v}\) multiplies a degree-\(d\) monomial by \((-1)^d\), so the even monomials survive and the odd ones change sign.

The Reynolds operator averages over all elements of the group

\[ R(f) = \frac{1}{|G|} \sum_{g\in G} g\cdot f = \tfrac12\big(e\cdot f+i\cdot f\big) = \tfrac12\big(f(\mathbf{v})+f(-\mathbf{v})\big) \]

This average is a general recipe for projecting a polynomial onto the invariants. On a few monomials

\[ R(v_x) = \tfrac12(v_x - v_x) = 0 \] \[ R(v_x^2) = \tfrac12(v_x^2 + v_x^2) = v_x^2 \] \[ R(v_x v_y) = \tfrac12(v_x v_y + v_x v_y) = v_x v_y \]

The odd monomial averages to zero, while the even ones survive. For example

\[ R(3 v_x + v_x^2 - 2 v_y v_z) = v_x^2 - 2 v_y v_z \]

The invariant ring is therefore the ring of even polynomials. It is generated by the six quadratics \(v_x^2,\ v_y^2,\ v_z^2,\ v_x v_y,\ v_x v_z,\ v_y v_z\), subject to relations among them.

Continuous Groups: Infinitesimal Constraints

A continuous group looks harder at first blush because it has infinitely many elements. However, for compact groups, we can just replace the average with an integral. Similarly, for a connected Lie group, we can differentiate the action at the identity and use its infinitesimal generators to obtain the invariants.

If \(X\) is an infinitesimal generator acting linearly on \(V\), its induced action on a polynomial is

\[ \widehat{X}f(x) = \nabla f(x)\cdot Xx \]

Invariance under the one-parameter group generated by \(X\) requires

\[ \widehat{X}f=0 \]

At a fixed degree, \(\widehat{X}\) is itself a matrix acting on the monomial basis. If \(X_1,\dots,X_r\) generate the Lie algebra, the invariant space is the common kernel

\[ \mathbb{R}[V]^G_d = \bigcap_{a=1}^r \ker \widehat{X}_a \]

Again, the calculation is just a stack of linear equations followed by computing a null space.

The rotation group \(\mathrm{SO}(3)\) imposes more conditions than inversion alone. On \(\mathbb{R}^3\), only the squared length remains

\[ q=v_x^2+v_y^2+v_z^2 \]

and

\[ \mathbb{R}[V]^{\mathrm{SO}(3)} = \mathbb{R}[q] \]

There is one invariant at each even degree and none at odd degree. Why did the answer shrink so much? The rotation generators add more rows to the constraint matrix, leaving a smaller common null space.

This is the same basic procedure I used in the computational invariant theory posts. There, the invariants separated orbits and classified a space of objects. Here I want to use them for something else: determining which mechanics a system can have.

Mechanics

An invariant polynomial is not yet mechanics. We still need a rule that picks a path through configuration space. Following the geometric controls post, I’ll assign a cost (the Lagrangian) to each possible path. The Lagrangian is the infinitesimal cost per unit time, and

\[ S = \int L\,dt \]

is the cost of the whole path. The system follows a stationary path of this action².

The Free Particle

How much of the free-particle Lagrangian can symmetry determine? Start with a particle moving through space and an unconstrained Lagrangian

\[ S = \int L(\mathbf{x}, \mathbf{v}, t)\, dt \qquad \mathbf{v} = \dot{\mathbf{x}} \]

with \(L\) so far unconstrained³.

Now impose the Galilean symmetries one at a time and see what is left (we should get kinetic energy).

Step 1: Spatial and Time Translations

Start with translations. A shift of the spatial origin

\[ \mathbf{x}\mapsto\mathbf{x}+\mathbf{a} \]

does not affect the velocity, so invariance demands

\[ L(\mathbf{x}+\mathbf{a},\mathbf{v},t) = L(\mathbf{x},\mathbf{v},t) \]

for all \(\mathbf{a}\).

Differentiating at \(\mathbf{a}=0\) gives

\[ \frac{\partial L}{\partial\mathbf{x}}=0 \]

Likewise, invariance under the time shift

\[ t\mapsto t+s \]

gives

\[ \frac{\partial L}{\partial t}=0 \]

The first constraint block removes \(\mathbf{x}\) and \(t\)

\[ L(\mathbf{x},\mathbf{v},t) \;\longrightarrow\; L(\mathbf{v}) \]

Step 2: Rotations

Next, impose rotations

\[ \mathbf{v}\mapsto R\mathbf{v} \]

An infinitesimal rotation about an axis \(\vec{\omega}\) sends

\[ \mathbf{v} \mapsto \mathbf{v}+\vec{\omega}\times\mathbf{v} \]

The determining equation is

\[ \frac{\partial L}{\partial\mathbf{v}} \cdot (\vec{\omega}\times\mathbf{v}) = 0 \qquad \forall\vec{\omega} \]

Using the scalar triple product, this becomes

\[ \vec{\omega}\cdot \left( \mathbf{v} \times \frac{\partial L}{\partial\mathbf{v}} \right) = 0 \qquad \text{for all }\vec{\omega} \]

so

\[ \mathbf{v} \times \frac{\partial L}{\partial\mathbf{v}} = 0 \]

The gradient of \(L\) must point along \(\mathbf{v}\). Therefore, \(L\) can depend on the velocity only through its squared length

\[ q=|\mathbf{v}|^2, \qquad L=\varphi(q) \]

The second constraint block gives

\[ L(\mathbf{v}) \;\longrightarrow\; \varphi(|\mathbf{v}|^2) \]

The costs left standing are the functions on the orbit space of the rotation group.

Step 3: Galilean Boosts

Finally, impose boosts

\[ \mathbf{v}\mapsto\mathbf{v}+\mathbf{u} \]

This condition is weaker, as a boost may change \(L\) by only a total time derivative, because that changes the action only at its endpoints

\[ \int_{t_i}^{t_f}\frac{dF}{dt}\,dt = F(t_f)-F(t_i) \]

The endpoints are fixed in the variational problem, so this term cannot move the stationary path.

For an infinitesimal boost, the determining equation is therefore

\[ \frac{\partial L}{\partial\mathbf{v}} \cdot\mathbf{u} = \frac{dF_{\mathbf{u}}}{dt} \]

for some boundary term \(F_{\mathbf{u}}\).

With

\[ L=\varphi(|\mathbf{v}|^2) \]

the left-hand side is

\[ 2\varphi'(|\mathbf{v}|^2)\, \mathbf{u}\cdot\mathbf{v} \]

Since

\[ \mathbf{u}\cdot\mathbf{v} = \frac{d}{dt} (\mathbf{u}\cdot\mathbf{x}) \]

this is a total derivative for arbitrary trajectories when \(\varphi'\) is constant. Thus

\[ \varphi(q)=c_0+c_1q \]

The constant \(c_0\) contributes only an additive constant to the action, so writing

\[ c_1=\tfrac12m \]

leaves

\[ L=\tfrac12m|\mathbf{v}|^2 \]

Which is the kinetic energy. The only physical coefficient left is the mass \(m\).

Free Particle as Linear Algebra

Let’s rerun the rotation step at fixed degree and to show how there is an associated finite linear-algebra problem (that the engine will solve).

Work at a fixed degree, say degree two⁴. Write the most general cost as a quadratic form in the velocity

\[ L = \mathbf{v}^{\top}C\mathbf{v} = \sum_{i,j}c_{ij}v_iv_j \]

where \(C\) is a symmetric \(3\times3\) matrix with six independent coefficients

\[ C = \begin{pmatrix} c_{xx} & c_{xy} & c_{xz}\\ c_{xy} & c_{yy} & c_{yz}\\ c_{xz} & c_{yz} & c_{zz} \end{pmatrix} \]

So the candidate space has six dimensions.

The rotation now acts on the coefficients. Under

\[ \mathbf{v}\mapsto R\mathbf{v} \]

the cost becomes

\[ L \mapsto \mathbf{v}^{\top} (R^{\top}CR) \mathbf{v} \]

Invariance is the linear condition

\[ R^{\top}CR=C \qquad \forall R\in\mathrm{SO}(3) \]

To turn this into a finite system, take the three infinitesimal rotation generators \(X_x,X_y,X_z\). Writing

\[ R(\theta)=I+\theta X+O(\theta^2) \]

gives

\[ R(\theta)^\top C R(\theta) = C+\theta(X^\top C+CX)+O(\theta^2) \]

The determining equation for each generator is therefore

\[ X^\top C+CX=0 \]

Each generator contributes a block of linear equations in

\[ (c_{xx},c_{yy},c_{zz},c_{xy},c_{xz},c_{yz}) \]

Stacking the three blocks gives one matrix system. Its solution kills every off-diagonal entry and forces the three diagonal entries to agree

\[ c_{xy}=c_{xz}=c_{yz}=0 \]

\[ c_{xx}=c_{yy}=c_{zz}\equiv\lambda \]

Thus

\[ C=\lambda I \]

and

\[ L = \lambda(v_x^2+v_y^2+v_z^2) = \lambda|\mathbf{v}|^2 \]

Stacking Multiple Symmetries

So far, every constraint came from one symmetry group, but a physical system usually has several.

Suppose we have a field space \(V\) acted upon by groups \(G_1, \dots, G_r\) through a representation \(\rho_i \colon G_i \to \mathrm{GL}(V)\). Identifying each group with its image, the \(G_i\) are subgroups of \(\mathrm{GL}(V)\).

Let \(G = \langle G_1, \dots, G_r \rangle \subseteq \mathrm{GL}(V)\) be the smallest subgroup containing them all. A theory respects every \(G_i\) exactly when it is invariant under \(G\). Thus

\[ \mathbb{R}[V]^{G} = \bigcap_{i=1}^{r} \mathbb{R}[V]^{G_i} \]

The invariant ring of \(G\) is the intersection of the individual invariant rings. Imposing the symmetries one at a time gives a descending chain

\[ \mathbb{R}[V]^{G_1} \supseteq \mathbb{R}[V]^{G_1} \cap \mathbb{R}[V]^{G_2} \supseteq \cdots \supseteq \mathbb{R}[V]^{G} \]

Each symmetry adds another block of equations and cuts the solution space. Because the endpoint is an intersection, the final answer does not depend on the order of the blocks.

General Mechanical Procedure

We can generalize our free particle example to obtain a general recipe.

Given the fields, the symmetry group, and a truncation rule:

1. List the field variables, together with any derivatives that are being retained.
2. List the monomials allowed by the degree and derivative bounds.
3. Compute the induced action of each symmetry generator on those monomials.
4. Convert invariance into linear equations on the monomial coefficients.
5. Stack the equations from all of the generators.
6. Compute the common null space.
7. Use a basis \(I_1,\dots,I_m\) of that null space to write the most general allowed expression \(F=\sum_{i=1}^m c_iI_i\)
8. Apply any additional quotients or quasi-symmetry conditions, such as equivalence up to a total derivative.
9. Derive the physical consequences of the surviving \(F\).

Across all degrees, the answer is a polynomial, or more generally a function, of the invariant generators, subject to the relations among them. At low order this gives a minimal theory. If we keep every local term order by order, we get the usual symmetry-compatible expansion. This is the Landau theory for a free energy (Landau 1937), or the logic of Weinberg’s Folk Theorem for a Lagrangian (Weinberg 1995).⁵

Symmetry decides which terms may appear. It normally does not determine their coefficients; those come from the material, a microscopic calculation, or a measurement.

Once we have

\[ F=\sum_i c_iI_i \]

the engine can generate the first layer of consequences symbolically: Euler-Lagrange equations, formal Noether currents, Hessians, and linearized operators. These are not the whole physics. Stability may still require spectral or nonlinear analysis, dispersion requires a background and treatment of zero modes, and gauge systems require constraint reduction. Coefficients still have to come from a microscopic calculation or measurement.

At a fixed truncation, finding the most general theory compatible with a declared set of linear or affine symmetries is therefore a finite linear-algebra problem.

Staged Reduction

What do we gain by imposing the groups one at a time?

General Procedure

The order cannot change the answer, as the allowed theory space is still the intersection of the individual invariant spaces, but it can change how much work we have to do.

At a fixed degree, suppose the first symmetry contributes a constraint matrix \(A_1\) and the second contributes \(A_2\). Imposing both at once gives

\[ \begin{pmatrix} A_1\\ A_2 \end{pmatrix} c = 0 \]

We can solve this system in the original coefficient space, or reduce it in stages.

First solve

\[ A_1 c = 0 \]

Let the columns of a matrix \(B_1\) form a basis of its null space. Every coefficient vector that survives the first symmetry can then be written as

\[ c = B_1 a \]

where \(a\) contains only the remaining free coefficients. Substituting this into the second block gives

\[ A_2 B_1 a = 0 \]

The second symmetry now acts only on what survived the first. This coefficient-space reduction is always valid, as it is just substitution into the remaining equations.

A separate question is whether \(G_2\) acts autonomously on the first invariant ring \(\mathbb{R}[V]^{G_1}\). For that interpretation, \(G_2\) must preserve the first invariant space. Normalization is one sufficient condition.

Suppose \(G_2\) normalizes \(G_1\), meaning

\[ g G_1 g^{-1} = G_1 \]

for every \(g\in G_2\). Then \(\mathbb{R}[V]^{G_1}\) is stable under the action of \(G_2\).

Indeed, let \(f\) be \(G_1\)-invariant, with \(g\in G_2\) and \(h\in G_1\). Then

\[ h\cdot(g\cdot f) = g\cdot\big((g^{-1}hg)\cdot f\big) \]

Since \(G_2\) normalizes \(G_1\), the element \(g^{-1}hg\) lies in \(G_1\) and therefore fixes \(f\). Hence

\[ h\cdot(g\cdot f) = g\cdot f \]

so \(g\cdot f\) is again \(G_1\)-invariant.

Then we can take the invariants in stages

\[ \mathbb{R}[V]^{\langle G_1,G_2\rangle} = \big(\mathbb{R}[V]^{G_1}\big)^{G_2} \]

When \(G_1\cap G_2=\{e\}\) as well, the generated group is the semidirect product

\[ G_1 \rtimes G_2 \]

We first find the invariants of the normal factor \(G_1\), use them as the new candidate space, and impose \(G_2\) there.

If the factors commute elementwise, the product is direct and either factor may be reduced first. If neither group normalizes the other, the coefficient-space substitution still works, we just lose the interpretation of the second step as an induced \(G_2\) action on \(\mathbb{R}[V]^{G_1}\).

Example: The Euclidean Group

Suppose we want to construct an interaction energy for a collection of particles. Their coordinates depend on where we place the origin and how we orient the axes, but the energy should not depend on those choices. What information survives them?

Translation invariance removes absolute position, leaving relative displacement vectors. Orthogonal invariance then removes absolute orientation and handedness, leaving distances and angles. The Euclidean group lets us perform those two reductions separately and see the physical variables emerge.

Consider the Euclidean group acting on \(k\) points

\[ \mathbf{x}_1,\dots,\mathbf{x}_k\in\mathbb{R}^n \]

with interaction energy

\[ F(\mathbf{x}_1,\dots,\mathbf{x}_k) \]

An orthogonal transformation carries a translation by \(\mathbf{t}\) to a translation by \(R\mathbf{t}\), so the translation group is normal in the Euclidean group

\[ E(n) = \mathbb{R}^n \rtimes O(n) \]

So we can compute its invariants in two stages: translations first, then the orthogonal group.

A translation moves every point together

\[ \mathbf{x}_i \mapsto \mathbf{x}_i + \mathbf{t} \]

The constraint is

\[ F(\mathbf{x}_1+\mathbf{t},\dots,\mathbf{x}_k+\mathbf{t}) = F(\mathbf{x}_1,\dots,\mathbf{x}_k) \]

for every \(\mathbf{t}\in\mathbb{R}^n\). Differentiating at \(\mathbf{t}=0\) gives

\[ \sum_{i=1}^k \frac{\partial F}{\partial\mathbf{x}_i} = 0 \]

The solutions depend only on relative positions. Taking

\[ \mathbf{d}_i = \mathbf{x}_i - \mathbf{x}_k \qquad i=1,\dots,k-1 \]

as the new variables reduces the original \(nk\) coordinates to \(n(k-1)\).

The orthogonal group now acts on this smaller space. Each difference transforms as

\[ \mathbf{d}_i \mapsto R\mathbf{d}_i \]

so the constraint is

\[ F(R\mathbf{d}_1,\dots,R\mathbf{d}_{k-1}) = F(\mathbf{d}_1,\dots,\mathbf{d}_{k-1}) \]

for every \(R\in O(n)\).

The \(O(n)\)-invariants are generated by the pairwise inner products

\[ \mathbf{d}_i \cdot \mathbf{d}_j \]

These are the entries of the Gram matrix, so they encode the lengths and angles of the configuration. There are

\[ \frac{k(k-1)}{2} \]

such entries. When \(k-1>n\), the Gram matrix has rank at most \(n\), so its entries also satisfy determinantal relations.

For the smallest nontrivial instance, take two points

\[ \mathbf{p}=(p_1,p_2) \qquad \mathbf{q}=(q_1,q_2) \]

in the plane, and a homogeneous quadratic interaction energy in their four coordinates

\[ \begin{aligned} F={}&c_1p_1^2+c_2p_2^2+c_3q_1^2+c_4q_2^2 +c_5p_1p_2+c_6p_1q_1\\ &+c_7p_1q_2+c_8p_2q_1+c_9p_2q_2+c_{10}q_1q_2 \end{aligned} \]

There are ten candidate coefficients.

The translation stage imposes

\[ \frac{\partial F}{\partial p_1} + \frac{\partial F}{\partial q_1} = 0 \]

and

\[ \frac{\partial F}{\partial p_2} + \frac{\partial F}{\partial q_2} = 0 \]

Each equation is an identity in the four coordinates, so the coefficient of each coordinate must vanish separately. Collecting the resulting equations gives

\[ \begin{aligned} 2c_1+c_6 &=0, & c_5+c_8 &=0, & 2c_3+c_6 &=0, & c_7+c_{10} &=0, \\ c_5+c_7 &=0, & 2c_2+c_9 &=0, & c_8+c_{10} &=0, & 2c_4+c_9 &=0 \end{aligned} \]

Stacked on the coefficient vector, these become the matrix system

\[ \begin{pmatrix} 2 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6\\ c_7\\ c_8\\ c_9\\ c_{10} \end{pmatrix} = 0 \]

The matrix has rank seven, so its null space is three-dimensional. The surviving combinations assemble into the difference

\[ \mathbf{d} = \mathbf{p} - \mathbf{q} \]

leaving

\[ F = a\,d_1^2+b\,d_1d_2+c\,d_2^2 \]

Translation invariance has replaced the ten coefficients of a general quadratic in the two positions with the three coefficients of a quadratic form in their separation.

The continuous rotation part of \(O(2)\) now acts entirely on the two difference variables. An infinitesimal rotation sends

\[ d_1 \mapsto d_1-\theta d_2 \qquad d_2 \mapsto d_2+\theta d_1 \]

Differentiating at \(\theta=0\) gives the determining equation

\[ \frac{\partial F}{\partial d_1}(-d_2) + \frac{\partial F}{\partial d_2}d_1 = 0 \]

Substituting the quadratic form gives

\[ b\,d_1^2 + 2(c-a)d_1d_2 - b\,d_2^2 = 0 \]

Since this is an identity in \(d_1\) and \(d_2\)

\[ b=0 \qquad c=a \]

The surviving quadratic is the squared distance

\[ F = a(d_1^2+d_2^2) = a|\mathbf{p}-\mathbf{q}|^2 \]

This is already invariant under reflections, so this is the most general homogeneous quadratic interaction energy for two points with Euclidean symmetry.

For \(k\) points, without the quadratic truncation, the complete Euclidean-invariant configurational energy has the form

\[ F(\mathbf{x}_1,\ldots,\mathbf{x}_k) = \Phi(G) \]

where \(G\) is the Gram matrix of the relative displacement vectors

\[ G_{ij} = \mathbf{d}_i\cdot\mathbf{d}_j \qquad i,j=1,\ldots,k-1 \]

The function \(\Phi\) can contain pair, angular, and genuine many-body dependence. No additive decomposition has been assumed. Staging the calculation has only removed the arbitrary origin and orientation, leaving the relational geometry on which the energy can depend.

The same symmetry problem can therefore be reduced one subgroup at a time without changing the final invariant space.

Fields with Derivatives

Next, we can extend this setup to include derivatives. Many systems involve derivatives: for example, elastic energy depends on displacement gradients, and magnetic energy depends on how the magnetization changes from point to point.

Collect the field and its derivatives through order \(k\) into one list, called its \(k\)-jet, and temporarily treat each entry as an independent variable. For a scalar field \(\phi(x_1,x_2)\) on the plane, the first jet is \((\phi,\phi_1,\phi_2)\), where \(\phi_i = \partial\phi/\partial x_i\). At fixed degree, a cost is again a polynomial with finitely many coefficients.

How does the symmetry act on the derivatives? We can determine the action by passing the symmetry action through the chain rule. This extension of the action is called prolongation.⁶ An internal transformation \(\phi \mapsto g\phi\) carries every derivative with it, so \(\phi_i \mapsto g\phi_i\). A spatial rotation \(x \mapsto Rx\) leaves a scalar \(\phi\) fixed and rotates its gradient, \(\nabla\phi \mapsto R\nabla\phi\).

For the linear field transformations and linear spatial actions used here, the prolonged action remains linear in the selected jet variables.

Example: Scalar Field on the Plane

Let’s try an example out: a scalar field on the plane. Suppose its energy density depends on the value of the field and its first spatial derivatives. If the plane is homogeneous and isotropic, the physics should not change when we translate the entire field configuration or rotate it. The question is which quadratic energy densities are compatible with these requirements.

Write the coordinates on the plane as \(\mathbf{x}=(x_1,x_2)\). The two first-order spatial derivatives of the field form its gradient

\[ \nabla\phi = \left( \frac{\partial\phi}{\partial x_1}, \frac{\partial\phi}{\partial x_2} \right) = (\phi_1,\phi_2) \]

Together with the field itself, these variables form the first jet \((\phi,\phi_1,\phi_2)\).

Homogeneity and isotropy mean invariance under two group actions. The translation group \(\mathbb{R}^2\) acts by shifting the argument of the field, while the rotation group \(\mathrm{SO}(2)\) acts by rotating it. For a scalar field, these actions are

\[ (T_{\mathbf a}\phi)(\mathbf x) = \phi(\mathbf x-\mathbf a) \qquad (R\phi)(\mathbf x) = \phi(R^{-1}\mathbf x) \]

Homogeneity means that the energy is invariant under every \(T_{\mathbf a}\), and isotropy means that it is invariant under every \(R\). On the first jet, translations forbid explicit dependence on \(\mathbf x\), while rotations leave \(\phi\) fixed and carry \(\nabla\phi\) as a vector.

Once dependence on \(\mathbf x\) has been removed, translations should impose no further constraint on the jet variables. The general quadratic energy density is

\[ F = c_1 \phi^2 + c_2\, \phi\, \phi_1 + c_3\, \phi\, \phi_2 + c_4\, \phi_1^2 + c_5\, \phi_1 \phi_2 + c_6\, \phi_2^2 \]

with six coefficients. An infinitesimal rotation fixes \(\phi\) and sends \(\phi_1 \mapsto \phi_1 - \theta\, \phi_2\) and \(\phi_2 \mapsto \phi_2 + \theta\, \phi_1\), the same rotation as in the previous example, so differentiating at \(\theta = 0\) gives

\[ \frac{\partial F}{\partial \phi_1}\,(-\phi_2) + \frac{\partial F}{\partial \phi_2}\, \phi_1 = c_3\, \phi\, \phi_1 - c_2\, \phi\, \phi_2 + c_5\, \phi_1^2 + 2(c_6 - c_4)\, \phi_1 \phi_2 - c_5\, \phi_2^2 = 0 \]

This identity must hold for every value of \(\phi\), \(\phi_1\), and \(\phi_2\). The coefficient of each monomial must therefore vanish, which gives \(c_2 = c_3 = c_5 = 0\) and \(c_6 = c_4\). The last three terms consequently reduce to \(c_4(\phi_1^2+\phi_2^2)\).

The sum \(\phi_1^2+\phi_2^2\) is the squared length of this gradient

\[ |\nabla\phi|^2 = \left(\frac{\partial\phi}{\partial x_1}\right)^2 + \left(\frac{\partial\phi}{\partial x_2}\right)^2 = \phi_1^2+\phi_2^2 \]

This is the same rotationally invariant quadratic that appeared for the pair of points, with the gradient now taking the place of the difference vector. The six coefficients reduce to two

\[ F = a\, \phi^2 + b\, |\nabla \phi|^2 \]

a potential term and a gradient term. These are the first two terms we would expect in a Landau free energy.

Note that we didn’t assume the gradient term, it just emergently popped out due to rotation invariance on the first jet at degree two.

Within this class, derivatives make the representation larger but don’t really change the calculations. Once we have prolonged the action, we list the monomials, impose one block per generator, and take the common null space.

After prolongation to a bounded jet, a local field theory becomes the same finite invariant problem.

Extracting the Physics

Once the engine has found an allowed polynomial theory

\[ F = \sum_i c_i I_i \]

we can ask what it predicts. Variational differentiation gives the equations of motion. The symmetry generators give conserved charges. Derivatives of \(F\) give stability and response. These are again finite symbolic calculations: differentiation, substitution, equation solving, and null spaces.

If \(F\) is a Lagrangian density \(L\), making the action stationary gives the Euler-Lagrange equations

\[ \frac{\partial L}{\partial \phi_a} - \frac{d}{dt} \frac{\partial L}{\partial \dot{\phi}_a} - \sum_i \frac{d}{dx_i} \frac{\partial L}{\partial \phi_{a,i}} = 0 \]

There is one equation for each field component. Total derivatives drop out because they change only the boundary contribution to the action. For example, a term linear in the time derivatives

\[ \sum_a I_a(\phi)\, \dot{\phi}_a \]

is locally a total derivative if

\[ \frac{\partial I_a}{\partial \phi_b} = \frac{\partial I_b}{\partial \phi_a} \]

for all \(a, b\).

What about conservation laws? For a vertical transformation in mechanics, meaning one that changes the configuration variables while leaving the time coordinate fixed, Noether’s theorem gives

\[ Q = \sum_a \frac{\partial L}{\partial \dot{\phi}_a}\, \delta \phi_a - F_X \]

Here \(\delta\phi_a\) is the generator’s action on the fields, and \(F_X\) is the boundary term in \(\operatorname{pr}(X)L = dF_X/dt\), which vanishes for a strict symmetry. Transformations of the spacetime coordinates require the full prolonged-current formula, while gauge symmetries require constraint analysis as well. For the vertical transformations used in the free-particle calculation, translations give momentum \(\mathbf{p}=m\mathbf{v}\), rotations give angular momentum \(\mathbf{q}\times\mathbf{p}\), and boosts give \(\mathbf{G}=m\mathbf{q}-\mathbf{p}t\).

If \(F\) is an energy or free energy, its stationary points satisfy

\[ \frac{\partial F}{\partial \phi_a} = 0 \]

The Hessian

\[ K_{ab} = \frac{\partial^2 F}{\partial \phi_a\, \partial \phi_b} \]

determines local stability. For finitely many variables, a positive-definite Hessian means a local minimum; a negative eigenvalue identifies an unstable direction. For fields, the same expansion produces a differential operator. The eigenvalues, or related Fourier transforms of this, give the stability conditions and linearized dispersion relations.

For a regular mechanical Lagrangian, we can also pass to Hamiltonian form. Define \(p_a = \partial L/\partial\dot q_a\). If the velocity Hessian \(\partial^2L/\partial\dot q_a\partial\dot q_b\) is nonsingular, we can solve for the velocities and set

\[ H = \sum_a p_a\, \dot{q}_a - L \]

Hamilton’s equations, \(\dot q_a=\partial H/\partial p_a\) and \(\dot p_a=-\partial H/\partial q_a\), are then equivalent to the Euler-Lagrange equations. If the velocity Hessian is singular, as it is for a Lagrangian linear in time derivatives, the regular Legendre transform fails. The Euler-Lagrange equations still work, but the Hamiltonian treatment needs constraint analysis.

The invariant basis can now be reduced modulo declared boundary terms and turned into equations of motion, charges, and local stability data.

Pipeline

Now let’s turn the hand calculation into code.

Data Structures

The data structures mostly follow the earlier computational invariant theory post. I will repeat the relevant pieces here.

Polynomials

Polynomials are stored as a dictionary from exponent tuples to exact rational coefficients:

Poly = dict[tuple[int, ...], Fraction]

The keys are exponent tuples \((\alpha_0,\dots,\alpha_{n-1})\), and the values are nonzero coefficients. The empty dictionary is the zero polynomial. For velocity variables \((v_0,v_1,v_2)\), the invariant \(|\mathbf{v}|^2\) is {(2,0,0): 1, (0,2,0): 1, (0,0,2): 1}. Addition, multiplication, and differentiation become dictionary operations on Fraction coefficients. The arithmetic is exact: two polynomials are equal exactly when their dictionaries are equal.

Group Actions

The engine sees a group action through an Action, which is just a bundle of capabilities. I have omitted a few optional orbit-computation methods:

@dataclass(frozen=True)
class Action:
    n_vars: int
    invariants_of_degree: Callable[[int], list[Poly]] | None = None
    is_invariant: Callable[[Poly], bool] | None = None
    hilbert_coeffs: Callable[[int], list[int]] | None = None

The different backends do different things. Finite-group and Lie-algebra actions construct explicit bases. Molien and Molien-Weyl backends may return only their dimensions. Classical-group backends can read generators directly from the fundamental theorems.

Group Backends

Where possible, every group encoding below produces the same Action interface. The rest of the pipeline does not need to know which backend produced it.

Finite Groups

A finite group is a list of rational matrices, one for each element. We may instead provide generators; the engine closes them under multiplication and checks the group axioms. These matrices describe the geometric action on physical space. For example, here is the inversion group \(\{e,i\}\):

import numpy as np
from symmetry_engine.groups import finite_group

inv = finite_group([-np.eye(3, dtype=int)])    # closure: {I, -I}

The field representation is supplied separately. A rep maps the eigenvalues of the defining action to the eigenvalues of the induced action on the field. The built-in constructors are scalar(), vector(), pseudovector(), and voigt(); voigt() is the representation for strain:

from symmetry_engine.reps import vector, pseudovector

eigs = inv['defining_eigs'](-np.eye(3, dtype=int))   # [-1, -1, -1]
vector()(eigs)           # [-1, -1, -1]: a polar vector flips
pseudovector()(eigs)     # [+1, +1, +1]: an axial vector does not

New reps can be added in the same form. They compose under direct sums, tensor products, and jets, so we can build the field content of a physical problem from a few primitives. The counting function takes the group and rep together and checks that their dimensions match:

from symmetry_engine.molien import molien

molien(inv, vector(), 4)         # [1, 0, 6, 0, 15]
molien(inv, pseudovector(), 4)   # [1, 3, 6, 10, 15]

The group is the same, but the counts differ. For a polar vector, all odd-degree invariants vanish and six quadratics remain, exactly as in the hand calculation. Inversion acts trivially on a pseudovector, so every monomial survives. The same field types work with the continuous encodings below.

Constructing the basis takes the same inputs. The rep carries its matrix realization along with its eigenvalues, so molien(inv, vector(), 4) and reynolds_matrix(inv, vector()) receive the same group and rep. One counts the invariants; the other projects onto them. invariants_of_degree applies the Reynolds operator

\[ R(f) = \frac{1}{|G|} \sum_{g \in G} g \cdot f \]

to every monomial in the basis. It writes the projected polynomials as coefficient vectors and row-reduces them to obtain a basis for the image of \(R\). When the group merely permutes monomials, this reduces to keeping representatives of monomial orbits. A general matrix action can mix several monomials, so the row-reduction step is essential.

The Molien series packages the number of independent invariants at every degree into one rational function

\[ H(t) = \sum_{d \ge 0} h_d\, t^d = \frac{1}{|G|} \sum_{g \in G} \frac{1}{\det\!\left(I - t\,\rho(g)\right)} \]

Here \(h_d\) is the dimension at degree \(d\). Since the eigenvalues of \(\rho(g)\) are roots of unity, the sum lies in \(\mathbb{Q}(t)\) and can be manipulated exactly.

Lie Groups

A connected Lie group is encoded by matrices for its infinitesimal generators. We stack their induced actions on the monomial basis and take the null space:

from symmetry_engine.lie_algebra_invariants import lie_algebra_action

Lx = [[0, 0, 0], [0, 0, -1], [0, 1, 0]]   # the three infinitesimal
Ly = [[0, 0, 1], [0, 0, 0], [-1, 0, 0]]   # rotations of R^3
Lz = [[0, -1, 0], [1, 0, 0], [0, 0, 0]]

so3 = lie_algebra_action([Lx, Ly, Lz], 3)
so3.invariants_of_degree(1)   # []
so3.invariants_of_degree(2)   # [v_0^2 + v_1^2 + v_2^2]

We get the same quadratic invariant as before

\[ |\mathbf{v}|^2 \]

This time it came from the null space of the infinitesimal generators, not an average. The same construction works for noncompact groups, where no Reynolds average exists. Disconnected components are added separately as finite constraints.

Compact Groups

For compact Lie groups, we can instead count by exact Haar integration. The engine stores an integration rule rather than a list of group elements.

The Molien integrand depends only on the representation’s eigenvalues, hence only on the conjugacy class. Weyl integration reduces the group integral to an integral over a maximal torus.

For \(\mathrm{SO}(3)\), write \(z = e^{i\theta}\). A rotation has eigenvalues

\[ z, \quad 1, \quad z^{-1} \]

and the Weyl measure is

\[ 2 - z - z^{-1} \]

On the torus, everything is a Laurent polynomial in \(z\). Its Haar integral is the constant term, since \(\oint z^k\,dz/z\) vanishes for \(k\neq0\). So integration becomes constant-term extraction on exponent dictionaries with Fraction coefficients. There is no numerical quadrature:

from symmetry_engine.molien_weyl import molien_series

weights = [(1,), (0,), (-1,)]                  # eigenvalues z, 1, 1/z
measure = {(1,): Fraction(-1), (0,): Fraction(2), (-1,): Fraction(-1)}
molien_series(weights, measure, 6, weyl_order=2)
# 1, 0, 1, 0, 1, 0, 1

The coefficients are

\[ 1, 0, 1, 0, 1, 0, 1 \]

There is one invariant at every even degree and none at odd degree. They are the powers of

\[ q = |\mathbf{v}|^2 \]

in agreement with \(\mathbb{R}[V]^{\mathrm{SO}(3)} = \mathbb{R}[q]\).

Classical Groups

For the classical groups, the First Fundamental Theorems already tell us the generators. There is no reason to reconstruct the ring from scratch.

For \(\mathrm{O}(n)\) acting on \(k\) vectors, the generators are the pairwise inner products. For \(\mathrm{SL}(n)\) they are the \(n \times n\) minors, and for \(\mathrm{Sp}(2n)\) they are the symplectic pairings. The theorems themselves, with proofs, are covered in the invariant theory post. As input, a classical group is just the pair \((n, k)\):

from invariants.classical import orthogonal_action

act = orthogonal_action(3, 2)      # O(3) on two vectors of R^3
act.invariants_of_degree(2)        # |v_1|^2,  v_1 . v_2,  |v_2|^2

These are the three inner products of two vectors: the Gram matrix from the Euclidean example. Higher-degree invariants are products of these generators. The method is_invariant checks whether a polynomial can be expressed in them.

The Physical System

The physical input says which fields exist and how each symmetry acts on them. A field may be a scalar, polar vector, axial vector, Voigt strain, gradient, or tensor product of these. Its type selects the representation.

This distinction matters under improper transformations and time reversal. A polar vector changes sign under inversion; an axial vector does not. Getting this choice wrong changes which crystal couplings are forbidden, which require chirality, and which require broken time-reversal symmetry.

The declaration is encoded by two dataclasses:

@dataclass(frozen=True)
class SymmetryDeclaration:
    name: str
    action: Action
    infinitesimal_generators: list[InfinitesimalGenerator]

@dataclass(frozen=True)
class PhysicalSystem:
    name: str
    n_vars: int
    symmetries: list[SymmetryDeclaration]

Each InfinitesimalGenerator stores the polynomial variation \(\delta q^i\). For a translation, \(\delta q=e_j\). For a rotation in the \(ij\)-plane

\[ \delta q^i = -q_j \qquad \delta q^j = q_i \]

The same generator imposes the symmetry constraint and later produces the Noether charge.

Driver

The function discover(system, max_degree, boost=None) -> Certificate runs the calculation. It first checks what the action can do. A constructive run needs invariants_of_degree; a counting-only run can stop at hilbert_coeffs. The constructive path is:

1. Compute the invariant basis for each declared symmetry through the working degree.
2. Intersect those spaces degree by degree, implementing the symmetry-stacking calculation above.
3. Construct the parameterized expression \(F = \sum_i c_i I_i\)
4. Impose any quasi-symmetry conditions, such as invariance up to a total derivative under Galilean boosts.
5. Run the calculations appropriate to the declared functional.⁷ A Lagrangian gives Euler-Lagrange equations and Noether charges. An energy or free energy gives Hessians and stability conditions. If the required terms are present, the engine can also extract dispersion relations, characteristic scales, and response tensors.

The declaration does not yet record which kind of functional \(F\) is, so the caller currently chooses this branch.

A separate module performs the regular Legendre transform and checks its identities symbolically. The certificate does not yet include the Hamiltonian data.

Certificates

A run returns a Certificate with the invariant basis, parameterized expression, derived quantities, and intermediate steps. Depending on the backend, this may include eigenvalue classes, Molien functions, series coefficients, null-space matrices, and invariant bases.

Why keep the intermediate steps? I want the calculation to be reproducible line by line, not just return a final formula. A separate Lean 4 layer can check selected arithmetic identities from some runs. That layer is still a work in progress, and any verification claim needs to state its exact coverage.

The Input Tuple

For the rest of the post, I will collect the whole declaration into one record. A catalog cell is

\[ \Pi = (\mathcal{G},\ V,\ \tau,\ \sim,\ A) \]

What goes in each slot?

• \(\mathcal{G}\) is a finite list of symmetry actions. They may be imposed jointly or, for a comparison, evaluated separately. An action may be given by exact matrices, infinitesimal generators, or an eigenvalue rule. A quasi-symmetry also carries its total-derivative condition.

• \(V\) is the carrier: either a finite list of variables, including derivatives through a stated jet order, or a family of representations.

• \(\tau\) is the truncation. It may bound polynomial degree and derivative order separately in different sectors, or cut off a representation family.

• \(\sim\) is the equivalence on candidates. In mechanics, \(\mathrm{im}\,d_t\) abbreviates the identification \(L \sim L + \tfrac{dF}{dt}\).

• \(A\) contains anything supplied after the invariant calculation: coefficient identifications, ansatz families, filling rules, or other auxiliary inputs.

The symbol \(\varnothing\) means there is no extra choice in that slot.

From one cell, the engine returns the allowed space, the general model on that space, and whatever consequences make sense for the carrier. These might be equations of motion and charges, representation dimensions and branchings, or a generating function. The next four examples are four choices of \(\Pi\).

All of these group encodings now feed the same physical declaration and return the same kind of certificate.

Example 1: The Free Particle

The free particle instantiates the input tuple as

\[ \Pi_1 = \big(\, \{\mathrm{E}(3),\ \mathbb{R}^3_{\mathbf{u}}\},\ (\mathbf{q}, \mathbf{v}),\ \deg \le 4,\ \mathrm{im}\, d_t,\ \{c_1 = \tfrac12 m\} \,\big) \longrightarrow L = c_0 + c_1\, |\mathbf{v}|^2 \]

The symmetry actions are the Euclidean group and the Galilean boosts \(\mathbf{v}\mapsto\mathbf{v}+\mathbf{u}\). The carrier is \((\mathbf{q},\mathbf{v})\), the truncation is \(\deg\le4\), and Lagrangians are identified modulo total time derivatives. The final entry, \(c_1=\tfrac12m\), is only a name for a coefficient after the calculation.

---

We already know the answer by hand, which makes the free particle a useful test. Can the same declaration recover it mechanically and show how?

Translations and rotations act on \((\mathbf{q},\mathbf{v})\). The Galilean boost is a quasi-symmetry, and the degree bound is four. We do not omit position dependence in advance. Translation invariance has to remove it. The engine then intersects the translation and rotation invariants, applies the boost condition, and computes the Noether charges.

from fractions import Fraction

from invariants.classical import orthogonal_action
from invariants.poly import mono, const, neg
from symmetry_engine.lie_algebra_invariants import translation_action
from symmetry_engine.system import (
    PhysicalSystem, SymmetryDeclaration, InfinitesimalGenerator,
    BoostConstraint,
)
from symmetry_engine.discover import discover

n = 6  # configuration-velocity variables (q_0, q_1, q_2, v_0, v_1, v_2)

translations = SymmetryDeclaration(
    name="Translation",
    action=translation_action((0, 1, 2), n),
    infinitesimal_generators=[
        InfinitesimalGenerator(f"translation_{k}", {k: const(1, n)})
        for k in range(3)
    ],
)

rotation_gens = []
for i in range(3):
    for j in range(i + 1, 3):
        rotation_gens.append(InfinitesimalGenerator(
            name=f"angular_momentum_{i}{j}",
            components={
                i: neg(mono(tuple(1 if k == j else 0 for k in range(n)))),
                j: mono(tuple(1 if k == i else 0 for k in range(n))),
            },
        ))
rotations = SymmetryDeclaration(
    name="O(3) rotation",
    action=orthogonal_action(3, 2),  # O(3) acting on the pair (q, v)
    infinitesimal_generators=rotation_gens,
)

system = PhysicalSystem(
    name="Free particle in R^3",
    n_vars=n,
    symmetries=[translations, rotations],
    velocity_indices=(3, 4, 5),
    var_names=("q_0", "q_1", "q_2", "v_0", "v_1", "v_2"),
)

cert = discover(system, max_degree=4,
                boost=BoostConstraint(name="Galilean boost", space_dim=3))
print(cert.summary(identify={"c_1": ("m", Fraction(1, 2))}))

=== Certificate: Free particle in R^3 ===
Configuration-velocity space: R^6
Symmetries: Translation, O(3) rotation

Invariant generators:
  joint_inv_0 = v_2^2 + v_1^2 + v_0^2
  (Exact intersection of Translation and O(3) rotation invariants)

Forced Lagrangian (most general):
  L = c_0 * (1) + c_1 * (v_2^2 + v_1^2 + v_0^2)

Noether rule:
  Given: L is invariant under the declared symmetry (or quasi-invariant
         up to a total derivative)
  Given: q(t) satisfies the Euler-Lagrange equations for L
  Then: the generator's charge J is conserved, d(J)/dt = 0
Applications:
  J_translation_0 = c_1 * (2*v_0)    [L is Translation-invariant]
  ...
  J_angular_momentum_01 = c_1 * (-2*q_1*v_0 + 2*q_0*v_1)    [L is O(3) rotation-invariant]
  ...

Parameter identification (supplied, not computed):
  c_1 := (1/2) * m
Identified charges:
  J_translation_0 = m * (v_0)
  ...
  J_angular_momentum_01 = m * (-q_1*v_0 + q_0*v_1)
  ...

Derivation:
  Step 1:
    Given: Translation acts on R^6
    Claim: The degree-2 invariant subspace has dimension 6
      Basis:
        v_2^2
        v_1*v_2
        ...
      Monomial dimension: 21
      Constraint rank: 15
      Invariant dimension: 6
    By: Kernel of the translation derivations d/dq_i (monomials free of
        the position variables)
  Step 2:
    Given: O(3) rotation acts on R^6
    Claim: The degree-2 invariant subspace has dimension 3
      Basis:
        q_2^2 + q_1^2 + q_0^2
        q_2*v_2 + q_1*v_1 + q_0*v_0
        v_2^2 + v_1^2 + v_0^2
      Monomial dimension: 21
      Constraint rank: 18
      Invariant dimension: 3
    By: Classical-group backend (generators from the First Fundamental Theorem)
  ...
  Step 4:
    Given: The 2 degree-2 invariant subspaces above
    Claim: The joint degree-2 invariant subspace has dimension 1
      Basis:
        v_2^2 + v_1^2 + v_0^2
      Per-symmetry dimensions: 6, 3
    By: Exact intersection of the invariant subspaces (kernel of the
        stacked coefficient system)
  Step 5:
    Given: Invariant ring generated by: v_2^2 + v_1^2 + v_0^2
    Given: First Fundamental Theorem for O(3) rotation
    Claim: Most general invariant Lagrangian of degree <= 4 has 3 free
           parameter(s): c_0, c_1, c_2
      Completeness check:
        degree 1: joint dimension 0, spanned by products 0
        degree 2: joint dimension 1, spanned by products 1
        degree 3: joint dimension 0, spanned by products 0
        degree 4: joint dimension 1, spanned by products 1
    By: Products of the joint quadratic invariants, with completeness
        verified degreewise against the exact joint dimensions
  Step 6:
    Given: Galilean boost: L(v + u) - L(v) is at most linear in v for all u
    Claim: The boost constraint removes 1 basis element(s), leaving 2 free
           parameter(s): c_0, c_1
      Input basis:
        1
        v_2^2 + v_1^2 + v_0^2
        v_2^4 + 2*v_1^2*v_2^2 + v_1^4 + 2*v_0^2*v_2^2 + 2*v_0^2*v_1^2 + v_0^4
      Removed by boost constraint:
        v_2^4 + 2*v_1^2*v_2^2 + v_1^4 + 2*v_0^2*v_2^2 + 2*v_0^2*v_1^2 + v_0^4
      Surviving basis:
        1
        v_2^2 + v_1^2 + v_0^2
    By: Galilean boost requires L to be at most quadratic in velocity
        (degree > 2 terms produce irreducible dependence on u that cannot
        be a total derivative)
  ...

Translation and rotation invariance together allow

\[ 1 \qquad |\mathbf{v}|^2 \qquad |\mathbf{v}|^4 \]

and the boost removes the quartic term. The two survivors are therefore

\[ L = c_0 + c_1 |\mathbf{v}|^2 \]

The constant remains for bookkeeping, although it has no dynamics. The substitution

\[ c_1 = \tfrac12 m \]

is supplied after the calculation. It is not derived. With this identification, the translation charges become the momenta \(mv_i\), and the rotation charges become angular momentum. For \(m\neq0\), the Euler-Lagrange equations give

\[ \dot{\mathbf{v}} = 0 \]

and hence uniform motion. If only the constant remains, there is no dynamics.

What happens if we drop translation invariance? Position dependence returns, and rotations alone allow three quadratic terms

\[ |\mathbf{q}|^2 \qquad \mathbf{q}\cdot\mathbf{v} \qquad |\mathbf{v}|^2 \]

The mixed term is a total derivative

\[ \mathbf{q}\cdot\mathbf{v} = \tfrac12\frac{d}{dt}|\mathbf{q}|^2 \]

so it cannot affect the equations of motion. Keeping the free kinetic term while allowing a rotationally symmetric potential gives

\[ L = \tfrac12 m |\mathbf{v}|^2 - V(|\mathbf{q}|) \]

the central-force problem. The harmonic oscillator and Kepler problem correspond to \(V \propto |\mathbf{q}|^2\) and \(V \propto 1/|\mathbf{q}|\), respectively.

Symmetry has fixed the form of the free kinetic term and the central-force extension, but has not fixed the coefficients or the function \(V\).

Example 2: The Nonlinear Schrödinger Equation in a First-Order \(U(1)\) Branch

The Schrödinger equation is first order in time, second order in space, and invariant under a constant change of phase. How much of its form follows after we declare those features?

Take a complex scalar field \(\psi(t,x)\) in one spatial dimension. I allow local terms built from the field, its first time derivative, and its first spatial derivative. Derivative terms stop at quadratic degree, while the potential continues through quartic degree. I impose spatial parity, keeping the branch with at most one factor of \(\dot{\psi}\), and identifying Lagrangians that differ by a total time derivative.

In the tuple notation, the declaration is

\[ \Pi_2 = \big(\, \{U(1),\ \mathbb{Z}_2^{x \mapsto -x}\},\ (\phi, \dot{\phi}, \phi'),\ \tau_2,\ \mathrm{im}\, d_t,\ \{\tfrac{1}{2m} = \tfrac{c_2}{c_1},\ U_0 = \tfrac{c_3}{c_1},\ G = \tfrac{2c_4}{c_1}\} \,\big) \] \[ \longrightarrow i\, \dot{\psi} = -\tfrac{1}{2m}\, \psi'' + U_0\, \psi + G\, |\psi|^2 \psi \]

Here \(\phi\) denotes the two real components of \(\psi\), and \(\tau_2\) is the truncation just described. I use units with \(\hbar=1\).

Phase Symmetry and Truncation

Why \(U(1)\)? The overall phase of a quantum state is unobservable. Multiplying the field by a constant phase should change nothing

\[ \psi \mapsto e^{i\theta} \psi \]

These transformations form \(U(1)\). The associated Noether charge is proportional to

\[ N = \int |\psi|^2\, dx \]

the total norm, or particle number.

The engine does its invariant calculation over the reals, so split the complex field into two real fields

\[ \psi = \phi_1 - i\phi_2 \]

and write

\[ \phi = (\phi_1, \phi_2) \]

The phase now acts as an ordinary rotation

\[ R(\theta) = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \]

So \(U(1)\) becomes \(\mathrm{SO}(2)\) on the pair \((\phi_1,\phi_2)\).

There is one spatial symmetry as well

\[ x \mapsto -x \]

I take \(\psi\) to be a spatial scalar. Parity leaves \(\phi\) and \(\dot\phi\) alone, but reverses \(\phi'\)

\[ (\phi, \dot{\phi}, \phi') \mapsto (\phi, \dot{\phi}, -\phi') \]

The carrier is

\[ V = (\phi, \dot{\phi}, \phi') \]

This is the first jet in one spatial dimension. Since the phase rotation is global, the same matrix acts on all three pairs

\[ R(\theta) \oplus R(\theta) \oplus R(\theta) \]

Notice what is missing: \(t\) and \(x\) themselves. The candidate Lagrangians therefore have no explicit coordinate dependence. Homogeneity in space and time is part of the declaration.

Next, decide which monomials to keep. Label each one by

\[ (a, b, c) = (\deg_{\phi},\ \deg_{\dot{\phi}},\ \deg_{\phi'}) \]

The certificate prints this triple next to each invariant. The truncation is

\[ \tau_2(a, b, c) \colon \quad \begin{cases} a + b + c \le 2 & \text{if } b + c > 0 \\ a \le 4 & \text{if } b = c = 0 \\ b \le 1 & \text{in the first-order branch} \end{cases} \]

Derivative terms stop at quadratic degree, but the potential continues through quartic degree. The branch we want contains at most one factor of \(\dot\phi\).

Why impose that last condition? Because phase symmetry permits both

\[ |\dot{\phi}|^2 \]

which gives second-order equations in time, and the antisymmetric pairing

\[ \varepsilon(\phi, \dot{\phi}) = \phi_1 \dot{\phi}_2 - \phi_2 \dot{\phi}_1 \]

which is linear in \(\dot\phi\) and gives first-order equations. \(U(1)\) does not choose between them. I choose the first-order branch and keep the second-order term as a control.

There is one more piece of bookkeeping. Identify

\[ L \sim L + d_t F \]

where \(F\) is a polynomial in the fields. A total time derivative changes the action at the endpoints but leaves the Euler-Lagrange equations alone. We could also quotient by total spatial derivatives, although that would not change the parity-even terms here.

The final slot contains the coefficient names

\[ \frac{1}{2m} = \frac{c_2}{c_1} \qquad U_0 = \frac{c_3}{c_1} \qquad G = \frac{2 c_4}{c_1} \]

These just rename ratios of free coefficients after the calculation. I assume \(c_1\neq0\) in the first-order branch.

What survives?

Deriving the Schrödinger Equation

At quadratic order, there are only two ways to pair real two-component vectors without breaking \(\mathrm{SO}(2)\). The first is the dot product

\[ \mathbf{a} \cdot \mathbf{b} \]

The second is the antisymmetric pairing

\[ \varepsilon(\mathbf{a}, \mathbf{b}) = a_1 b_2 - a_2 b_1 \]

Apply both to \(\phi\), \(\dot\phi\), and \(\phi'\). This gives nine terms

\[ \begin{aligned} L_2 = {}& a_1\, |\phi|^2 + a_2\, |\dot{\phi}|^2 + a_3\, |\phi'|^2 + a_4\, \phi \cdot \dot{\phi} + a_5\, \phi \cdot \phi' + a_6\, \dot{\phi} \cdot \phi' \\ &+ a_7\, \varepsilon(\phi, \dot{\phi}) + a_8\, \varepsilon(\phi, \phi') + a_9\, \varepsilon(\dot{\phi}, \phi') \end{aligned} \]

Parity kills the four terms with one factor of \(\phi'\)

\[ \phi \cdot \phi' \qquad \dot{\phi} \cdot \phi' \qquad \varepsilon(\phi, \phi') \qquad \varepsilon(\dot{\phi}, \phi') \]

One more term looks dynamical but is not

\[ \phi \cdot \dot{\phi} = \phi_1 \dot{\phi}_1 + \phi_2 \dot{\phi}_2 = \frac{d}{dt} \frac12 |\phi|^2 \]

The boundary term disappears in the quotient.

The antisymmetric pairing is different

\[ \varepsilon(\phi, \dot{\phi}) = \phi_1 \dot{\phi}_2 - \phi_2 \dot{\phi}_1 \]

If it were \(dF/dt\) for some \(F(\phi_1,\phi_2)\), then

\[ \frac{\partial F}{\partial \phi_1} = -\phi_2 \qquad \frac{\partial F}{\partial \phi_2} = \phi_1 \]

But the mixed partials would be

\[ \frac{\partial^2 F}{\partial \phi_2\, \partial \phi_1} = -1 \qquad \frac{\partial^2 F}{\partial \phi_1\, \partial \phi_2} = 1 \]

They disagree. No such \(F\) exists.

Now remove \(|\dot\phi|^2\) to take the first-order branch. Three quadratic structures remain

\[ \varepsilon(\phi, \dot{\phi}) \qquad |\phi'|^2 \qquad |\phi|^2 \]

The potential is even simpler. Every derivative-free \(U(1)\) invariant is built from

\[ |\phi|^2 = \phi_1^2 + \phi_2^2 \]

so quartic order adds

\[ |\phi|^4 = (\phi_1^2 + \phi_2^2)^2 \]

An additive constant has no effect on the dynamics. After relabeling the coefficients, the answer is

\[ L = c_1\, \varepsilon(\phi, \dot{\phi}) - c_2\, |\phi'|^2 - c_3\, |\phi|^2 - c_4\, |\phi|^4 \]

The important term is the first one. It is linear in \(\dot\phi\), but unlike \(\phi\cdot\dot\phi\), it is not a boundary term. It can generate first-order evolution.

Varying the two real components gives

\[ \begin{aligned} 2 c_1 \dot{\phi}_2 + 2 c_2 \phi_1'' - 2 c_3 \phi_1 - 4 c_4 |\phi|^2 \phi_1 &= 0 \\ -2 c_1 \dot{\phi}_1 + 2 c_2 \phi_2'' - 2 c_3 \phi_2 - 4 c_4 |\phi|^2 \phi_2 &= 0 \end{aligned} \]

Now put the two real fields back together

\[ \psi = \phi_1 - i\phi_2 \]

The two real equations become one complex equation

\[ i\, \dot{\psi} = -\frac{c_2}{c_1}\, \psi'' + \frac{c_3}{c_1}\, \psi + \frac{2 c_4}{c_1}\, |\psi|^2 \psi \]

Introducing

\[ \frac{1}{2m} = \frac{c_2}{c_1} \qquad U_0 = \frac{c_3}{c_1} \qquad G = \frac{2 c_4}{c_1} \]

gives

\[ i\, \frac{\partial \psi}{\partial t} = -\frac{1}{2m}\, \frac{\partial^2 \psi}{\partial x^2} + U_0\, \psi + G\, |\psi|^2 \psi \]

This is the cubic nonlinear Schrödinger equation, or Gross-Pitaevskii equation. Set \(c_4=0\) and it becomes linear. Set \(c_3=c_4=0\) and it becomes free. A constant \(U_0\) can also be removed by a time-dependent phase.

The calculation has fixed the form of the Lagrangian, not its coefficients.⁸

Engine Output

Does the engine find the same thing? Name the real variables

\[ (p_1, p_2) \leftrightarrow \phi \qquad (d_1, d_2) \leftrightarrow \dot{\phi} \qquad (g_1, g_2) \leftrightarrow \phi' \]

and act on them with three copies of the same rotation.

For this run, I replace the continuous rotation group by the cyclic subgroup \(C_{12}\)

\[ R\!\left(\frac{2\pi j}{12}\right) \qquad j = 0, \ldots, 11 \]

Why is this exact? After complexification, each variable has charge \(+1\) or \(-1\). A degree-four monomial has total charge between \(-4\) and \(4\). A \(C_{12}\) invariant must have charge zero modulo twelve. The only possibility in that range is zero. So \(C_{12}\) and \(U(1)\) have the same invariants through degree four.

The \(\mathrm{O}(2)\) control uses the dihedral group \(D_{12}\). It adds \(\psi\mapsto\bar\psi\), or \(\operatorname{diag}(1,-1)\) on each real pair.

The internal-symmetry bases are constructed as follows:

from symmetry_engine.so2_surrogates import cn_on_copies, dn_on_copies
from symmetry_engine.reynolds_rational import invariant_polynomials_rational

fields = ["p1", "p2", "d1", "d2", "g1", "g2"]
basis_so2 = invariant_polynomials_rational(cn_on_copies(12, 3), 2,
                                           n=6, var_names=fields)
basis_o2 = invariant_polynomials_rational(dn_on_copies(12, 3), 2,
                                          n=6, var_names=fields)
len(basis_so2), len(basis_o2)   # (9, 6)

The full experiment applies spatial parity, the truncation \(\tau_2\), and the total-derivative equivalence. Its degree-two output is

degree-2 bases constructed: SO(2) 9, O(2) 6
   (2, 0, 0)  p1**2 + p2**2
   (1, 1, 0)  d1*p1 + d2*p2
   (1, 1, 0)  -d1*p2 + d2*p1
   (1, 0, 1)  g1*p1 + g2*p2
   (1, 0, 1)  -g1*p2 + g2*p1
   (0, 2, 0)  d1**2 + d2**2
   (0, 1, 1)  d1*g1 + d2*g2
   (0, 1, 1)  d1*g2 - d2*g1
   (0, 0, 2)  g1**2 + g2**2
  per-sector counts match exact charge counting
  SO(2)-only terms = the 3 epsilon invariants: PASS

There are the same nine terms as before: three norms, three dot products, and three antisymmetric pairings. The \(\mathrm{O}(2)\) control loses all three antisymmetric terms, including the one needed for first-order evolution.

The parity and boundary-term stages report

parity filter: kept 5, dropped ['g1*p1 + g2*p2', '-g1*p2 + g2*p1', 'd1*g1 + d2*g2', 'd1*g2 - d2*g1']
total-derivative quotient: dropped d1*p1 + d2*p2 = d/dt(p1**2/2 + p2**2/2)  PASS

Parity removes the four terms with one spatial derivative. The boundary-term test then finds

\[ d_1 p_1 + d_2 p_2 \]

as the time derivative of

\[ \frac12 (p_1^2 + p_2^2) \]

This is the same boundary term we removed by hand.

In the derivative-free quartic sector, it finds

degree-4 potential sector (4,0,0): p1**4 + 2*p1**2*p2**2 + p2**4  PASS

This is

\[ (p_1^2 + p_2^2)^2 \]

or \(|\phi|^4\).

The first-order branch now assembles the Lagrangian and varies it

BRANCH B (first-order): L = c1*(-d1*p2 + d2*p1) - c2*(g1**2 + g2**2) - c3*(p1**2 + p2**2) - c4*(p1**4 + 2*p1**2*p2**2 + p2**4)
  EL_1: Eq(2*c1*Derivative(phi2(t, x), t) + 2*c2*Derivative(phi1(t, x), (x, 2)) - 2*c3*phi1(t, x) - 4*c4*phi1(t, x)**3 - 4*c4*phi1(t, x)*phi2(t, x)**2, 0)
  EL_2: Eq(-2*c1*Derivative(phi1(t, x), t) + 2*c2*Derivative(phi2(t, x), (x, 2)) - 2*c3*phi2(t, x) - 4*c4*phi1(t, x)**2*phi2(t, x) - 4*c4*phi2(t, x)**3, 0)
  GP match (psi = phi1 - i phi2):  PASS
    i dpsi/dt = -(1/2m) psi'' + U_0 psi + G |psi|^2 psi
    m = c1/(2*c2),  U_0 = c3/c1,  G = 2*c4/c1
    m = c1/(2 c2): PASS;  c4 = 0 -> linear Schrodinger (G = 0): PASS

The controls give the other two branches

BRANCH A (retains |dphi/dt|^2): EL is second order in time  PASS

O(2) control: surviving first-order kinetic terms = 0 (epsilon forbidden, phi.phidot is a boundary term) -> no Schrodinger-type equation  PASS

Keeping \(|\dot\phi|^2\) gives second-order equations. Enlarging \(U(1)\) to \(\mathrm{O}(2)\) removes the useful antisymmetric term, and the only other term linear in \(\dot\phi\) is a boundary term. The Schrödinger branch needs both choices.

The added \(\mathrm{O}(2)\) reflection acts as complex conjugation at fixed \(t\), which is an internal control on the representation, not physical antiunitary time reversal, which also reverses time.

The hand calculation and certificate agree. Four coefficients remain, and one ratio is called the mass

\[ m = \frac{c_1}{2 c_2} \]

Phase symmetry has not determined its value. Galilean invariance would connect the same ratio to the mass in the projective boost law, but it would not determine the value either.

Example 3: The Lowest-Order Spin Splitting in CrSb

(AI disclosure: ChatGPT assisted with the background exposition on altermagnets. I am not a specialist in this area.)

Magnetism is one of the main ways we store and control information. A magnetic region can hold a stable state, like a tiny switch that stays pointed one way or the other, which is why magnetic effects are used in memory, sensors, and many electronic devices.

Ordinary ferromagnets are comparatively easy to read and control because their magnetic moments add to a net magnetization. But that magnetization also produces stray fields, which can interfere with nearby components as devices become smaller. Antiferromagnets avoid these stray fields because their opposing magnetic moments cancel. The same cancellation makes their magnetic state harder to detect and manipulate.

Altermagnets have recently attracted attention because they may offer a middle ground between ferromagnets and antiferromagnets. Like antiferromagnets, their opposing magnetic moments cancel, leaving no net magnetization. But like ferromagnets, they can still separate up- and down-spin electronic states in energy. This combination could be useful for dense, fast memory and other devices that control information through electron spin, although most such applications remain prospective (Song et al. 2025; Jungwirth et al. 2026).

CrSb, or chromium antimonide, is one of the clearest experimental examples. Measurements have directly observed its spin-split electronic bands and mapped how the splitting changes across three-dimensional momentum space (Reimers et al. 2024; Yang et al. 2025).

In this example, we use the symmetry engine to determine the form of that splitting. The previous examples recovered familiar theories from familiar symmetries. Here, the target is a detailed pattern measured in a real material. We expect the splitting to first appear at a particular polynomial degree, vanish on several planes, and vary with direction in a specific way. We will give the engine the relevant crystal symmetries, but not the observed formula, and ask it to recover that structure. This is a partial reproduction of the recent experimental result rather than a microscopic model of CrSb.

Engine Input

Electron states in a crystal are labeled by crystal momentum

\[ \mathbf{k} = (k_x,k_y,k_z) \]

Crystal momentum plays the role of momentum while accounting for the repeating structure of the crystal. Several electronic states may share the same \(\mathbf{k}\), and their energies trace out bands as \(\mathbf{k}\) varies.

We do not need a general theory of electronic bands for this calculation. The engine uses the three components of \(\mathbf{k}\) as polynomial variables and treats the splitting within one band at a time.

We consider the two states whose spins point up and down along the magnetic axis. Write their energies as

\[ E_\uparrow(\mathbf{k}) = \varepsilon_0(\mathbf{k})+\Delta(\mathbf{k}) \qquad E_\downarrow(\mathbf{k}) = \varepsilon_0(\mathbf{k})-\Delta(\mathbf{k}) \]

The function \(\varepsilon_0\) is their average energy, while \(\Delta\) is half the spin splitting. The full energy difference has magnitude \(2|\Delta(\mathbf{k})|\).

We work near \(\Gamma\), which means the point \(\mathbf{k}=0\). Near this point, a smooth function such as \(\Delta(\mathbf{k})\) can be approximated by a polynomial in \(k_x\), \(k_y\), and \(k_z\). We ask the engine to search through degree four and determine the first degree at which a nonzero splitting is allowed.

The calculation uses three local symmetry rules

Operation on momentum	Required action on \(\Delta\)
Rotate the \(x\)-\(y\) plane by \(60^\circ\)	Change sign
Reflect \(k_z\mapsto-k_z\)	Change sign
Reflect \(k_y\mapsto-k_y\)	Change sign

The chromium sites form two interpenetrating magnetic sublattices, meaning two sets of sites with oppositely directed magnetic moments. Each of these operations exchanges the two sublattices. This in turn exchanges the positive and negative spin-energy shifts, forcing \(\Delta\) to change sign.

The invariant backend normally searches for expressions that remain unchanged. We introduce a formal variable \(s\) that changes sign under the same three operations. The product \(s f(\mathbf{k})\) is then invariant precisely when \(f(\mathbf{k})\) has the required sign changes. The variable \(s\) is bookkeeping rather than a physical coordinate.

The input tuple is

\[ \Pi_3 = \left( \mathcal{G}_{\mathrm{CrSb}}^{(\Gamma)}, (\mathbf{k},s), \left\{\deg_{\mathbf{k}}\leq4,\ \deg_s=1\right\}, \varnothing, \varnothing \right) \]

Here \(\mathcal{G}_{\mathrm{CrSb}}^{(\Gamma)}\) is the group generated by the three local actions in the table. The carrier \((\mathbf{k},s)\) contains three physical momentum coordinates and one formal sign coordinate. The truncation keeps one factor of \(s\) and momentum degree at most four. There is no equivalence relation on the candidate splittings and no auxiliary input after the invariant calculation.

Derivation

First impose the \(60^\circ\) rotation. Write the in-plane momentum as

\[ k_x = k_\perp\cos\theta \qquad k_y = k_\perp\sin\theta \]

The rotation sends \(\theta\mapsto\theta+\pi/3\). The first polynomial angular patterns that reverse sign under this shift are

\[ k_x^3 - 3k_xk_y^2 \]

and

\[ k_y(3k_x^2 - k_y^2) \]

They are proportional to \(\cos3\theta\) and \(\sin3\theta\). When \(\theta\) increases by \(\pi/3\), their arguments increase by \(\pi\), so both expressions reverse sign. The rotation therefore requires the in-plane part of the splitting to begin at cubic degree.

Next impose the horizontal reflection \(k_z\mapsto-k_z\). Since this reflection must also reverse \(\Delta\), the splitting must be odd in \(k_z\). The lowest possible expression contains one factor of \(k_z\)

\[ \Delta(\mathbf{k}) = k_z \left[ a\left(k_x^3-3k_xk_y^2\right) + b\,k_y\left(3k_x^2-k_y^2\right) \right] \]

The constants \(a\) and \(b\) are the two coefficients left after the first two symmetry rules. The cubic in-plane pattern and the additional factor of \(k_z\) make the first possible term quartic.

Finally, impose the vertical reflection \(k_y\mapsto-k_y\). The first cubic is even in \(k_y\)

\[ k_x^3-3k_xk_y^2 \mapsto k_x^3-3k_xk_y^2 \]

This cubic cannot produce the required sign change, so \(a=0\). The second cubic is odd in \(k_y\)

\[ k_y(3k_x^2-k_y^2) \mapsto -k_y(3k_x^2-k_y^2) \]

This cubic has the required sign change and therefore survives. The unique lowest-order splitting is therefore

\[ \boxed{ \Delta(\mathbf{k}) = \lambda k_zk_y(3k_x^2-k_y^2) } \]

Symmetry fixes the polynomial up to the coefficient \(\lambda\). This degree-four angular pattern is called \(g\)-wave. The name describes the shape of the momentum-space splitting and does not mean that the electron occupies an atomic \(g\) orbital.

Consequences

Using the polar coordinates above

\[ k_y(3k_x^2 - k_y^2) = k_\perp^3\sin 3\theta \]

The magnitude of the spin splitting is therefore

\[ \left| E_\uparrow(\mathbf{k}) - E_\downarrow(\mathbf{k}) \right| = 2|\lambda|\,|k_z|\,k_\perp^3|\sin 3\theta| \]

The signed splitting has a \(\sin3\theta\) angular form: rotating the momentum by \(60^\circ\) reverses which spin state has the higher energy. Its magnitude, proportional to \(|\sin3\theta|\), repeats every \(60^\circ\).

The polynomial factors as

\[ \Delta(\mathbf{k}) = \lambda k_zk_y \left(\sqrt{3}k_x-k_y\right) \left(\sqrt{3}k_x+k_y\right) \]

The factored polynomial therefore vanishes on four planes through \(\Gamma\)

\[ k_z=0 \]

and

\[ k_y=0 \qquad k_y=\sqrt{3}\,k_x \qquad k_y=-\sqrt{3}\,k_x \]

These are the four nodal planes. On each plane, \(\Delta=0\), so the up and down spin states have the same energy.

Engine Output

target: invariants linear in s
  momentum degree 0: dimension 0
  momentum degree 1: dimension 0
  momentum degree 2: dimension 0
  momentum degree 3: dimension 0
  momentum degree 4: dimension 1
basis in Cartesian coordinates, up to scale:
  s*k_z*k_y*(3*k_x**2 - k_y**2)

This is the first example in which the engine reproduces a detailed, material-specific result rather than a standard textbook theory. The important output is not only the final polynomial. The engine identifies degree four as the first allowed order and computes that the allowed space at this degree is one-dimensional. The symmetry-fixed angular dependence and the four nodal planes then follow from the single returned basis element.

The calculation uses the same representation and invariant-space operations as the earlier examples. It can therefore be repeated without inventing a separate derivation for every material. Replacing the group, carrier, or degree bound produces another input tuple and another cell in the same catalog.

The calculation determines the shape of the splitting, but not its scale. The coefficient \(\lambda\) must come from experiment or a microscopic calculation. The result is also a local expansion around \(\Gamma\), or \(\mathbf{k}=0\). A model valid throughout the crystal’s momentum space would have to respect the periodicity of the lattice rather than use a polynomial centered at one point.

Example 4: A Mechanical Domain Coupling in CrSb

(AI disclosure: ChatGPT assisted with the background exposition on altermagnets. I am not a specialist in this area.)

If altermagnets are to be useful for memory, observing a spin splitting is not enough. A memory element needs at least two stable states and a way to choose between them. In a magnetic material, those states often appear as domains, which are regions that realize different symmetry-related arrangements of the magnetic moments.

CrSb has two such domains related by reversing every ordered moment. The reversal leaves the net magnetization at zero, but it reverses the momentum-space spin splitting. After choosing which arrangement to call positive, let \(\eta=+1\) label that domain and \(\eta=-1\) label the reversed domain. Their local splittings are

\[ \Delta_\eta(\mathbf{k}) = \eta|\lambda|\, k_z k_y(3k_x^2 - k_y^2) \]

Time reversal reverses every magnetic moment, so it sends \(\eta\mapsto-\eta\) and exchanges the two splitting patterns. Choosing the domain therefore chooses the sign of the splitting at every momentum.

The absence of net magnetization makes the domains difficult to manipulate in the same way as ordinary ferromagnetic domains. Their order is tied to the crystal symmetry, however, which makes a mechanical perturbation a natural possibility. Experiments have already shown that changing the crystal symmetry with strain can manipulate the altermagnetic order in CrSb (Zhou et al. 2025).

Here we ask a narrower question. Which local mechanical fields can distinguish the two domains under the CrSb symmetries? A perturbation can favor one domain only if it lowers the energy of one arrangement relative to the other. The contribution to an effective energy density must therefore be odd in \(\eta\)

\[ F_{\mathrm{bias}} = -h_\eta\eta \]

A positive \(h_\eta\) favors the \(\eta=+1\) domain, while a negative \(h_\eta\) favors the \(\eta=-1\) domain. The engine must now search for a mechanical expression that can play the role of \(h_\eta\). Unlike Example 3, this is a candidate coupling rather than a reproduction of an established formula. We will construct the mechanical carrier before choosing a particular drive.

Engine Input

The carrier needs to describe the domain and the mechanical disturbance. Let \(\mathbf{u}(\mathbf{r},t)\) be the displacement of the lattice from its equilibrium position. Moving the entire crystal by the same amount does not stretch or rotate it, so a uniform displacement cannot distinguish the domains. The first useful information comes from spatial derivatives of \(\mathbf{u}\).

The symmetric part of the displacement gradient is the strain tensor

\[ \varepsilon_{ij} = \frac12 \left( \partial_i u_j+\partial_j u_i \right) \]

Strain measures local stretching and shear. The antisymmetric part describes local rotation. For a time-dependent disturbance, its angular velocity is

\[ \mathbf{\Omega} = \frac12\nabla\times\dot{\mathbf{u}} \]

The domain label, strain, and angular velocity transform differently. Let \(R_g\) be the spatial matrix for a symmetry \(g\), and let \(\chi_g\) be \(+1\) if the symmetry preserves the magnetic domain and \(-1\) if it exchanges the two domains. The three fields transform as

\[ \eta\mapsto\chi_g\eta \qquad \varepsilon\mapsto R_g\varepsilon R_g^{\mathsf T} \qquad \mathbf{\Omega}\mapsto\det(R_g)R_g\mathbf{\Omega} \]

The determinant appears because angular velocity is an axial vector: reflecting space changes a rotation axis differently from an ordinary displacement vector. Time reversal leaves strain unchanged but reverses both \(\eta\) and \(\mathbf{\Omega}\). Therefore, a term containing one factor each of \(\eta\), \(\varepsilon\), and \(\mathbf{\Omega}\) is even under time reversal and is not excluded on that ground.

The local CrSb model uses a \(60^\circ\) rotation in the \(x\)-\(y\) plane, the reflection \(z\mapsto-z\), and the reflection \(y\mapsto-y\). Each operation exchanges the two magnetic sublattices and therefore sends \(\eta\mapsto-\eta\). Let \(\mathcal{G}_{\mathrm{CrSb}}^{(\mathrm{local})}\) denote the group generated by these three actions.

The carrier contains one scalar domain label, six strain components, and three angular-velocity components. In representation notation, the input tuple is

\[ \Pi_4 = \left( \mathcal{G}_{\mathrm{CrSb}}^{(\mathrm{local})}, \mathbb{R}_\eta \oplus \mathrm{Sym}^2(\mathbb{R}^3)_\varepsilon \oplus (\det\otimes\mathbb{R}^3)_\Omega, \tau_{1,1,1}, \varnothing, \varnothing \right) \]

The terms \(\mathbb{R}_\eta\), \(\mathrm{Sym}^2(\mathbb{R}^3)_\varepsilon\), and \((\det\otimes\mathbb{R}^3)_\Omega\) compactly encode the scalar, strain, and axial-vector pieces. The truncation selects terms containing one factor of each

\[ \tau_{1,1,1} \colon \quad \left( \deg_\eta, \deg_\varepsilon, \deg_{\mathbf{\Omega}} \right) = (1,1,1) \]

The two empty slots mean that no equivalence relation or auxiliary identification is imposed in this run. Two control runs select the lower sectors \((1,1,0)\) and \((1,0,1)\) to test strain and angular velocity separately.

Writing \(\mathcal{C}\) for the invariant space, the two lower-order searches return no invariant

\[ \dim\mathcal{C}_{\eta\varepsilon}=0 \qquad \dim\mathcal{C}_{\eta\Omega}=0 \]

The trilinear sector is one-dimensional

\[ \dim\mathcal{C}_{\eta\varepsilon\Omega}=1 \]

The result is not that strain can couple or that rotation can couple, as neither of those can distinguish the domains by itself under these symmetries. The first allowed local mechanical signal is their product.

Derivation

The engine now searches terms containing one factor each of the domain sign \(\eta\), the strain \(\varepsilon\), and the angular velocity \(\mathbf{\Omega}\). The relevant in-plane strain combinations are

\[ Q_1 = \varepsilon_{xx} - \varepsilon_{yy}, \qquad Q_2 = 2\varepsilon_{xy} \]

These combinations record directional stretching and shear. The in-plane angular velocity rotates once with the crystal, while this trace-free in-plane strain rotates twice. Their product therefore rotates three times. Under a \(60^\circ\) crystal rotation, the product rotates through \(180^\circ\) and acquires a minus sign. This matches the sign change of \(\eta\). Two orientations survive the rotation

\[ C_3 = \Omega_xQ_1 - \Omega_yQ_2 \qquad S_3 = \Omega_xQ_2 + \Omega_yQ_1 \]

The horizontal reflection changes the signs of both \(C_3\) and \(S_3\). The vertical reflection changes the sign of \(C_3\) but leaves \(S_3\) unchanged. Since \(\eta\) changes sign under both reflections, only \(\eta C_3\) remains invariant. The unique result is

\[ \boxed{ I_{\mathrm{mech}} = \eta \left[ \Omega_x(\varepsilon_{xx} - \varepsilon_{yy}) - 2\Omega_y\varepsilon_{xy} \right] } \]

The most general coupling in this sector is therefore

\[ \mathcal{L}_{\mathrm{mech}} = \gamma I_{\mathrm{mech}} \]

The symmetry engine finds no coupling between the domain and strain alone, and no coupling between the domain and angular velocity alone. The first allowed mechanical coupling requires both fields simultaneously. Since \(\mathbf{\Omega}\) contains a time derivative, this is a driven coupling rather than a static strain energy.

A surface acoustic wave provides one possible drive. Such a wave carries sound along the surface of a solid. In a Rayleigh-type wave, the atoms move both along the direction of travel and perpendicular to the surface. A phase difference between those motions makes each atom trace an ellipse.

Let the wave travel in the direction

\[ \mathbf{q} = q(\cos\theta,\sin\theta,0) = q\hat{\mathbf{q}} \]

The angle \(\theta\) is measured from the \(x\)-axis, and \(2\pi/q\) is the wavelength. At the surface, use the displacement

\[ \mathbf{u}(\mathbf{r},t) = \operatorname{Re} \left[ \left( U_L\hat{\mathbf{q}} + U_z\hat{\mathbf{z}} \right) e^{i(\mathbf{q}\cdot\mathbf{r} - \omega t)} \right] \]

The complex amplitudes \(U_L\) and \(U_z\) describe motion along the direction of travel and normal to the surface. Their relative phase records the direction and shape of the ellipse. Substituting this displacement into the strain and angular-velocity definitions and averaging over one oscillation gives

\[ \boxed{ \overline{\mathcal{L}}_{\mathrm{mech}} = \frac{\gamma\eta\omega q^2}{4} \operatorname{Im}(U_zU_L^*) \sin 3\theta } \]

The factor \(\operatorname{Im}(U_zU_L^*)\) is nonzero only when the two motions are out of phase, so the atoms trace ellipses rather than moving back and forth along a line. The factor \(\sin3\theta\) fixes the dependence on the direction of travel. If this cycle-averaged term acts as an effective domain bias, its sign depends on the domain, the direction of the ellipse, and the direction in which the wave travels.

Change or limit	Consequence
Put the two motions in phase	The average vanishes
Reverse the direction of the ellipse	The coupling changes sign
Exchange the two magnetic domains	The coupling changes sign
Set \(\theta=n\pi/3\)	The coupling vanishes
Remove the time dependence	\(\mathbf{\Omega}=0\) and the coupling vanishes

The signed coupling has the same \(\sin3\theta\) angular form that appeared in the spin splitting: changing the propagation angle by \(60^\circ\) reverses the bias, while its magnitude repeats every \(60^\circ\). The present expression describes a mechanical coupling rather than an electronic energy. The wave provides a concrete drive only after the engine has found the allowed field combination.

Engine Output and Limits

target sectors:
  linear in eta and strain:                    dimension 0
  linear in eta and angular velocity:          dimension 0
  linear in eta, strain, and angular velocity: dimension 1

returned basis, up to scale:
  eta*(Omega_x*eps_xx - Omega_x*eps_yy - 2*Omega_y*eps_xy)

surface-wave substitution:
  cycle-averaged L_mech =
    (gamma*eta*omega*q**2/4)*Im(Uz*conj(UL))*sin(3*theta)

The certificate records both vanishing control sectors and the one-dimensional target sector. Its single basis element reproduces the coupling derived above, while the surface-wave substitution confirms the angular dependence, zeros, and sign reversals. The later sweeps automate this comparison across larger field libraries.

The returned term is a symmetry-allowed candidate, not yet a demonstrated switching mechanism. Symmetry does not determine whether \(\gamma\) is nonzero or large enough to move a domain wall. A full physical treatment would also need the actual surface symmetry, the wave’s depth profile, integration by parts in the complete action, and a model of the domain dynamics.

The literature status remains provisional. Nearby work studies dynamic strain coupling in altermagnets (Steward et al. 2023), surface-acoustic-wave-driven spin currents in altermagnetic films (Gunnink et al. 2025), magneto-rotation coupling in ferromagnets (Xu et al. 2020), and chiral phonons in CrSb (Rieger et al. 2026). I have not found this exact CrSb coupling in that literature.

Sweeps

Given the above, how do we run the engine such that we can just plug in new symmetries and representations and try to get out new theories?

The engine evaluates one input tuple at a time

\[ \Pi = (\mathcal{G},\ V,\ \tau,\ \sim,\ A) \]

We need some search procedure to decide which of these to try.

Group and Field Libraries

There is no finite library of all symmetry groups. There are complete lists within a declared class: the 32 crystallographic point groups, the 122 magnetic point groups, the 230 space groups, and the 1,651 magnetic space groups. Finite groups can also be enumerated up to a chosen order, while compact Lie groups fall into classified families. In both cases, removing the bound gives an infinite collection.

A group sweep must therefore begin by declaring its universe

\[ \mathfrak{G} = \left\{ \text{groups included in the search} \right\} \]

The crystallographic point groups form one possible \(\mathfrak{G}\). The magnetic point groups form another. A new group can be added once its generators and their actions have been supplied. The engine does not require the group to belong to one of the built-in crystallographic lists.

There is also no library of all physical fields. A symmetry group determines possible representations, but it does not tell us which of them corresponds to an electric field, a strain, an order parameter, or any other quantity in a particular system. That physical interpretation is input.

We can still generate a bounded library of representation types. Let \(V_0=\mathbb{R}^3\) denote the polar-vector representation. For spatial symmetries, the natural primitives include

\[ \mathbf{1}, \qquad \det, \qquad V_0, \qquad \det\otimes V_0 \]

These are the scalar, pseudoscalar, polar-vector, and axial-vector representations. The engine combines them through direct sums and tensor products. Symmetric tensor powers give fields such as strain

\[ S = \operatorname{Sym}^2(V_0) \]

while prolongation adds derivatives through a chosen jet order. Magnetic groups also require a time-reversal parity for each field.

The resulting field library is finite only after we impose bounds. We might limit the tensor rank, the number of fields in a coupling, the polynomial degree, and the derivative order. Let

\[ \mathfrak{F}_{\mathcal{G}}^{(r,d,k)} \]

denote the carriers generated for \(\mathcal{G}\) through tensor rank \(r\), polynomial degree \(d\), and jet order \(k\). The search space is then

\[ \mathcal{S} = \left[ (\mathcal{G},\ V,\ \tau,\ \sim,\ A) \right]_{\substack{ \mathcal{G}\in\mathfrak{G}\\ V\in\mathfrak{F}_{\mathcal{G}}^{(r,d,k)} }} \]

This is not a search over every imaginable theory. It is a complete search within the declared group library, representation grammar, and cutoffs.

The familiar crystal-tensor tables (Nye 1985; Le Quang and He 2011) provide a check. Taking \(\mathfrak{G}\) to be the 32 crystallographic point groups and using the carriers

\[ W_{\mathrm{piezo}} = V_0\otimes\operatorname{Sym}^2(V_0) \]

\[ W_{\mathrm{elastic}} = \operatorname{Sym}^2\left(\operatorname{Sym}^2(V_0)\right) \]

\[ W_{\mathrm{flexo}} = V_0\otimes\operatorname{Sym}^2(V_0)\otimes V_0 \]

reproduces the 106 benchmark entries used by the engine. This sweep validates the library and the induced actions. It is not a discovery.

The search begins when the same construction is applied to response spaces that have not yet been checked against a complete table. Running over the 122 magnetic point groups requires time-reversal-labeled representations and produces candidate classifications for altermagnetic harmonics, magneto-odd elasticity, flexomagnetism, and magnetic second-harmonic generation. These calculations are complete within their declared bounds, while their comparison with the specialty literature remains pending.

We can instead fix one group and vary the field library. For CrSb, the implemented fields include the domain sign \(\eta\), strain \(\varepsilon\), lattice angular velocity \(\mathbf{\Omega}\), electric and magnetic fields, current, and the temperature gradient

\[ \mathfrak{F}_{\mathrm{CrSb}} = \left\{ \eta,\varepsilon,\mathbf{\Omega},E,B,j,\nabla T \right\} \]

If the search requires exactly one factor of \(\eta\) and allows total degree through three, it generates sectors of the form

\[ \eta \qquad \eta F_i \qquad \eta F_iF_j \]

with \(F_i\) and \(F_j\) drawn from the remaining field library. The vanishing of the \(\eta\varepsilon\) and \(\eta\mathbf{\Omega}\) sectors, and the one-dimensional \(\eta\varepsilon\mathbf{\Omega}\) sector from Example 4, are then outputs of the search rather than cases chosen afterward.

The group and field libraries define what the engine can search mechanically. New primitives or new group actions enlarge that search space, while the cutoffs keep it finite.

LLM-Based Proposals

Another option is to use LLMs to propose groups, field sectors, or truncations. The invariant calculation then proceeds, and the output can be passed to other agents for literature search.

Conclusion

Based on symmetries acting on a given representation, we can compute the allowable theories with a mechanical framework. This lets us put a large amount of results from mathematical physics into a common, automatable language. Furthermore, this lets us search the space of theories; in fact, I’ve already run a fair number of sweeps already using this system⁹.

This does not “solve physics,” nor is it meant to, as symmetry cannot supply coefficients, microscopic physics, filling rules, or novelty. That being said, this could be a good tool to reduce the need to rely on the literature and to work in conjunction with LLMs.

Beyond these examples, I have already implemented extensions over symmetry-breaking lattices, Landau transitions, empirical strata, indistinguishable group pairs, Lie-algebra landscapes, conservation-law drift, stability, persistence, and graded spurion expansions, together with a cohomology layer for extensions and anomalies. These need their own checks and will be the subjects of future posts, but the point is that this system can be extended to catalog and search a wider swathe of ideas in mathematical physics.

Appendix: Sweep Inventory

Forthcoming.

References

Gergely, László Á. 2002. “On Hamiltonian Formulations of the Schrödinger System.” Annals of Physics 298 (2): 394–402. https://doi.org/10.1006/aphy.2002.6262.

Gunnink, Pieter M., Jairo Sinova, and Alexander Mook. 2025. “Surface Acoustic Wave Driven Acoustic Spin Splitter in \(d\)-Wave Altermagnetic Thin Films.” arXiv Preprint arXiv:2502.18007, ahead of print. https://doi.org/10.48550/arXiv.2502.18007.

Jungwirth, Tomáš, Jairo Sinova, Peter Wadley, et al. 2026. “Altermagnetic Spintronics.” Nature Physics, ahead of print. https://doi.org/10.1038/s41567-026-03337-w.

Landau, L. D. 1937. “On the Theory of Phase Transitions.” Zh. Eksp. Teor. Fiz. 7: 19–32.

Le Quang, H., and Q.-C. He. 2011. “The Number and Types of All Possible Rotational Symmetries for Flexoelectric Tensors.” Proceedings of the Royal Society A 467: 2369–86. https://doi.org/10.1098/rspa.2010.0521.

Leutwyler, H. 1994. “Nonrelativistic Effective Lagrangians.” Physical Review D 49 (6): 3033–43. https://doi.org/10.1103/PhysRevD.49.3033.

Nye, J. F. 1985. Physical Properties of Crystals: Their Representation by Tensors and Matrices. Oxford University Press.

Olver, Peter J. 1993. Applications of Lie Groups to Differential Equations. 2nd ed. Vol. 107. Graduate Texts in Mathematics. Springer. https://doi.org/10.1007/978-1-4612-4350-2.

Reimers, Sonka, Lukas Odenbreit, Libor Šmejkal, et al. 2024. “Direct Observation of Altermagnetic Band Splitting in CrSb Thin Films.” Nature Communications 15: 2116. https://doi.org/10.1038/s41467-024-46476-5.

Rieger, Philipp, Markus Weißenhofer, Sergiy Mankovsky, and Peter M. Oppeneer. 2026. “Symmetry-Enforced Chiral Phonons in Altermagnets via Magnon-Phonon Coupling.” arXiv Preprint arXiv:2607.06792, ahead of print. https://doi.org/10.48550/arXiv.2607.06792.

Song, Cheng, Hua Bai, Zhiyuan Zhou, et al. 2025. “Altermagnets as a New Class of Functional Materials.” Nature Reviews Materials 10: 473–85. https://doi.org/10.1038/s41578-025-00779-1.

Steward, Charles R. W., Rafael M. Fernandes, and Jörg Schmalian. 2023. “Dynamic Paramagnon-Polarons in Altermagnets.” Physical Review B 108: 144418. https://doi.org/10.1103/PhysRevB.108.144418.

Weinberg, Steven. 1995. The Quantum Theory of Fields, Volume 1: Foundations. Cambridge University Press. https://doi.org/10.1017/CBO9781139644167.

Xu, Mingran, Kei Yamamoto, Jorge Puebla, et al. 2020. “Nonreciprocal Surface Acoustic Wave Propagation via Magneto-Rotation Coupling.” Science Advances 6: eabb1724. https://doi.org/10.1126/sciadv.abb1724.

Yang, Guowei, Zhanghuan Li, Sai Yang, et al. 2025. “Three-Dimensional Mapping of the Altermagnetic Spin Splitting in CrSb.” Nature Communications 16: 1442. https://doi.org/10.1038/s41467-025-56647-7.

Zhou, Zhiyuan, Xingkai Cheng, Mengli Hu, et al. 2025. “Manipulation of the Altermagnetic Order in CrSb via Crystal Symmetry.” Nature 638: 645–50. https://doi.org/10.1038/s41586-024-08436-3.

AI Disclosure

I used AI to compile this post from notes and edit it. AI helped draft Examples 3,4.

Footnotes

1. Essentially the smallest nontrivial example.

2. Often described as the path of least action, although technically the physical path is an extremal, or stationary point, of the action. For general systems (not necessarily physics) I like to think of the Lagrangian as the cost density over a system trajectory.

3. The invariance steps below are pure invariant theory and do not themselves assume a variational principle. A symmetry that leaves \(L\) unchanged forces it to be a function of the invariants, whatever \(L\) means. Least action enters when we allow a symmetry to change \(L\) by a total derivative, since a boundary term shifts every path’s action by an endpoint contribution and so cannot move its extremals.

4. Nothing is lost by fixing the degree. A linear group action sends degree-\(d\) monomials to degree-\(d\) polynomials, so the invariance conditions never mix degrees, and the invariants of an arbitrary polynomial are the invariants of its homogeneous pieces. Each degree is its own finite linear-algebra problem.

5. Both of these use the same rule, which is to write every local term allowed by symmetry, then truncate by order. Landau theory applies the rule to the free energy near a phase transition (Landau 1937). Weinberg’s Folk Theorem applies it to a Lagrangian consistent with the assumed symmetries, unitarity, and cluster decomposition (Weinberg 1995). It is “folk” because it is a motivated principle rather than a proved theorem, although rejecting it would overturn the usual practice of effective field theory. The claims we have here are local. Global and topological terms need additional data and belong to a separate cohomological layer in this framework, which is developed but out of scope for this write-up.

6. Prolongation and differential invariants are treated systematically by Olver (1993).

7. The Lagrangian-to-Noether half of this machinery is developed in a different setting elsewhere on this site, in a Lagrangian for games and approximate Noether theorems.

8. The derivation follows the logic of nonrelativistic effective field theory (Leutwyler 1994) and the standard variational formulation of the Schrödinger field (Gergely 2002). The claim recorded in the catalog is a reproduction.

9. Results pending. I’ve run a number of sweeps and there are some candidates but of unknown (likely low) impact. It seems that the field of physics is “efficient” in that most of the low hanging fruit has been done. That being said, a lot of the techniques here are reasonably old, so likely extensions to this engine are required to put us in novel territory.