Matsubara Basis

Finite-temperature GF2 and GW calculations store functions on the imaginary-time interval $0\leq \tau\leq\beta$ or on the Matsubara-frequency axis. These functions are smooth compared with real-frequency spectra, but they still carry information from a large real-energy window. A useful numerical representation must therefore resolve low-energy features near the chemical potential, high-energy structure from core or virtual states, and the high-frequency moments needed for accurate densities and energies.

The connection to real frequencies can be written as a spectral representation. For a fermionic Green’s function, a common form is

Here $\rho(\omega)$ is the spectral density, $\beta$ is the inverse temperature, and $\omega_{\max}$ is an energy cutoff large enough to contain the relevant spectral weight. Bosonic response functions have a closely related kernel with Bose statistics. The dimensionless difficulty of the representation is controlled mainly by $\Lambda=\beta\omega_{\max}$: low temperatures and wide energy windows require more information.¹

Uniform imaginary-time or Matsubara-frequency grids are inefficient for this problem. A uniform Matsubara grid spends many points at high frequencies where the function is smooth, while still needing enough points near zero frequency and enough tail information to preserve moments. Early realistic finite-temperature calculations therefore used more compact grids and transforms: Legendre-polynomial representations, orthogonal-polynomial transform grids, and cubic-spline interpolation on sparse Matsubara grids all reduced the number of stored points by exploiting smoothness and known high-frequency behavior.²

Chebyshev Basis

The Chebyshev representation starts from a simple map from imaginary time to the interval $[-1,1]$, $x(\tau)=2\tau/\beta-1$. A matrix element of a Green’s function or self-energy is then expanded in Chebyshev polynomials of the first kind,

The prime on the sum denotes the conventional half weight of the zeroth coefficient. In practice the coefficients can be obtained from values at Chebyshev nodes by a discrete cosine transform, and the function can be evaluated stably by recurrence relations.²

The main advantage of the Chebyshev basis is mathematical control. Chebyshev approximation theory gives exponential convergence for analytic functions on the interval, and the operations needed in many-body calculations - evaluation, convolution, Fourier transformation, and even solving Dyson-like equations - can be formulated in coefficient space.² The method is also transparent: if the high-order coefficients stop decreasing, the basis is not yet converged.

The disadvantage is compactness. The Chebyshev basis is general and robust, but it is not adapted to the analytic-continuation kernel. Large $\Lambda$ values can require hundreds of basis functions in realistic calculations.

Intermediate Representation

The intermediate representation (IR), introduced by Shinaoka and collaborators, builds the basis from the spectral kernel itself rather than from a fixed family of polynomials. With $x=2\tau/\beta-1$, $y=\omega/\omega_{\max}$, and $\Lambda=\beta\omega_{\max}$, the fermionic kernel becomes

The IR basis is defined by a singular-value decomposition of this kernel,

The functions $u_l$ form an orthonormal imaginary-time basis, $v_l$ form the corresponding real-frequency basis, and the singular values $s_l$ decay rapidly. An imaginary-time function and the corresponding spectral density can be expanded as $G(\tau)=\sum_l g_l u_l(x(\tau))$ and $\rho(\omega)=\sum_l \rho_l v_l(\omega/\omega_{\max})$, with coefficients related by the singular values of the kernel.¹

This construction is numerical: the IR basis functions are not elementary polynomials and are usually generated by libraries such as irbasis or sparse-ir. The payoff is compactness. The basis is adapted to the map between imaginary time and real frequency, so it generally needs fewer coefficients than Chebyshev for the same $\Lambda$ and target accuracy. The original IR paper also shows that the Legendre basis appears as a limiting case at small $\Lambda$, while the IR becomes increasingly advantageous as $\Lambda$ grows.¹, ³

The numerical nature of IR changes the workflow. One chooses $\Lambda$ and an accuracy, obtains precomputed or generated basis functions, and checks convergence in the number of retained singular functions. This is less elementary than a polynomial expansion, but it is often a better match to the Matsubara problem encountered in finite-temperature ab initio calculations.⁴

Sparse Sampling and Transforms

A compact basis is only useful if one can move accurately between coefficient space, imaginary-time values, and Matsubara-frequency values. The sparse-sampling approach addresses this directly. For a basis $F_l^\alpha(\tau)$, where $\alpha$ denotes fermionic or bosonic statistics, one writes

Sparse time points $\bar\tau_k$ and sparse Matsubara frequencies $i\bar\omega_k$ are chosen so that the basis matrices are well conditioned. Then the same coefficients can be recovered from either domain,

For Chebyshev, the time nodes are tied to roots of Chebyshev polynomials. For IR, the time and frequency sampling points are chosen from the structure of the numerical basis functions, and the Matsubara points are often distributed almost logarithmically from low to high frequency.⁵

This sparse-sampling layer is what makes the basis practical in GF2 and GW calculations. Diagrams are naturally evaluated in imaginary time, Dyson-like equations are naturally solved in Matsubara frequency, and compact coefficients provide the bridge between those two domains. The same idea also handles the change between fermionic quantities such as $G$ and $\Sigma$ and bosonic quantities such as polarizations and screened interactions.⁵

DLR Basis

The discrete Lehmann representation (DLR) starts from the same spectral-kernel viewpoint but chooses a different basis. Instead of using orthonormal singular functions of the kernel, DLR selects a small set of real-frequency nodes $\omega_l$ and uses the corresponding kernel columns as basis functions:

The selected frequencies, imaginary-time nodes, and Matsubara-frequency nodes are obtained by rank-revealing linear algebra such as pivoted QR or interpolative decompositions. The rank $r$ depends mainly on $\Lambda=\beta\omega_{\max}$ and a requested tolerance $\epsilon$.⁶

DLR is less orthogonal than IR, but it is very explicit. The basis functions are exponentials in imaginary time and simple rational functions on the Matsubara axis, so transforms between time and frequency can be written analytically once the DLR coefficients are known. Library work such as libdlr and cppdlr packages these grids, interpolation problems, and transforms for practical use.⁶, ⁷, ⁸

How to Choose

All of these bases solve the same problem: represent smooth imaginary-axis functions with far fewer degrees of freedom than a uniform grid. Legendre and spline approaches are simple, controlled, and historically important. Chebyshev is especially attractive when one wants direct polynomial approximation with clear convergence diagnostics. IR is usually more compact because it is adapted to the spectral kernel. Sparse sampling turns Chebyshev or IR coefficients into practical time and Matsubara grids. DLR uses the same kernel information in an interpolation-oriented form that gives simple basis functions and direct transforms.

In Green calculations the best basis is therefore a numerical choice, not a change in physics. The Hamiltonian, density matrix, and static term $\Sigma^\infty$ are unchanged. The basis controls how accurately and efficiently the dynamical functions on the imaginary axis are stored, transformed, and iterated.