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>;

public:
// 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;
elementary_space_gamma&
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;

public:
// 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
index_sequence.push_back(std::get<0>(it->indices()));
// 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] https://en.wikipedia.org/wiki/Generalized_permutation_matrix
//
// '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];
break;
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];
break;
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];
break;
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];
break;
}
}
// 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;
}