Discrete Lehmann Representation

Discrete Lehmann Representation#

Theory#

We explain how to transform expansion coefficients in IR to the discrete Lehmann representation (DLR) [1]. We model the spectral function as

\[ \rho(\omega) = \sum_{p=1}^L c_p \delta(\omega - \bar{\omega}_p), \]

where sampling frequencies \(\{\bar{\omega}_1, \cdots, \bar{\omega}_{L}\}\) are chosen to the extrema of \(V'_{L-1}(\omega)\) [2]. In the original paper [1], a more and different systematic way of choosing poles was used. We can transform data between \(\rho_l\) and \(c_p\) through the relation

\[ \rho_l = \sum_{p=1}^L \boldsymbol{V}_{lp} c_p, \]

where the matrix \(\boldsymbol{V}_{lp}~[\equiv V_l(\bar{\omega}_p)]\) is well-conditioned. As a result, in DLR, the Green’s function is represented as

\[ G(\mathrm{i}\omega_n) = \sum_{p=1}^L \frac{c_p}{\mathrm{i}\omega_n - \bar{\omega}_p}. \]

[1] J. Kaye, K. Chen, O. Parcollet, Phys. Rev. B 105, 235115 (2022).
[2] H. Shinaoka et al., arXiv:2106.12685v2.

We consider fermions and the semi circular DOS

\[ \rho(\omega) = \sqrt{\frac{2}{\pi}} \sqrt{1-\omega^2}. \]

The corresponding Green’s function is given by

\[ G(\mathrm{i}\omega_n) = \int_{-1}^1 \mathrm{d}\omega \frac{\rho(\omega)}{\mathrm{i}\omega_n - \omega}. \]

The Green’s function is expanded in IR as

\[ G_l = - S_l \int_{-1}^1 \mathrm{d}\omega \rho(\omega) V_l(\omega). \]

Below, we demonstrate how to transform \(G_l\) to the DLR coefficients \(c_p\).

using SparseIR
using Plots
gr() # USE GR backend
using LaTeXStrings

Implementation#

Create basis object#

wmax = 1.0
lambda_ = 1e+4
beta = lambda_/wmax

basis = FiniteTempBasis(Fermionic(), beta, wmax, 1e-15)
print(length(basis))

Setup model#

rho(omega) = sqrt(1-omega^2)/sqrt(0.5*π)

omega = LinRange(-wmax, wmax, 1000)
plot(omega, rho.(omega), xlabel=latexstring("\\omega"), ylabel=latexstring("\\rho(\\omega)"))

rhol = overlap(basis.v, rho)
ls = collect(1:length(basis))
plot(ls[1:2:end], abs.(rhol)[1:2:end], marker=:cross, yaxis=:log, ylims=(1e-5,1))

gl = - basis.s .* rhol
plot(ls[1:2:end], abs.(gl)[1:2:end], marker=:cross, ylims=(1e-5,10), yaxis=:log)

Create a DLR object and perform transformation#

dlr = DiscreteLehmannRepresentation(basis)

# To DLR
g_dlr = SparseIR.from_IR(dlr, gl)

plot(dlr.poles, g_dlr, marker=:cross)

# Transform back to IR from DLR
gl_reconst = SparseIR.to_IR(dlr, g_dlr)

plot(
    [abs.(gl), abs.(gl_reconst), abs.(gl-gl_reconst)],
    label=["Exact" "Reconstructed from DLR" "error"],
    marker=[:cross :x :circle], line=(nothing,nothing,nothing), yaxis=:log,
    ylims=(1e-18,10)
)

┌ Warning: Invalid negative or zero value 0.0 found at series index 4 for log10 based yscale
└ @ Plots ~/.julia/packages/Plots/rz1WP/src/utils.jl:106
┌ Warning: Invalid negative or zero value 0.0 found at series index 4 for log10 based yscale
└ @ Plots ~/.julia/packages/Plots/rz1WP/src/utils.jl:106
┌ Warning: Invalid negative or zero value 0.0 found at series index 4 for log10 based yscale
└ @ Plots ~/.julia/packages/Plots/rz1WP/src/utils.jl:106

Evaluation on Matsubara frequencies#

#v = 2 .* collect(-1000:10:1000) .+ 1
v = FermionicFreq.(2 .* collect(-1000:10:1000) .+ 1)
iv = SparseIR.valueim.(v, beta)

newaxis = [CartesianIndex()]
transmat = 1 ./ (iv[:,newaxis] .- dlr.poles[newaxis,:])
giv = transmat * g_dlr

smpl = MatsubaraSampling(basis; sampling_points=v)
giv_exact = evaluate(smpl, gl)

plot(
    imag.(iv), [imag.(giv_exact), imag.(giv)], marker=[:cross :x], line=(nothing,nothing),
    xlabel=latexstring("\\nu"),
    ylabel=latexstring("G(\\mathrm{i}\\omega_n)")
    )

┌ Warning: Sampling matrix is poorly conditioned (cond = 4.001560882675462e15).
└ @ SparseIR ~/.julia/packages/SparseIR/KrYmQ/src/sampling.jl:81