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 -{\overline{\omega }}_p),$$

where sampling frequencies $\{{\overline{\omega }}_1,\cdots ,{\overline{\omega }}_L\}$ are chosen to the extrema of $V_{L-1}^{\prime }(\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{\mathit{V}}_{lp}{c}_p,$$

where the matrix ${\mathit{V}}_{lp}[\equiv V_l({\overline{\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-{\overline{\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))

104

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