Operators

Besides Hilbert spaces, the second key objects in quantum mechanics are operators. In a many-body setting, we consider operators of the form, $$ O = \sum_{A\subseteq \mathcal{L}} c_A O_A, $$ where \(O_A\) denotes a local operator acting on sites \(A=\{a_1, \ldots, a_{l_A}\}\) and \(\mathcal{L}\) denotes the lattice and \(c_{A}\) are coupling constants. In the case of the Heisenberg model, we would thus have \(\mathcal{O}_{A} = \mathbf{S}_i\cdot\mathbf{S}_j\) and \(c_A = J\). In XDiag, the local operators are represented via an Op object while the sum is represented by an OpSum object. These values of the coupling constants \(c_A\) can either be a real or complex number, or a string that later needs to be replaced. The Hamiltonian of a spin \(S=1 / 2\) Heisenberg chain is created in the following way:

JuliaC++

ops = OpSum()
for i in 1:N
    s1 = i
    s2 = mod1(i+1, N)
    ops += "J" * Op("SdotS", [s1, s2])
end
ops["J"] = 1.0

auto ops = OpSum();
for (int i=0; i<N; ++i) {
    int s1 = i;
    int s2 = (i+1) % N;
    ops += "J" * Op("SdotS", {s1, s2});
}
ops["J"] = 1.0;

We first create an empty OpSum and then add additional terms to it. The first part of the product denotes the coupling constant, here given as a string. Alternatively, one could have directly used real / complex numbers here. The second part of the product is a single Op object. It is created with two inputs:

• The type, here SdotS denoting an operator of the form \(\mathbf{S}_{i} \cdot \mathbf{S}_{j}\). XDiag features a wide variety of different operator types.
• An array defining which site the operator lives on. Notice, that in Julia we start counting the sites from 1, while in C++ we start counting the sites from 0.

Which operator types are available, and on which blocks they may be used, is documented in two places:

• On the documentation page of each block type (for example Spinhalf, tJ, Electron, Boson, or Fermion), a table lists all operator types that can be applied to that particular block.
• A central register of every operator type, together with its definition, its required number of sites, and the blocks it is defined for, is collected on the operator types page.

Coupling constants

The coefficient \(c_A\) multiplying an operator can be given in two ways. It can be a plain real or complex number, which is the most direct choice. Alternatively, it can be a named coupling given as a string, such as "J" in the example above. Named couplings are convenient because they let us define the structure of a model once and only fix the numerical values later. This is especially useful when the same Hamiltonian is diagonalized repeatedly for different parameter values, or when parameters are read from an input file.

The value of a named coupling is assigned using the [] operator, and the currently assigned values can be queried the same way. Once all named couplings have been assigned a value, the plain() function returns a copy of the OpSum in which every named coupling has been substituted by its numerical value.

JuliaC++

# "J" is a named coupling. Its value is set with the [] operator, ...
ops["J"] = 1.0
# ... and plain() substitutes every named coupling by its numerical value.
ops_plain = plain(ops)

// "J" is a named coupling. Its value is set with the [] operator, ...
ops["J"] = 1.0;
// ... and plain() substitutes every named coupling by its numerical value.
auto ops_plain = ops.plain();

Most computational routines call plain() internally, so in practice it is rarely necessary to invoke it by hand. An OpSum that still contains unresolved named couplings cannot be turned into a matrix or applied to a state.

The operator algebra

Beyond forming linear combinations, XDiag now supports the full set of algebraic operations on operators. In addition to addition, subtraction, and multiplication by a scalar, two OpSum objects can be multiplied with one another. This product is the ordinary composition of operators and is distributed over the sum,

\[ \Big(\sum_i c_i M_i\Big)\Big(\sum_j d_j N_j\Big) = \sum_{i,j} (c_i d_j)\, M_i N_j. \]

Together with the vector space operations, this multiplication turns the set of OpSum objects into an algebra. Note that the product is generally non-commutative, reflecting the fact that quantum mechanical operators need not commute. A simple example is the product of two single-site spin operators,

JuliaC++

# OpSums can be multiplied, forming the (non-commutative) operator algebra product
sz0 = 1.0 * Op("Sz", 1)
sz1 = 1.0 * Op("Sz", 2)
szsz = sz0 * sz1

// OpSums can be multiplied, forming the (non-commutative) operator algebra product
auto sz0 = 1.0 * Op("Sz", 0);
auto sz1 = 1.0 * Op("Sz", 1);
auto szsz = sz0 * sz1;

This makes it possible to build composite operators, such as correlation functions or powers of a Hamiltonian, directly from elementary building blocks.

Hermitian conjugation

The algebra is complemented by an operation of hermitian conjugation, computed with the hc function. It returns the hermitian conjugate \(\mathcal{O}^\dagger\) of an Op or OpSum, complex-conjugating the coefficients and conjugating each elementary operator (for example \(\mathrm{hc}(S^+) = S^-\)).

JuliaC++

# hermitian conjugation via hc(), e.g. hc(S+) = S-
sp = 1.0 * Op("S+", 1)
sm = hc(sp)
sx = sp + hc(sp)   # a hermitian combination

// hermitian conjugation via hc(), e.g. hc(S+) = S-
auto sp = 1.0 * Op("S+", 0);
auto sm = hc(sp);
auto sx = sp + hc(sp);   // a hermitian combination

Hermitian conjugation is an involution: applying it twice returns the original operator, \((\mathcal{O}^\dagger)^\dagger = \mathcal{O}\), and it reverses the order of products, \((\mathcal{A}\mathcal{B})^\dagger = \mathcal{B}^\dagger \mathcal{A}^\dagger\). An algebra equipped with such an operation is called an involutive algebra (or \(*\)-algebra). The OpSum objects of XDiag therefore realize precisely the mathematical structure that operators in quantum mechanics are expected to have, and hc is the natural tool for building hermitian Hamiltonians, for instance by symmetrizing a non-hermitian term as A + hc(A).