Lattice Momentum Tower of States for the Heisenberg Chain

This example demonstrates how to set up a tower of states analysis with respect to translation symmetries using the example of a spin-\(1/2\) chain

\[ H = J \sum_{\langle i, j \rangle} \bm{S}_i \cdot \bm{S}_j, \]

with periodic boundary conditions and \(0 < J\) for ferromagnetic interactions.

The code below demonstrates how the lowest energy levels for each value of the lattice momentum can be computed using the Lanczos algorithm. The Julia version also contains a simple plotting method, leading to the following tower of states with respect to the lattice momentum.

As can be seen in the resulting figure, the tower of states nicely captures the dispersion relation of spinon excitations (e.g. flipping one spin relative to the ground state creates two spinons both carrying spin-\(1/2\)).

This examples also demonstrates that the ground state is not necessarily found in the trivial representation of the translation symmetry, as the ground state resides at lattice momentum \(k=\pi\) in the figure. In fact, the ground state alternates between \(k=0\) and \(k=\pi\) for \(N\in 2(2\mathbb{N}+1)\) and \(N \in 4 \mathbb{N}\), respectively.

C++Julia

#include <vector>
#include <complex>
#include <math.h>
#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{

        // specify system parameters
        int N = 18;
        double J = 1.0;

        // fix number of up spins
        int Nup = N/2;

        // build Hamiltonian
        auto H = OpSum();
        for (int i=0; i<N; i++){
            H += "J" * Op("SdotS", {i, (i+1)%N});
        }
        H["J"] = J;

        // set up translation symmetry 
        // step 1: 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);

        // step 2: 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); 

        // step 3: 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);
        };

        // sort states in Sz_tot = 0 sector by translation irrep
        int lvl_per_block = 5;
        std::vector<double> TOS_x_vals(N*lvl_per_block); 
        std::vector<double> TOS_y_vals(N*lvl_per_block); 
        for (int k = 0; k < N; k++){
            auto sym_block = Spinhalf(N, Nup, irreps[k]);
            // set neigvals = ... for reasonable convergence of lowest eigenvalues
            int neigvals = 3;
            auto lanczos_result = eigvals_lanczos(H, sym_block, neigvals);
            for (int i=0; i<lvl_per_block; i++){
                TOS_x_vals[k*lvl_per_block + i] = (2 * M_PI * k) / N;
                TOS_y_vals[k*lvl_per_block + i] = lanczos_result.eigenvalues[i];
            }
        };

        // use the vectors TOS_x_vals and TOS_y_vals to create e.g. the TOS plot ...
        // ... for a simple plotting method check the Julia version of this example

        return 0;

    }catch(Error e){
        error_trace(e);
    }

}

using XDiag
using Plots 

function main()

    # specify system parameters
    N = 18
    J = 1.0

    # fix number of up spins
    Nup = N÷2

    # build Hamiltonian
    H = OpSum()
    for i in 1:N
        H += "J" * Op("SdotS", [i,mod1(i+1, N)])
    end
    H["J"] = J

    # set up translation symmetry
    T = Permutation(circshift(1:N, -1))
    C_N = PermutationGroup([T^k for k in 0:(N-1)])
    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 ]   

    # sort states in Sz_tot = 0 sector by translation irrep
    lvl_per_block = 5
    TOS_x_vals = Vector{Float64}(undef, N*lvl_per_block)
    TOS_y_vals  = Vector{Float64}(undef, N*lvl_per_block)
    for k in eachindex(irreps)
        sym_block = Spinhalf(N, Nup, irreps[k])
        # set neigvals = ... for reasonable convergence of lowest eigenvalues
        lanczos_result = eigvals_lanczos(H, sym_block, neigvals=3) 
        for i in 1:lvl_per_block
            TOS_x_vals[(k-1)*lvl_per_block + i] = 2*pi*(k-1)/N
            TOS_y_vals[(k-1)*lvl_per_block + i] = lanczos_result.eigenvalues[i]
        end
    end

    # create TOS plot
    scatter(
        TOS_x_vals,
        TOS_y_vals,
        title="Tower of States, AFM Heisenberg Chain, N=$N\n",
        xlabel="lattice momentum k",
        ylabel="lowest eigenstates",
        legend=false,
        xticks = ([0:π/2:2*π;], ["0", "\\pi/2", "\\pi", "3\\pi/2", "2\\pi"]),
        right_margin=5Plots.mm,
        dpi=300)
end


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


main()