\(S=\frac12\) chain: Symmetries and Ground State Correlators

Author: Paul Ebert

This example demonstrates how ground state expectation values can be computed for the spin-\(\frac{1}{2}\) XXX antiferromagnetic chain

\[ \mathcal{H} = \sum_{\langle i, j \rangle} \bm{S}_i \cdot \bm{S}_j = \sum_{\langle i, j \rangle} (S^x_i S^x_j + S^y_i S^y_j + S^z_i S^z_j) \]

with periodic boundary conditions. The Hamiltonian is invariant under translations and decomposes into \(N\) irreducible representations (irreps) labelled by lattice momentum \(2\pi/N \times k\) where \(k=0, 1, \ldots, N-1\). The Lanczos algorithm is run on all momentum sub-blocks, effectively exchanging a longer runtime for less memory consumption. An important caveat when working with representations is that the operators inside expectation values must be symmetrized with respect to the symmetry group.

C++Julia

#include <xdiag/all.hpp>

using namespace xdiag;

std::complex<double> C_N_character(int N, int k, int p) {
  return std::exp(((2 * k * p) * M_PI / (double)N) *
                  std::complex<double>(0, 1.));
}

int main() try {
  // periodic N-site Heisenberg antiferromagnet
  int N = 14;
  auto H = OpSum();
  for (int i = 0; i < N; i++) {
    H += "J" * Op("SdotS", {i, (i + 1) % N});
  }
  H["J"] = 1.0;

  // define shift by one site
  std::vector<int64_t> T_perm(N, 0);
  for (int i = 0; i < N - 1; i++) {
    T_perm[i] = i + 1;
  }
  auto T = Permutation(T_perm);

  // define cyclic group
  std::vector<Permutation> perms(N);
  for (int i = 0; i < N; i++) {
    perms[i] = pow(T, i);
  }
  auto C_N = PermutationGroup(perms);

  // define irreps from character tables, labelled by momentum (2pi/N ×) k
  std::vector<Representation> irreps(N);
  std::vector<std::complex<double>> characters(N, 0.0);
  for (int k = 0; k < N; k++) {
    for (int p = 0; p < N; p++) {
      characters[p] = C_N_character(N, k, p);
    };
    irreps[k] = Representation(C_N, characters);
  }

  // check all translation symmetry blocks
  double e0 = 0.0;
  State psi0;
  for (int k = 0; k < N; k++) {
    auto block = Spinhalf(N, irreps[k]);
    auto [e0k, psik] = eig0(H, block);
    if (e0k < e0) {
      e0 = e0k;
      psi0 = psik;
    }
  }

  // we measure the correlation (which needs to besymmetrized)
  auto corr_op = Op("SdotS", {0, N / 2 - 1});
  auto corr = inner(symmetrize(corr_op, C_N), psi0);

  Log("Ground state energy: {:.12f}", e0);
  Log("Ground state correlator: {:.12f}", corr);
} catch (Error e) {
  error_trace(e);
}

using XDiag

function main()

    # periodic N-site Heisenberg antiferromagnet
    N = 14
    H = OpSum()
    for i in 1:N
        H += "J" * Op("SdotS", [i, mod1(i+1, N)])
    end
    H["J"] = 1.0 

    # define spin correlator at "half-chain" distance
    corr_op = OpSum()
    corr_op += Op("SdotS", [1, N ÷ 2])

    # define shift by one site
    T = Permutation(circshift(1:N, -1))

    # define cyclic group
    C_N = PermutationGroup([T^k for k in 0:(N-1)])

    # define irreps from character tables, labelled by momentum (2pi/N ×) k
    character_table = [ [C_N_character(N, k, p) for p in 0:(N-1)] for k in 0:(N-1)]
    irreps = [ Representation(C_N, character_table[k]) for k in 1:N ]   

    # check all translation symmetry blocks
    e0 = Inf
    psi0 = nothing
    for k in eachindex(irreps)
        block = Spinhalf(N, irreps[k]) 
        e0k, psik = eig0(H, block)
        if  e0k < e0
            e0 = e0k
            psi0 = psik
        end
    end

    # now the operator must be symmetrized!
    corr = inner(symmetrize(corr_op, C_N), psi0)

    println("Ground state energy: $e0")
    println("Ground state correlator: $corr")
    println("State vector length = ", size(psi0))

end

function C_N_character(N::Int, k::Int, p::Int)
    exp(im * 2 * pi * k * p / N)
end

main()