# Minimax Isometry Method

We present a compressive sensing approach for the long standing problem of Matsubara summation in many-body perturbation theory. By constructing low-dimensional, almost isometric subspaces of the Hilbert space we obtain optimum imaginary time and frequency grids that allow for extreme data compression of fermionic and bosonic functions in a broad temperature regime. The method is applied to the random phase and self-consistent GW approximation of the grand potential and integration and transformation errors are investigated for Si and SrVO_3 .

There are no comments yet.

## Authors

• 1 publication
• 1 publication
04/26/2016

### Compressive phase-only filtering at extreme compression rates

We introduce an efficient method for the reconstruction of the correlati...
11/27/2015

### Gradient Estimation with Simultaneous Perturbation and Compressive Sensing

This paper aims at achieving a "good" estimator for the gradient of a fu...
04/28/2020

### Phase reconstruction with iterated Hilbert transforms

We present a study dealing with a novel phase reconstruction method base...
03/06/2017

### Generalizing CoSaMP to Signals from a Union of Low Dimensional Linear Subspaces

The idea that signals reside in a union of low dimensional subspaces sub...
10/14/2021

### Learning a Compressive Sensing Matrix with Structural Constraints via Maximum Mean Discrepancy Optimization

We introduce a learning-based algorithm to obtain a measurement matrix f...
08/26/2020

### Fast Bayesian Force Fields from Active Learning: Study of Inter-Dimensional Transformation of Stanene

We present a way to dramatically accelerate Gaussian process models for ...
05/21/2020

### Blind Two-Dimensional Super Resolution in Multiple Input Single Output Linear Systems

In this paper, we consider a multiple-input single-output (MISO) linear ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction: Matsubara Technique

The Matsubara technique is a way to formulate quantum field theory at finite-temperature. More precisely, it makes use of the Wick rotationWick (1954), which transforms the real time axis of Minkowski spacetime to the imaginary time axis . Because real space remains unchanged by this transformation, spacetime becomes essentially euclidean, so that this approach is also known as euclidean quantum field theory (QFT). A summary on euclidean QFT is given in the appendix. Here we focus on the Matsubara method.

As Matsubara has shown, the imaginary time-integrals in finite-temperature perturbation theory are restricted to the interval .Matsubara (1955) This has the advantage that one can expand the imaginary time-dependence of the corresponding integrands into a Fourier series, such that imaginary-time integration becomes essentially an (infinite) series over discrete Fourier coefficients. The corresponding discrete frequencies are known as Matsubara frequencies and it is important for us to distinguish between fermionic, denoted by in the following, and bosonic Matsubara frequencies, denoted by in the remainder of this paper.

Fermionic Matsubara frequencies represent the non-zero Fourier modes of fermionic functions, while bosonic frequencies are the non-zero modes of bosonic functions. This is explained in more detail below by means of the free-electron Green’s function (Feynman propagator) and the irreducible polarizability, the building blocks of finite-temperature perturbation theory.

The free Feynman propagator in imaginary time represents a prototype of a fermionic function. In a one-electron basis the free propagator is diagonal and the entries readNegele and Orland (1988); Fetter and Walecka (2003)

 gα(−iτ)=e−(ϵα−μ)τ[(1−fα)Θ(τ)−fαΘ(−τ)], (1)

where , and are the one electron energy, the chemical potential and the Fermi function (see Equ. 69), respectively. Here is the Heaviside step function,Olver et al. and it is the reason why changes sign at . Also, the presence of the step functions implies

 gα(−iβ+iτ)=−gα(iτ),0<τ<β. (2)

This anti-symmetric property has an important effect on the Fourier series in the interval

 gα(−iτ) =1β∞∑m=−∞~gα(iωm)e−iωmτ (3) ~gα(iωm) =∫β2−β2dτgα(−iτ)eiωmτ, (4)

because it contains only fermionic frequencies

 ωm=2m+1βπ,m∈Z. (5)

This holds true for all fermionic functions on the imaginary time axis, including the self-energy (see appendix B).

An example for a bosonic function is the independent particle polarizability, which is diagonal and has the entries

 χαγ(−iτ)=−gα(−iτ)gγ(+iτ). (6)

However, in contrast to Equ. (2), bosonic functions do not change sign in imaginary time, but are symmetric

 χαγ(−iβ+iτ)=χαγ(+iτ),0<τ<β. (7)

Consequently, the Fourier expansion

 χαγ(−iτ) =1β∞∑n=−∞~χαγ(iνn)e−iνnτ (8) ~χαγ(−iνn) =∫β2−β2dτχαγ(−iτ)eiνnτ (9)

contains only bosonic frequencies

 νn=2nβπ,n∈Z. (10)

It is often argued that bosonic (fermionic) functions are periodic (anti-periodic) in . This is strictly speaking not correct. The free propagator, for instance, is defined a priori only in the fundamental imaginary time interval , because (1) grows or decays exponentially for arguments outside (green lines in Fig. 1). The same holds true for the irreducible polarizability. In fact, only the Fourier expansions (8) and (3) define periodic and anti-periodic functions in with (anti-) period . This behavior is illustrated in Fig. 1 showing a typical fermionic and bosonic function. In practice, it is quite important to recall that exponentially growing terms in propagators are not present, because the -integrations are performed over or, equivalently, over . Due to consistency with our previous papersKaltak et al. (2014a, b); Liu et al. (2016) we work in the interval. Furthermore, we mean by correlation function either a bosonic or fermionic function.

The Matsubara summation has one major drawback. Matsubara series converge very slowly with the cutoff frequency, as shown in II.4. One thus requires compressed Matsubara representations of bosonic and fermionic correlation functions. For this purpose we studied isometries in metric spaces and developed the following formalism.

## Ii Minimax Isometry

To keep the notation simple, we consider the case eV in the following. The general case follows from scaling relations that are discussed in II.2 and II.3.

Deriving an alternative to the Matsubara technique requires one additional abstraction level. For this purpose we write for the energies and assume that ( is allowed in the following). Every energy

is represented by a vector

in a Hilbert space and it is assumed that and are two complete basis sets in imaginary time and frequency, such that the identity operator can be expressed as

 11= (11) 11= ∑n∈Z|n⟩⟨n| (12)

for all with representing the set of integers. From a functional analysis perspective, one says that the two spaces and are isometric w.r.t. to the scalar product induced norm . This isometry is effectively a simple basis transformation that does not change the induced norm, since

 (13)

If is the time and the frequency basis, then and are the matrix elements of the forward and backward basis transformation, respectively as for instance given in Eqs. (3),(4),(8) and (9).111 The isometry (13) is known as Parseval theoremvon Querenburg (2013) or, if

is continuous, Plancherel theorem for Fourier transforms.

Plancherel (1910). Consequently, the two spaces and are equivalent and span the same Hilbert space . This equivalence holds true only if infinitely many basis vectors are considered; for finite dimensional subspaces the perfect isometry (13) is violated.

One may illustrate the violation of the isometry (13) with the discrete Fourier transform (DFT) having the bases with and (truncated bosonic Matsubara grid). The corresponding completeness relations (11), (12) become projectors onto finite dimensional subspaces and have the form

 P= 12NN∑k=1|τk⟩⟨τk| (14) ~P= N∑n=−N|n⟩⟨n|. (15)

Only in the limit the projectors approach the identity operator . For finite , the isometry (13) is violated, but can be replaced by a so-called -isometryFleming and Jamison (2002); Ding (1988)

 ∥P−~P∥∞:=maxa≤x≤b|⟨x|P−~P|x⟩|≤ε. (16)

In the following, we will drop the subscript and write for the maximum norm .

Of interest to us is the magnitude of and especially how it decreases with increasing . For instance, in the case of the DFT is the Riemann sum of the integral in (13) of order and is known to be a poor method to evaluate integrals. As a consequence, of the Matsubara grid is a weakly decaying function in and cannot be used for our purposes, as shown in section II.5.

The following question naturally arises: how can one determine -isometric subspaces , and , such that the completeness relations (11) and (12) are approximated as good as possible for all vectors with ?

Using the notation in (16), the answer to this question are the solutions of following minimax problems:

 minσk>0,τk∈[0,1/2] ∥∥ ∥∥11−N∑k=1σk|τk⟩⟨τk|∥∥ ∥∥ (17) minλk>0,ωk>0 ∥∥ ∥∥11−N∑k=1λk|ωk⟩⟨ωk|∥∥ ∥∥. (18)

Provided the solutions exist, they are known to yield errors that decay exponentially with .Braess (1986) In the following we prove that (17) and (18) satisfy our requirements.

To prove the assertion above it suffices to show that the minimax errors are an upper bound for the isometry violation in (16). Therefore, assume and are the solutions of (17) and (18) and and the corresponding projectors, respectively. Then a positive number exists (for every given ) as an upper bound for (17) and (18) and one can write

 ∥∥11−P∗∥∥≤12ε∥∥11−~P∗∥∥=∥∥~P∗−11∥∥≤12ε. (19)

Adding both inequalities in (19) and using the triangle inequality (satisfied by every normvon Querenburg (2013)) one obtains

 (20)

Last inequality implies (16) for the projectors and concludes our proof. ∎

This is a quite remarkable result, because it means that the projectors and converge to the identity operator and, thus, define -isometric topological vector spaces that have the approximation property.Axler et al. (1999) A summary of -isometric bases is given in Tab. 1 and its selection is motivated below.

Note that the discussion above does not give a prescription how to determine the transformation ; a corresponding method is presented in section II.5.

The proof above contains only an upper bound for the transformation error in (20). This upper bound is inherited from the sum of the convergence rate of the minimax solutions in the - and -domain. This convergence rate has been studied by Braess and Hackbusch for the minimax problem IC in the -domain listed in Tab. 1. They obtained for , where belongs to the largest possible error for a given order .Braess and Hackbusch (2005) Our numerical experiments discussed in section IV indicate similar convergence rates for all other minimax problems in Tab. 1. In contrast, the DFT or Matsubara grid has only a linear rate of convergence .

### ii.1 Motivation of basis functions

In this subsection, we motivate our choice of -isometric basis functions tabulated in Tab. 1. We start with the IC bases, which has been used in previous publications by the authors to construct optimized minimax grids for low scaling random phase and algorithms at zero temperature.Kaltak et al. (2014a, b); Liu et al. (2016) Specifically, of IC describes the imaginary time dependence of the independent particle polarizability at zero temperature, while the corresponding -basis functions describe its imaginary frequency dependence.Kaltak et al. (2014a) Then isometric bases are the ideal choice for the zero-temperature formalism of quantum field theory, because the corresponding conserved -norm (forth column of Tab. 1) is the key quantity for the second order contributions to the correlation energy (see appendix (76) or Ref. Kaltak et al., 2014a; Takatsuka et al., 2008). These contributions involve energy denominators of the form and are considered to be bound, i.e. virtual states with energies are separated by a band gap from occupied states with energy . The induced norm is important and essentially tells the optimization in the minimax problem, which contributions to the energy are most relevant. In our case contributions from small energy differences dominate over contributions from large energy differences.

The -basis function of IC has the same time dependence (apart from the opposite sign), while the imaginary frequency dependence of the cosine transformation differs considerably from the one obtained from the sine transform (compare of IC and IC). It comes with no surprise that the corresponding minimax grids for and differ too.Liu et al. (2016) However, the minimax isometry guarantees that one can map in time between IC and IC with high precision; a fact that has been exploited in a low scaling self-consistent approach at zero temperature.Grumet et al. (2018)

Next, we consider the four basis functions of group IA and IB. They can be grouped into bosonic (IA and IB) and fermionic (IA and IB) pairs, where the -bases for bosonic (fermionic) functions are defined for bosonic (fermionic) frequencies a priori (see next section for additional motivation of these functions). The symbol denotes the Fourier transformations (cosine or sine) of these functions. When optimizing the frequency grids using the minimax algorithm, we allow

in IA and IB to deviate from the corresponding Matsubara grid. Indeed, the corresponding Minimax solutions are non-uniformly distributed, but nevertheless closely match Matsubara frequencies at small

. It turns out, as shown in section IV, that this freedom allows us to describe the high frequency behavior of the correlation functions with high precision even in low dimensional subspaces . The corresponding -isometric time basis functions have the fermionic anti-symmetry [Equ. (2)] and bosonic symmetry [Equ. (7)] for , respectively, and are illustrated in Fig. 2.

Similar to the zero-temperature case, the conserved -norm describes the second order contribution to the correlation part of the grand-canonical potential (see appendix 76). Thus, it can be employed for data compression when calculating the correlation part of the grand canonical potential in the random phase approximation (see appendix B). The corresponding imaginary time and frequency grid are discussed in II.2 and II.3, respectively.

There is one further important point to note here. Some bosonic functions, such as the polarizability or screened potential can be entirely presented by IA basis functions since (green lines in Fig. 2 right). This means that the second order part of the grand-canonical potential or the RPA correlation energy can rely on grids constructed for IA, only. This is not necessarily so for bosonic functions constructed e.g. by the product of the self-energy and Green’s functions, as is the case in the Galitskii-Migdal (GM) formula.

In fact, the minimax isometry method is not straightforwardly applicable to the GM formula of the grand canonical potential; here both, even and odd basis functions in the frequency domain (IB

and IA) contribute to the grand potential and have different -norms (compare third column of IA and IB). We, therefore, propose an alternative approach in this case that is based on the minimization of the -quadrature error instead (see section II.4).

### ii.2 Imaginary time grid

To construct an imaginary time grid, we make use of the scaling properties

 τj→βτj,σj→βσj (21)

that allow to recover the time quadrature for an arbitrary interval from the unscaled solution determined for . Furthermore, we allow for infinitesimal small excitation energies and accordingly set . The minimization interval is then characterized by one parameter , as for the IC quadrature.Kaltak et al. (2014a)

We now discuss how to choose the imaginary time grid. This remains a somewhat involved issue even with the yet gained insight. The main problem is that we want to use one and only one time grid, since this will allow us to represent the Green’s functions and the bosonic quantities on the same time grid permitting one to calculate for instance the polarizability as on that single grid. The even and odd basis functions of the IA and IB -isometry in Tab. 1 are obviously identical for fermions and bosons for positive ,

 ϕ(x,τ)= 12coshx2(1−2τ)coshx2,τ>0 (22) ψ(y,τ)= 12sinhy2(1−2τ)coshy2,τ>0. (23)

Even without considering the related norm in Tab. 1, these functions seem to be a handy and intuitive choice, since they can be used to express the diagonal elements of the free propagator of (1)

 (24)

as well as the matrix elements of the polarizability (6). Optimization of the time grid for even and odd functions, however, yields different time grids. As already motivated in the previous section, we have decided to use the optimal even time grid (IA) as a common grid for both fermionic and bosonic functions and repeat the arguments here. (i) The second order and RPA correlation energy depends only on bosonic functions, e.g. the polarizability. Hence the fermionic functions are only used at an intermediate stage. (ii) The polarizability possesses a special symmetry matching the IA -isometry in Tab. 1. (iii) Finally, the imaginary time grid for the even functions yields a small minimax error also for the odd basis functions for the entire interval , with negligible deviations even for . This follows from the fact that the corresponding norms differ only for very small arguments (compare third column of IA and IB in Tab. 1). Consequently, we solve the minimax problem (17) only for and use the same time grid points for the odd basis functions.

To obtain the time grid points , it is convenient to rewrite the minimax problem (17) into the following form

 minσj>0,τj∈(0,1/2)max0≤x≤b∣∣ ∣∣∥x∥22−N∑j=1σjϕ2j(x)∣∣ ∣∣ (25)

with and the maximum one-electron energy considered. Then it becomes evident that (25) is a non-linear fitting problem of separable type,Golub and Pereyra (2003) which in general has only a solution, if every basis function is linearly independent and has less than zeros. The alternant theorem then impliesBraess (1986); Hammerlin and Hoffmann (1994) a set of points (alternant) and a set of non-linear equations

 ∥x∗j∥22−N∑k=1σ∗k|⟨τ∗k∣∣x∗j⟩|2=EN(−1)j (26)

with

 EN=±maxa≤x≤b∣∣ ∣∣|⟨ϕ|x⟩|2−N∑k=1σkϕk(x)∣∣ ∣∣ (27)

being positive (negative) if the l.h.s. of (26) is positive (negative) at .

The alternant theorem provides the basis for the non-linear Remez algorithm that has been used successfully in other papers and yields the minimax solution .Braess and Hackbusch (2005); Takatsuka et al. (2008); Kaltak et al. (2014a) The minimax solution yields abscissas in the unscaled interval . Furthermore, the corresponding weights are positive and holds true. This is important for the application in many-body theory, since the conservation of particles is guaranteed with increasing including particles with energy .

Furthermore, the same minimax quadrature can be used to evaluate ”off-diagonal terms”, since errors of the form

 η(x,y)=⟨x|y⟩−N∑k=1σ∗jϕ∗k(x)ϕ∗k(y) (28)

are exponentially suppressed with increasing quadrature order . An example is given in Fig. 3 that demonstrates how the error is uniformly distributed in the entire region and bounded for .

Before we discuss the construction of the frequency grids, we note that a similar basis has been recently used to compress Green’s functions on the imaginary time axis in quantum Monte Carlo algorithms.Shinaoka et al. (2017) We believe that there is a close connection to our method, since we found that the quadrature obtained from the solution of (25) is also a good approximation to the solution for the corresponding problem for the odd basis function (23) and both, and linear combined yield Shinoaka’s basis. However, Shinoaka et al. determine the grid as the solution of an integral equation and the connection to -isometric subspaces is not immediately evident. We, therefore, prefer to solve the non-linear optimization problem (25) instead.

### ii.3 Bosonic Frequency Grid

To construct the bosonic Matsubara grid we use the -isometric basis of the even time basis (22), specifically again the IA basis

 ~ϕn(x)=xx2+ν2ntanhx2. (29)

The motivation behind this choice is three-fold. Firstly, it is obtained from the cosine transformation of the even time basis (22) evaluated for bosonic Matsubara frequencies . Thus it describes the imaginary frequency dependence of the polarizability (6) that is of bosonic nature; the IA basis is obtained from the sine transform of these functions and is evaluated at fermionic frequencies and hence irrelevant for the evaluation of bosonic integrals. Secondly, we can use the minimax isometry method to switch between the frequency and time representation of the polarizability with high precision. This follows from the theorem proved in section II Equ. (16). Lastly, the infinite bosonic Matsubara series of the RPA grand potential (79

) can be evaluated with high precision without using any interpolation technique.

In practice, the unscaled bosonic frequency quadrature for is determined first and following scaling relations are used to obtain the result for arbitrary inverse temperatures

 νk→νkβ,λk→λkβ. (30)

The corresponding minimax problem reads

 minλk>0,νk∈(0,∞)max0≤x≤b∣∣ ∣∣∥x∥22−N∑k=1λk~ϕ2k(x)∣∣ ∣∣ (31)

where the -norm is given in Tab. 1. The solution is called IA-quadrature in the following, in agreement with the notation used in Tab. 1.

Before we investigate the convergence of the IA-quadrature, we discuss the construction of the fermionic grid.

### ii.4 Fermionic Frequency Grid

In this section, we discuss the construction of a compressed fermionic Matsubara quadrature. More precisely, we look for a quadrature that is converging exponentially with the number of grid points for the GM expression for the grand potential (80). In contrast to the polarization function and second order correlation energies, this requires an accurate handling of fermionic functions of the type IA and IB in the frequency domain, which is an intricate problem.

First, we consider the frequency dependence of the free propagator (1). The cosine and sine transformations of the odd and even time basis functions (23) and (22) for fermionic frequencies (5) are determined as:

 IB1:(???)→um(x)= ⟨um|x⟩=xx2+ω2m (32) IA2:(???)→vm(x)= ⟨vm|x⟩=ωmx2+ω2m. (33)

Then the non-interacting propagator (1) on the fermionic Matsubara axis reads

 ~gα(iωm)=um(ϵα−μ)+ivm(ϵα−μ). (34)

Second, we observe that every fermionic function can be decomposed into terms that are even and odd in , including the product of the propagator and self-energy as it appears in the GM grand potential (80). It is, obvious, that only the real part of the product contributes to the total energy. Thus the most general matrix element, which gives a non-zero contribution to the GM grand potential has the form222The imaginary part is proportional to odd terms in the frequency and vanishes when the sum over all (positive and negative) fermionic Matsubara frequencies is carried out.

 (35)

where and are the poles of the Green’s function and the self-energy on the real-frequency axis and the spectral densities, respectively. Without loss of generality, we set and assume that the magnitudes of the poles are smaller than a positive number, i.e. .

Third, we note that the analogue of the IA-quadrature of bosonic functions (31) for fermionic ones

 minσk>0,ωk∈(0,∞)max0≤x≤b∣∣ ∣∣∥x∥22−N∑k=1σku2k(x)∣∣ ∣∣ (36)

yields the IB-quadrature, see Tab.1. Unfortunately, the IB-quadrature only allows to evaluate the first term on the r.h.s. of (35) accurately, but fails for the product of two odd functions . Similarly, the IA-quadrature obtained from the minimax problem

 minσk>0,ωk∈(0,∞)max0≤x≤b∣∣ ∣∣∥x∥22−N∑k=1σkv2k(x)∣∣ ∣∣ (37)

that approximates the same norm as the time and IA-quadrature, describes only the second term in (35). Consequently, neither the IB- nor the IA-quadrature can be used for our purposes.

A solution to this dilemma can be found by approximating the -norm instead

 (38)

A proof of this identity is found in the appendix C. That is, we determine the solution of

 minγk,ωk>0max0≤x≤b∣∣ ∣∣∥x∥1−N∑k=1γkuk(x)∣∣ ∣∣. (39)

Because holds true for any ,Fleming and Jamison (2002) the solution , called F-quadrature in the following, yields linearly independent basis functions that span a larger vector space than the basis obtained from (36). Our numerical experiments discussed below show that the F-quadrature evaluates the infinite sum over both terms in (35) with high precision for increasing .

A similar technique was employed by OzakiOzaki (2007), who used the same basis functions , but determined the abscissas from the partial fraction decomposition of the hyperbolic tangent in combination with a continued fraction representation of the hypergeometric function to derive a compressed form of (38).

Both, the F-quadrature as well as Ozaki’s hypergeometric quadrature (OHQ) use essentially a rational polynomial approximation to the hyperbolic tangent. In the following, we show why this approach provides also a good approximation of fermionic Matsubara series, such as the the density matrix for holes (upper sign) and electrons (lower sign)

 Γ=±limη→0±G(−iη)=±limη→0±1β∞∑n=−∞~G(iωn)e−iωnη (40)

where and is the interacting Green’s function in imaginary time and on the Matsubara axis, respectively. Specifically, we show that the last expression on the r.h.s. of (40) can be approximated with the following quadrature formula

 Γ≈sgn(η)211+N∑k=1σk2[~G(iωk)+~G(−iωk)], (41)

where

is the identity matrix in the considered basis and

are either the OHQ- or F-quadrature points.

To understand Equ. (41) and, in general, why an approximation to the hyperbolic tangent provides an excellent approach to compress fermionic Matsubara series, we consider a general correlation function that is analytic in the complex plane with a branch cut on the real axis and decays with or faster to zero for ; for instance or . Following Fetter and Walecka,Fetter and Walecka (2003) one introduces an auxiliary function with an infinitesimal to force the complex contour integral over the infinite large outer circle in Fig. 4 to vanish, that is

 ∮Cdz2πi~A(z)hη(z)=0. (42)

Regardless of the specific choice of (discussed below), one can easily show using the residue theorem and the contours depicted in Fig. 4 following identities:

 ∑n∈ZResz=iωn[~A(z)hη(z)]=∮Fdz2πi~A(z)hη(z)=−∮Bdz2πi~A(z)hη(z)=∫∞−∞dω1πIm[~A(ω)hη(ω)] (43)

Apart from condition (42), the auxiliary function has to be chosen such that the l.h.s. in (43) gives the fermionic Matsubara series , which imposes two conditions on . Firstly, must have an infinite number of poles located at (crosses in Fig. 4). Secondly, the corresponding residue has to be for .

If is of order for the contour integral (42) is zero also for constant functions in and one can choose

 hη(z)=12tanhz2=12∑n∈Z[1z−iωn+1z+iωn], (44)

where the last line follows from (38) and shows the locations of the poles of . Thus for the fermionic Matsubara series is independent of the directional limit . The replacement of the hyperbolic tangent by a rational polynomial with poles on the imaginary axis in (43) shows that both terms in the GM energy (35) can be well approximated using the F-quadrature.

In contrast, for correlation functions that decay only with at the sign of the infinitesimal matters. This includes as well as any mean field terms, e. g. Green’s function times the mean field Hamiltonian . For instance, the limit in (40) gives the electron density matrix, while gives the density matrix of holes. For functions of order at one, therefore, has to add a term to the hyperbolic tangent. As can be shown easily, the form for for which (42) holds true is333 Note, this expression is obtained from the free propagator (1) by replacing , and using the identity .

 hη(z)=[sgn(η)2+12tanhz2]e−zη. (45)

Inserting Equ. (45) into the r.h.s. of (43) yields

 ∑n∈Z~A(iωn)e−iωnη=sgn(η)2∫∞−∞dω1πIm[~A(z)e−zη]+∫∞−∞dω1πIm[~A(z)12tanhz2e−zη], (46)

In the last term on the r.h.s. the evaluation of the limit can be performed before integration, because the integrand is of order for [see Equ. (44)]. The corresponding integral over the arch vanishes, so that the last term in (46) on the r.h.s. can be rewritten into the Matsubara series of that is independent of the sign of . This is the convergent part of the Matsubara series and the term that can be evaluated using quadratures. Note that the first term on the r.h.s. of (46) cannot be written into a Matsubara series, because the integrand diverges for prohibiting the closure of the integration contour at infinity. However, for one has

 limη→0±sgn(η)2∫∞−∞dω1πIm[~G(z)e−zη]=±1211. (47)

This concludes our proof of Equ. (41). We call this term, therefore, the divergent part of the Matsubara series, although the term is finite in any practical calculation (number of electrons/holes is finite in practice). We use (41) for the evaluation of the density matrix in self-consistent calculations at finite temperature (see section IV).

In summary, one can say that the approximation of the hyperbolic tangent by rational polynomials with poles only on the imaginary axis gives rise to fermionic Matsubara quadratures that describe the convergent part of the Matsubara representation. To obtain the directional limits of weakly decaying correlation functions, such as the propagator of electrons or holes, the integral over the spectral function has to be added or subtracted, respectively. The evaluation of the GM energy does not require this term, because decays with . Analogous bosonic quadratures can be obtained by approximating the -norm of the hyperbolic cotangent, but this was not further investigated.

We have compared our F-quadrature with Ozaki’s hypergeometric quadrature (OHQ) by means of calculating the GM factor (35) for a model that includes 50 randomly sampled poles in and 50 poles in . The results for and eV are shown in Fig. 5 and are contrasted to the grid convergence for the ordinary fermionic Matsubara quadrature . It can be seen that the F-quadrature outperforms the OHQ in all cases, especially for (low temperatures). This can be explained by the fact that the F-grid minimizes the quadrature error for all energies uniformly in the interval . The corresponding OHQ-quadrature error is non-uniformly distributed in the same interval and has the effect that at high values the convergence is almost linear with the number of grid points for small . The same figure, also shows the linear convergence of the conventional Matsubara grid and demonstrates its pathology in practice.

### ii.5 Minimax Isometry Transformation

We have seen how different basis functions for the time and frequency domain give rise to different grids. In this section we study the error made by transforming an object represented on the time grid to the frequency axis. As a measure for the transformation error we use

 ~E(ω)=mintωk∈R∥∥ ∥∥⟨ω|x⟩−N∑k=1tωk⟨τ∗k∣∣x⟩∥∥ ∥∥22, (48)

where is either the even or odd time basis function (22),(23) evaluated at the minimax time grid obtained from (25) and acts as a placeholder for one of the basis functions in the frequency domain listed in Tab. 1 with being a positive real number. The -norm is evaluated by sampling the values with 100 points determined from the alternant of the minimax time problem in (25) and additional uniformly distributed points in each of the sub-intervals .

The solution of the ordinary least square problem for the frequency

is then given by the corresponding normal equationPress et al. (2007)

 100∑i=1⟨ω∣∣X∗i⟩⟨X∗i∣∣τ∗k⟩=N∑j=1tωj100∑i=0⟨τ∗j∣∣X∗i⟩⟨X∗i∣∣τ∗k⟩,k=1,⋯,N. (49)

The transformation error is rather insensitive to changes of the number of sampling points ; the alternant points of the time grid also often suffice in practice.

Furthermore, Equ. (48) allows to plot as a function of the frequency. This gives independent insight, on which frequencies one is supposed to use in combination with a certain set of time basis functions, independent of the previous considerations (see Fig. 6).

Transformation to the IA frequency basis functions (blue line), clearly shows that the error is minimal at the previously determined IA frequency points (blue triangles), and transformation to the IA frequency basis functions (green line) shows that the error is smallest at the previously determined IA frequency points (green diamonds). The reason for this behavior is due to the fact that the IA, IA and the time quadrature for the even time basis [Equ. (22)] possess the same approximation property and span -dimensional, almost isometric subspaces of the Hilbert space as proven in section II. The good agreement is a numerical confirmation that the previously determined frequency grids are optimal.

Thus, for polarizabilities and second order correlation energies, the optimal frequency points are clearly the IA-quadrature points (triangles). They approach the conventional bosonic frequencies [Equ. (10)] at small . The corresponding weights (not shown) approach (except for the first frequency ). Higher quadrature frequency points (as well as weights) deviate considerably from the conventional bosonic Matsubara points . This behavior is very similar to the bosonic grid presented by Hu et al. that is based on the continued fraction decomposition of the hyperbolic cotangent, the analogue of Ozaki’s method for bosons.Hu et al. (2010) However, the IA-quadrature has the advantage that the error is minimized uniformly for all transition energies , while the continued fraction method yields non-uniformly distributed errors in general.

Transformation from the even time basis [Equ. (22)] to the fermionic frequency basis IA yields further important insight. As already emphasized the minimax IA frequency points match exactly those frequency points where the error for transformation into the IA basis functions is minimal. On the other hand, the IB and F minimax grid points are chosen to optimally present odd time basis functions [Equ. (23)] using the corresponding frequency basis [Equ. (32)]. At small frequencies, these points are slightly shifted away from the optimal IA frequency points, causing somewhat larger transformation errors for even functions to IA basis functions. This is to be expected, since the points have been chosen to approximate a different scalar product than for IA (and IA), see fourth column in Tab. (1). Specifically, the IB minimax frequencies are by construction optimal to present odd time basis functions. Although, IB and IA minimax points are close at low frequencies, they progressively move away at higher frequencies, which prohibits the construction of a common frequency grid that can present the GM energy well (Sec. II.4). From this figure, it is somewhat unclear, why the F frequency grid works well, although it is noteworthy that the corresponding frequency points lie roughly at the positions where IA and IA errors intersect. This might implying an equally acceptable representation of odd and even functions.

#### ii.5.1 ε-isometric time grids of the F-quadrature

To avoid ambiguity, we call the time grid defined in section II.2 the IA time grid in the following, while the IB time grid is the same minimax solution (25), but with norm and basis function of the IB isometry listed in Tab. 1.

Recapitulating the previous section, a natural question arises: Is there an optimum time grid for the F-quadrature? In analogy, to Equ. (48) this grid may be defined by the minima of the inverse transformation error

 (50)

where are the abscissa of the F-quadrature. Table 1 shows that there are two possible choices for the time and frequency basis functions ; one function describing the transformation error from odd to even basis functions (Equ. (33) to Equ. (22)) and one from Equ. (32) to Equ. (23). Both functions are plotted in Fig. 7 (blue and green line), respectively.

The figure clearly shows that the minima of both error functions differ and implies two -isometric time grids for the frequency F-grid. This is analogous to the forward transformation errors discussed in the previous section, where the A- and A-frequency grids are -isometric to the IA time grid, respectively. For the F-grid, however, the IA and IB minimax solutions in time coincide with the minima of the error functions only for , for larger values of the transformation error minima (-isometric grids) deviate from the IA and IB minimax grid points, respectively (compare minima of green and blue curve with triangles and diamonds in Fig. 7).

The small IB transformation error for follows from the fact that for small the time basis becomes . Per construction (see (38)), the hyperbolic tangent function is approximated well by the basis (32) using the F-grid. The IA transformation error (blue line), in contrast, is several orders of magnitude larger at , since the time basis function is constant and the frequency basis functions represent constants only poorly. The deviation of the IA and IB time grids from the -isometric time grids at higher values is not surprising, since the F-grid deviates from the -isometric frequency grids, that is the A and B-grid discussed in section II.3 and II.4.

In summary, we recommend using the IA time grid presented in section II.2 in combination with the A frequency grid for bosonic functions (see II.3) and the F-grid for fermionic functions. First, the exact -isometric time points of the F-grid is only known numerically from inspection of the transformation error; an analogue to the minimax isometry method is not known to us. Second, Green’s function in imaginary time can be contracted without error, whilst at the same time transformation errors to the imaginary frequency axis are controlled. Before we demonstrate these advantages in section IV, we discuss the following details about our implementation in the Vienna ab initio software package (VASP).Kresse and Joubert (1999)

## Iii Technical Details

The implementation of the finite temperature RPA and algorithms is the same as the zero-temperature ones,Kaltak et al. (2014b); Liu et al. (2016) with three exceptions.

• The zero-temperature frequency grid is replaced by the IA-grid discussed in II.3 for the bosonic correlation functions , while the F-quadrature from II.4 replaces the grid for the fermionic functions and .

• All correlation functions are evaluated on the same imaginary time grid presented in sec. II.2.

• The occupied and unoccupied Green’s function need to be set up carefully considering the partial occupancies in each system.

The last point requires some clarification. The Green’s function for positive and negative times can be combined to a full Green’s function using Heaviside theta functions

 G(−iτ)=Θ(τ)¯¯¯¯G(τ)−Θ(−τ)G––(τ). (51)

At zero temperature the occupied and unoccupied imaginary time Green’s function readRojas et al. (1995); Kaltak et al. (2014b)

 G––(τ)|β=∞= ∑αΘ(ϵF−ϵα)e−(ϵα−ϵF)τ (52) ¯¯¯¯G(τ)|β=∞= ∑αΘ(ϵα−ϵF)e−(ϵα−ϵF)τ. (53)

Here the Fermi energy makes sure that and contains only unoccupied (occupied) one-electron states. This changes as temperature increases, which follows from and the limit

 Θ(±ϵF∓ϵα)=limβ→∞1e±β(ϵα−μ)+1, (54)

implying the form given in (1) for the full Green’s function (51). Consequently, the Green’s function needs to include also partially occupied states at finite temperature and vise versa, so that the positive and negative imaginary time Green’s functions

 G––(τ)= ∑αfαe−(ϵα−μ)τ,τ<0 (55) ¯¯¯¯G(τ)= ∑α(1−f