Advanced: A user-defined algebra
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 the commutation relations
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) {}
// Standard constructors and destructor
generator_gamma(generator_gamma const&) = default;
generator_gamma(generator_gamma&&) noexcept = default;
~generator_gamma() override = default;
// Generator objects must be immutable
generator_gamma& operator=(generator_gamma const&) = delete;
generator_gamma& operator=(generator_gamma&&) noexcept = delete;
// 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)
var_number swap_with(base const& g2, linear_function_t& f) const override {
// Do g1 and g2 have the same index?
bool diag = base::equal(g2);
// Minkowski metric
int index = std::get<0>(base::indices());
int eta = diag ? (index == 0 ? 1 : -1) : 0;
// Set f(g) to be the constant 2\eta(g1, g2)
f.set(2 * eta);
// Return the 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());
// Here, we call base::shared_from_this() to obtain a valid std::shared_ptr
// that co-owns *this. This approach allows to save memory, and is strongly
// preferred to making a copy of *this wrapped in a new shared pointer and
// passing it to f.set().
// For more details on shared_from_this() see
// https://en.cppreference.com/w/cpp/memory/enable_shared_from_this.html
if(index == 0) {
// f(g) = 0 + 1*(*this)
f.set(0, shared_from_this(), 1);
} else {
// f(g) = 0 + (-1)*(*this)
f.set(0, shared_from_this(), -1);
}
}
// String representation of this \gamma-matrix
std::string to_string() const override {
int index = std::get<0>(base::indices());
return "γ^" + std::to_string(index);
}
};
} // namespace libcommute
Note
If the constants \(c\), \(F_{\alpha\beta}\),
\(f_{\alpha\beta}^\gamma\) belong to a certain category of numbers
(integers, rationals or general real values), then the generator is compatible
with all ScalarType’s belonging to
the same or wider category.
libcommute::linear_function internally stores its coefficients as
instances of the variant type libcommute::var_number. In the
particular example presented here, all coefficients are integers, which are
automatically converted to libcommute::var_number in the expressions
similar to f.set(2 * eta). In the case of the rational coefficients,
it is recommended to explicitly construct the values by calling
var_number(numerator, denominator) instead of just passing floating-point
numbers. This way, the defined generators can be used with a rational
ScalarType.
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η(" << mu << ", " << nu << ") = "
<< (gamma_mu * gamma_nu + gamma_nu * gamma_mu - 2 * eta(mu, nu))
<< '\n';
}
}
Let us also define the fifth gamma matrix
and check that
// \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 << "γ^5 - conj(γ^5) = " << (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 << "{γ^5, " << gamma_mu
<< "} = " << (gamma5 * gamma_mu + gamma_mu * gamma5) << '\n';
}