Analytic continuation with Pad\'e decomposition

TL;DR

This paper introduces a Padé decomposition method to improve the accuracy of analytic continuation of Green's functions by enhancing the convergence of Matsubara frequency summations, enabling more precise input data for the Padé approximation.

Contribution

The paper proposes a novel Padé decomposition approach that significantly accelerates Matsubara frequency summation convergence, improving the analytic continuation process.

Findings

Enhanced convergence of Matsubara sums using Padé decomposition

More accurate high-precision Green's function data achievable

Improved reliability of Padé analytic continuation results

Abstract

The ill-posed analytic continuation problem for Green's functions or self-energies can be done using the Pad\'e rational polynomial approximation. However, to extract accurate results from this approximation, high precision input data of the Matsubara Green's function are needed. The calculation of the Matsubara Green's function generally involves a Matsubara frequency summation which cannot be evaluated analytically. Numerical summation is requisite but it converges slowly with the increase of the Matsubara frequency. Here we show that this slow convergence problem can be significantly improved by utilizing the Pad\'e decomposition approach to replace the Matsubara frequency summation by a Pad\'e frequency summation, and high precision input data can be obtained to successfully perform the Pad\'e analytic continuation.

Tables

Table 1. Table 1: Values of the first 20 pairs of ξ j subscript 𝜉 𝑗 \xi_{j} and η j subscript 𝜂 𝑗 \eta_{j} defined in the Padé spectral decomposition of the Fermi function.

$j$	${\eta }_j$	${\xi }_j/\pi$	$j$	${\eta }_j$	${\xi }_j/\pi$
1	1.00	1.00	11	1.00	21.0
2	1.00	3.00	12	1.03	23.0
3	1.00	5.00	13	1.22	25.2
4	1.00	7.00	14	1.70	28.1
5	1.00	9.00	15	2.52	32.2
6	1.00	11.00	16	3.90	38.5
7	1.00	13.00	17	6.59	48.7
8	1.00	15.00	18	13.1	67.3
9	1.00	17.00	19	36.8	111
10	1.00	19.00	20	332	332

Equations

$$K^R(\omega )=K(i{\omega }_n\to \omega +i{0}^{+}),$$

$$A(\omega )=-\frac{1}{\pi }\mathrm{Im}{K}^R\left(\omega +i{0}^{+}\right).$$

$$K\left(z\right)=\frac{p_0+p_{1}z+\cdots +p_{r-1}{z}^{r-1}}{q_0+q_{1}z+\cdots +q_{r-1}{z}^{r-1}+z^r},$$

$$K(q,i{\omega }_n)=\frac{1}{\beta N}\sum _{k,i{\omega }_m}G(k,i{\omega }_m)G(q+k,i{\omega }_n+i{\omega }_m),$$

$$K\left(q,i{\omega }_n\right)=\frac{1}{N}\sum _k\frac{f\left({\xi }_k\right)-f\left({\xi }_{k+q}\right)}{i{\omega }_n+{\xi }_k-{\xi }_{k+q}},$$

$$f\left(x\right)=\frac{1}{e^x+1}=\frac{1}{2}-\sum _{j=1}^{N\to \infty }\frac{2{\eta }_{j}x}{x^2+{\xi }_j^2},$$

$$K\left(q,i{\omega }_n\right)$$

Full text

The theory of the Matsubara Green function provides a convenient and powerful tool to calculate physical properties at finite temperatures Mahan-2000 . In this theory, a Matsubara correlation function corresponding to a dynamic response function measured by experiments is evaluated at a discrete set of imaginary Matsubara frequencies. An analytic continuation is then performed from this Matsubara correlation function to extract the dynamic response function in real frequency. This analytic continuation can be carried out straightforwardly if the Matsubara correlation function can be expressed in a simple analytic formula. However, in many cases, the Matsubara correlation functions can be evaluated only approximately. In these cases, the analytic continuation becomes a numerically ill-posed inverse problem, and remains one of the most challenging issues in computational physics.

Different methods have been proposed to perform analytic continuation according to the accuracy of the input Matsubara correlation functions. The data of the Matsubara correlation function generated by quantum Monte Carlo simulations (QMC)Blankenbecler-1981 ; Hirsch-1983 contain stochastic noise induced by random sampling. For this kind of data, various numerical algorithms, such as the maximum entropy Jarrell-1996 ; Dominic-2016 ; Goulko-2017 , the singular value decomposition Creffield-1995 ; Gunnarsson-2010 and the stochastic analytical continuation Beach-2004 ; Sandvik-2016 ; Qin-2016 , have been developed and successfully applied to various quantum systems. However, these methods have difficulties in reproducing fine structures in dynamic response functions. If, on the other hand, high precision input data are available, another method called the Padé analytic continuation, Vidberg-1977 ; Beach-2000 ; Ostlin-2012 ; Osolin-2013 ; Schott-2016 can be used. In this method, the analytic continuation is carried out by interpolating the input date defined in the Matsubara frequency domain with a rational function that is assumed to represent approximately the Matsubara correlation function. However, the Padé analytic continuation is sensitive to the accuracy of the input data. A tiny noise in the input may lead to a large error in the final result, and high precision calculation of the Matsubara correlation function is needed in order to use the Padé formula Beach-2000 .

In the Feynman diagram calculation of quantum field theory, a Matsubara correlation function is often determined by a summation over a set of Matsubara frequencies. In many cases, this summation cannot be carried out analytically and one has to resort to numerical calculations. A problem often encountered is that the summation converges very slowly with the Matsubara frequency, rendering high precision data required in the Padé analytic continuation difficult to obtain. It is well known that the summation over Matsubara frequencies is equivalent to the contour integration around the poles of Fermi or Bose distribution functions. The poor convergence of the summation can be overcome by finding new poles of the Fermi or Bose distribution functionsTaisuke-2007 ; Croy-2009 ; Hu-2010 ; Jie-2011 . An efficient approach along this line is the Padé decomposition first introduced by OzakiTaisuke-2007 , and the poles are called Padé frequencies. In this letter, we apply this approach to solve the slow converging problem encountered in the analytic continuation using the Padé rational polynomial approximation.

Let us assume $K(i\omega_{n})$ to be the Matsubara Green function corresponding to a dynamic correlation function $K^{R}(\omega)$, where $\omega_{n}$ is the Matsubara frequency, and $\omega$ is the real frequency. As $K(z)$ is analytic in the whole upper complex plan of $z$ excluding the real axis, the retarded Green function $K^{R}(\omega)$ can be obtained from $K(i\omega_{n})$ by analytic continuation

[TABLE]

provided that the analytic expression of $K\left(i\omega_{n}\right)$ is known. The spectral function $A\left(\omega\right)$ is determined by the imaginary part of $K^{R}(\omega+i0^{+})$ by the formula

[TABLE]

However, as already mentioned, the analytic expression of $K\left(i\omega_{n}\right)$ is not always available, and we can only calculate the correlation function $K\left(i\omega_{n}\right)$ at some Matsubara frequencies $i\omega_{n}$ numerically. According to a theorem proven by Baym and Mermin Baym-1961 , there exists a unique analytic continuation of the Matsubara correlation function, provided the values of $K\left(i\omega_{n}\right)$ for an infinite set of points, including the points at infinity. This is apparently impossible and we have to resort to some approximations.

Padé analytic continuation is based on the assumption that the Matsubara correlation function can be approximately represented by a rational function of degree $r$

[TABLE]

where $p_{r}$ and $q_{r}$ are complex coefficients, which can be determined by solving $2r$ linear equations from $2r$ arbitrary but different input points $\left\{i\omega_{n},K\left(i\omega_{n}\right)\right\}$ Beach-2000 . After the coefficients $\{p_{r},q_{r}\}$ are determined, we can replace $z$ by $\omega+i0^{+}$ to obtain the retarded correlation function $K^{R}(\omega)$. The procedure of this analytic continuation is illustrated in Fig. 1.

When the degree of the rational function $r$ is large, the ratio between the largest and the smallest input imaginary frequencies, $\omega_{n}$, could be very large. In this case, high precision calculation is needed to accurately determine the coefficient matrix of the $2r$ linear equations Beach-2000 . To understand this more clearly, let us take the correlation function determined by the particle-hole bubble diagram in one dimension,

[TABLE]

as an example to examine how the results of Padé analytic continuation are affected by the accuracy of calculations. In Eq. (4), $G\left(k,i\omega_{m}\right)=1/(i\omega_{m}-\xi_{k})$ is the single-particle Green function with $\xi_{k}=-2\cos k$. $\omega_{m}$ and $\omega_{n}$ are the fermionic and bosonic Matsubara frequencies. In the following we use $D$ to denote the decimal digits of the precision, and set $q=\pi$, temperature $T=0.1$, and the lattice size $N=50$. In the analytic continuation, we replace $z$ by $\omega+i\delta$ with $\delta=0.01$ in the Padé function (3).

For this simple correlation function, the Matsubara frequency summation in Eq. (4) can be carried out exactly. This gives a rigorous expression

[TABLE]

that can be used to benchmark the results obtained numerically with the Padé analytic continuation. Here $f\left(x\right)$ in Eq. (5) means Fermi distribution function. The analytic continuation of Eq. (5) can be carried out by replacing $i\omega_{n}$ with $\omega+i\delta$, and the corresponding spectral functions, which would serve as the exact results, can be obtained. For the exact results, we also set $\delta=0.01$, $T=0.1$ and the lattice size $N=50$.

To perform the Padé analytic continuation, we use Eq. (5) to generate $2r$ ($r=30$) data points $\left\{i\omega_{n},K\left(\pi,i\omega_{n}\right)\right\}$ with precision $D=16$ and $D=56$, respectively. The Padé coefficients are then determined by solving the $2r$ linear equations numerically. Fig. 2 compares the spectral function obtained by the Padé analytic continuation with the exact result obtained from Eq. (5). With the double precision calculation ($D=16$), we find that the Padé analytic continuation fails to reproduce all the fine structures of the spectral function except the peak structure at the highest frequency. By increasing the precision to $D=56$, the spectral function calculated by the Padé approximation begins to show the multi-peak structures of the spectral function. However, it still deviates significantly from the exact one. To reproduce accurately the exact result with a difference less than $10^{-15}$, we find that the precision has to be increased to $D=82$.

In the case that the Matsubara frequency summation cannot be evaluated analytically, we have to carry out the summation numerically. For a direct sum of the Matsubara frequency, one has to introduce a frequency cutoff. As the summation converges slowly, it is almost impossible to obtain the input data with sufficiently high precision that can be used in the Padé analytic continuation. This problem, however, can be solved by replacing the Matsubara frequency summation with a Padé frequency summation Taisuke-2007 ; Croy-2009 ; Hu-2010 ; Jie-2011 .

The Padé frequencies are determined by the Padé decomposition of the Fermi function

[TABLE]

where $\pm i\xi_{j}$ and $\eta_{j}$ represent the poles and the corresponding residues of the Fermi function in this decomposition scheme Taisuke-2007 . Table 1 lists the values of $\xi_{j}$ and $\eta_{j}$ up to $N=20$. Unlike the Matsubara decomposition, the poles of the Padé decomposition are unevenly distributed.

Using this Padé decomposition formula, we can rewrite Eq. (4) as

[TABLE]

which differs from Eq. (4) in two respects: (1) the Matsubara frequency summation is changed to the Padé frequency summation, and (2) the constant Matsubara residues are replaced by the Padé residues.

The Padé decomposition can greatly improve the accuracy in the calculation of the Matsubara frequency summation. For example, just using $N_{P}=70$ Padé frequencies, we can already calculate the value $K(\pi,i2\pi T)$ with a precision $D=100$. By contrast, the precision one can obtain for this quantity by directly summing up $N_{M}=10000$ Matsubara frequencies is just $D=3$.

The upper panel of Fig. 3 shows the spectral function of $A\left(\pi,\omega\right)$ obtained by the Padé analytic continuation with the input data calculated by directly summing over $N_{M}=600$ and $N_{M}=10000$ Matsubara frequencies using Eq. (4), respectively. The exact result is also shown in this figure for comparison. The difference between the results obtained by the Padé formula and the exact one is apparently very large. The improvement by increasing $N_{M}$ from 600 to 10000 is very small. The spectral function obtained with $N_{M}=10000$ even becomes negative in some regimes which is clearly unphysical. For comparison, the lower panel of Fig. 3 shows the spectral function obtained by the Padé analytic continuation with the input data calculated by summing over $N_{P}=600$ Padé frequencies using Eq. (7), which agrees perfectly with the exact one.

The accuracy of the Padé analytic continuation also depends on the choice of the power $r$ of the rational function. One should use a sufficiently large $r$ to carry out the analytic continuation so that the results are converged. As the denominator in the rational function (3) is a polynomial of order $r$, the number of its roots cannot exceed $r$. To reproduce accurately a Green’s function, $r$ should be equal to or larger than the number of the poles of that function. However, if $r$ is larger than the number of poles, the presence of the extra poles would also introduce errors in the calculation since it is difficult to make the contribution of these redundant poles cancel each other. In the case that the exact number of poles is unknown, it is desired to use a larger $r$ and to carry out the calculation with high precision. Some other schemes have also been proposed to deal with this so-called pole-zero pairs problem Beach-2000 ; Ostlin-2012 ; Osolin-2013 ; Schott-2016 .

In summary, we have shown that high precision input of the Matsubara Green function is needed to carry out the analytic continuation using the Padé rational polynomial approximation. When the summation over Matsubara frequencies cannot be carried out analytically, the replacement with summation over Padé frequencies can significantly improve the convergent behavior. This approach can be used to calculate the Matsubara Green function with high precision. It can be also used in the calculation of quantum Monte Carlo and dynamic mean-field theory, whenever a Matsubara frequency summation is present.

This work is supported by the National Natural Science Foundation of China (Grants No. 11474331 and No. 11190024)