Compressing Green's function using intermediate representation between imaginary-time and real-frequency domains

TL;DR

This paper introduces an intermediate representation (IR) for compressing imaginary-time data, significantly improving the efficiency of quantum Monte Carlo calculations for many-body systems.

Contribution

The authors develop a model-independent IR method that compactly represents correlation functions, enhancing the analysis of quantum Monte Carlo data.

Findings

IR provides a significantly more compact form of correlation functions.

IR improves the efficiency of quantum Monte Carlo calculations.

The framework can be broadly applied to many-body system analyses.

Abstract

New model-independent compact representations of imaginary-time data are presented in terms of the intermediate representation (IR) of analytical continuation. This is motivated by a recent numerical finding by the authors [J. Otsuki et al., arXiv:1702.03056]. We demonstrate the efficiency of the IR through continuous-time quantum Monte Carlo calculations of an Anderson impurity model. We find that the IR yields a significantly compact form of various types of correlation functions. The present framework will provide general ways to boost the power of cutting-edge diagrammatic/quantum Monte Carlo treatments of many-body systems.

Equations

$$\mathbf{G}$$

$$G(\tau )$$

$${\rho }_{\mathrm{F}}(\omega )$$

$${\rho }_{\mathrm{B}}(\omega )$$

$$K_{\mathrm{F}}(\tau ,\omega )$$

$$K_{\mathrm{B}}(\tau ,\omega )$$

$$K_{\mathrm{F}}(x,y)$$

$$K_{\mathrm{B}}(x,y)$$

$$K(x,y)$$

$$G(\tau )$$

$$\rho (\omega )$$

$$g_l=-s_l{\rho }_l.$$

$$\mathcal{H}$$

$$G_{\sigma }(\tau )$$

$${\chi }^{\mathrm{ch}}(\tau )$$

$$\tilde{\chi }({\tau }_{14},{\tau }_{24},{\tau }_{34})\equiv ⟨T_{\tau }{c}_{{\sigma }_1}^{†}({\tau }_1)c_{{\sigma }_2}({\tau }_2)c_{{\sigma }_3}^{†}({\tau }_3)c_{{\sigma }_4}({\tau }_4)⟩$$

$$-⟨T_{\tau }{c}_{{\sigma }_1}^{†}({\tau }_1)c_{{\sigma }_2}({\tau }_2)⟩⟨T_{\tau }{c}_{{\sigma }_3}^{†}({\tau }_3)c_{{\sigma }_4}({\tau }_4)⟩,$$

$${\tilde{\chi }}_{{\sigma }_1{\sigma }_2{\sigma }_3{\sigma }_4}({\tau }_{12},{\tau }_{34},{\tau }_{14})\equiv \sum _{l{l}^{\prime }\ge 0}\sum _{m\in \mathcal{Z}}\sqrt{\frac{2}{N_l}}\sqrt{\frac{2}{N_{l^{\prime }}}}{\beta }^{-3}(-1{)}^{l^{\prime }+1}$$

$$u_l[x({\tau }_{12})]u_{l^{\prime }}[x({\tau }_{34})]e^{\mathrm{i}{\omega }_m{\tau }_{14}}{\tilde{\chi }}_{l{l}^{\prime }}^{{\sigma }_1{\sigma }_2{\sigma }_3{\sigma }_4}(\mathrm{i}{\omega }_m).$$

$$K(x,y)$$

$$s_l$$

$$u_l(x)$$

$$v_l(y)$$

$$s_l$$

$$\mathbf{K}$$

$$s_l$$

$$=\int \mathrm{d}\tilde{x}\mathrm{d}\tilde{y}{\tilde{u}}_l^{\ast }(\tilde{x})\tilde{K}(\tilde{x},\tilde{y}){\tilde{v}}_l(\tilde{y}),$$

$${\tilde{u}}_l(\tilde{x})$$

$${\tilde{v}}_l(\tilde{y})$$

$$\tilde{K}(\tilde{x},\tilde{y})$$

$$x$$

$$y$$

$$K_{\mathrm{F}}(x,y)$$

$$K_{l{l}^{\prime }}$$

Full text

I Introduction

Many-body theories based on Matsubara Green’s function are powerful tools to study correlated systems. Elaborate diagrammatic methods have been widely used for investigating static and dynamic responses of the systems Baym and Kadanoff (1961); Baym (1962); Bickers and Scalapino (1989). Modern quantum Monte Carlo (QMC) methods even provide access to numerically exact ground-state and dynamical properties of lattice models and impurity models Werner et al. (2006); Gull et al. (2008, 2011); Sandvik and Kurkijärvi (1991); Rubtsov et al. (2005); Rombouts et al. (1999); Prokof’Ev et al. (1998); Beard and Wiese (1996); Sandvik and Kurkijärvi (1991); Iazzi and Troyer (2015); Wang et al. (2015); Blankenbecler et al. (1981); Kawashima and Harada (2004); Prokof’ev and Svistunov (1998). In these numerical calculations, however, one frequently faces two problems: (1) storage size and postprocessing cost of imaginary-time objects and (2) analytical continuation to the real-frequency axis.

The first issue becomes problematic in solving low-energy lattice models. For instance, one needs to treat two-particle quantities for computing lattice susceptibilities. Two-particle quantities also play a central role in some diagrammatic extensions of dynamical mean-field theory (DMFT) Georges et al. (1996) for describing non-local spatial correlations Toschi et al. (2007); Rubtsov et al. (2008). A recent technical advance is the compact representation of the imaginary-time dependence in terms of Legendre polynomials Boehnke et al. (2011). Efforts have been also made to describe the high-frequency asymptotic behavior of two-particle objects Li et al. (2016); Wentzell et al. (2016); Kuneš (2011). However, the application of these elaborate methods to realistic models is still too computationally expensive. A similar problem appears in quantum chemistry calculations based on a single-particle-level perturbative approach Kananenka et al. (2016); Rusakov and Zgid (2016). In this case, one needs to treat a much wider energy range than the low-energy models. Thus, there is a high demand for a more compact representation as a key ingredient in cutting-edge simulations of many-body systems.

The second problem is ill-conditioned analytical continuation from imaginary-time data to real-frequency axis. One example is estimating the spectral function from imaginary-time Green’s function computed in QMC simulations. The problem can be formulated as the linear equation

[TABLE]

where $\boldsymbol{G}$ and $\boldsymbol{\rho}$ are vectors representing imaginary-time and real-frequency data and the matrix $\boldsymbol{K}$ a kernel. Since $\boldsymbol{K}$ is ill-conditioned, the singular values of $\boldsymbol{K}$ decay very fast. As a result, most of independent components in $\boldsymbol{\rho}$ give almost no contribution to $\boldsymbol{G}$. Thus, if one simply minimizes $|\boldsymbol{G}+\boldsymbol{K}\boldsymbol{\rho}|^{2}$ with respect to $\boldsymbol{\rho}$, any errors in $\boldsymbol{G}$ are enormously amplified in $\boldsymbol{\rho}$.

The authors have recently developed a new method for analytical continuation of QMC data Otsuki et al. (2017). We demonstrated that, using a modern information theory called “sparse modeling”, relevant information can be successfully extracted from imaginary-time data with statistical errors. One of the key steps in this method is to transform the original data into a basis obtained by the singular value decomposition (SVD) of the matrix $\boldsymbol{K}$. As a result, after the errors are properly removed, the original data are expressed with only a few components. A similar observation was made in previous studies where SVD was employed in the context of analytical continuation CREFFIELD et al. (1995); Jarrell and Gubernatis (1996); Gunnarsson et al. (2010). These strongly suggest a possibility that this basis, which plays a key role in the analytical continuation, settles the first issue on the storage size and computational cost.

In this paper, we show that this model-independent basis can be used to compress various types of imaginary-time objects. We coin the term “intermediate representation (IR)” for this basis as it is defined between real-frequency and imaginary-time domains (Fig. 1). After investigating the properties of the IR in detail, we assess its efficiency for a single-site Anderson impurity model through continuous-time QMC simulations. We thus demonstrate that the IR provides significantly compact representations of the single-particle Green’s function, the charge susceptibility and the generalized susceptibility.

II Properties of basis functions

We start our discussion by considering the spectral (Lehmann) representation of a single-particle Green’s function $G$

[TABLE]

where we take $\hbar=1$ and $0\leq\tau\leq\beta$. This equation is reduced to Eq. (1) when the variables $\tau$ and $\omega$ are discretized. The spectra function $\rho(\omega)$ is given by

[TABLE]

or

[TABLE]

in the fermionic/bosonic case, respectively. The kernel is defined correspondingly

[TABLE]

or

[TABLE]

The extra $\omega$’s in Eqs. (4) and (6) are introduced to avoid a singularity of the kernel at $\omega=0$. We assume that the spectral function is bounded in the interval $[-\omega_{\mathrm{max}},~{}\omega_{\mathrm{max}}]$. Note that there is a similar spectral representation for the self-energy of a system of fermions Luttinger (1961).

For convenience, we transform the variables $\tau$ and $\omega$ into dimensionless variables $x\equiv 2\tau/\beta-1\in[-1,1]$ and $y\equiv\omega/\omega_{\mathrm{max}}\in[-1,1]$. The kernels are then rewritten as

[TABLE]

Here we introduced a dimensionless parameter $\Lambda\equiv\beta\omega_{\mathrm{max}}$. The IR changes its form depending on the value of $\Lambda$ as we will see below.

The IR is now defined through the decomposition of the kernels as

[TABLE]

This decomposition can be performed by SVD of a kernel matrix $\boldsymbol{K}$ defined on a dense uniform mesh: $\boldsymbol{K}=\boldsymbol{U}\boldsymbol{S}\boldsymbol{V}^{\dagger}$. In the continuous limit, column vectors of $\boldsymbol{U}$ ($\boldsymbol{V}$) yield orthonormal basis set $\{u_{l}(x)\}$ in the $\tau$ domain ($\{v_{l}(y)\}$ in the $\omega$ domain)  111 For a large $\Lambda$, some of zeros of $\{u_{l}(x)\}$ and $\{v_{l}(y)\}$ are distributed close to $x=\pm 1$ and $y=0$, respectively. For uniform meshes, this leads to a slow convergence of results with respect to the number of mesh points. To improve the convergence, we actually use a non-uniform mesh. Note that special care must be taken in the transformation of $\boldsymbol{K}$ to make sure that $\boldsymbol{U}$ and $\boldsymbol{V}$ converge to the same $\{u_{l}(x)\}$ and $\{v_{l}(y)\}$ obtained by uniform meshes. A technical note is given in the Appendix.. $s_{l}$ ($>0$) are singular values in non-ascending order. In the literature, SVD was utilized in the context of analytical continuation CREFFIELD et al. (1995); Jarrell and Gubernatis (1996); Gunnarsson et al. (2010). Our idea is to represent imaginary-time dependence by $u_{l}(x)$ to acquire compact forms of correlation functions.

We expand a given imaginary-time object $G(\tau)$ and the corresponding spectral function $\rho(\omega)$, respectively, in terms of $\{u_{l}(x)\}$ and $\{v_{l}(y)\}$ as

[TABLE]

Using Eqs. (2) and (9), we can demonstrate that the coefficients $g_{l}$ and $\rho_{l}$ have one-to-one correspondence

[TABLE]

Note that the singular values $\{s_{l}\}$ decay at least exponentially as shown in the upper panel of Fig. 2. It leads to an exponential decay of $g_{l}$, provided that $\rho_{l}$ does not grow for large $l$. In practical cases, we have confirmed that $\rho_{l}$ vanishes as $l$ increases, and $g_{l}$ decays even faster than $s_{l}$. The expansion in Eq. (10), therefore, may be truncated at a certain order, which will be demonstrated later using QMC simulations.

Here, we investigate the properties of the IR basis to get intuitive understanding why the expansion converges fast. Figure 2 shows the basis functions computed for the fermionic case. The real-frequency basis functions $v_{l}(y)$ have fine structure around $\omega=0$, which becomes shaper as $\Lambda$ is increased. This is consistent with that the kernel does not filter out the fine structure of a spectral function at small $\omega$. For $u_{l}(x)$, two notable features are clearly discernible: $u_{l}(x)$ is an even/odd function for even/odd $l$, and there are $l$ zeros 222A previous report on the number of zeros can be found in Ref. CREFFIELD et al., 1995.. More importantly, we found that $u_{l}(x)$ [and $v_{l}(y)$] converges to the $l$-th Legendre polynomial $P_{l}(x)$ up to a normalization factor as $\Lambda\rightarrow 0$ 333A note is given in the Appendix. In the note, we numerically and algebraically show that the first few basis functions converge to the Legendre polynomials in the high temperature limit.. It means that our representation using the IR basis includes the Legendre representation as a special limit. However, since this limit corresponds to the high-$T$ limit, the Legendre expansion may not be efficient for low $T$. A difference between the IR basis and the Legendre polynomials becomes clear as $\Lambda$ is increased: The values of $u_{0}(x)$ and $u_{1}(x)$ change more rapidly around $x=\pm 1$, which resemble the behavior of diagonal and off-diagonal elements of the Green’s function, respectively. Moreover $u_{0}(0)$ becomes suppressed, similarly to the low-$T$ behavior of the diagonal elements $G(\tau=\beta/2)\propto T$. Therefore, an efficient descriptions in terms of the IR basis is expected especially at low $T$.

III Results of quantum Monte Carlo simulations

III.1 Model

Now we demonstrate the efficiency of the IR for describing various types of imaginary-time objects. As a simple example, we consider the particle-hole symmetric single-site Anderson impurity model defined by the Hamiltonian

[TABLE]

with $\mu=U/2$ and $\sigma$ is spin index. $c_{\sigma}$ and $c^{\dagger}_{\sigma}$ are annihilation and creation operators at the impurity site, while $a_{k\sigma}$ and $a^{\dagger}_{k\sigma}$ are those of the bath sites ($k$ is the internal degree of freedom of the bath). The distribution of $\epsilon_{k}$ is a semicircular density of states of width $4$. We solve the model and compute correlation functions by means of the hybridization expansion continuous-time Monte Carlo technique Werner et al. (2006).

III.2 Single-particle Green’s function

First, we discuss the impurity single-particle Green’s function defined as $G_{\sigma}(\tau)=-\braket{c_{\sigma}(\tau)c^{\dagger}_{\sigma}(0)}$ ($0\leq\tau\leq\beta$). We expand $G_{\sigma}(\tau)$ in terms of an orthogonal basis set $\{f_{l}(x)\}$ [$P_{l}$ or $u_{l}$] as

[TABLE]

where $x(\tau)=2\tau/\beta-1$ and $\int_{-1}^{1}f_{l}(x)f_{l^{\prime}}(x)dx=N_{l}\delta_{ll^{\prime}}$. We directly measure the coefficients $G_{l}^{\sigma}$ in QMC simulations as described in Ref. Boehnke et al., 2011.

In Fig. 3, we show the coefficients $G_{l}$ obtained for $U=4$ and $\beta=100$. The large-$l$ asymptomatic behavior of the Legendre representation is known to be exponential Boehnke et al. (2011), while the Matsubara-frequency representation has a $1/\mathrm{i}\omega_{n}$ tail. As expected, the IR yields coefficients decaying even faster than the Legendre basis. One can expect that the most compact representation is obtained when $\Lambda/\beta$ matches the actual width of the spectrum. This suggests a practical way to choose an appropriate value of $\Lambda$. Actually, the optimal value obtained is $\Lambda\simeq 1000$ for $\beta=100$, being consistent with the largest dimensionless energy scale of the system, i.e., $\beta U$, $\beta W=400$. As $\Lambda$ exceeds the optimal value, the efficiency gets worse only slowly. In particular, we observed the non-monotonic behavior of $G_{l}$ around $l=5$ for $\Lambda>500$, which signals that $\Lambda$ exceeds an optimal value.

In Fig. 3, we also show $G(\tau)$ reconstructed from the coefficients for $l\leq 6$. The data obtained by the IR ($\Lambda=500$) shows a perfect agreement with the numerically exact data, while the truncation in the Legendre representation results in large Gibbs oscillations.

III.3 Charge susceptibility

Second, as a typical bosonic quantity, we analyze the charge susceptibility $\chi^{\mathrm{ch}}(\tau)$ defined by

[TABLE]

on the interval $[0,\beta]$ ($n_{\sigma}\equiv\sum_{\sigma}c^{\dagger}_{\sigma}c_{\sigma}$). We expand the $\tau$ dependence using Eq. (14) in terms of the bosonic IR or the fermionic IR. Strange as the latter may sound, it is possible since the basis functions $v_{l}(x)$ always form a complete basis set on the interval $[-1,1]$. Figure 4 shows the results obtained for $\beta=100$. Remarkably, the bosonic IR requires only few coefficients beyond statistical errors. Again, the most compact representation is obtained when $\Lambda/\beta$ matches the spectral width. On the other hand, surprisingly, the fermionic IR is better than the Legendre basis [the lower panel of Fig. 4]. However, the most compact representation is obtained with the bosonic IR with $\Lambda\simeq 500$. This shows the importance of using the bosonic IR for bosonic quantities.

III.4 Two-particle Green’s function

Finally, we demonstrate that the IR’s for single-particle Green’s functions can be used for expanding objects with multiple time indices. As an example, we consider the generalized susceptibility defined by

[TABLE]

where $\tau_{ab}\equiv\tau_{a}-\tau_{b}$ and the second term subtracts the trivial contribution of the bubble diagram. In Ref. Boehnke et al., 2011, Boehnke et al. introduced the mixed representation of Legendre polynomials and bosonic Matsubara frequencies, in which $\tau_{12}$ and $\tau_{34}$ dependence of a two-particle object is expanded in terms of the Legendre polynomials, while the $\tau_{14}$ dependence is described through Fourier modes $e^{\mathrm{i}\omega_{m}\tau_{14}}$ Boehnke et al. (2011). This is motivated by that fact that the Bethe-Salpeter equation is diagonal in the bosonic frequency $\mathrm{i}\omega_{m}$ that is connected to $\tau_{14}$ through a Fourier transformation.

We now simply replace the Legendre polynomials by the IR for the fermionic kernel. This leads to

[TABLE]

Hereafter, we consider only the spin diagonal components $\tilde{\chi}_{\uparrow\uparrow\uparrow\uparrow}~{}(=\tilde{\chi}_{\downarrow\downarrow\downarrow\downarrow})$ and drop the spin indices. In practice, we measure the coefficients of the first term in Eq. (16) in QMC simulations and subsequently subtract the bubble-diagram contribution (second term) computed from the data of the single-particle Green’s function. The accumulation of the coefficients for a single bosonic frequency requires $O(k^{4}N^{2})$ operations, where $k$ is the expansion order of QMC and $N$ is the number of the IR basis functions or Legendre polynomials for fermionic frequencies. Thus, any reduction of $N$ will significantly speed-up QMC measurement. We refer the readers to Ref. Boehnke et al., 2011 for more technical details.

We show the results at the zero bosonic frequency $\mathrm{i}\omega_{m}=0$ for $\beta=25$ and $100$ in Fig. 5. In the mixed representation of the IR, the coefficients for even $l$ and $l^{\prime}$ show a much faster decay than those for the Legendre basis. The other coefficients for odd $l$ or odd $l^{\prime}$ take on a very small value, compatible with a vanishing value within statistical errors. As $\beta$ increases, the decay becomes slower for the Legendre basis. On the other hand, for the IR basis, the rate of the decay depends on the value of $\Lambda$ slightly. But one does not observe a noticeable slowdown in the decay if the value of $\Lambda$ is chosen appropriately. As a result, the new basis becomes more superior as $\beta$ increases.

An interesting observation is that the most compact representation is obtained for a value of $\Lambda$ close to the optimal one for the single-particle Green’s function for $\beta=100$.

To demonstrate advantages of the new basis functions in practical QMC calculations, we measured the computational time of the measurements of two-particle Green’s function. The results for $\beta=100$ are shown in Fig. 6. We employed 50 (Legendre), 10 ($\Lambda=500$), 20 ($\Lambda=10000$) basis functions in the measurement so that the data beyond the noise level are accumulated (see Fig. 5). The simulations were performed on a 2.5GHz Intel Xeon CPU (E5-2680 v3) without parallelization. One can see that the measurement for $\Lambda=500$ is faster than the case of the Legendre polynomials approximately by 24 times, being consistent with the expected scaling.

IV Summary

In summary, we proposed the new compact representations of imaginary-time data, which was named IR, through the lenses of analytical continuation. The new basis does not depend on the details of the systems. In particular, we studied the properties of the IR’s for fermionic and bosonic Green’s functions. We found that the conventional Legendre basis corresponds to the high-$T$ limit of the IR. The IR was applied to QMC simulations of the single-site quantum impurity model. We confirmed that the present method yields significantly compact form of various imaginary-time correlation functions than the conventional ones. This allows the direct measurement of Green’s functions in a compressed form, which reduces the computational cost and memory usage.

An optimal value of $\Lambda$ depends on temperature and the width of spectral function. The numerical tests indicate that the data remain compact even when $\Lambda$ exceeds the optimal value. Thus, one does not have to tune the value of $\Lambda$ very precisely. It may be practically efficient enough to let the value of $\Lambda$ be on the large side.

The present scheme provides a new approach to solve technical issues in a variety of state-of-the-art treatments of many-body quantum systems. For instance, one may be able to perform diagrammatic calculations with Bethe-Salpeter/parquet equations in the IR. Promising applications are the diagrammatic extensions of DMFT (dual fermions Toschi et al. (2007) and dynamical vertex approximation Rubtsov et al. (2008)) and the computation of lattice susceptibilities within DMFT. On the other hand, the kernel for the Keldysh Green’s function is also known to be ill-conditioned Dirks et al. (2013). The present scheme may be easily applied to non-equilibrium cases. Furthermore, a modern regularization technique will enable to separate relevant information from statistical noise in the IR Otsuki et al. (2017).

Acknowledgements.

We are grateful to Emanuel Gull, Kristjan Haule, Yusuke Nomura and Philipp Werner for their useful comments on the manuscript. HS thanks fruitful discussions with Lewin Boehnke on the Legendre basis.. HS was supported by JSPS KAKENHI Grant No. 16H01064 (J-Physics), 16K17735. JO was supported by JSPS KAKENHI Grant No. 26800172, 16H01059 (J-Physics). MO was supported by MEXT KAKENHI Grant No. 25120008, JST CREST and JSPS KAKENHI No. 16H04382. KY was supported by Building of Consortia for the Development of Human Resources in Science and Technology, MEXT, Japan. Part of the calculations were run on the ISSP supercomputing system with codes based the ALPSCore libraries Gaenko et al. (2016) and ALPSCore/CT-HYB Shinaoka et al. (2017).

Appendix A Decomposition of the kernel

A.1 Singular value decomposition by means of a linear mesh

We start our discussion with

[TABLE]

This equation is recast into

[TABLE]

Now, we introduce equally-spaced $N$ points $x_{n}$ and $y_{m}$ on the interval $[-1,1]$ ($n,m=0,\cdots,N-1$). We approximate $u_{l}(x)$ and $v_{l}(y)$ as

[TABLE]

Substituting these into Eq. (19) leads to

[TABLE]

where $\boldsymbol{u}_{l}$ and $\boldsymbol{v}_{l}$ are the column vectors whose elements are $u_{l,n}$ and $v_{l,m}$, respectively. For $N\gg 1$, the orthonormal condition of $u_{l}(x)$ and $v_{l}(y)$ is equivalent to that of the column vectors. The matrix element of $\boldsymbol{K}$ at $(n,m)$ is given by $\frac{4}{N^{2}}K(x_{n},y_{m})$. Such vectors can be computed by a singular value decomposition (SVD) of the matrix $\boldsymbol{K}$ as

[TABLE]

where $\boldsymbol{S}$ is the diagonal matrix whose elements are given by $s_{l}$, $\boldsymbol{U}=(\boldsymbol{u}_{0},\cdots,\boldsymbol{u}_{N-1})$ and $\boldsymbol{V}=(\boldsymbol{v}_{0},\cdots,\boldsymbol{v}_{N-1})$.

A.2 Double-exponential mesh

The linear mesh is not optimal because the values of the basis functions change very rapidly around $x=\pm 1$ and $y=0$. Let us consider the change of variables $x=f(\tilde{x})$ and $y=g(\tilde{y})$. Then, Eq. (18) reads

[TABLE]

where

[TABLE]

Here, we assume $f^{\prime}(\tilde{x})\geq 0$ and $g^{\prime}(\tilde{y})\geq 0$. Note that $\tilde{u}_{l}(\tilde{x})$ and $\tilde{v}_{l}(\tilde{y})$ are orthonormal functions with respect to $\tilde{x}$ and $\tilde{y}$, respectively. One can also solve Eq. (25) in a way analogous to the solution of Eq. (25).

This provides the possibility of choosing appropriate transformations so that we have more dense points in the regions where the values of the basis functions change rapidly, i.e., $x=\pm 1$ and $y=0$. In particular, we adopt the double-exponential transformations Takahasi and Mori (1974)

[TABLE]

with the cutoff $|\tilde{x}|\leq 4$ and $|\tilde{y}|\leq 4$. This transformation maps $x\in[-1,1]$ to $\tilde{x}\in[-\infty,\infty]$. The cutoff can be introduced very safely because the derivative of $f(\tilde{x})$ and $g(\tilde{y})$ show a double-exponential decay. We found that $N=1001$ gives sufficiently accurate solutions for our purpose. For more technical details, please study the Python scripts provided in the supplemental material.

Appendix B Asymptotic behavior of basis functions of Intermediate representation (IR)

In Fig. 7, we show that the basis functions $u_{l}(x)$ and $v_{l}(y)$ converge to the Legendre polynomials in the limit of $\Lambda\rightarrow 0$ 444 We provide a python script for generating Fig. 7 at a public git repository Shinaoka and in the Suppelmental Material.. This is clearly seen in the numerical data shown in Fig. 7. Due to the fast decay of the singular values $s_{l}$, it is numerically difficult to compute the basis functions for large $l$ accurately when $\Lambda\ll 1$.

To study the asymptotic behavior for $\Lambda\rightarrow 0$ more precisely, we algebraically expand the fermionic kernel [Eq. (7)] in terms of the Legendre polynomials as

[TABLE]

In practical, we first expand $K_{\mathrm{F}}(x,y)$ in powers of $x$ and $y$. Then, we compute the expansion coefficients in Eq. (32) by performing the integration over $x$ and $y$. For $N=4$, we obtained

[TABLE]

The basis functions of the IR can be computed by diagonalizing the matrix $\boldsymbol{K}^{T}\boldsymbol{K}$ ($y$ space) and the matrix $\boldsymbol{K}\boldsymbol{K}^{T}$ ($x$ space), respectively. The singular values of $\boldsymbol{K}$ are given by

[TABLE]

This suggests that the $l$-th singular value is $O(\Lambda^{l})$.

The matrix representation of the basis functions for $y$ reads

[TABLE]

where the $l$-th column corresponds to the $l$-th basis function (up to sign factors). Similarly, the basis functions for $x$ reads

[TABLE]

One can clearly see that these basis functions converge to the Legendre polynomials as $\Lambda\rightarrow 0$.