Usage examples

Kanamori interaction Hamiltonian: Basic usage of expressions

The Kanamori multiorbital Hamiltonian describes Coulomb repulsion of \(t_{2g}\) electronic states as relevant to a transition-metal ion in a cubic crystal field with an octahedral environment [GdMM2013]. In its rotationally invariant form, the Hamiltonian reads

\[\begin{split}\begin{multline} \hat H_K = U \sum_{m} n_{m\uparrow} n_{m\downarrow} + (U-2J)\sum_{m\neq m'} n_{m\uparrow} n_{m'\downarrow} +\\+ (U-3J)\sum_{m<m',\sigma} n_{m\sigma} n_{m'\sigma} -J\sum_{m\neq m'} c^\dagger_{m\uparrow} c_{m\downarrow} c^\dagger_{m'\uparrow} c_{m'\uparrow} +J\sum_{m\neq m'} c^\dagger_{m\uparrow} c^\dagger_{m\downarrow} c_{m'\downarrow} c_{m'\uparrow}. \end{multline}\end{split}\]

The orbital indices \(m\), \(m'\) run over the \(t_{2g}\) triplet \(p_x, p_y, p_z\). \(U\) and \(J\) are Coulomb integrals (independent parameters of the model). The rotationally invariant Hamiltonian can be expressed in terms of a few integrals of motion,

\[\hat H_{t_{2g}} = (U-3J)\frac{\hat N (\hat N-1)}{2} - 2J\mathbf{S}^2 - \frac{J}{2}\mathbf{L}^2 + \frac{5}{2}J \hat N,\]

where \(\hat N\) is the total number of electrons on the shell,

\[\hat N = \sum_{m\sigma} n_{m\sigma},\]

\(\mathbf{S}\) is the total spin vector,

\[\mathbf{S} = \frac{1}{2} \sum_m \sum_{\sigma\sigma'} c^\dagger_{m\sigma} \boldsymbol{\tau}_{\sigma\sigma'} c_{m\sigma'},\]

and \(\mathbf{L}\) is the total orbital isospin,

\[L_m = i \sum_{mm'} \sum_{\sigma} \epsilon_{mm'm''} c^\dagger_{m'\sigma} c_{m''\sigma}.\]

The following script verifies that \(\hat H_K = \hat H_{t_{2g}}\) indeed.

  1#
  2# ## Kanamori interaction Hamiltonian
  3# ### and its expression in terms of integrals of motion N, S^2 and L^2.
  4#
  5#   "Strong Correlations from Hund's Coupling",
  6#   A. Georges, L. de' Medici and J. Mravlje,
  7#   Annu. Rev. Condens. Matter Phys. 2013. 4:137–78,
  8#   https://doi.org/10.1146/annurev-conmatphys-020911-125045
  9#
 10
 11from itertools import product
 12
 13# Polynomial expression with real and complex coefficients
 14from pycommute.expression import ExpressionR, ExpressionC
 15
 16# Factory functions returning expressions for fermionic creation, annihilation
 17# and occupation number operators.
 18from pycommute.expression import c_dag, c, n
 19
 20# Orbital degeneracy of the shell (t_{2g} triplet)
 21n_orbs = 3
 22# Coulomb integrals U and J
 23U = 4.0
 24J = 0.2
 25
 26# Kanamori Hamiltonian (Eq. (2))
 27H_K = ExpressionR()
 28
 29# Intraorbital density-density interaction terms
 30for m in range(n_orbs):
 31    H_K += U * n(m, "up") * n(m, "down")
 32
 33# Interorbital density-density interaction terms (different spins)
 34for m1, m2 in product(range(n_orbs), range(n_orbs)):
 35    if m1 == m2:
 36        continue
 37    H_K += (U - 2 * J) * n(m1, "up") * n(m2, "down")
 38
 39# Interorbital density-density interaction terms (same spin)
 40for m1, m2 in product(range(n_orbs), range(n_orbs)):
 41    if m1 >= m2:
 42        continue
 43    H_K += (U - 3 * J) * n(m1, "up") * n(m2, "up")
 44    H_K += (U - 3 * J) * n(m1, "down") * n(m2, "down")
 45
 46# Spin-flip terms
 47for m1, m2 in product(range(n_orbs), range(n_orbs)):
 48    if m1 == m2:
 49        continue
 50    H_K += \
 51        -J * c_dag(m1, "up") * c(m1, "down") * c_dag(m2, "down") * c(m2, "up")
 52
 53# Pair-hopping terms
 54for m1, m2 in product(range(n_orbs), range(n_orbs)):
 55    if m1 == m2:
 56        continue
 57    H_K += J * c_dag(m1, "up") * c_dag(m1, "down") * c(m2, "down") * c(m2, "up")
 58
 59# Print the Hamiltonian
 60print("H_K =", H_K)
 61
 62#
 63# ### Integrals of motion N, S^2 and L^2 (Eq. (4))
 64#
 65
 66# Total number of particles.
 67N = ExpressionR()
 68for m in range(n_orbs):
 69    N += n(m, "up") + n(m, "down")
 70
 71# Total spin operators S_x, S_y, S_z.
 72Sx = ExpressionC()
 73Sy = ExpressionC()
 74Sz = ExpressionC()
 75for m in range(n_orbs):
 76    Sx += 0.5 * (c_dag(m, "up") * c(m, "down") + c_dag(m, "down") * c(m, "up"))
 77    Sy += 0.5j * (c_dag(m, "down") * c(m, "up") - c_dag(m, "up") * c(m, "down"))
 78    Sz += 0.5 * (n(m, "up") - n(m, "down"))
 79
 80# Operator S^2 = S_x S_x + S_y S_y + S_z S_z.
 81S2 = Sx * Sx + Sy * Sy + Sz * Sz
 82
 83
 84# Levi-Civita symbol \epsilon_{ijk}.
 85def eps(i, j, k):
 86    return (j - i) * (k - j) * (k - i) / 2
 87
 88
 89# Orbital isospin generators L_x, L_y, L_z.
 90Lx = ExpressionC()
 91Ly = ExpressionC()
 92Lz = ExpressionC()
 93
 94for spin in ("up", "down"):
 95    for m1, m2 in product(range(n_orbs), range(n_orbs)):
 96        Lx += 1j * eps(0, m1, m2) * c_dag(m1, spin) * c(m2, spin)
 97        Ly += 1j * eps(1, m1, m2) * c_dag(m1, spin) * c(m2, spin)
 98        Lz += 1j * eps(2, m1, m2) * c_dag(m1, spin) * c(m2, spin)
 99
100# Operator L^2 = L_x L_x + L_y L_y + L_z L_z.
101L2 = Lx * Lx + Ly * Ly + Lz * Lz
102
103#
104# ### Hamiltonian as a function of N, S^2 and L^2 (Eq. (7))
105#
106
107H_t2g = (U - 3 * J) / 2 * N * (N - 1.0) - 2 * J * S2 - J / 2 * L2 \
108    + (5.0 / 2) * J * N
109
110# Must be zero
111print("H_K - H_t2g =", H_K - H_t2g)

Holstein model

This example shows how to construct Hamiltonian of the Holstein model on a square lattice and to check that the total number of electrons is a conserved quantity of the model.

The Hamiltonian considered here is a sum of three terms.

  • An electronic tight-binding model on a square 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}.\]

  • A harmonic oscillator at each lattice site (a localized phonon),

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

  • A coupling between the electrons and the phonons.

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

Instead of writing the sums over lattice sites explicitly, we call library functions tight_binding(), dispersion() and holstein_int(). We also make use of the NetworkX package to easily generate the adjacency matrix of the periodic square lattice.

  1#
  2# ## Holstein model on a square lattice with periodic boundary conditions
  3#
  4
  5from itertools import product
  6import numpy as np
  7
  8from pycommute.expression import ExpressionR, conj, FERMION, BOSON
  9from pycommute.expression import n
 10from pycommute.models import tight_binding, dispersion, holstein_int
 11
 12from networkx.generators.lattice import grid_2d_graph
 13from networkx.linalg.graphmatrix import adjacency_matrix
 14
 15#
 16# ### Hamiltonian of the electronic tight-binding model on a square lattice
 17#
 18
 19# Number of lattice sites in each direction
 20# (the total number of sites is N*N).
 21N = 10
 22
 23# Electron hopping constant - energy parameter of the TB model.
 24t = 2.0
 25
 26# Use NetworkX to construct the periodic square lattice (graph)
 27lat = grid_2d_graph(N, N, periodic=True)
 28
 29# Create lists of indices for electronic spin-up and spin-down operators.
 30# lat.nodes() returns a list of N^2 pairs of site indices (x, y).
 31indices_up = [(x, y, "up") for x, y in lat.nodes()]
 32indices_dn = [(x, y, "down") for x, y in lat.nodes()]
 33
 34# A sum of tight-binding Hamiltonians for both spins.
 35# The hopping matrix passed to tight_binding() is proportional to the
 36# adjacency matrix of the lattice.
 37hopping_matrix = -t * adjacency_matrix(lat).todense()
 38H_e = tight_binding(hopping_matrix, indices=indices_up) \
 39    + tight_binding(hopping_matrix, indices=indices_dn)
 40
 41#
 42# ### Hamiltonian of phonons localized at lattice sites
 43#
 44
 45# Frequency of the localized phonon.
 46w0 = 0.5
 47
 48# A lists of indices for bosonic operators, simply the (x, y) pairs
 49indices_phonon = lat.nodes()
 50
 51# All N^2 phonons have the same frequency
 52phonon_freqs = w0 * np.ones(N ** 2)
 53H_ph = dispersion(phonon_freqs, indices=indices_phonon, statistics=BOSON)
 54
 55#
 56# ### Hamiltonian of electron-phonon coupling
 57#
 58
 59# Electron-phonon coupling constant.
 60g = 0.1
 61
 62H_e_ph = holstein_int(g * np.ones(N ** 2),
 63                      indices_up=indices_up,
 64                      indices_dn=indices_dn,
 65                      indices_boson=indices_phonon)
 66
 67# Complete Holstein Hamiltonian.
 68H_H = H_e + H_ph + H_e_ph
 69
 70# Print H_H. There will be quite a lot of terms for the 100-site lattice!
 71print("H_H =", H_H)
 72
 73# Check hermiticity of H_H.
 74print(r"H_H - H_H^\dagger =", H_H - conj(H_H))
 75
 76# Check that H_H commutes with the total number of electrons N_e.
 77N_e = ExpressionR()
 78for spin in ("up", "down"):
 79    for x, y in product(range(N), range(N)):
 80        N_e += n(x, y, spin)
 81
 82print("[H_H, N_e] =", H_H * N_e - N_e * H_H)
 83
 84#
 85# ### Iteration interface
 86#
 87
 88# Iterate over monomial-coefficient pairs in polynomial expression H_H
 89for monomial, coeff in H_H:
 90    print("Coefficient:", coeff)
 91    # Iterate over algebra generators (creation/annihilation operators)
 92    # in the monomial
 93    for generator in monomial:
 94        # Detect what algebra this generator belongs to
 95        if generator.algebra_id == FERMION:
 96            print("\tFermionic", end=' ')
 97        elif generator.algebra_id == BOSON:
 98            print("\tBosonic", end=' ')
 99        # Creation or annihilation operator?
100        if generator.dagger:
101            print("creation operator", end=' ')
102        else:
103            print("annihilation operator", end=' ')
104        # Extract indices carried by the generator
105        print("with indices", list(generator.indices))
106        # N.B. generator.indices is an instance of a special tuple-like type
107        # `pycommute.expression.Indices`. Types of its elements are restricted
108        # to integer and string, and its comparison operators behave differently
109        # from those of the Python tuple. Otherwise, it supports `len()`,
110        # indexed element access and iteration protocol.

Spin-1/2 Heisenberg chain

The spin-1/2 Heisenberg chain is a textbook example of an integrable quantum system. Its Hamiltonian

\[\hat H = g \sum_i \mathbf{S}_i \cdot \mathbf{S}_{i+1}\]

conserves three projections of the total spin

\[\mathbf{S} = \sum_i \mathbf{S}_i\]

as well as a series of higher order charges \(Q_n\). Existence of these charges can be derived from the transfer matrix theory. Explicit expressions for \(Q_n\) were obtained in [GM94]. The following script constructs Hamiltonian of the Heisenberg chain with periodic boundary conditions (see heisenberg()) and checks that \([\hat H, \mathbf{S}] = 0\), \([\hat H, Q_n] = 0\) and \([Q_n, Q_m] = 0\) for \(m,n = 3,4,5\).

 1#
 2# ## Periodic spin-1/2 Heisenberg chain and its integrals of motion
 3#
 4# Expressions for the integrals of motion are taken from
 5#
 6#   "Quantum Integrals of Motion for the Heisenberg Spin Chain",
 7#   M. P. Grabowski and P. Mathieu
 8#   Mod. Phys. Lett. A, Vol. 09, No. 24, pp. 2197-2206 (1994),
 9#   https://doi.org/10.1142/S0217732394002057
10#
11
12from numpy import array, zeros, dot, cross
13from pycommute.expression import S_x, S_y, S_z
14from pycommute.models import heisenberg
15
16# Number of spins in the chain
17N = 20
18# Heisenberg exchange constant
19g = 2
20
21# List of 3-component spin vectors {S_0, S_1, ..., S_{N-1}}
22S = [array([S_x(i), S_y(i), S_z(i)]) for i in range(N)]
23
24# Matrix of exchange constants between spins i and j
25exchange_matrix = zeros((N, N))
26# Set elements corresponding to the nearest neighbors to -g
27# (index shift modulo N ensures periodic boundary conditions).
28for i in range(N):
29    exchange_matrix[i, (i + 1) % N] = -g
30
31# Hamiltonian of the spin-1/2 Heisenberg chain.
32H = heisenberg(exchange_matrix)
33
34# Total spin of the chain.
35S_tot = array(sum(S))
36
37# All three components of S_tot commute with the Hamiltonian.
38print("[H, S_x] =", (H * S_tot[0] - S_tot[0] * H))
39print("[H, S_y] =", (H * S_tot[1] - S_tot[1] * H))
40print("[H, S_z] =", (H * S_tot[2] - S_tot[2] * H))
41
42# Higher charge Q_3 (1st line of Eq. (10)).
43Q3 = sum(dot(cross(S[i], S[(i + 1) % N]), S[(i + 2) % N]) for i in range(N))
44print("[H, Q3] =", (H * Q3 - Q3 * H))
45
46# Higher charge Q_4 (2nd line of Eq. (10)).
47Q4 = sum(4.0 * dot(
48    cross(
49        cross(S[i], S[(i + 1) % N]), S[(i + 2) % N]
50    ), S[(i + 3) % N]) for i in range(N))
51Q4 += sum(dot(S[i], S[(i + 2) % N]) for i in range(N))
52print("[H, Q4] =", (H * Q4 - Q4 * H))
53
54# Higher charge Q_5 (3rd line of Eq. (10)).
55Q5 = sum(4.0 * dot(
56    cross(
57        cross(
58            cross(S[i], S[(i + 1) % N]),
59            S[(i + 2) % N]
60        ), S[(i + 3) % N]),
61    S[(i + 4) % N]) for i in range(N))
62Q5 += sum(dot(cross(S[i], S[(i + 2) % N]), S[(i + 3) % N]) for i in range(N))
63Q5 += sum(dot(cross(S[i], S[(i + 1) % N]), S[(i + 3) % N]) for i in range(N))
64print("[H, Q5] =", (H * Q5 - Q5 * H))
65
66# Check that the higher charges pairwise commute.
67print("[Q3, Q4] =", (Q3 * Q4 - Q4 * Q3))
68print("[Q3, Q5] =", (Q3 * Q5 - Q5 * Q3))
69print("[Q4, Q5] =", (Q4 * Q5 - Q5 * Q4))

Spectrum of Tavis-Cummings model

In this example we demonstrate how to use tools from the loperator module along with NumPy’s eigensolvers to diagonalize Hamiltonians of simple quantum models.

We will consider a two-atom Jaynes-Cummings (Tavis-Cummings) Hamiltonian [TC68], which models a system of two qubits coherently coupled via a bosonic degree of freedom (quantum oscillator),

\[\hat H = \epsilon (\hat S_{z,0} + \hat S_{z,1}) + \omega \hat a^\dagger \hat a + g_0 (\hat a^\dagger \hat S_{-,0} + \hat a \hat S_{+,0}) + g_1 (\hat a^\dagger \hat S_{-,1} + \hat a \hat S_{+,1}).\]

An expression of this form can be conveniently generated by a single call to jaynes_cummings().

Since the occupation number of the bosonic mode can formally grow to \(+\infty\), we have to artificially restrict the dimension of the bosonic Hilbert space in order to be able to construct a finite matrix representation of our problem.

Note

Performance of repeated calls to LOperatorR (or LOperatorC) in order to construct a matrix representation of a linear operator is limited by Python method call overhead. For large-scale problems it is advised to call the utility function make_matrix().

 1#
 2# ## Diagonalization of a two-qubit Tavis-Cummings model
 3#
 4
 5import numpy as np
 6
 7# Utility function used to construct the generalized Jaynes-Cummings model.
 8from pycommute.models import jaynes_cummings
 9
10# Hilbert spaces.
11from pycommute.loperator import HilbertSpace, make_space_boson, make_space_spin
12
13# Real-valued linear operator object.
14from pycommute.loperator import LOperatorR
15
16# Build the matrix form of a linear operator.
17from pycommute.loperator import make_matrix
18
19# Transition frequencies of qubits.
20eps = np.array([1.0, 1.0])
21
22# Oscillator frequency.
23# (In a more general case of M oscillators, omega would be a length-M array).
24omega = np.array([0.8])
25
26# Qubit-oscillator coupling constants as a 2x1 array.
27g = np.array([[0.5], [0.6]])
28
29# Create the Tavis-Cummings Hamiltonian.
30H = jaynes_cummings(eps, omega, g)
31
32# Construct state space of our problem as a direct product of two
33# two-dimensional Hilbert spaces (qubits) and one truncated bosonic Hilbert
34# space.
35# make_space_boson(16, ...) returns the truncated bosonic space with allowed
36# occupation numbers N = 0, 1, ..., 15.
37hs = HilbertSpace([
38    make_space_spin(1 / 2, 0),  # Qubit 1: spin-1/2, index 0
39    make_space_spin(1 / 2, 1),  # Qubit 2: spin-1/2, index 1
40    make_space_boson(16, 0)     # Oscillator, index 0
41])
42
43# Construct a linear operator corresponding to 'H' and acting in the Hilbert
44# space 'hs'.
45H_op = LOperatorR(H, hs)
46
47#
48# ### Prepare a matrix representation of `H_op`
49#
50
51# Method I (manual).
52H_mat1 = np.zeros((hs.vec_size, hs.vec_size))
53for i in range(hs.vec_size):
54    # A column vector psi = {0, 0, ..., 1, ..., 0}
55    psi = np.zeros(hs.vec_size)
56    psi[i] = 1.0
57    # Act with H_op on psi and store the result the i-th column of H_mat
58    H_mat1[:, i] = H_op * psi
59
60# Method II (automatic and faster).
61H_mat2 = make_matrix(H_op, hs)
62print("Max difference between H_mat1 and H_mat2:",
63      np.max(np.abs(H_mat1 - H_mat2)))
64
65#
66# ### Use NumPy to compute eigenvalues of `H_mat1`
67#
68
69E = np.linalg.eigvals(H_mat1)
70
71print("Energies:", np.sort(E))

Sectors of a multi-orbital interaction Hamiltonian

When implementing an exact diagonalization algorithm, it is usually advantageous to take an extra preparatory step and split the Hilbert space of the problem into subspaces (sectors) characterized by a set of quantum numbers. In many situations, it is possible to find the sectors even without knowing the quantities conserved by the Hamiltonian (see, for instance, [SKFP16]).

In the script shown below we construct an electron-electron interaction Hamiltonian of an atomic \(d\)-shell using Slater parameterization. Then, we employ the SpacePartition and BasisMapper utility classes (and optionally the make_matrix() function) to independently diagonalize the Hamiltonian within each sector.

 1#
 2# ## Sector-wise diagonalization of the Slater interaction Hamiltonian
 3#
 4
 5import numpy as np
 6
 7# Utility function used to construct the Slater interaction Hamiltonian.
 8from pycommute.models import slater_int
 9
10# Hilbert space.
11from pycommute.loperator import HilbertSpace
12
13# Real-valued linear operator object.
14from pycommute.loperator import LOperatorR
15# Automatic space partition and basis state mapper
16from pycommute.loperator import make_space_partition, BasisMapper
17# Build the matrix form of a linear operator
18from pycommute.loperator import make_matrix
19
20# Values of Slater radial integrals F^0, F^2, F^4
21F = np.array([4.0, 1.0, 0.2])
22
23# Build an expression for the interaction Hamiltonian.
24H = slater_int(F)
25
26# Analyze structure of 'H' and construct a suitable Hilbert space.
27hs = HilbertSpace(H)
28
29# Construct a linear operator corresponding to H and acting in the Hilbert
30# space 'hs'.
31H_op = LOperatorR(H, hs)
32
33# The first value returned by 'make_space_partition()' is a
34# SpacePartition object. It represents a partition of the full Hilbert space
35# into sectors, i.e. subspaces invariant under the action of 'H_op'.
36#
37# The second returned value is a dictionary {(i, j): value} of all non-vanishing
38# matrix elements of H_{ij}. By definition, all matrix elements of 'H_op'
39# between different sectors vanish.
40sp, matrix_elements = make_space_partition(H_op, hs)
41
42# Print out information about the revealed sectors.
43print("Dimension of full Hilbert space:", sp.dim)
44print("Total number of sectors:", sp.n_subspaces)
45for i in range(sp.dim):
46    print("Many-body basis state %d belongs to sector %d" % (i, sp[i]))
47
48# Compile lists of basis states spanning each sector.
49sectors = sp.subspace_bases()
50
51# Diagonalize 'H_op' within each sector
52for n, sector in enumerate(sectors):
53    sector_dim = len(sector)
54    print("Diagonalizing sector %d of H_op" % n)
55    print("Sector size is %d" % sector_dim)
56
57    # A BasisMapper object translates indices of basis states from the full
58    # Hilbert space to a given sector.
59    basis_mapper = BasisMapper(sector)
60
61    #
62    # Prepare a matrix representation of `H_op` within current sector
63    #
64
65    # Method I (manual).
66    H_mat1 = np.zeros((sector_dim, sector_dim))
67    for i in range(sector_dim):
68        # A column vector psi = {0, 0, ..., 1, ..., 0}
69        psi = np.zeros(sector_dim)
70        psi[i] = 1.0
71        # This vector will receive the result of H_op * psi
72        phi = np.zeros(sector_dim)
73
74        # Since both 'psi' and 'phi' are 1D NumPy arrays with size of the chosen
75        # sector, 'H_op' cannot act on them directly as it expects vectors of
76        # size hs.dim.
77        # Instead, we are going to use our basis mapper object to construct
78        # special views of psi and phi.
79        psi_view = basis_mapper(psi)
80        phi_view = basis_mapper(phi)
81
82        # Now, 'H_op' can act on the mapped views.
83        H_op(psi_view, phi_view)
84
85        # Store H_op * psi in the i-th column of 'H_mat1'.
86        H_mat1[:, i] = phi
87
88    # Method II (automatic and faster).
89    H_mat2 = make_matrix(H_op, sector)
90    print("Max difference between H_mat1 and H_mat2:",
91          np.max(np.abs(H_mat1 - H_mat2)))
92
93    # Use NumPy to compute eigenvalues of 'H_mat1'.
94    E = np.linalg.eigvals(H_mat1)
95    print("Energies:", np.sort(E))

N-fermion sector views

This example shows how to partially diagonalize the single band Fermi-Hubbard model on a large 2D square lattice with periodic boundary conditions.

\[\hat H = -t\sum_{\langle i,j \rangle, \sigma} c^\dagger_{i,\sigma} c_{j,\sigma} -\mu \sum_{i, \sigma} n_{i,\sigma} + U \sum_i n_{i,\uparrow} n_{i,\downarrow}.\]

The number of atoms in the considered lattice (cluster) is set to 16, which makes for an intractably large Hilbert space of dimension \(2^{32}\). Since the Hamiltonian of the Hubbard model preserves total numbers of spin-up and spin-down electrons, one can diagonalize the Hamiltonian within a single \(N\)-fermion sector or within a \((N_\uparrow, N_\downarrow)\)-multisector.

An \(N\)-fermion sector is a subspace of a full Hilbert space, which is spanned by all basis states with a fixed total occupation of fermionic degrees of freedom (FDOF). Similarly, a multisector is spanned by all basis states with a fixed occupation \(N_1\) of a subset of the FDOF \(\{S_1\}\), occupation \(N_2\) of another subset \(\{S_2\}\) and so on. There can be any number of pairs \((\{S_i\}, N_i)\) (sectors contributing to the multisector) as long as all the subsets \(\{S_i\}\) are disjoint.

In our example we consider a moderately sized sector with \(N = 2\) (\(\dim = {32 \choose 2} = 496\)) and a multisector with \(N_\uparrow = 1, N_\downarrow = 1\) (\(\dim = {16 \choose 1}{16 \choose 1} = 256\)).

Note

In general, the N-fermion (multi)sector views and functions do not require a purely fermionic system. The definitions of a sector and a multisector stand in presence of bosons or spins.

  1#
  2# ## Fermi-Hubbard model on a large square lattice
  3# ### and diagonalization of its moderately sized sectors
  4#
  5
  6import numpy as np
  7
  8# Utility functions used to construct the Fermi-Hubbard Hamiltonian.
  9from pycommute.models import tight_binding, dispersion, hubbard_int
 10
 11# Hilbert space.
 12from pycommute.loperator import HilbertSpace
 13
 14# Real-valued linear operator object.
 15from pycommute.loperator import LOperatorR
 16# Real-valued N-fermion (multi)sector views and related functions
 17from pycommute.loperator import (
 18    NFermionSectorViewR,
 19    NFermionMultiSectorViewR,
 20    n_fermion_sector_size,
 21    n_fermion_multisector_size
 22)
 23
 24from networkx.generators.lattice import grid_2d_graph
 25from networkx.linalg.graphmatrix import adjacency_matrix
 26
 27#
 28# ### Hamiltonian of a tight-binding model on a square lattice
 29#
 30
 31# Number of lattice sites in each direction
 32# (the total number of sites is Nx * Ny).
 33Nx = 4
 34Ny = 4
 35
 36# Hopping constant
 37t = 0.5
 38# Chemical potential
 39mu = 1.0
 40# Coulomb repulsion
 41U = 2.0
 42
 43# Use NetworkX to construct the periodic square lattice (graph)
 44lat = grid_2d_graph(Nx, Ny, periodic=True)
 45
 46# Create lists of indices for spin-up and spin-down operators.
 47# lat.nodes() returns a list of Nx * Ny pairs of site indices (x, y).
 48indices_up = [(x, y, "up") for x, y in lat.nodes()]
 49indices_dn = [(x, y, "down") for x, y in lat.nodes()]
 50
 51# A sum of tight-binding Hamiltonians for both spins.
 52# The hopping matrix passed to tight_binding() is proportional to the
 53# adjacency matrix of the lattice.
 54hopping_matrix = -t * adjacency_matrix(lat).todense()
 55H = tight_binding(hopping_matrix, indices=indices_up) \
 56    + tight_binding(hopping_matrix, indices=indices_dn)
 57
 58# Add the chemical potential terms
 59H += dispersion(-mu * np.ones(len(indices_up)), indices=indices_up)
 60H += dispersion(-mu * np.ones(len(indices_dn)), indices=indices_dn)
 61
 62# Add the Hubbard interaction term
 63H += hubbard_int(U * np.ones(len(indices_up)),
 64                 indices_up=indices_up,
 65                 indices_dn=indices_dn)
 66
 67# Analyze structure of H and construct a suitable Hilbert space.
 68hs = HilbertSpace(H)
 69print("Full Hilbert space dimension:", hs.dim)
 70
 71# Construct a linear operator corresponding to H and acting in the Hilbert
 72# space 'hs'.
 73H_op = LOperatorR(H, hs)
 74
 75#
 76# ### Diagonalize the N = 2 sector of the model
 77#
 78
 79N = 2
 80
 81sector_size = n_fermion_sector_size(hs, N)
 82print("Size of the N = 2 sector:", sector_size)
 83
 84# Prepare a matrix representation of H_op within the N = 2 sector.
 85H_mat = np.zeros((sector_size, sector_size))
 86for i in range(sector_size):
 87    # A column vector psi = {0, 0, ..., 1, ..., 0}
 88    psi = np.zeros(sector_size)
 89    psi[i] = 1.0
 90    # This vector will receive the result of H_op * psi
 91    phi = np.zeros(sector_size)
 92
 93    # Since both psi and phi are 1D NumPy arrays with size of the chosen
 94    # sector, H_op cannot act on them directly as it expects vectors of
 95    # size hs.dim.
 96    # Instead, we are going to use real-valued 2-fermion sector views of
 97    # psi and phi.
 98    psi_view = NFermionSectorViewR(psi, hs, N)
 99    phi_view = NFermionSectorViewR(phi, hs, N)
100
101    # Now, H_op can act on the mapped views.
102    H_op(psi_view, phi_view)
103
104    # Store H_op * psi in the i-th column of H_mat.
105    H_mat[:, i] = phi
106
107# Use NumPy to compute eigenvalues of H_mat.
108E = np.linalg.eigvals(H_mat).real
109print("10 lowest eigenvalues of the N = 2 sector:", np.sort(E)[:10])
110
111#
112# ### Diagonalize the N_up = 1, N_down = 1 multisector of the model
113#
114
115N_up = 1
116N_dn = 1
117
118# Define sectors contributing to the multisector
119sector_up = (indices_up, N_up)
120sector_dn = (indices_dn, N_dn)
121
122sectors = [sector_up, sector_dn]
123
124multisector_size = n_fermion_multisector_size(hs, sectors)
125print("Size of the N_up = 1, N_down = 1 multisector:", multisector_size)
126
127# Prepare a matrix representation of H_op within the multisector.
128H_mat = np.zeros((multisector_size, multisector_size))
129for i in range(multisector_size):
130    # A column vector psi = {0, 0, ..., 1, ..., 0}
131    psi = np.zeros(multisector_size)
132    psi[i] = 1.0
133    # This vector will receive the result of H_op * psi
134    phi = np.zeros(multisector_size)
135
136    # Since both psi and phi are 1D NumPy arrays with size of the chosen
137    # sector, H_op cannot act on them directly as it expects vectors of
138    # size hs.dim.
139    # Instead, we are going to use real-valued multisector views of psi and phi.
140    psi_view = NFermionMultiSectorViewR(psi, hs, sectors)
141    phi_view = NFermionMultiSectorViewR(phi, hs, sectors)
142
143    # Now, H_op can act on the mapped views.
144    H_op(psi_view, phi_view)
145
146    # Store H_op * psi in the i-th column of H_mat.
147    H_mat[:, i] = phi
148
149# Use NumPy to compute eigenvalues of H_mat.
150E = np.linalg.eigvals(H_mat).real
151print("10 lowest eigenvalues of the multisector:", np.sort(E)[:10])

Advanced: A user-defined algebra

It is possible to define generators of a new algebra in addition to the fermionic, bosonic and spin operators to be used in expressions. A basic introduction into the topic can be found on the respective documentation page of libcommute. Here, we show how to define the algebra of Dirac \(\gamma\)-matrices in a Python script by extending the base class pycommute.expression.Generator.

  1#
  2# ## User-defined algebra of Dirac gamma matrices
  3#
  4
  5# Algebra generators and expressions
  6from pycommute.expression import (Indices,
  7                                  Generator,
  8                                  MIN_USER_DEFINED_ALGEBRA_ID,
  9                                  Monomial,
 10                                  ExpressionC,
 11                                  conj)
 12
 13
 14#
 15# Our generator type: a gamma matrix with one integer index \mu = 0, ..., 3
 16#
 17class GeneratorGamma(Generator):
 18
 19    # Initialize the base class by passing the index \mu to it.
 20    def __init__(self, index):
 21        Generator.__init__(self, Indices(index))
 22
 23    # Algebra ID of this generator. In this case, the lowest algebra ID
 24    # available to user-defined algebras. Further algebras can use IDs
 25    # MIN_USER_DEFINED_ALGEBRA_ID + 1, ...
 26    def algebra_id(self):
 27        return MIN_USER_DEFINED_ALGEBRA_ID
 28
 29    # This method will be called for g1 = self and g2 such that g1 > g2.
 30    # We must transform the product g1 * g2 and put it into the
 31    # canonical order,
 32    #
 33    # g1 * g2 -> -g2 * g1 + 2\eta(g1, g2) = -g2 * g1 + f.
 34    def swap_with(self, g2, f):
 35        # Update the LinearFunctionGen object 'f'
 36        #
 37        # Set the constant term 2\eta(g1, g2), where \eta is the diagonal
 38        # Minkowski metric tensor.
 39        mu = self.indices[0]
 40        nu = g2.indices[0]
 41        f.const_term = (mu == nu) * (2 if mu == 0 else -2)
 42        # 'f' gets no contributions linear in the generators.
 43        f.terms = []
 44        return -1     # Return the coefficient in front of g2 * g1.
 45
 46    # This method will be called for g1 = self and g2 that are already
 47    # canonically ordered, g1 <= g2.
 48    #
 49    # It tries to simplify squares of gamma matrices,
 50    # f = (\gamma^0)^2 = I_4
 51    # f = (\gamma^k)^2 = -I_4 for k = 1,2,3.
 52    def simplify_prod(self, g2, f):
 53        if self == g2:
 54            # Update the LinearFunctionGen object 'f'.
 55            #
 56            # The constant term of 'f' is equal to either 1 or -1, depending
 57            # on which gamma matrix is being squared.
 58            mu = self.indices[0]
 59            f.const_term = 1 if (mu == 0) else -1
 60            # 'f' gets no contributions linear in the generators.
 61            f.terms = []
 62            return True   # A square of a gamma matrix has been simplified.
 63        else:
 64            return False  # No simplification has been applied.
 65
 66    # Compute Hermitian conjugate of a gamma matrix g = self,
 67    # f = (\gamma_0)^+ = \gamma_0
 68    # f = (\gamma_k)^+ = -\gamma_k, k=1,2,3.
 69    #
 70    # N.B.: If this method is not overridden in the subclass, all generators
 71    # of the new algebra are assumed to be Hermitian.
 72    def conj(self, f):
 73        # Update the LinearFunctionGen object 'f'
 74        #
 75        # Hermitian conjugate of 'g' contains no constant contribution.
 76        f.const_term = 0
 77        # The only linear in generators term is proportional to the argument
 78        # g = self, and the coefficient in front of it depends on the index.
 79        mu = self.indices[0]
 80        f.terms = [(self, 1 if (mu == 0) else -1)]
 81
 82    # String representation of this generator
 83    def __repr__(self):
 84        mu = self.indices[0]
 85        return f"γ^{mu}"
 86
 87
 88# Convenience function: Make an expression out of one gamma matrix.
 89def gamma(mu):
 90    return ExpressionC(1.0, Monomial([GeneratorGamma(mu)]))
 91
 92
 93#
 94# ### Check that expressions with gamma matrices behave as expected
 95#
 96
 97# Minkowski metric tensor
 98def eta(mu, nu):
 99    return (mu == nu) * (1.0 if mu == 0 else -1.0)
100
101
102# Check the commutation relations
103for mu in range(4):
104    for nu in range(4):
105        gamma_mu = gamma(mu)
106        gamma_nu = gamma(nu)
107
108        print(f"{{{gamma_mu}, {gamma_nu}}} - 2η({mu}, {nu}) =",
109              (gamma_mu * gamma_nu + gamma_nu * gamma_mu - 2 * eta(mu, nu)))
110
111# \gamma^5
112gamma5 = 1j * gamma(0) * gamma(1) * gamma(2) * gamma(3)
113
114# \gamma^5 is Hermitian ...
115print("γ^5 - conj(γ^5) =", gamma5 - conj(gamma5))
116# ... and anti-commutes with \gamma_\mu.
117for mu in range(4):
118    gamma_mu = gamma(mu)
119    print(f"{{γ^5, {gamma_mu}}} =", gamma5 * gamma_mu + gamma_mu * gamma5)

References

[GdMM2013]

“Strong Correlations from Hund’s Coupling”, A. Georges, L. de’ Medici and J. Mravlje, Annu. Rev. Condens. Matter Phys. 2013. 4:137–78, https://doi.org/10.1146/annurev-conmatphys-020911-125045

[GM94]

“Quantum Integrals of Motion for the Heisenberg Spin Chain”, M. P. Grabowski and P. Mathieu, Mod. Phys. Lett. A, Vol. 09, No. 24, pp. 2197-2206 (1994), https://doi.org/10.1142/S0217732394002057

[TC68]

“Exact Solution for an N-Molecule-Radiation-Field Hamiltonian”, M. Tavis and F. W. Cummings Phys. Rev. 170, 379 (1968) https://doi.org/10.1103/PhysRev.170.379

[SKFP16]

“TRIQS/CTHYB: A continuous-time quantum Monte Carlo hybridisation expansion solver for quantum impurity problems”, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet, Comp. Phys. Comm. 200, March 2016, 274-284, http://dx.doi.org/10.1016/j.cpc.2015.10.023 (section 4.2)