Discrete Lehmann representation of imaginary time Green's functions

We present an efficient basis for imaginary time Green's functions based on a low rank decomposition of the spectral Lehmann representation. The basis functions are simply a set of well-chosen exponentials, so the corresponding expansion may be thought of as a discrete form of the Lehmann representation using an effective spectral density which is a sum of $\delta$ functions. The basis is determined only by an upper bound on the product $\beta \omega_{\max}$, with $\beta$ the inverse temperature and $\omega_{\max}$ an energy cutoff, and a user-defined error tolerance $\epsilon$. The number $r$ of basis functions scales as $\mathcal{O}\left(\log(\beta \omega_{\max}) \log (1/\epsilon)\right)$. The discrete Lehmann representation of a particular imaginary time Green's function can be recovered by interpolation at a set of $r$ imaginary time nodes. Both the basis functions and the interpolation nodes can be obtained rapidly using standard numerical linear algebra routines. Due to the simple form of the basis, the discrete Lehmann representation of a Green's function can be explicitly transformed to the Matsubara frequency domain, or obtained directly by interpolation on a Matsubara frequency grid. We benchmark the efficiency of the representation on simple cases, and with a high precision solution of the Sachdev-Ye-Kitaev equation at low temperature. We compare our approach with the related intermediate representation method, and introduce an improved algorithm to build the intermediate representation basis and a corresponding sampling grid.