Introducing a new algebra is as easy as deriving a class from the abstract base libcommute::generator. libcommute’s DSL can work with algebras, whose generators $$g_\alpha$$ obey commutation relations

$g_\alpha g_\beta - c g_\beta g_\alpha = F_{\alpha\beta} + \sum_\gamma f_{\alpha\beta}^\gamma g_\gamma$

with real constants $$c$$, $$F_{\alpha\beta}$$ and $$f_{\alpha\beta}^\gamma$$. Such algebraic structures include the Lie and Clifford algebras.

In the following example we will define the algebra of Dirac gamma matrices $${\gamma^0,\gamma^1,\gamma^2,\gamma^3}$$. The generators will carry one integer index (0, 1, 2 or 3), which will also fix their canonical ordering. The identities we will be using in the code below are

• Anti-commutation relation $$\{\gamma^\mu, \gamma^\nu\} = 2\eta^{\mu\nu} \hat I_4$$, where $$\eta^{\mu\nu}$$ is the Minkowski metric tensor with signature $$(+,-,-,-)$$ and $$\hat I_4$$ is the four-dimensional identity matrix.

• Squares of gamma matrices,

$\begin{split}(\gamma^0)^2 &= \hat I_4,\\ (\gamma^k)^2 &= -\hat I_4 \text{ for } k=1,2,3,\\\end{split}$
• Hermitian conjugates

$\begin{split}(\gamma^0)^\dagger &= \gamma^0,\\ (\gamma^k)^\dagger &= -\gamma^k.\end{split}$
#include <libcommute/expression/generator.hpp>

// For min_user_defined_algebra_id
#include <libcommute/algebra_ids.hpp>

#include <complex>
#include <iostream>

namespace libcommute {

//
// First, we define an ID for our new algebra
//

// Use the lowest algebra ID available to user-defined algebras
static constexpr int gamma = min_user_defined_algebra_id;

//
// Our generator type: a gamma matrix with one integer index
//
class generator_gamma : public generator<int> {

using base = generator<int>;
using linear_function_t = typename base::linear_function_t;

public:
// Algebra ID of this generator
int algebra_id() const override { return libcommute::gamma; }

// Constructor: Just pass the index to the base class
explicit generator_gamma(int index) : base(index) {}

// Virtual copy-constructor.
// Make a smart pointer that manages a copy of this generator
std::unique_ptr<base> clone() const override {
// With C++14 or newer, libcommute::make_unique() will be resolved to
// std::make_unique(). Otherwise, a custom implementation will be used.
return make_unique<generator_gamma>(*this);
}

// This function will be called for g1 = *this and g2 such that g1 > g2.
// We must transform the product g1 * g2 and put it into the
// canonical order,
//
// g1 * g2 -> -g2 * g1 + 2\eta(g1, g2)
double swap_with(base const& g2, linear_function_t& f) const override {

// Do g1 and g2 have the same indices?
bool diag = base::equal(g2);

// Minkowski metric
int index = std::get<0>(base::indices());
double eta = diag ? (index == 0 ? 1 : -1) : 0;

// Set f(g) to be the constant 2\eta(g1, g2)
f.set(2 * eta);

// Return coefficient in front of g2 * g1 in the transformed expression
return -1;
}

// This function will be called for g1 = *this and g2 that are already
// canonically ordered, g1 <= g2.
//
// It tries to simplify squares of gamma matrices,
// (\gamma^0)^2 = I_4
// (\gamma^k)^2 = -I_4 for k = 1,2,3
bool simplify_prod(base const& g2, linear_function_t& f) const override {
if(*this == g2) {
// Replace the square with a constant
int index = std::get<0>(base::indices());
f.set(index == 0 ? 1 : -1);
return true;
} else
// Not a square, cannot simplify
return false;
}

// Hermitian conjugate of this generator as a linear function of generators.
// \gamma^0 is Hermitian and \gamma^k are anti-Hermitian
void conj(linear_function_t& f) const override {
int index = std::get<0>(base::indices());
if(index == 0) {
// f(g) = 0 + 1*(*this)
f.set(0, clone(), 1);
} else {
// f(g) = 0 + (-1)*(*this)
f.set(0, clone(), -1);
}
}

// Stream output function
std::ostream& print(std::ostream& os) const override {
int index = std::get<0>(base::indices());
return os << "\\gamma^" << index;
}
};

} // namespace libcommute


It is usually worth defining a factory function that creates an expression containing one generator with a unity prefactor.

#include <libcommute/expression/expression.hpp>

namespace libcommute {

//
// Define a factory function to simplify construction of expressions
//

expression<std::complex<double>, int> make_gamma(int index) {
using expr_t = expression<std::complex<double>, int>;
// Return an expression made of one monomial with coefficient = 1.
return expr_t(1.0, expr_t::monomial_t(generator_gamma(index)));
}

} // namespace libcommute


Now we can check that generators of our algebra actually fulfil the canonical anti-commutation relations.

  using namespace libcommute;

// Minkowski metric tensor
auto eta = [](int mu, int nu) -> double {
if(mu != nu) return 0;
return mu == 0 ? 1.0 : -1.0;
};

// Check the commutation relations
for(int mu = 0; mu < 4; ++mu) {
for(int nu = 0; nu < 4; ++nu) {
auto gamma_mu = make_gamma(mu);
auto gamma_nu = make_gamma(nu);

std::cout << "{" << gamma_mu << ", " << gamma_nu << "}"
<< " - 2\\eta(" << mu << ", " << nu << "} = "
<< (gamma_mu * gamma_nu + gamma_nu * gamma_mu - 2 * eta(mu, nu))
<< '\n';
}
}


Let us also define the fifth gamma matrix

$\gamma^5 = i \gamma^0 \gamma^1 \gamma^2 \gamma^3$

and check that

$\begin{split}(\gamma^5)^\dagger = \gamma^5,\\ \{\gamma^5, \gamma^\mu\} = 0.\end{split}$
  // \gamma^5
std::complex<double> const I(0, 1);
auto gamma5 =
I * make_gamma(0) * make_gamma(1) * make_gamma(2) * make_gamma(3);

// \gamma^5 is Hermitian ...
std::cout << "gamma5 - conj(gamma5) = " << (gamma5 - conj(gamma5)) << '\n';

// ... and anti-commutes with \gamma^\mu.
for(int mu = 0; mu < 4; ++mu) {
auto gamma_mu = make_gamma(mu);
std::cout << "{gamma5, " << gamma_mu
<< "} = " << (gamma5 * gamma_mu + gamma_mu * gamma5) << '\n';
}
}