# 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