Normal Ordering: Making Physics Legal

Here’s the central problem: quantum chemistry gives us operator expressions in arbitrary order, but the physics only works if we put them in normal order — all creation operators ($a^\dagger$) before all annihilation operators ($a$).

Reordering isn’t free. Fermions obey the canonical anti-commutation relations (CAR):

This means:

1. Swapping two operators flips the sign (the “anti” in anti-commutation)
2. Swapping $a_i$ past $a^\dagger_i$ (same index) generates an identity term

Let’s see this in action with a simple example — $a_0 a^\dagger_1$ (annihilate before create, wrong order):

let disordered =
    P<IxOp<uint32, LadderOperatorUnit>>.Apply [|
        LadderOperatorUnit.FromUnit(false, 0u)   // a₀  (Lower first — wrong!)
        LadderOperatorUnit.FromUnit(true, 1u)     // a†₁ (Raise second)
    |]
    |> S<IxOp<uint32, LadderOperatorUnit>>.Apply

// Normal-order it:
let ordered =
    LadderOperatorSumExpr<FermionicAlgebra>.ConstructNormalOrdered disordered

// Result: −1 · a†₁ a₀
// The sign flipped because we swapped a fermion past another!
// No identity term because the indices differ (0 ≠ 1).

Now watch what happens with same-index operators — $a_0 a^\dagger_0$:

let sameIndex =
    P<IxOp<uint32, LadderOperatorUnit>>.Apply [|
        LadderOperatorUnit.FromUnit(false, 0u)    // a₀
        LadderOperatorUnit.FromUnit(true, 0u)     // a†₀
    |]
    |> S<IxOp<uint32, LadderOperatorUnit>>.Apply

let result = LadderOperatorSumExpr<FermionicAlgebra>.ConstructNormalOrdered sameIndex
// Result: 1 − a†₀ a₀
// TWO terms! The identity (δ₀₀ = 1) plus the reordered product with a sign flip.

This is the CAR in action: $a_0 a^\dagger_0 = \delta_{00} - a^\dagger_0 a_0 = 1 - a^\dagger_0 a_0$.

The algebra is pluggable

The key design insight: commutation relations are not hard-coded. They live behind an interface:

type ICombiningAlgebra<'op> =
    abstract Combine :
        P<IxOp<uint32,'op>> -> C<IxOp<uint32,'op>> -> P<IxOp<uint32,'op>>[]

The Combine method takes a partial product and the next operator, returning all terms that result from appending — including any terms generated by commutation relations.

FockMap ships two implementations:

The generic LadderOperatorSumExpr<'algebra> type accepts either. The compiler enforces that you don’t accidentally mix fermionic and bosonic algebra within the same expression.

How the pipeline works

The normal ordering pipeline proceeds in stages:

1. Take the input S<IxOp<uint32, LadderOperatorUnit>> expression.
2. For each product term, repeatedly combine adjacent operators via the algebra’s Combine method, collecting all generated terms.
3. Sort indices within raise/lower groups for canonical form.
4. Return the normalised S<...>.

Because S<'T> keys its terms by string representation, like terms are automatically combined at every step — coefficients sum, duplicates vanish.

Bosonic normal ordering

Compare the fermionic result above with the bosonic version of $b_0 b^\dagger_0$:

let bosonicExpr =
    P<IxOp<uint32, LadderOperatorUnit>>.Apply [|
        LadderOperatorUnit.FromUnit(false, 0u)   // b₀
        LadderOperatorUnit.FromUnit(true, 0u)     // b†₀
    |]
    |> S<IxOp<uint32, LadderOperatorUnit>>.Apply

let bosonicResult = constructBosonicNormalOrdered bosonicExpr
// Result: 1 + b†₀ b₀
// Note the PLUS sign! Bosons commute — no sign flip.

Fermionic: $a_0 a^\dagger_0 = 1 - a^\dagger_0 a_0$. Bosonic: $b_0 b^\dagger_0 = 1 + b^\dagger_0 b_0$. That single sign difference is the entire distinction between matter and light.

Next: Your First Encoding — mapping ladder operators to Pauli strings