Skip to content

Dynamical structure factor

Uses the Lanczos algorithm[1] to calculate the dynamical spin spectral function of the spin-1/2 Heisenberg chain. The structure factor is defined as

$$ S^{zz}({\bf k},\omega) = \int dt e^{-i\omega t}\langle S^z_{\bf k}(t)S^z_{\bf -k}\rangle, $$ where \(S^z_{\bf k}=\frac{1}{\sqrt{N}}\sum_i e^{-i\bf{k\cdot r}_i}S^z_i\) and \(S^z_i\) is a spin operator on site \(i\). The calculation is done using the Lehmann representation

$$ S^{zz}({\bf k},\omega)=\frac{1}{N}\sum_{m=1}^M |\langle\Psi_0|S^z_{\bf k}|\psi_m\rangle|^2\delta(\omega-\epsilon_m+E_0), $$ where \(|\Psi_0\rangle\) is the ground state which needs to be computed. The algorithm proceeds in 3 steps:

  1. Compute the ground state \(|\Psi_0\rangle\) (is the appropriate symmetry sector, e.g., usually \({\bf k}=0\)).

  2. Find the operator \(S^z_{\bf k}\), e.g. symmetrize with respect to the appropriate irrep, and calculate \(|\tilde{\Psi}_0\rangle=S^z_{\bf k}|\Psi_0\rangle\).

  3. Rerun the Lanczos algorigthm using the normalized state \(|\Psi_1\rangle=|\tilde{\Psi}_0\rangle/\sqrt{\langle\tilde{\Psi}_0|\tilde{\Psi}_0\rangle}\).

main example code

    using Pkg;
    Pkg.instantiate();
    using XDiag
    using HDF5

    function main()
      say_hello()
      filename = "../../misc/data/examples_output/spinhalf_chain_structure_factor.h5"
      outfile = h5open(filename, "w")

      # define model
      N = 16
      nup = div(N, 2)
      ops = OpSum()
      for i in 1:N
        ops += "J" * Op("SdotS", [i, mod1(i + 1, N)])
      end
      ops["J"] = 1.0
      @show ops

      # Create the permutation group
      @show perm = Permutation([mod1(s + 1, N) for s in 1:N])
      @show group = PermutationGroup([perm^s for s in 0:N-1])

      # Create the irreps at momenta k
      irreps = Representation[]
      for k in 0:N-1
        phase = exp(2im * pi * k / N)
        characters = [phase^s for s in 0:N-1]
        irrep = Representation(group, characters)
        # @show irrep
        push!(irreps, irrep)
      end

      # compute groundstate (known to be at k=0)
      println("Computing ground state ...")
      block = Spinhalf(N, nup, irreps[1])
      e0, gs = eig0(ops, block)
      make_complex!(gs)
      println("done.")
      println("Ground state energy: $e0")
      outfile["e0"] = e0

      # loop through momenta
      for k in 1:N
        println("coumputing structure factor at q=$k")
        S_q = symmetrize(Op("Sz", 1), irreps[k])
        Av = apply(S_q, gs)
        nrm = norm(Av)
        Av /= nrm

        res = eigvals_lanczos_inplace(ops, Av)
        outfile["$k/norm"] = nrm
        outfile["$k/alphas"] = res.alphas
        outfile["$k/betas"] = res.betas
        outfile["$k/eigs"] = res.eigenvalues
      end
    end

    main()

    #include <xdiag/all.hpp>

    using namespace xdiag;
    using namespace std::complex_literals;
    using fmt::format;

    constexpr double XDIAG_PI = 3.14159265358979323846;

    int main() try {
      say_hello();

      // IO
      std::string filename = XDIAG_DIRECTORY
          "/misc/data/examples_output/spinhalf_chain_structure_factor.h5";
      auto outfile = FileH5(filename, "w!");

      // Define the nearest-neighbor Heisenberg model
      int N = 16;
      int nup = N / 2;
      OpSum ops;
      for (int i = 0; i < N; ++i) {
        ops += "J" * Op("SdotS", {i, (i + 1) % N});
      }
      ops["J"] = 1.0;

      // Create the permutation group
      std::vector<int> translation;
      for (int s = 0; s < N; ++s) {
        translation.push_back((s + 1) % N);
      }
      Permutation perm(translation);
      std::vector<Permutation> perms;
      for (int s = 0; s < N; ++s) {
        perms.push_back(pow(perm, s));
      }
      auto group = PermutationGroup(perms);

      // Create the irreps at momenta k
      std::vector<Representation> irreps;
      for (int k = 0; k < N; ++k) {
        complex phase = exp(2i * XDIAG_PI * (double)k / (double)N);
        std::vector<complex> characters;
        for (int s = 0; s < N; ++s) {
          characters.push_back(pow(phase, s));
        }
        auto irrep = Representation(group, characters);
        irreps.push_back(irrep);
      }

      // compute groundstate (known to be at k=0)
      Log("Computing ground state ...");
      auto block = Spinhalf(N, nup, irreps[0]);
      auto [e0, gs] = eig0(ops, block);
      gs.make_complex();
      Log("done.");
      Log("Ground state energy: {:.12f}", e0);
      outfile["e0"] = e0;

      // loop through momenta
      for (int k = 0; k < N; ++k) {
        Log("S^zz(k,w) at momentum {}", k);
        auto S_q = symmetrize(Op("Sz", 0), irreps[k]);
        auto Av = apply(S_q, gs);
        auto nrm = norm(Av);
        Av /= nrm;
        XDIAG_SHOW(nrm);

        auto res = eigvals_lanczos_inplace(ops, Av);
        outfile[format("{}/norm", k)] = nrm;
        outfile[format("{}/alphas", k)] = res.alphas;
        outfile[format("{}/betas", k)] = res.betas;
        outfile[format("{}/eigs", k)] = res.eigenvalues;
      }

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

visualization script

Postprocessing is here split into two parts:

  1. Compute the poles/weights of the spectral function. The \(m\)-th pole appears at frequency \(\omega = \epsilon_m -E_0\) with the associated weights given by the norm of \(|\tilde{\Psi}_0\rangle\) and the projection \(w_m=\langle \tilde{\Psi}_1|\psi_m\rangle\). The latter is obtained from the first eigenvector of the tridiagonal matrix.

  2. Spread out the poles and weights using a Gaussian kernel of width \(\eta\) onto some frequency interval, artificially broadening the \(\delta\)-functions.

   using XDiag
   using LinearAlgebra
   using HDF5
   using Plots

   function plot_sf(outfile::String, omegas::Vector{Float64}, eta::Float64 = 0.15)
       alls = Matrix{Float64}(undef, 16, length(omegas))
       for k in 0:15
           p, w = poles_weights(outfile, k)
           s = spectrum(p, w, omegas, eta)
           alls[k+1, :] = s
       end
       return heatmap(collect(range(0, 15, step=1)), omegas, alls',
           xlabel="\$k\$", ylabel="\$\\omega\$",
           interpolation=:false)
   end

   function poles_weights(outfile::String, k::Int64)
       e0 = h5read(outfile, "/e0")
       nrm, alphas, betas, eigs = h5read.((outfile,), ["/$k/norm", "/$k/alphas", "/$k/betas", "/$k/eigs"])
       tmat = SymTridiagonal(alphas, betas)
       es, evecs = eigen(tmat)
       # es == eigs

       poles = es .- e0
       weights = evecs[1, :] .^ 2 .* nrm^2
       return poles, weights
   end

   function spectrum(poles::Vector{Float64}, weights::Vector{Float64},
       omegas::Vector{Float64}, eta::Float64=0.1; cutoff::Float64=0.0, filter::Float64=1e-6)
       fpoles = [poles[i] for i in 1:length(weights) if abs(weights[i]) > filter]
       fweights = [weights[i] for i in 1:length(weights) if abs(weights[i]) > filter]
       diffs = omegas .- fpoles'
       gaussians = exp.(-(diffs ./ (2.0 * eta)) .^ 2.0) ./ (eta * sqrt(2.0 * pi))
       if abs(cutoff) > 1e-12
           idxs = findall(abs.(fpoles) .> cutoff)
           s = zeros(length(omegas))
           for idx in idxs
               s += gaussians[:, idx] * fweights[idx]
           end
           return s
       else
           return gaussians * fweights
       end
   end

   let
       source = "../../misc/data/examples_output/spinhalf_chain_structure_factor.h5"
       omegas = collect(range(0.0, 3.0, length=1000))
       eta = 10 * (omegas[2] - omegas[1])
       fig = plot_sf(source, omegas, eta)
       savefig(fig, "spekt.pdf")
   end

references

[1] Prelovšek, P., & Bonča, J. (2013). Ground state and finite temperature Lanczos methods. Strongly Correlated Systems: Numerical Methods, 1-30.