⇦Chinese Physics Letters, 2017, Vol. 34, No. 7, Article code 077102 Analytic Continuation with Padé Decomposition * Xing-Jie Han(韩兴杰)1,2, Hai-Jun Liao(廖海军)1,2, Hai-Dong Xie(谢海东)1,2, Rui-Zhen Huang(黄瑞珍)1,2, Zi-Yang Meng(孟子杨)1,2, Tao Xiang(向涛)1,2,3** Affiliations 1Institute of Physics, Chinese Academy of Sciences, Beijing 100190 2University of Chinese Academy of Sciences, Beijing 100049 3Collaborative Innovation Center of Quantum Matter, Beijing 100190 Received 11 April 2017 *Supported by the National Natural Science Foundation of China under Grant Nos 11474331 and 11190024.
**Corresponding author. Email: txiang@iphy.ac.cn Citation Text: Han X J, Liao H J, Xie H D, Huang R Z and Meng Z Y et al 2017 Chin. Phys. Lett. 34 077102 Abstract The ill-posed analytic continuation problem for Green's functions or self-energies can be carried out using the Padé rational polynomial approximation. However, to extract accurate results from this approximation, high precision input data of the Matsubara Green function are needed. The calculation of the Matsubara Green 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é decomposition approach to replace the Matsubara frequency summation by a Padé frequency summation, and high precision input data can be obtained to successfully perform the Padé analytic continuation. DOI:10.1088/0256-307X/34/7/077102 PACS:71.15.Dx, 02.70.Hm © 2017 Chinese Physics Society Article Text The theory of the Matsubara Green function provides a convenient and powerful tool to calculate physical properties at finite temperatures.[1] 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)[2,3] contain stochastic noise induced by random sampling. For this kind of data, various numerical algorithms, such as the maximum entropy,[4-6] the singular value decomposition[7,8] and the stochastic analytical continuation,[9-11] 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[12-16] can be used. In this method, the analytic continuation is carried out by interpolating the input data 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 to use the Padé formula.[13] 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 functions.[17-20] An efficient approach along this line is the Padé decomposition first introduced by Ozaki,[17] 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^{\rm 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^{\rm R}(\omega)$ can be obtained from $K(i\omega_n)$ by analytic continuation $$ K^{\rm R}(\omega)=K(i \omega_n \rightarrow \omega + i 0^+),~~ \tag {1} $$ provided that the analytic expression of $K(i\omega_{n})$ is known. The spectral function $A(\omega)$ is determined by the imaginary part of $K^{\rm R}(\omega +i0^+)$ by the formula $$ A(\omega) =- \frac{1}{\pi} {\rm Im} K^{\rm R}(\omega +i0^+).~~ \tag {2} $$ However, as already mentioned, the analytic expression of $K(i\omega_{n})$ is not always available, and we can only calculate the correlation function $K(i\omega_{n})$ at some Matsubara frequencies $i\omega_{n}$ numerically. According to a theorem proven by Baym and Mermin,[21] there exists a unique analytic continuation of the Matsubara correlation function, provided the values of $K(i\omega_{n})$ for an infinite set of points, including the points at infinity. This is apparently impossible and we have to resort to some approximations.
cpl-34-7-077102-fig1.png
Fig. 1. Procedure of the Padé analytic continuation: the Padé rational function $K(z)$, defined by Eq. (3), is first determined from a set of input data $\{i\omega_{n},K(i\omega_{n})\}$, and then the retarded Green function $K^{\rm R}(\omega +i0^+)$ is obtained from $K(z)$ by substituting $z$ with $\omega +i0^+$.
The Padé analytic continuation is based on the assumption that the Matsubara correlation function can be approximately represented by a rational function of degree $r$, $$ K(z) =\frac{p_{0}+p_{1}z+\ldots +p_{r-1}z^{r-1}}{q_{0}+q_{1}z+\ldots +q_{r-1}z^{r-1}+z^{r}},~~ \tag {3} $$ 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 $\{i\omega_{n}, K(i\omega_{n})\}$.[13] After the coefficients $\{p_{r},q_{r}\}$ are determined, we can replace $z$ by $\omega+i0^+$ to obtain the retarded correlation function $K^{\rm 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.[13] To understand this more clearly, let us take the correlation function determined by the particle-hole bubble diagram in one dimension, $$ 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}),~~ \tag {4} $$ as an example to examine how the results of Padé analytic continuation are affected by the accuracy of calculations. In Eq. (4), $G(k,i\omega_{m})=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 $$ K(q,i\omega_{n}) =\frac{1}{N}\sum_{k}\frac{f(\xi _{k}) -f(\xi _{k+q}) }{i\omega_{n}+\xi _{k}-\xi _{k+q}},~~ \tag {5} $$ which can be used to benchmark the results obtained numerically with the Padé analytic continuation. Here $f(x)$ in Eq. (5) means the 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$.
cpl-34-7-077102-fig2.png
Fig. 2. The spectral function of the particle-hole correlation function, $A(\pi, \omega)$, obtained by the Padé analytic continuation with $D=16$ (red dashed line), $D=56$ (blue dashed line) and $D=82$ (green dashed line), respectively. The exact result (black solid line) is also shown for comparison.
To perform the Padé analytic continuation, we use Eq. (5) to generate $2r$ ($r=30$) data points $\{i\omega_{n},K(\pi,i\omega_{n})\}$ with precision $D=16$ and $D=56$, respectively. The Padé coefficients are then determined by solving the $2r$ linear equations numerically. Figure 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.[17-20] The Padé frequencies are determined by the Padé decomposition of the Fermi function $$ f(x) =\frac{1}{e^{x}+1}=\frac{1}{2}-\sum_{j=1}^{N\rightarrow \infty }\frac{2\eta _{j}x}{x^{2} + \xi _{j}^{2}},~~ \tag {6} $$ where $\pm i\xi _{j}$ and $\eta _{j}$ represent the poles and the corresponding residues of the Fermi function in this decomposition scheme.[17] 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.
Table 1. Values of the first 20 pairs of $\eta_j$ and $\xi_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
Using this Padé decomposition formula, we can rewrite Eq. (4) as $$\begin{align} K(q,i\omega_{n}) =\,&\frac{1}{\beta N}\sum_{k}\sum_{j=-N_{\rm P}}^{N_{\rm P}} \eta _{j}G^{(0) }(k,i\xi _{j}/\beta)\\ &\cdot G^{(0)}(q+k,i\omega_{n}+i\xi _{j}/\beta),~~ \tag {7} \end{align} $$ 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.
cpl-34-7-077102-fig3.png
Fig. 3. Comparison of the spectral function $A(\pi, \omega)$ obtained by the Padé analytic continuation with the exact one (black solid lines). The input data $\{i\omega_{n},K(\pi,i\omega_{n})\} $ used in the Padé analytic continuation are calculated by summing $N_{\rm M}=600$ (red solid line) and $N_{\rm M}=10000$ (blue solid line) Matsubara frequencies (upper panel), and $N_{\rm P}=600$ (red dashed line) Padé frequencies (lower panel), respectively.
The Padé decomposition can greatly improve the accuracy in the calculation of the Matsubara frequency summation. For example, just using $N_{\rm 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_{\rm M}=10000$ Matsubara frequencies is just $D=3$. The upper panel of Fig. 3 shows the spectral function of $A(\pi,\omega)$ obtained by the Padé analytic continuation with the input data calculated by directly summing over $N_{\rm M}=600$ and $N_{\rm 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_{\rm M}$ from 600 to 10000 is very small. The spectral function obtained with $N_{\rm 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_{\rm 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.[13-16] 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. References Monte Carlo calculations of coupled boson-fermion systems. IDiscrete Hubbard-Stratonovich transformation for fermion lattice modelsBayesian inference and the analytic continuation of imaginary-time quantum Monte Carlo dataAlgorithms for optimized maximum entropy and diagnostic tools for analytic continuationNumerical analytic continuation: Answers to well-posed questionsSpectral Weight Function for the Half-Filled Hubbard Model: A Singular Value Decomposition ApproachAnalytical continuation of imaginary axis data for optical conductivityIdentifying the maximum entropy method as a special limit of stochastic analytic continuationConstrained sampling method for analytic continuationAmplitude Mode in Three-Dimensional Dimerized AntiferromagnetsSolving the Eliashberg equations by means ofN-point Pad� approximantsReliable Padé analytical continuation method based on a high-accuracy symbolic computation algorithmOne-particle spectral function and analytic continuation for many-body implementation in the exact muffin-tin orbitals methodPad? approximant approach for obtaining finite-temperature spectral functions of quantum impurity models using the numerical renormalization group techniqueAnalytic continuation by averaging Pad? approximantsContinued fraction representation of the Fermi-Dirac function for large-scale electronic structure calculationsPartial fraction decomposition of the Fermi functionCommunication: Pad? spectrum decomposition of Fermi function and Bose functionPad? spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systemsDetermination of Thermodynamic Green's Functions
[1]Mahan G D 2000 Many Particle Physics (New York: Kluwer Academic/Plenum Publishers)
[2] Blankenbecler R, Scalapino D J and Sugar R L 1981 Phys. Rev. D 24 2278
[3] Hirsch J E 1983 Phys. Rev. B 28 4059(R)
[4] Jarrell M and Gubernatis J E 1996 Phys. Rep. 269 133
[5] Bergeron D and Tremblay A M S 2016 Phys. Rev. E 94 023303
[6] Goulko O, Mishchenko A S, Pollet L, Prokofév N and Svistunov B 2017 Phys. Rev. B 95 014102
[7] Creffield C E, Klepfish E G, Pike E R and Sarkar S 1995 Phys. Rev. Lett. 75 517
[8] Gunnarsson O, Haverkort M W and Sangiovanni G 2010 Phys. Rev. B 82 165125
[9] Beach K S D 2004 arXiv:cond-mat/0403055[cond-mat.str-el]
[10] Sandvik A W 2016 Phys. Rev. E 94 063308
[11] Qin Y Q, Normand B, Sandvik A W and Meng Z Y 2017 Phys. Rev. Lett. 118 147207
[12] Vidberg H J and Serene J 1977 J. Low Temp. Phys. 29 179
[13] Beach K S D, Gooding R J and Marsiglio F 2000 Phys. Rev. B 61 5147
[14] Östlin A, Chioncel L and Vitos L 2012 Phys. Rev. B 86 235107
[15] Osolin Ž and Žitko R 2013 Phys. Rev. B 87 245135
[16] Schött J, Locht I L M, Lundin E, Grånäs O, Eriksson O and Marco I D 2016 Phys. Rev. B 93 075104
[17] Ozaki T 2007 Phys. Rev. B 75 035123
[18] Croy A and Saalmann U 2009 Phys. Rev. B 80 073102
[19] Hu J, Xu R X and Yan Y J 2010 J. Chem. Phys. 133 101106
[20] Hu J, Luo M, Jiang F, Xu R X and Yan Y J 2011 J. Chem. Phys. 134 244106
[21] Baym G and Mermin D 1961 J. Math. Phys. 2 232