lxvm
commented
1 month ago Thanks for your insights on this problem @jasonkaye . I agree that the "Fermi sea" causes difficulties for integration due to the infinite limits. In most cases this semi-infinite integral is handled automatically with a change of variables, which works reliably except when unless the self energy is limited to a frequency interval. (We can work around this by interpolating $\Sigma$ using AAA to get a representation on the real axis.) However, this doesn't work for the real part of the Green's function, whose integral doesn't converge since the decay is too slow. The consequence of this is that theElectronDensitySolver only integrates the imaginary part of G, which could also lead to some peak-missing issues. Additionally, truncation to a finite interval $[a, b]$ could also cause peak missing if the chosen interval is too large. For reference, the upper limit always gets truncated according to the decay of the Fermi function.
Adding an option to use the DLR using Lehmann.jl sounds great and can be implemented in addition to the existing code using the same API, or we can add something new if we want to be able to return the real part of the Green's function as well.
In the Fermi liquid regime, $\eta \sim T^2$, so the cost scaling of adaptive integration is $\sim \log(1/\eta) ~ \log(1/T)$ and the DLR rank scales as $\sim \log(1/T)$, so the computational complexities are comparable.
There is a potential problem in computing the electron density via the ElectronDensitySolver function: if the density of states was only computed on a domain $[a,b]$ which is not sufficient to cover its full support, the integral from $\int{-\infty}^{\infty}$ in the definition of the electron density will be truncated to $\int{a}^{b}$, yielding some error. This is particularly problematic for the $\int_{-\infty}^{a}$ part, which doesn't benefit from the exponential decay of the Fermi function. At high temperatures, this requires computing the density of states on a very large domain, which can be challenging and requires careful by-hand convergence with respect to the domain size.
I believe there is a more efficient and robust approach, which is to use the fact that the electron density is given, up to factors of $\beta$, by the sum of the Green's function $G$ over all Matsubara frequencies. Using the discrete Lehmann representation (DLR), it is sufficient to obtain the values of the Green's function at the DLR Matsubara frequency nodes. Then, it can be shown that the Matsubara sum is given, again up to factors of $\beta$, by the sum of the DLR coefficients, which can be obtained from its values at the DLR nodes. This can be done efficiently, to arbitrary accuracy. Now, instead of computing the Green's function at many real frequency points to obtain the density of states, one need only compute it at a small number of Matsubara frequency points. This yields a smaller number of integrals, each of which is easier to compute owing to some additional smearing coming from the replacement of real frequencies with imaginary freqencies. It replaces the issue of truncating the density of states by a DLR spectral cutoff parameter $\Lambda$, which can straightforwardly be increased until convergence is reached (resulting in a very slowly-growing increase in the number of DLR nodes).
In practice, I think the easiest way to do this is using Lehmann.jl. It should be quite straightforward, as it only requires the simplest functionality of the DLR. I'd be happy to collaborate on this to get it up and running.