Dynamical spin structure factor
Author Martin Ulaga
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:
-
Compute the ground state \(|\Psi_0\rangle\) (is the appropriate symmetry sector, e.g., usually \({\bf k}=0\)).
-
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\).
-
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:
-
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.
-
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.