Background

This page gives a brief review of the discrete Lehmann representation (DLR), and establishes the definitions and conventions used in cppdlr (which can vary from one reference to another). If you are already familiar with the DLR, you should probably still read the section on the relative imaginary time format below. For a more detailed description of the DLR, please see the references listed on the main page. For examples of the implementation of these concepts in cppdlr, please see the examples page.

Imaginary time Green’s functions and the Lehmann representation

The single-particle imaginary time Green’s function is defined in terms of the time-ordered expectation value as

\[G_{ab}(\tau) = - \langle \mathcal{T} c_a(\tau) c_b^\dagger(0) \rangle,\]

where \(c^\dagger_b(0)\) is the creation operator for a particle in state \(b\) at time \(0\) and \(c_a(\tau)\) is the annihilation operator for a particle in state \(a\) at time \(\tau\). The Green’s function is defined on the interval \(\tau \in (0, \beta)\), where \(\beta\) is the inverse temperature, but it can be extended to \(\tau \in (-\beta, 0)\) using the periodicity and anti-periodicity properties

\[G_{ab}(\tau) = \xi G_{ab}(\beta + \tau),\]

of bosonic (\(\xi = 1\)) and fermionic (\(\xi = -1\)) particles, respectively.

The spectral Lehmann representation of a Green’s function is given by

\[G(\tau) = \int_{-\infty}^\infty K(\tau,\omega) \rho(\omega) \, d\omega,\]

where \(\rho(\omega)\) is the spectral function corresponding to the Green’s function \(G\), and \(K(\tau,\omega)\) is the analytic continuation kernel, given by

\[K(\tau, \omega) = -\frac{e^{-\omega \tau}}{1 + e^{-\beta \omega}}.\]

Taking the Fourier transform to the imaginary (or Matsubara) frequency domain gives

\[G(i \nu_n) = \int_{-\infty}^\infty K(i \nu_n,\omega) \rho(\omega) \, d\omega,\]

with

\[K(i \nu_n, \omega) = \frac{1}{i\nu_n - \omega}\]

for fermionic Green’s functions, and

\[K(i \nu_n, \omega) = \frac{\tanh (\beta \omega/2)}{i\nu_n - \omega}\]

for bosonic Green’s functions. Here, the Matsubara frequencies are defined as \(i \nu_n = 2 n \pi i/\beta\) and \(i \nu_n = (2n+1) \pi i/\beta\) for bosonic and fermionic Green’s functions, respectively.

Discrete Lehmann representation

The discrete Lehmann representation (DLR) is constructed by making a low-rank approximation of the analytic continuation kernel (using the interpolative decomposition). Let us define the dimensionless cutoff paramter \(\Lambda \equiv \beta \omega_{\max}\), where \(\omega_{\max}\) is defined such that \(\rho(\omega) = 0\) outside of \([-\omega_{\max},\omega_{\max}]\). In practice, \(\beta\) is typically known, and \(\omega_{\max}\) can be estimated. \(\Lambda\) is a user-specified parameter, and in the typical case that \(\omega_{\max}\) is not known exactly, results can be converged with respect to \(\Lambda\). Given \(\Lambda\) and another user-specified error tolerance parameter \(\epsilon\), the DLR expansion of an imaginary time Green’s function \(G(\tau)\) is given by

\[\begin{equation} G(\tau) \approx \sum_{l=1}^r K(\tau,\omega_l) \widehat{g}_l, \label{dlrexp} \tag{1} \end{equation}\]

with equality to accuracy \(\epsilon\) in a suitable norm (discussed below). Here, the \(r\) DLR frequencies \(\omega_l\) determining the DLR basis functions \(K(\tau,\omega_l)\) are carefully chosen (by a pivoted Gram-Schmidt procedure) depending only on \(\Lambda\) and \(\epsilon\), but not on \(G\) itself. As for the closely related intermediate representation (implemented, for example, in sparse-ir), for which the basis functions are orthogonal but non-explicit, we have \(r = \mathcal{O}(\log(\Lambda) \log(\epsilon^{-1}))\). Thus the DLR enables a highly efficient and high-order accurate discretization of all imaginary time Green’s functions, with a number of degrees of freedom independent of the specific structure of a given Green’s function (beyond its cutoff parameter \(\Lambda\)).

Constructing a DLR expansion

In practice, the DLR coefficients \(\widehat{g}_l\) must be determined from some samples of \(G(\tau)\). This can be done by fitting to data, e.g. via ordinary least squares, or by interpolation at the DLR imaginary time nodes \(\tau_k\). These are \(r\) imaginary time points which, given the DLR frequencies \(\omega_l\), are also chosen by a pivoted Gram-Schmidt procedure. In particular, given the values \(G(\tau_k)\), we can solve the \(r \times r\) linear system (or interpolation problem)

\[G(\tau_k) = \sum_{l=1}^r K(\tau_k,\omega_l) \widehat{g}_l\]

to obtain the DLR coefficients. \(G(\tau)\) can then be evaluated using its DLR expansion \(\eqref{dlrexp}\).

DLR in the Matsubara frequency domain

Fourier transform of \(\eqref{dlrexp}\) yields

\[\begin{equation} G(i \nu_n) \approx \sum_{l=1}^r K(i \nu_n,\omega_l) \widehat{g}_l, \label{dlrexp_imfreq} \tag{2} \end{equation}\]

so we see that the DLR expansion can be evaluated directly in imaginary time or imaginary frequency, i.e. the Fourier transform is performed analytically. As in imaginary time, the DLR coefficients can be obtained by solving the \(r \times r\) interpolation problem

\[G(i \nu_{n_k}) = \sum_{l=1}^r K(i \nu_{n_k},\omega_l) \widehat{g}_l\]

at the \(r\) DLR imaginary frequency nodes \(i \nu_{n_k}\), whereupon \(G(i \nu_n)\) can be evaluated using \(\eqref{dlrexp_imfreq}\) (or \(G(\tau)\) can be evaluated using \(\eqref{dlrexp}\)).

Operations in the DLR basis

Since the DLR basis functions are known analytically, common linear operations can be straightforwardly performed by representing them in the DLR basis. These include

Fourier transform: as explained above, one can switch between imaginary time and imaginary frequency representations via the DLR expansion, with no additional Fourier transform operation
Products: in imaginary time or imaginary frequency, by simply multiplying the functions on the DLR grid, i.e. \(H(\tau_k) = F(\tau_k) G(\tau_k)\), whereupon the DLR expansion of the result can be recovered
Imaginary time convolution: this includes the full convolution \(H(\tau) = \int_0^\beta F(\tau-\tau') G(\tau') \, d\tau'\), which requires using the periodicity/anti-periodicity condition, or the time-ordered convolution \(H(\tau) = \int_0^\tau F(\tau-\tau') G(\tau') \, d\tau\)
Linear functionals: e.g. inner products with a given function, evaluation at a point, etc…

All such operations take the form of vectors/matrices/tensors acting on \(r \times 1\) vectors, which represent the DLR expansion of a Green’s function \(G\) (either the vector of DLR coefficients of \(G\), or the vector of values of \(G\) at the DLR nodes). Common operations are implemented in cppdlr in a user-friendly manner, and the implementation of new operations should be requested on the GitHub issues page.

Imaginary time point format

First, in cppdlr imaginary time points are scaled from the interval \([0,\beta]\) to the interval \([0,1]\). This is because cppdlr works with dimensionless variables whenever possible, so in many functions it is unnecessary to specify the inverse temperature \(\beta\) explicitly.

Second, cppdlr stores imaginary time points in a peculiar manner, called the relative time format. This is a subtle issue which ``cppdlr`` users should be aware of, in particular if one wants to supply imaginary time points at which to evaluate a DLR expansion. For the TLDR, skip to the guidelines below. For an even more detailed discussion of this issue than the one given here, see Appendix C of this paper.

The relative time format works as follows. Points \(\tau \in [0, 0.5]\) are represented normally. However, instead of representing points \(\tau \in (0.5,1)\) directly, we instead store the number \(\tau^* = \tau-1\). In other words, we store the negative distance of \(\tau\) to 1, rather than tau itself. Recovering \(\tau\) in the standard absolute time format is straightforward, and is implemented by the function rel2abs.

The reason for this has to do with maintaining full relative accuracy in floating point arithmetic. To evaluate the kernel \(K(\tau,\omega)\), we sometimes need to compute the value \(1-\tau\) for \(\tau\) very close to 1. If we work with tau directly, there is a loss of accuracy due to catastrophic cancellation, which begins to appear in extreme physical regimes and at very high requested accuracies. If we instead compute \(\tau^*\) to full relative accuracy and work with it directly rather than with \(\tau\), for example by exploiting symmetries of \(K(\tau,\omega)\) to avoid ever evaluating \(1-\tau\), we can maintain full relative accuracy.

This annoyance is the price of maintaining full accuracy in floating point arithmic. But it is largely ignoreable if the loss of accuracy is not noticeable in your application, as will be the case for many users.

Simply follow these guidelines:

Use functions provided by cppdlr to carry out all imaginary time operations whenever possible. This will usually hide this technical complication.
In a situation in which you want to provide a point \(\tau\) at which to evaluate a DLR, there are two options:
- (The power user option) Compute \(\tau^*\), defined above, to full relative accuracy, and provide this according to the instructions in the relevant functions, thereby maintaining full relative accuracy in calculations, or
- (The typical user option) If you don’t care about the (usually minor) digit loss which comes from ignoring this subtlety, simply convert your point \(\tau\) in the standard, absolute format (a point \(\tau \in [0,1]\)) to the relative format \(\tau^*\) defined above using the abs2rel function. Since the point will have started its life in the absolute format, converting it to relative format cannot recover full relative accuracy, but it still needs to be converted in order to be compatible with cppdlr subroutines.
If you happen to want to evaluate a Green’s function on an equispaced grid on \([0,1]\) in imaginary time, use the function eqpts_rel to generate this grid in the relative format.

Matsubara frequency point format

We define the Matsubara (or imaginary) frequency points as \(i \nu_n = (2 n + 1) \pi i/\beta\) for fermionic Green’s functions, and \(i \nu_n = 2 n \pi i/\beta\) for bosonic Green’s functions. In cppdlr, Matsubara frequency points are represented by specifying the integer n, the inverse temperature \(\beta\), and whether the point is a fermionic or bosonic Matsubara frequency using the statistic_t specifier.

Measuring error

In what sense is the error of the DLR guaranteed to be less than \(\epsilon\)? Strictly speaking, if one obtains the DLR expansion by interpolation at the DLR nodes or by fitting on a finite number of grid points, there is no guarantee. However, in practice we observe very well-controlled errors, and we can give a more useful and practical answer. We consider the case of scalar-valued Green’s functions for simplicity. Let us define the \(L^2(\tau)\) norm as follows:

\[||G||_{L^2(\tau)}^2 = \frac{1}{\beta} \int_0^\beta d\tau \, |G(\tau)|^2.\]

By Parseval’s theorem, this is equal to the \(l^2(i \nu_n)\) norm:

\[||G||_{l^2(i \nu_n)}^2 = \frac{1}{\beta^2} \sum_{n=-\infty}^\infty |G(i \nu_n)|^2.\]

Ignoring errors at the level of machine precision, it can be shown that given an imaginary time Green’s function \(G(\tau)\) with spectral function \(\rho(\omega)\), there exists a DLR expansion \(G_{\text{DLR}}(\tau) \approx \sum_{l=1}^r K(\tau,\omega_l) \widehat{g}_l\) such that

\[||G - G_{\text{DLR}}||_{L^2(\tau)} = ||G - G_{\text{DLR}}||_{l^2(i \nu_n)} \leq ||\rho||_{L^1(\omega)} \epsilon.\]

Here, \(||\rho||_{L^1(\omega)} = \int_{-\infty}^\infty d\omega \, |\rho(\omega)|\), which is typically \(1\). Note: The inequality holds as long as the DLR frequencies \(\omega_l\) are constructed using the pivoted Gram-Schmidt procedure on a suitably fine discretization of \(K(\tau, \omega)\) with rows weighted by appropriate quadrature weights. The latter point is not discussed in the original reference on the DLR, but will be in an upcoming publication.

This is an existence statement: can we actually find such a DLR expansion? Least squares fitting on a sufficiently fine imaginary time grid would yield a DLR expansion satisfying this error bound. However, in practice, we also find that interpolation at the DLR nodes is sufficient to produce an expansion of comparable quality, satisfying the bound, or very nearly so.

We therefore recommend measuring error in either the \(L^2(\tau)\) or the \(l^2(i \nu_n)\) norm, as defined above. If one wishes to obtain a given accuracy in the \(L^\infty\) norm—that is, pointwise accuracy—in imaginary time or frequency, we would recommend decreasing \(\epsilon\) further.

Symmetrized DLR grids

By default, the DLR frequencies \(\omega_l\) are not chosen to be symmetrized about \(\omega = 0\), nor are the DLR imaginary time nodes \(\tau_k\) chosen to be symmetrized about \(\tau = \beta/2\) or the imaginary frequency nodes about \(i \nu_n = 0\). Indeed, this would represent an additional constraint in the pivoted Gram-Schmidt procedure used to select the points, and is not necessary in most applications. However, in some cases, it might be desirable to have symmetric frequencies and grids, and cppdlr provides this functionality via symmetrization flags. Please see the list of other cppdlr capabilities section on the examples page for a list of cppdlr tests which showcase this functionality.

We make a note about symmetrization for bosonic Green’s functions. In this case, we always include \(i \nu_n = 0\) as a DLR imaginary frequency node. Since we do not include the point \(\omega = 0\) in the \(r\) symmetrized DLR frequencies, \(r\) is even in this case, and the extra imaginary frequency grid points makes the number of DLR imaginary frequency nodes \(r+1\). To obtain the DLR coefficients from samples of the Green’s function at the DLR imaginary frequency nodes, we therefore solve a slightly overdetermined (\((r+1) \times r\), rather than \(r \times r\)) linear system using the least squares method.