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\). 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.