Kanamori interaction Hamiltonian

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'} d^\dagger_{m\uparrow} d_{m\downarrow} d^\dagger_{m'\uparrow} d_{m'\uparrow} +J\sum_{m\neq m'} d^\dagger_{m\uparrow} d^\dagger_{m\downarrow} d_{m'\downarrow} d_{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 (Eq. (7) of [GdMM2013]),

\[\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 program shows that \(\hat H_K = \hat H_{t_{2g}}\) indeed.

  1#include <cmath>
  2#include <iostream>
  3#include <string>
  4
  5#include <libcommute/libcommute.hpp>
  6
  7using namespace libcommute;
  8
  9int main() {
 10
 11  // For functions c_dag(), c() and n().
 12  using namespace static_indices;
 13
 14  // Orbital degeneracy of the shell (t_{2g} triplet)
 15  int const n_orbs = 3;
 16  // Coulomb integrals U and J
 17  double const U = 4.0;
 18  double const J = 0.2;
 19
 20  // Kanamori Hamiltonian (Eq. (2))
 21  static_indices::expr_real<int /* orbital index */,
 22                            std::string /* spin index */>
 23      H_K;
 24
 25  // Intraorbital density-density interaction terms
 26  for(int m = 0; m < n_orbs; ++m) {
 27    H_K += U * n(m, "up") * n(m, "down");
 28  }
 29
 30  // Interorbital density-density interaction terms (different spins)
 31  for(int m1 = 0; m1 < n_orbs; ++m1) {
 32    for(int m2 = 0; m2 < n_orbs; ++m2) {
 33      if(m1 == m2) continue;
 34      H_K += (U - 2 * J) * n(m1, "up") * n(m2, "down");
 35    }
 36  }
 37
 38  // Interorbital density-density interaction terms (same spin)
 39  for(int m1 = 0; m1 < n_orbs; ++m1) {
 40    for(int m2 = 0; m2 < n_orbs; ++m2) {
 41      if(m1 >= m2) continue;
 42      H_K += (U - 3 * J) * n(m1, "up") * n(m2, "up");
 43      H_K += (U - 3 * J) * n(m1, "down") * n(m2, "down");
 44    }
 45  }
 46
 47  // Spin-flip terms
 48  for(int m1 = 0; m1 < n_orbs; ++m1) {
 49    for(int m2 = 0; m2 < n_orbs; ++m2) {
 50      if(m1 == m2) continue;
 51      H_K += -J * c_dag(m1, "up") * c(m1, "down") * c_dag(m2, "down") *
 52             c(m2, "up");
 53    }
 54  }
 55
 56  // Pair-hopping terms
 57  for(int m1 = 0; m1 < n_orbs; ++m1) {
 58    for(int m2 = 0; m2 < n_orbs; ++m2) {
 59      if(m1 == m2) continue;
 60      H_K +=
 61          J * c_dag(m1, "up") * c_dag(m1, "down") * c(m2, "down") * c(m2, "up");
 62    }
 63  }
 64
 65  // Print the Hamiltonian
 66  std::cout << "H_K = " << H_K << '\n';
 67
 68  //
 69  // Integrals of motion N, S^2 and L^2 (Eq. (4)).
 70  //
 71
 72  std::complex<double> const I(0, 1);
 73
 74  // Total number of particles
 75  decltype(H_K) N;
 76  for(int m = 0; m < n_orbs; ++m) {
 77    N += n(m, "up") + n(m, "down");
 78  }
 79
 80  // Total spin operators S_x, S_y, S_z.
 81  static_indices::expr_complex< // to allow for complex coefficients in S_x, S_y
 82      int,
 83      std::string>
 84      Sx, Sy, Sz;
 85  for(int m = 0; m < n_orbs; ++m) {
 86    Sx += 0.5 * (c_dag(m, "up") * c(m, "down") + c_dag(m, "down") * c(m, "up"));
 87    Sy += 0.5 * I *
 88          (c_dag(m, "down") * c(m, "up") - c_dag(m, "up") * c(m, "down"));
 89    Sz += 0.5 * (n(m, "up") - n(m, "down"));
 90  }
 91  // Operator S^2 = S_x S_x + S_y S_y + S_z S_z
 92  auto S2 = Sx * Sx + Sy * Sy + Sz * Sz;
 93
 94  // Levi-Civita symbol \epsilon_{ijk}
 95  auto eps = [](int i, int j, int k) -> double {
 96    return (j - i) * (k - j) * (k - i) / 2;
 97  };
 98
 99  // Orbital isospin generators L_x, L_y, L_z.
100  static_indices::expr_complex< // to allow for complex coefficients
101      int,
102      std::string>
103      Lx, Ly, Lz;
104
105  for(std::string spin : {"up", "down"}) {
106    for(int m1 = 0; m1 < n_orbs; ++m1) {
107      for(int m2 = 0; m2 < n_orbs; ++m2) {
108        Lx += I * eps(0, m1, m2) * c_dag(m1, spin) * c(m2, spin);
109        Ly += I * eps(1, m1, m2) * c_dag(m1, spin) * c(m2, spin);
110        Lz += I * eps(2, m1, m2) * c_dag(m1, spin) * c(m2, spin);
111      }
112    }
113  }
114  // Operator L^2 = L_x L_x + L_y L_y + L_z L_z
115  auto L2 = Lx * Lx + Ly * Ly + Lz * Lz;
116
117  // Hamiltonian as a function of N, S^2 and L^2 (Eq. (7)).
118  auto H_t2g = (U - 3 * J) / 2 * N * (N - 1.0) - 2 * J * S2 - J / 2 * L2 +
119               (5.0 / 2) * J * N;
120
121  // Must be zero
122  std::cout << "H_K - H_t2g = " << (H_K - H_t2g) << '\n';
123
124  return 0;
125}
[GdMM2013] (1,2)

“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