Advanced: Linear operator representation of a user-defined algebra

Having defined a new algebra one can take a step further and implement the linear operator representation of it. Let us build on the example where we introduced the algebra of 4-dimensional \(\gamma\)-matrices. This page will outline the extra steps needed to plug the new algebra into the linear operator framework.

  • Declare a new elementary space class. The elementary space must share algebra ID with generator_gamma. Its dimension is 4, so the number of occupied bits is \(b = 2\).

    #include "gamma.hpp" // Definition of 'gamma' algebra ID
    #include <libcommute/loperator/elementary_space.hpp>
    #include <libcommute/loperator/loperator.hpp>
    #include <libcommute/loperator/monomial_action.hpp>
    #include <algorithm>
    #include <array>
    #include <vector>
    namespace libcommute {
    // 4-dimensional elementary space gamma matrices act in
    class elementary_space_gamma : public elementary_space<int> {
      using base = elementary_space<int>;
      // Since all 4 gamma matrices act in the same elementary space,
      // we can initialize the base class with any number
      elementary_space_gamma() : base(0) {}
      elementary_space_gamma(elementary_space_gamma const&) = default;
      elementary_space_gamma(elementary_space_gamma&&) noexcept = default;
      elementary_space_gamma& operator=(elementary_space_gamma const&) = default;
      operator=(elementary_space_gamma&&) noexcept = default;
      ~elementary_space_gamma() override = default;
      // Virtual copy-constructor.
      // Make a smart pointer that manages a copy of this elementary space
      std::unique_ptr<base> clone() const override {
        return make_unique<elementary_space_gamma>(*this);
      // Algebra ID, must be the same with generator_gamma::algebra_id()
      int algebra_id() const override { return libcommute::gamma; }
      // We need 2 bits to enumerate all 4 = 2^2 basis vectors
      int n_bits() const override { return 2; }
    } // namespace libcommute
  • Specialize class template monomial_action for the new algebra. This specialization will describe how a product \(\gamma^\mu \gamma^\nu \gamma^\kappa \ldots\) acts on a single basis state \(|n\rangle, n = \{0,1,2,3\}\).

    // Action of a product of gamma matrices (a gamma-matrix monomial)
    // on a basis vector.
    namespace libcommute {
    template <> class monomial_action<libcommute::gamma> {
      // Indices of matrices in the product, right to left.
      // This order is chosen because the matrix on the right acts on a state first
      std::vector<int> index_sequence;
      // m_range is a pair of iterators over a list of generator_gamma objects.
      // This range represents the product of gamma matrices we want to act with.
      monomial_action(monomial<int>::range_type const& m_range,
                      hilbert_space<int> const& hs) {
        // Iterate over matrices in the product, left to right
        for(auto it = m_range.first; it != m_range.second; ++it)
          // Collect indices of the matrices
        // Reverse the order
        std::reverse(index_sequence.begin(), index_sequence.end());
      // Act on a basis state
      // In the present implementation, libcommute requires all algebra generators
      // to be represented by generalized permutation matrices [1], i.e. matrices
      // with exactly one non-zero element in each row and each column.
      // Equivalently, any generator acting on a basis state must give a single
      // basis state multiplied by a constant.
      // [1]
      // 'index' is the index (a 64-bit unsigned integer) of the basis state we are
      // acting upon. It must also receive the index of the resulting basis state.
      // 'coeff' must be multiplied by the overall constant factor acquired as
      // a result of monomial action.
      template <typename ScalarType>
      inline bool act(sv_index_type& index, ScalarType& coeff) const {
        std::complex<double> const I(0, 1);
        // Act with all matrices in the product, right to left
        for(int i : index_sequence) {
          switch(i) {
          case 0:
            // Action of \gamma^0
            // It is diagonal => 'index' does not change
            coeff *= std::array<double, 4>{1, 1, -1, -1}[index];
          case 1:
            // Action of \gamma^1
            coeff *= std::array<double, 4>{-1, -1, 1, 1}[index];
            index = std::array<sv_index_type, 4>{3, 2, 1, 0}[index];
          case 2:
            // Action of \gamma^2
            coeff *= std::array<std::complex<double>, 4>{-I, I, I, -I}[index];
            index = std::array<sv_index_type, 4>{3, 2, 1, 0}[index];
          case 3:
            // Action of \gamma^3
            coeff *= std::array<std::complex<double>, 4>{-1, 1, 1, -1}[index];
            index = std::array<sv_index_type, 4>{2, 3, 0, 1}[index];
        // This 'true' signals that the action result is not identically zero.
        // Returning 'false' in cases when it is zero would help improve
        // performance.
        return true;
    } // namespace libcommute

    In general, monomial_action is defined as follows.

    template<int AlgebraID>
    class monomial_action

    Defined in <libcommute/loperator/monomial_action.hpp>

    Action of a product of generators belonging to the same algebra AlgebraID on a basis vector.

    template<typename ...IndexTypes>
    monomial_action(monomial<IndexTypes...>::range_type const &m_range, hilbert_space<IndexTypes...> const &hs)

    m_range is a pair of monomial::const_iterator’s. This range represents the product of generators one wants to act with. hs is the Hilbert space to act in.

    template<typename ScalarType>
    bool act(sv_index_type &index, ScalarType &coeff) const

    Act on a basis state. It is assumed the action of a generator (and, therefore, of a monomial) maps a single basis state to a single basis state multiplied by a constant.

    index is the index of the input and output basis states.

    coeff must be multiplied by the overall constant factor acquired as a result of monomial action.

    This method should return false if the action result is the identical zero and true otherwise.

  • Make a linear operator object with the new algebra ID and apply it to 4-dimensional state vectors as usual.

    int main() {
      using namespace libcommute;
      // Hilbert space made of one elementary space for gamma matrices.
      hilbert_space<int> hs{elementary_space_gamma()};
      std::complex<double> const I(0, 1);
      // Expression for \gamma^5
      auto gamma5 =
          I * make_gamma(0) * make_gamma(1) * make_gamma(2) * make_gamma(3);
      // Linear operator representation of \gamma^5
      auto gamma5op =
          loperator<std::complex<double>, libcommute::gamma>(gamma5, hs);
      // Build the explicit matrix representation of \gamma^5
      // Preallocate state vectors.
      std::vector<std::complex<double>> psi(4);
      std::vector<std::complex<double>> phi(4);
      // Iterate over basis states
      for(int i = 0; i < 4; ++i) {
        // Reset vectors |\psi> and |\phi> to zero.
        std::fill(psi.begin(), psi.end(), 0);
        std::fill(phi.begin(), phi.end(), 0);
        psi[i] = 1; // 'psi' is i-th basis vector now
        phi = gamma5op(psi);
        // Print the result
        std::cout << "gamma5|" << i << "> = ";
        for(int j = 0; j < 4; ++j) {
          std::cout << "+" << phi[j] << "|" << j << ">";
        std::cout << '\n';
      return 0;