Introduction

Dealing with involved Hamiltonians and other operators from the quantum many-body theory in C++ can be tough. In some cases, quantum-mechanical operators can be uniformly represented by finite-dimensional matrices. This representation, however, is often impractical. Storing the matrices can quickly prove infeasible as the amount of required memory grows exponentially with the number of degrees of freedom. For this reason, many computational programs in the field use sparse matrices or hard-coded procedures that describe how said operators act on quantum states. These implementations usually accept a few Hamiltonian parameters as input, but switching to a more general form/adding more terms to the Hamiltonian requires a considerable code rewrite.

The goal of libcommute’s Domain-Specific Language (DSL) is to streamline this coding task. It introduces an abstraction of the polynomial quantum-mechanical operator expression, which can be manipulated as easily as an equation written on a piece of paper.

As a primer, let us consider the following simple program that constructs Hamiltonian of an electronic tight-binding model on a square \(10\times 10\) lattice with only nearest-neighbour hopping allowed,

\[\hat H_\text{e} = -t \sum_\sigma \sum_{\langle i,j\rangle} c^\dagger_{i,\sigma} c_{j,\sigma}.\]

#include <string>

int main() {

  //
  // Let us define Hamiltonian of an electronic tight-binding model
  // on a square lattice.
  //

  // Number of lattice sites in each direction
  // (the total number of sites is N * N)
  int const N = 10;

  // Electron hopping constant - energy parameter of the TB model
  double const t = 2.0;

  // TB Hamiltonian as an expression with real coefficients and
  // statically-typed indices. In this case, the indices are a pair of
  // integers - lattice site coordinates - and a spin label ("up" or "down").
  libcommute::static_indices::expr_real<int, int, std::string> H_e;

  // Use functions c_dag() and c() that return fermionic creation/annihilation
  // operators.
  using libcommute::static_indices::c_dag;
  using libcommute::static_indices::c;

  // Are two sites neighbors along an axis with periodicity?
  auto neighbors = [](int i, int j) {
    return std::abs(i - j) == 1 || std::abs(i - j) == N;
  };

  // Iterate over spin projections
  for(auto spin : {"up", "down"}) {
    // Iterate over all lattice sites with coordinates i = (ix, iy)
    for(int ix = 0; ix < N; ++ix) {
      for(int iy = 0; iy < N; ++iy) {
        // Iterate over all lattice sites with coordinates j = (jx, jy)
        for(int jx = 0; jx < N; ++jx) {
          for(int jy = 0; jy < N; ++jy) {
            // Skip all pairs of lattice sites i and j that are not
            // nearest-neighbors. The modulus operation accounts for
            // periodic boundary conditions on the lattice.
            if((neighbors(ix, jx) && iy == jy) ||
               (ix == jx && neighbors(iy, jy))) {
              // Add a hopping term
              H_e += -t * c_dag(ix, iy, spin) * c(jx, jy, spin);
            }
          }
        }
      }
    }
  }

Now, let us add a harmonic oscillator at each lattice site (a localized phonon),

\[\hat H_\text{ph} = \omega_0 \sum_i a^\dagger_i a_i.\]

  // Frequency of the localized phonon
  double const w0 = 0.5;

  //
  // Hamiltonian of phonons localized at lattice sites.
  // We want this object to have the same type as H_e.
  //
  decltype(H_e) H_ph;

  // Use functions a_dag() and a() that return bosonic creation/annihilation
  // operators.
  using libcommute::static_indices::a_dag;
  using libcommute::static_indices::a;

  // Iterate over all lattice sites
  for(int ix = 0; ix < N; ++ix) {
    for(int iy = 0; iy < N; ++iy) {
      // Energy of the localized phonon at site (ix, iy)
      H_ph += w0 * a_dag(ix, iy, "") * a(ix, iy, "");
    }
  }

Note

We had to assign an empty spin label “” to the bosons, because all operators in H_ph have to carry exactly three indices with the last one being a string. It is possible to overcome this limitation and put just two integer indices on \(a^\dagger\)/\(a\) by switching to the dynamically-typed indices. Be aware, however, that the dynamic indices require C++17 and may result in less type-safe code.

Finally, we are going to couple electrons with phonons and arrive at Hamiltonian of the Holstein model \(\hat H_H\),

\[\hat H_\text{e-ph} = g \sum_\sigma \sum_i n_{i,\sigma}(a^\dagger_i + a_i),\]

\[\hat H_H = \hat H_\text{e} + \hat H_\text{ph} + \hat H_\text{e-ph}.\]

  // Electron-phonon coupling constant
  double const g = 0.1;

  //
  // Hamiltonian of electron-phonon coupling.
  //
  decltype(H_e) H_e_ph;

  // Use function n() that returns the fermionic number operator n = c_dag * c
  using libcommute::static_indices::n;

  // Iterate over spin projections
  for(auto spin : {"up", "down"}) {
    // Iterate over all lattice sites
    for(int ix = 0; ix < N; ++ix) {
      for(int iy = 0; iy < N; ++iy) {
        // Electron-phonon coupling at site (ix, iy)
        H_e_ph += g * n(ix, iy, spin) * (a_dag(ix, iy, "") + a(ix, iy, ""));
      }
    }
  }

  // Holstein Hamiltonian.
  auto H_H = H_e + H_ph + H_e_ph;

  // Print H_H. There will be quite a lot of terms for the 100-site lattice!
  std::cout << "H_H = " << H_H << '\n';

  // Check hermiticity of H_H
  std::cout << "H_H - H_H^\\dagger = " << (H_H - conj(H_H)) << '\n';

  // Check that H_H commutes with the total number of electrons
  decltype(H_H) N_e;
  for(auto spin : {"up", "down"}) {
    for(int ix = 0; ix < N; ++ix) {
      for(int iy = 0; iy < N; ++iy) {
        N_e += n(ix, iy, spin);
      }
    }
  }

  std::cout << "[H_H, N_e] = " << (H_H * N_e - N_e * H_H) << '\n';

  return 0;
}

Now that we have the complete Hamiltonian object, we could proceed along one of the following routes.

Make a loperator object out of \(\hat H_H\) and use it to act on state vectors in a finite-dimensional Hilbert space. This is a common step in implementing Exact Diagonalization algorithms.
Use the iteration interface to analyze the structure of the Hamiltonian term by term.
Transform the Hamiltonian by applying a function to each term.