Green's function

Uses the Lanczos algorithm[1] to calculate the conductivity of the planar \(t\)-\(J\) model with one hole. See also the documentation page for the spin structure factor. The ground state conductivity is geven by the current autocorrelator

\[ T\sigma(\omega)=-i\int dt e^{-i\omega t}\langle J(t)J\rangle. \]

The plotting script includes a simple benchmark exploring the role of the number of Lanczos iterations on the result.

example code

JuliaC++

    using XDiag
    using HDF5

    function main()
      say_hello()

      # IO
      latticeInput = "../../misc/data/square.16.tJ.toml"
      lfile = FileToml(latticeInput)

      filename = "../../misc/data/examples_output/tJ_conductivity_jll.h5"
      outfile = h5open(filename, "w")

      # Lanczos parameters
      precision = 0.0 # turns off precision checking
      maxiters = 200

      # Define the model
      N = 16
      nup = 8
      ndn = 7 # one hole
      ops = read_opsum(lfile, "Interactions")
      t = 1.0
      J = 0.3
      ops["T"] = t
      ops["J"] = J
      @show(ops)
      irrep = read_representation(lfile, "Gamma.C1.A")

      # compute groundstate
      println("Computing ground state ...")
      block = tJ(N, nup, ndn, irrep)
      e0, gs = eig0(ops, block)
      println("done.")
      println("Ground state energy: $e0")
      outfile["e0"] = e0

      current = symmetrize(1.0 * im * Op("Hop", [1, 2]), irrep)
      Av = apply(current, gs)
      @show nrm = norm(Av)
      Av /= nrm

      res = eigvals_lanczos_inplace(ops, Av; neigvals=1, precision=precision, max_iterations=maxiters)
      @show res.eigenvalues[1]
      outfile["norm"] = nrm
      outfile["alphas"] = res.alphas
      outfile["betas"] = res.betas
      outfile["eigs"] = res.eigenvalues
    end

    main()

    #include <xdiag/all.hpp>

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

    int main() try {
      say_hello();
      set_verbosity(0);
      // IO
      std::string latticeInput =
          XDIAG_DIRECTORY "/misc/data/square.16.tJ.toml";
      auto lfile = FileToml(latticeInput);

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

      // Lanczos parameters
      auto precision = 0; // turns off convergence checking
      auto maxiters = 200;

      // Define the model
      int N = 16;
      int nup = 8;
      int ndn = 7; // one hole
      auto ops = read_opsum(lfile, "Interactions");
      auto t = 1.;
      auto J = 0.3;
      ops["T"] = t;
      ops["J"] = J;
      XDIAG_SHOW(ops);

      auto irrep = read_representation(lfile, "Gamma.C1.A");

      // compute groundstate
      Log("computing ground state");
      Log("Computing ground state ...");
      auto block = tJ(N, nup, ndn, irrep);
      XDIAG_SHOW(block);
      auto [e0, gs] = eig0(ops, block);
      Log("Ground state energy: {:.12f}", e0);
      outfile["e0"] = e0;

      // create & apply current operator
      Log("applying curr. operator");
      auto current = symmetrize(1i * Op("Hop", {0, 1}),
                                irrep);
      // XDIAG_SHOW(current);
      auto Av = apply(current, gs);
      auto nrm = norm(Av);
      Av /= nrm;
      XDIAG_SHOW(nrm);
      outfile["norm"] = nrm;

      // second run
      Log("second Laczos run on A|gs>");
      auto res = eigvals_lanczos_inplace(ops, Av, 1, precision, maxiters);
      XDIAG_SHOW(res.eigenvalues(0));
      outfile["alphas"] = res.alphas;
      outfile["betas"] = res.betas;
      outfile["eigs"] = res.eigenvalues;

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

plotting script

Julia

    using XDiag
    using LinearAlgebra
    using HDF5
    using Plots
    using IterTools

    # pythonplot()

    function poles_weights(outfile::String; M::Int64=-1)
        e0 = h5read(outfile, "e0")
        nrm, alphas, betas, eigs = h5read.((outfile,), ["norm", "alphas", "betas", "eigs"])
        if M > 0
            try
                alphas = alphas[1:M]
                betas = betas[1:M]
            catch e
                println(e)
            end
        end
        tmat = SymTridiagonal(alphas, betas)
        es, evecs = eigen(tmat)
        # @show es[1], eigs[1]

        poles = es .- e0
        weights = evecs[1, :] .^ 2 .* nrm^2
        sortdisp(poles, weights/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

    function sortdisp(p,w;max=10)
        upto=min(max,length(p))
        display([p[sortperm(w,rev=true)][1:upto] sort(w,rev=true)[1:upto]])
    end

    let
        source = "../../misc/data/examples_output/tJ_conductivity_jll.h5"
        omegas = collect(range(0., 10.0, length=40000))
        # eta = omegas[2]-omegas[1]
        eta = 0.1
        fig = plot()
        for M in [30, 50, 100, 200]
            p, w = poles_weights(source; M=M)
            s = spectrum(p, w, omegas, eta)
            plot!(omegas, s, label="M=$M")
        end
        savefig(fig, "cond_jll.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.