cppqed-multidiagonal

A header-only C++ library implementing sparse multidiagonal operator algebra for quantum systems of arbitrary arity, extracted from the C++QED framework.

What is a multidiagonal operator?

A multidiagonal operator over a Hilbert space of dimension $N$ is a sparse matrix storing an arbitrary set of diagonals at arbitrary offsets. In the context of quantum optics and cavity QED, Hamiltonians and jump operators naturally take this form. For example, the harmonic oscillator ladder operators $a$, $a^\dagger$ are single off-diagonals; number operators $a^\dagger a$ are main diagonals, and so is the Kerr term $(a^\dagger)^2 a^2 = a^\dagger a,(a^\dagger a - 1)$; two-photon terms like $a^2$ sit two diagonals off the main one.

For a composite system of $R$ subsystems, MultiDiagonal<RANK> stores an operator of the form

$$A = \sum_{\mathbf{o} \in \mathbb{K}} \sum_{\mathbf{n}} \alpha^{\mathbf{o}}_{\mathbf{n}} \left|\mathbf{n}\right\rangle\left\langle \mathbf{n}+\mathbf{o}\right|$$

with $\mathbb{K} \subset \mathbb{Z}^R$ an arbitrary set of signed offset vectors, one entry per subsystem axis, and $\mathbf{n}$ running over the multi-indices for which both $\mathbf{n}$ and $\mathbf{n}+\mathbf{o}$ lie in the truncated Hilbert space. A positive entry thus means an upper diagonal on that axis. The coefficients need not factorize over the subsystems; direct products, built with operator%, are the special case where they do.

The signed-offset representation

Each diagonal is keyed by a signed offset per axis — positive = upper, negative = lower, zero = main diagonal. The offset counts the quanta a diagonal moves along each axis, and the algebra reduces to arithmetic on these counts:

• composition: $\text{offset}_{AB} = \text{offset}_A + \text{offset}_B$ per axis, for arbitrary rank
• Hermitian conjugation: negation of all offsets (plus conjugation of the stored values)
• direct product: concatenation of the offset arrays, reflecting the tensor product structure exactly

The diagonal at signed offset $o$ (per axis, dimension $D$) has extent $D-|o|$, and entry $k$ holds the matrix element $A_{k+o^-,,k+o^+}$ with $o^\pm=\max(\pm o,0)$. This min(row,col)-anchored storage is symmetric under transposition, which is what makes conjugation a pure key negation. Application to a state vector costs $O(N \cdot d)$ for $d$ diagonals; composition runs over the $d_A \cdot d_B$ diagonal pairs, each in one pass over the window where the two diagonals overlap.

Why not Tridiagonal?

The predecessor class Tridiagonal<RANK> in C++QED v2 stored exactly three diagonals per axis — at positions $-K$, $0$, $+K$ for a single shared offset $K$ — and was not closed under composition: the product of two tridiagonal operators is in general pentadiagonal, and with different offsets is not tridiagonal at all. MultiDiagonal lifts all these restrictions — any number of diagonals at independent arbitrary offsets, and the operator algebra is closed under composition (operator*), direct product (operator%), addition, and Hermitian conjugation.

The truncated algebra is represented exactly, including its physical deviations from the untruncated one: e.g. $[a,a^\dagger]$ at cutoff $D$ is $\mathrm{diag}(1,\ldots,1,1-D)$, not $\mathbb{1}$ — as it must be for any finite-dimensional representation.

Usage

Operators are built from single diagonals with multidiagonal::diagonal(offset, values), and combined with +, * (composition $AB$, the offsets add up) and % (direct product $A\otimes B$, the left operand's axes first). * and % share the multiplicative precedence, group from the left, and bind more tightly than +; % is not the element-wise product it is in Armadillo. A Jaynes–Cummings coupling between a qubit and a mode with cutoff 5, applied to the state $|0\rangle|1\rangle$:

#include <multidiagonal/MultiDiagonal.h>

#include <cmath>

using namespace cppqedutils;                     // dcomp, MultiArray
namespace md = quantumoperator::multidiagonal;

int main()
{
  const size_t N=5;                              // cutoff of the mode
  std::vector<dcomp> ladder(N-1);
  for (size_t k=0; k<N-1; ++k) ladder[k]=std::sqrt(k+1.);

  auto a {md::diagonal(+1,std::move(ladder))};   // ⟨n|a|n+1⟩ = √(n+1): the diagonal at offset +1
  auto sm{md::diagonal(+1,{1.})};                // σ⁻ = |0⟩⟨1|

  auto n {hermitianConjugateOf(a) * a};                                  // composition: the offsets add up to 0
  auto jc{hermitianConjugateOf(sm) % a + sm % hermitianConjugateOf(a)};  // direct product: qubit on axis 0, mode on axis 1

  MultiArray<dcomp,2> psi{{2,N},multiarray::zeroInit<dcomp>}, dpsi{{2,N},multiarray::zeroInit<dcomp>};
  psi(0,1)=1.;                                   // |0⟩|1⟩
  jc(0.,psi,dpsi);                               // accumulates jc|0⟩|1⟩ = |1⟩|0⟩ into dpsi
}

Application accumulates into its last argument and applies no factor $1/i$ (C++QED stores $H/i$). A diagonal whose coefficients don't factorize over the axes goes into the operator's diagonals map directly. CI compiles and runs this example against the installed library: it is the first part of tests/consumer/consumer.cc.

This is the beginning of the README. Read the rest at the original. Read the full README ↗