Bayesian parametric analytic continuation of Green's functions

TL;DR

This paper introduces Bayesian parametric analytic continuation (BPAC), a method for inferring spectral functions from noisy Green's function data, providing a more reliable alternative to the maximum entropy method, especially for complex spectral structures.

Contribution

The paper presents BPAC, a novel Bayesian approach using model comparison and nested sampling for analytic continuation of Green's functions, improving reliability over traditional methods.

Findings

BPAC accurately reconstructs spectral functions from noisy data.

BPAC can assess the support for specific spectral structures.

Application to AIM demonstrates effectiveness in complex scenarios.

Abstract

Bayesian parametric analytic continuation (BPAC) is proposed for the analytic continuation of noisy imaginary-time Green's function data as, e.g., obtained by continuous-time quantum Monte Carlo simulations (CTQMC). Within BPAC, the spectral function is inferred from a suitable set of parametrized basis functions. Bayesian model comparison then allows to assess the reliability of different parametrizations. The required evidence integrals of such a model comparison are determined by nested sampling. Compared to the maximum entropy method (MEM), routinely used for the analytic continuation of CTQMC data, the presented approach allows to infer whether the data support specific structures of the spectral function. We demonstrate the capability of BPAC in terms of CTQMC data for an Anderson impurity model (AIM) that shows a generalized Kondo scenario and compare the BPAC reconstruction to the MEM as well as to the spectral function obtained from the real-time fork tensor product state impurity solver where no analytic continuation is required. Furthermore, we present a combination of MEM and BPAC and its application to an AIM arising from the ab initio treatment of SrVO$_3$.

Tables

Table 1. Table 1: Prior ranges and estimated parameters for model M 4 subscript 𝑀 4 M_{4} for the spin-down spectral function.

	${\mu }_n$	$a_n$	${\gamma }_n^{\text{l}}$	${\gamma }_n^{\text{r}}$
prior range	(-0.4, -0.2)	(0, 1)	(0.001, 0.3)	(0.001, 0.3)
estimator	-0.32 $\pm$ 0.04	0.7 $\pm$ 0.1	0.20 $\pm$ 0.06	0.18 $\pm$ 0.06
prior range	(-0.07, -0.03)	(0, 0.5)	(0.001, 0.1)	(0.001, 0.1)
estimator	-0.05 $\pm$ 0.01	0.08 $\pm$ 0.04	0.06 $\pm$ 0.03	0.04 $\pm$ 0.02
prior range	(-0.02, 0.02)	(0, 0.5)	(0.001, 0.1)	(0.001, 0.1)
estimator	0.001 $\pm$ 0.003	0.022 $\pm$ 0.005	0.008 $\pm$ 0.004	0.006 $\pm$ 0.003
prior range	(0.1, 0.25)	(0, 1)	(0.001, 0.3)	(0.001, 0.3)
estimator	0.22 $\pm$ 0.02	0.6 $\pm$ 0.1	0.06 $\pm$ 0.01	0.09 $\pm$ 0.05

Table 2. Table 2: Prior ranges and estimated parameters for model M 2 subscript 𝑀 2 M_{2} for the spin-up spectral function.

	${\mu }_n$	$a_n$	${\gamma }_n^{\text{l}}$	${\gamma }_n^{\text{r}}$
prior range	( 0.03, 0.07)	( 0, 0.1)	(0.001, 0.1)	(0.001, 0.1)
estimator	0.063 $\pm$ 0.001	0.078 $\pm$ 0.002	0.0116 $\pm$ 0.0004	0.016 $\pm$ 0.002
prior range	( 0.2, 0.3)	( 0, 1.5)	(0.001, 0.5)	(0.001, 0.5)
estimator	0.218 $\pm$ 0.002	1.06$\pm$0.03	0.0320 $\pm$ 0.0007	0.100 $\pm$ 0.007

Equations

$$G(\tau )=\frac{1}{\beta }\sum _{{\omega }_n}{e}^{-\mathrm{i}{\omega }_n\tau }\mathcal{G}(\mathrm{i}{\omega }_n).$$

$$\mathcal{G}(z)=\mp \int d\omega \frac{A(\omega )}{z-\omega }.$$

$$G(\tau )=\int d\omega \frac{e^{-\omega \tau }}{e^{-\beta \omega }+1}A(\omega )=\int d\omega K(\tau ,\omega )A(\omega ).$$

$$S(A)=\int d\omega \left(A(\omega )-D(\omega )-A(\omega )\ln \frac{A(\omega )}{D(\omega )}\right)$$

$$f_n(\omega |{\tilde{\alpha }}_n)=\{\begin{array}{ll}\sum _{i=1}^{C_n}{a}_n^i\frac{2{\gamma }_n^{i,\text{l}}}{{\gamma }_n^{i,\text{l}}+{\gamma }_n^{i,\text{r}}}L(\omega |{\mu }_n,{\gamma }_n^{i,\text{l}})& \text{forleavevmodenobreak }\omega <{\mu }_n\\ \sum _{i=1}^{C_n}{a}_n^i\frac{2{\gamma }_n^{i,\text{r}}}{{\gamma }_n^{i,\text{l}}+{\gamma }_n^{i,\text{r}}}L(\omega |{\mu }_n,{\gamma }_n^{i,\text{r}})& \text{forleavevmodenobreak }{\mu }_n\le \omega \end{array}.$$

$$L(\omega |\mu ,\gamma )$$

$$A(\omega |\mathbf{\alpha })=\{\begin{array}{ll}{a}^{\text{l}}|\omega -{\mu }^{\text{l}}{|}^{-{\nu }^{\text{l}}}& \text{forleavevmodenobreak }\omega <{\Omega }^{\text{l}}\\ {\sum }_{n=1}^{N_{\text{p}}}{f}_n(\omega |{\tilde{\alpha }}_n)& \text{forleavevmodenobreak }{\Omega }^{\text{l}}\le \omega \le {\Omega }^{\text{r}}\\ a^{\text{r}}|\omega -{\mu }^{\text{r}}{|}^{-{\nu }^{\text{r}}}& \text{forleavevmodenobreak }{\Omega }^{\text{r}}<\omega \end{array}.$$

$$\underset{\text{posterior}}{\underbrace{p(\mathbf{\alpha }|{\mathbf{G}}_{\text{d}},\Sigma ,M,\mathcal{I})}}\underset{\text{data evidence}}{\underbrace{p({\mathbf{G}}_{\text{d}}|M,\Sigma ,\mathcal{I})}}$$

$$L(\mathbf{\alpha }):$$

$${\chi }^2(\mathbf{\alpha })$$

$$\pi (\mathbf{\alpha })$$

$$=\prod _{i=1}^{N_{\alpha }}\frac{1}{{\alpha }_i^{\text{max}}-{\alpha }_i^{\text{min}}}\Theta ({\alpha }_i^{\text{min}}<{\alpha }_i<{\alpha }_i^{\text{max}}).$$

$$p(\mathbf{A}|{\mathbf{G}}_{\text{d}},\Sigma ,M,\mathcal{I})=\int d\mathbf{\alpha }\underset{\delta (\mathbf{A}-\mathbf{A}(\mathbf{\alpha }))}{\underbrace{p(\mathbf{A}|\mathbf{\alpha },M,\mathcal{I})}}p(\mathbf{\alpha }|{\mathbf{G}}_{\text{d}},\Sigma ,M,\mathcal{I}).$$

$$P(M|{\mathbf{G}}_{\text{d}},\Sigma ,\mathcal{I})=\frac{1}{p({\mathbf{G}}_{\text{d}}|\mathcal{I})}p({\mathbf{G}}_{\text{d}}|\Sigma ,M,\mathcal{I})P(M|\mathcal{I}).$$

$$\mathcal{O}=\frac{P(M_1|{\mathbf{G}}_{\text{d}},\Sigma ,\mathcal{I})}{P(M_2|{\mathbf{G}}_{\text{d}},\Sigma ,\mathcal{I})}=\underset{\text{Bayes factor}}{\underbrace{\frac{p({\mathbf{G}}_{\text{d}}|M_1,\Sigma ,\mathcal{I})}{p({\mathbf{G}}_{\text{d}}|M_2,\Sigma ,\mathcal{I})}}}\underset{\text{prior odds}}{\underbrace{\frac{P(M_1|\mathcal{I})}{P(M_2|\mathcal{I})}}}.$$

$$p({\mathbf{G}}_{\text{d}}|M,\Sigma ,\mathcal{I})$$

$$=\int d\lambda X(\lambda ),$$

$$X(\lambda )=\int d\mathbf{\alpha }\pi (\mathbf{\alpha })\Theta \left(L(\mathbf{\alpha })>\lambda \right)$$

$$p({\mathbf{G}}_{\text{d}}|M,\Sigma ,\mathcal{I})$$

$$\approx \sum _{n=0}^{\infty }\Delta X_n\mathcal{L}(X_n)$$

$$\approx \sum _{n=0}^{n_{\text{max}}}\Delta X_n{\lambda }_n^{\ast },$$

$$\pi (\mathbf{\alpha }|{\lambda }_n^{\ast })=\frac{\pi (\mathbf{\alpha })}{X({\lambda }_n^{\ast })}\Theta (L(\mathbf{\alpha })>{\lambda }_n^{\ast }).$$

$${\hat{H}}_{\text{AI}}=\sum _{i\sigma }{ϵ}_{i\sigma }{\hat{n}}_{i\sigma }+\frac{1}{2}\sum _{ij\sigma }{U}_{ij}{\hat{n}}_{i\sigma }{\hat{n}}_{j\bar{\sigma }}+\frac{1}{2}\sum _{i\ne j,\sigma }{V}_{ij}{\hat{n}}_{i\sigma }{\hat{n}}_{j\sigma }.$$

$${\hat{H}}_{\text{AIM}}={\hat{H}}_{\text{AI}}+\sum _{ik\sigma }{\tilde{V}}_{ik}\left({\hat{a}}_{i\sigma }^{†}{\hat{c}}_{ik\sigma }^{}+h.c.\right)+\sum _{ik\sigma }{ϵ}_{ik}{\hat{n}}_{ik\sigma }.$$

$$U=\left(\begin{array}{cc}\tilde{U}& \tilde{U}\\ \tilde{U}& \tilde{U}\end{array}\right)\text{leavevmodenobreak leavevmodenobreak leavevmodenobreak andleavevmodenobreak leavevmodenobreak leavevmodenobreak }V=\left(\begin{array}{cc}0& \tilde{U}\\ \tilde{U}& 0\end{array}\right)$$

$${\Delta }_i(\omega )$$

$$p(\mathbf{A}|{\mathbf{G}}_{\text{d}},\Sigma ,\mathcal{I})$$

$$K$$

$$\mathbf{G}$$

$${\chi }^2$$

$$\displaystyle{\bf A}_{\text{ML}}=\sum_{i}\frac{\big{(}{\bf G}_{\text{d}}\cdot{\bf u}_{i}\big{)}}{E_{i}}\;{\bf v}_{i}\;.$$

$$\mathbf{A}$$

Full text

I Introduction

In quantum many-body physics, Green’s functions are often calculated in the imaginary-time domain, e.g., by numerical approaches, such as quantum Monte Carlo Hirsch and Fye (1986); Foulkes et al. (2001); Werner et al. (2006); Gull et al. (2011). The imaginary-time Green’s function $G(\tau)$ is related to the spectral function $A(\omega)$ by a Laplace transform, the so-called analytic continuation (AC). Obtaining $A(\omega)$ from $G(\tau)$ corresponds to inverting a Fredholm integral of the first kind and small changes in $G(\tau)$ correspond to large differences in $A(\omega)$. Therefore, the inversion of this problem is highly ill posed and very unstable against numerical noise, and even errors at the level of machine precision can lead to unphysical results in practice.

Many different methods to perform the AC have been proposed, e.g., series expansions such as the Padé method Ferris-Prabhu and Withers (1973), machine learning Arsenault et al. (2017), stochastic methods White et al. (1991); Sandvik (1998); Mishchenko et al. (2000); Beach (2004); Fuchs et al. (2010); Sandvik (2016); Ghanem (2017), and the maximum entropy method (MEM) Silver et al. (1990); Gubernatis et al. (1991); Gull and Skilling (1984); Bao et al. (2016); Kraberger et al. (2017). The latter is a consistent approach as it is based on Bayesian probability theory; however, a highly ignorant entropic prior is used, which merely accounts for positivity and additivity of the reconstructed signal.

Here, we propose a physically motivated prior that takes into account the knowledge of typical structures of a spectral density, which results in a parametric instead of a form-free reconstruction. The Bayesian parametric analytic continuation (BPAC) is based on Bayesian parameter estimation von der Linden. Volker Dose. Udo von Toussaint (2014); Gregory (2005); Jaynes et al. (2003); Sivia and Skilling (2006) to obtain a parametrized spectral function. To be precise, we use asymmetric Lorentzians as well as suitable tails to build up the spectral function. To validate parametrizations, we use Bayesian model comparison. The required evidence integrals are computed by nested sampling (NESA) Skilling (2004). With this approach, we can compare parametrizations, e.g., with a different number of asymmetric Lorentzians. Compared to the other methods, this allows to ask specific questions about the spectral function, e.g., about the reliability of peaks in the spectral function.

We demonstrate the capability of BPAC in terms of the spectral function of an Anderson impurity model (AIM), which exhibits a generalized Kondo scenario. We calculate the imaginary-time Green’s function of the AIM with continuous-time quantum Monte Carlo (CTQMC) Werner et al. (2006); Parcollet et al. (2015); Seth et al. (2016) and compare the spectral functions, obtained from BPAC with that of a MEM reconstruction. The MEM spectrum shows a peak close to the Abrikosov-Suhl resonance, but it is unclear whether this feature is physical or an artifact of the AC. BPAC can answer this question, showing that it is, in fact, an artifact. Additionally, we successfully compare the BPAC result to the solution obtained with the recently developed real-time fork tensor product state (FTPS) impurity solver Bauernfeind et al. (2017); Bauernfeind (2018), which directly computes the spectral density, without any AC.

In addition, we present a combination of MEM and BPAC and its application to an AIM arising from the dynamical mean-field theory (DMFT) Metzner and Vollhardt (1989); Georges and Kotliar (1992) treatment of SrVO3. The spectral function obtained with the real-time FTPS solver shows a three-peak structure in the upper Hubbard band, which is absent in the CTQMC + MEM spectral function. We investigate the question whether the absence of this structure is due to a failure of MEM or due to the ill-posed nature of the AC.

The applications of BPAC presented in this paper demonstrate that BPAC is a valuable addition to nonparametric reconstruction methods, such as MEM, e.g., to assess whether the data support specific features found in the MEM reconstruction.

The paper is organized as follows: We first introduce the AC problem in Sec. II. In Sec. III, we define and evaluate our parametrizations of the spectral function. Bayesian parameter estimation and model comparison are discussed in Sec. IV. Finally, in Secs. V and VI, we demonstrate the capability of BPAC, first, on an AIM that exhibits a generalized Kondo scenario and, second, on an AIM stemming from the ab initio treatment of SrVO3.

II Analytic continuation and maximum entropy methods

Dynamical correlation functions in imaginary time as obtained from CTQMC, obey the (anti-) periodicity relation $G(\tau+\beta)=\mp G(\tau)$. The upper sign (-) holds for fermions, and the lower sign (+) holds for bosons. Due to (anti-) periodicity, $G(\tau)$ is uniquely determined by its values in the interval $\tau\in[0,\beta)$, and its discrete Fourier representation is

[TABLE]

The sum is over the Matsubara frequencies $\omega_{n}=(2n+1)\pi/\beta$ for fermions and $\omega_{n}=2n\pi/\beta$ for bosons, where $n\in\mathbb{Z}$. The retarded Green’s function $\mathcal{G}(\omega+{\mathrm{i}}0^{+})$ and Matsubara Green’s function $\mathcal{G}({\mathrm{i}}\omega_{n})$ are related through the analyticity of $\mathcal{G}(z)$. The spectral function $A(\omega)=-\frac{1}{\pi}\Im\mathcal{G}(\omega+i0^{+})$ determines

[TABLE]

Merging Eqs. 1 and 2 produces the relation between the imaginary-time Green’s function $G(\tau)$ and the spectral function $A(\omega)$, e.g., for fermions,

[TABLE]

To handle the problem numerically, we discretize the functions $G(\tau)$ and $A(\omega)$, i.e., $({\bf G})_{n}=G_{n}=G(\tau_{n})$ and $({\bf A})_{m}=A_{m}=A(\omega_{m})$. Consequently, discretizing the kernel $K_{nm}=K(\tau_{n},\omega_{m})$ produces the matrix equation ${\bf G}=K{\bf A}$. Note that, as shown by Ref. Ghanem (2017), choosing the discretization grid already includes prior information and is equivalent to imposing a default model. Here, we restrict ourselves to linear discretization grids.

The determination of ${\bf G}$ from ${\bf A}$ is straight forward, but the inversion ${\bf A}=K^{-1}{\bf G}$ is a highly ill-posed inversion problem, which is impossible to tackle without taking the noise statistics and reliable prior knowledge consistently into account. In assuming a multivariate normal distribution with the covariance matrix $\Sigma$ of the QMC data vector ${\bf G}_{\text{d}}$, the maximum likelihood (ML) estimator ${\bf A}_{\text{ML}}$ is obtained by minimization of $\chi^{2}({\bf A})=(K{\bf A}-{\bf G}_{\text{d}})^{T}\Sigma^{-1}(K{\bf A}-{\bf G}_{\text{d}})$. Due to the ill posedness of the problem, ${\bf A}_{\text{ML}}$ is, in general, not a satisfying solution, e.g., negative, spiky, and unnormalized. Additional information, e.g., positivity, smoothness, etc., can be incorporated to regularize the problem. As shown by Skilling Skilling (1989) on a rigorous probabilistic footing, introducing an entropy term,

[TABLE]

and maximizing $-\frac{1}{2}\chi^{2}({\bf A})+\alpha S({\bf A})$, where $S({\bf A})$ is the discretized version of the entropy $S(A)$, regularizes the problem. In this so-called MEM, the standard model $D(\omega)$ determines the prior information about the spectral function, and the hyperparameter $\alpha$, roughly speaking, determines the mixing ratio between the ML solution and the standard model $D(\omega)$. A small $\alpha$ produces the ML solution, whereas for large $\alpha$, the spectral function approaches the standard model $D(\omega)$. The hyperparameter $\alpha$ can be adjusted in various ways, e.g., historic MEM Gull S.F. and Daniell G.J. (1978); Gull and Skilling (1984), classical MEM Gull (1989); Fischer et al. (1996), and Bryan MEM Bryan (1990).

III Parametrization of spectral functions

In the present paper, we propose BPAC. This approach circumvents the ill-posed problem to a large extend by representing the spectrum using only a few parameters.

We want to build up the spectral function as a sum of peaks where each peak is supposed to correspond to a real peak in the spectrum, e.g., the Abrikosov-Suhl resonance, the left and right Hubbard band, etc. Due to the natural line width of spectral lines, the obvious choice is to use Lorentzian functions. General peaks in spectral functions are not single Lorentzians and can show shoulders or plateaus, e.g., between a Hubbard band and the Abrikosov-Suhl resonance. Therefore, we introduce a sum of asymmetric Lorentzian functions for each peak and add additional tails to describe the decay of the spectrum at higher energies. It depends on the desired accuracy $|A(\omega)-A_{0}(\omega)|$ of the reconstructed spectrum $A(\omega)$ to the true spectral function $A_{0}(\omega)$, whether using one or more Lorentzians per peak is more appropriate. In our parametrization, the $n$th peak is located at position $\mu_{n}$ and is built up by $C_{n}$ asymmetric Lorentzians, each $i\in\{1,2,...,C_{n}\}$ having its own individual amplitude $a^{i}_{n}$, left width $\gamma^{i,\text{l}}_{n}$, and right width $\gamma^{i,\text{r}}_{n}$, resulting in $(3C_{n}+1)$ parameters ${\tilde{\alpha}}_{n}=\{\mu_{n},\{a^{i}_{n},\gamma^{i,\text{l}}_{n},\gamma^{i,\text{r}}_{n}\}\}$,

[TABLE]

$L(\omega|\mu,\gamma)$ denotes the normalized Lorentzian with center $\mu$ and width $\gamma$,

[TABLE]

We refer to the case of $C_{n}>1$ as split Lorentzian. The prefactors ${2\gamma^{i,\text{l}}_{n}}{(\gamma^{i,\text{l}}_{n}+\gamma^{i,\text{r}}_{n})^{-1}}$ and ${2\gamma^{i,\text{r}}_{n}}{(\gamma^{i,\text{l}}_{n}+\gamma^{i,\text{r}}_{n})^{-1}}$ in Eq. 5 ensure continuity and normalization of the asymmetric Lorentzians. The spectrum in the frequency interval $I_{\Omega}=[\Omega^{\text{l}},\Omega^{\text{r}}]$ is described by a superposition of the Lorentzians $f_{n}(\omega|\tilde{\alpha}_{n})$. For the tails of the spectrum outside the interval $I_{\Omega}$, a power-law decay is more appropriate. Then, the total spectral function becomes

[TABLE]

$N_{\text{p}}$ is the total number of peaks and $\bm{\alpha}$ includes the parameters of the Lorentzians $\bm{\tilde{\alpha}}$ and the six parameters of the tails, namely, $\{\mu^{\text{l}},\nu^{\text{l}},\Omega^{\text{l}},\mu^{\text{r}},\nu^{\text{r}},\Omega^{\text{r}}\}$. The parameters $a^{\text{l}}$ and $a^{\text{r}}$ are determined by forcing the spectral function to be continuous. Hence, we end up with a spectral function described by $N_{\alpha}=\left(\sum_{n}^{N_{\text{p}}}(3C_{n}+1)\right)+6$ parameters.

In a first test, we analyze how well parametrized spectral functions using Lorentzians of increasing complexity as defined in Eq. 7 can represent typical physical spectra. As a test case, we use the spin-down and spin-up spectra $A(\omega)$ of the AIM discussed in Sec. V.

Fig. LABEL:fig:fitA shows the (reference) spectra obtained by the FTPS impurity solver ($A_{0}(\omega)$, solid black lines) for the (a) spin-down and the (b) spin-up electrons and compares $A_{0}(\omega)$ to approximations of increasing complexity with the generalized Lorentzian ansatz (colored lines). For the (a) spin-down spectrum, we used a four-peak ($N_{\text{p}}=4$) spectral function whereas the (b) spin-up spectrum is approximated by a two-peak ($N_{\text{p}}=2$) spectral function. We determine the parameters in Eq. 7 by the least-squares approach. We find that using asymmetric Lorentzians instead of symmetric ones, allows to describe the peaks much better, whereas including tails leads to visible improvements in the high-energy regions, see the lower plots in Figs. LABEL:fig:fitA (a) and LABEL:fig:fitA (b). Using split Lorentzians further decreases the deviation to the reference spectrum. Note that the small oscillations in the lower plots of Fig. LABEL:fig:fitA (a) and (b) are artifacts caused by the Fourier transformation of the finite-time solution of the FTPS impurity solver.

To summarize, we find that the parametrization of Eq. 7 is highly flexible and allows to represent reliably the entire structure of the spectrum.

IV Parameter estimation and model comparison

IV.1 Bayesian data analysis

In this section, we discuss how to determine the parameters $\bm{\alpha}$ of the spectral function ${\bf A}(\bm{\alpha})$ and how to judge which parametrization is supported best by the data. From now on, we call a parametrization of the spectral function the model $M=M({\bf C})$, depending implicitly on the peak complexity ${\bf C}$, which also includes the number of peaks $N_{\text{p}}=\dim({\bf C})$. Note that $M$ does not define the values of the model parameters. Bayes’ theorem gives

[TABLE]

where ${\bf G}_{\text{d}}$ stands for the $N_{\text{d}}$ imaginary-time data points from CTQMC. We assume that they have a multivariate Gaussian stochastic error. The corresponding covariance matrix is denoted by $\Sigma$. The kernel $K$ is included in the conditional complex ${\cal I}$. The likelihood is therefore,

[TABLE]

Since correlations are negligible in the data sets used in the present paper we take $\Sigma_{ij}=\sigma_{i}^{2}\delta_{ij}$ in the following.

We use the prior probability to restrict the parameter space for two reasons: First, to only obtain physical results, e.g., by forcing the spectral function to be positive ($a_{n}^{i}>0$, $\forall n,i$); second, to build in and test additional knowledge about the spectral function, e.g., by forcing a peak to appear in a chosen energy interval to analyze whether the data supports this peak. For example, a question that can be answered by BPAC could be: Is there a side peak left to the Abrikosov-Suhl resonance? Bayesian model comparison allows to judge whether the model including the additional peak is more probable than the model without this peak. Apart from the above restrictions, we used a flat prior,

[TABLE]

In the present paper, we restrict ourselves to this prior, which is simple to implement and computationally inexpensive since sampling from uniform distributions is cheap. More advanced priors are possible although, e.g., including transformation invariance, smoothness, using testable information, such as the normalization of the spectral function, or even using the entropic prior, may help to improve the results.

Primarily, we are interested in the probability density for the spectral function $p({\bf A}|{\bf G}_{\text{d}},\Sigma,M,{\cal I})$, which we easily obtain from the posterior distribution $p(\bm{\alpha}|{\bf G}_{\text{d}},\Sigma,M,{\cal I})$ by using the marginalization rule,

[TABLE]

Additionally, we want to calculate the data evidence $p({\bf G}_{\text{d}}|\Sigma,M,{\cal I})$, which allows to assess the probability of different models relative to each other since the probability for model $M$ is proportional to the data evidence,

[TABLE]

In the so-called odds ratio, the ratio between the probabilities for models $M_{1}$ and $M_{2}$, the unknown probability $p({\bf G}_{\text{d}}|{\cal I})$ cancels out, and we get

[TABLE]

In the present paper, we set the prior odds to one since we do not want to favor any model.

IV.2 Evaluating posterior and data evidence

In this section, we want to give a very brief introduction to NESA, which is a method providing both the data evidence and the samples from the posterior. We provide the basic equations in this section but refer to Refs. Sivia and Skilling (2006); Skilling (2006); von der Linden. Volker Dose. Udo von Toussaint (2014) for the detailed derivation of NESA.

Skilling Skilling (2004) proposed to write the data evidence integral as the Lebesgue integral,

[TABLE]

where the integral over the prior mass,

[TABLE]

runs over the likelihood values $\lambda$. Equivalently, the data evidence can be written as

[TABLE]

where the integral is approximated by the Riemann sum and $\Delta X_{n}=X_{n}-X_{n+1}$. The likelihood $\mathcal{L}(X)$ is a monotonically decreasing function of the prior mass, and the computation is as complicated as the original evaluation of the data evidence. Skilling proposed a stochastic approach to sample $\mathcal{L}(X)$ based on order statistics providing the likelihood minima $\{\lambda^{*}_{n}\}$. The pseudo code is given in Algorithm IV.2.

{pseudocode}

[ruled]NESA algorithm{λ^_n},{α^n},n_max **input parameters: ** N_w, ϵλ

**initialize ** λ^*_0=0,  n=0,

draw $N_{\text{w}}$ configurations  {α_i}  at random from π(α | λ^*_0)  (Eq. 17)

take the smallest likelihood value $\lambda^{*}=\min\left\{\lambda_{i}={L}(\bm{\alpha}_{i})\right\}$ and its configuration $\bm{\alpha}^{*}$

**set ** λ^_n=1:= λ^ and α^_n=1:= α^

\WHILE|(

λ^_n+1-λ^_n)/λ^n+1|   >   ϵλ* \DO\BEGINn \GETSn+1

replace configuration α^* with a new configuration drawn from   π(α | λ^*_n)

determine the smallest likelihood $\lambda^{*}=\min\left\{\lambda_{i}={L}(\bm{\alpha}_{i})\right\}$ and its configuration $\bm{\alpha}^{*}$

**set ** λ^_n+1=λ^ and α^_n+1=α^ \END

**set ** n = n_max

\RETURN

{λ^_n},{α^_n},n_max

The nested sampling moves in configuration space ensure that even well-separated peaks of the likelihood function are sampled correctly. The crucial step for the NESA algorithm is to draw from the constrained prior probability,

[TABLE]

This probability density represents the normalized prior restricted to areas, where $L(\bm{\alpha})$ exceeds the $\lambda^{*}_{n}$ threshold. In the present paper, $\pi(\bm{\alpha})$ is constant within the prior constraints according to Eq. 10. Therefore, we need to draw samples from the uniform distribution constrained by both the likelihood and the prior. A simple way to draw a sample from Eq. 17 is to clone an existing configuration, which obviously fulfills all constraints, and perform an ordinary Markov chain Monte Carlo update obeying $\pi(\mathbf{x}\,\bm{|}\,\lambda^{*}_{n})$. We implemented local updates in the parameters and monitored autocorrelations, which can become considerable depending on the problem.

The prior masses can be derived using order statistics as shown in detail in Ref. von der Linden. Volker Dose. Udo von Toussaint (2014). We can write the $n$th prior mass as $X_{n}=\prod_{\nu=1}^{n}\theta_{\nu}$, where the shrinking factors $\theta_{\nu}$ are independent and identically distributed random variables and obey the first-order statistic of the uniform probability density function (PDF), the $\beta$ distribution $p(\theta)=\theta^{N_{\text{w}}-1}/N_{\text{w}}$. Knowing the distribution of $X_{n}$ and, therefore, of $\Delta X_{n}$, allows to calculate the expectation value and variance of the Riemann sum in Eq. 16 given the set of likelihood minima $\{\lambda^{*}_{n}\}$ obtained from Algorithm IV.2.

Also posterior samples can be generated from a single NESA run by reusing the samples $\{\bm{\alpha}^{*}_{n}\}$ obtained from Algorithm IV.2. Eq. 16 shows that the $n$th NESA step contributes with weight $\Delta X_{n}\lambda^{*}_{n}$ to the Riemann sum for calculating the data evidence. Therefore, samples from the posterior PDF can be provided by choosing $n$ with the corresponding $\bm{\alpha}^{*}_{n}$ according to its weight $p(n)\propto\Delta X_{n}\lambda^{*}_{n}$, e.g., by inversion sampling. With such posterior samples $\{\bm{\alpha}_{\nu}\}$, the expectation value of the spectral function can be obtained as $\langle{\bf A}\rangle=\frac{1}{N_{\nu}}\sum_{\nu=1}^{N_{\nu}}{\bf A}(\bm{\alpha}_{\nu})$.

There are different improvements of NESA going beyond the algorithm we presented within this section, which may increase the performance, e.g., updating more configurations at once Sivia and Skilling (2006); von der Linden. Volker Dose. Udo von Toussaint (2014), using a parallel version of NESA Henderson and Goggans (2014), extending the update method in the prior sampling Skilling (2012), or using the knowledge of the position of the minima obtained by optimization algorithms Martiniani et al. (2014).

V Application I: BPAC

In this section, we apply BPAC to an impurity problem closely related the one studied in Ref. Rumetshofer et al. (2019). In Sec. V.1, we define the AIM, which we solve, subsequently, using three different methods: FTPS, CTQMC + MEM, and CTQMC + BPAC. Technical details of these methods are given in Sec. V.2, whereas the comparison of the results is given in Sec. V.3.

V.1 The impurity problem

The Hamiltonian of the isolated multiorbital Anderson impurity with on-site energies $\epsilon_{i\sigma}$ and interaction parameters $U_{ij}$ for electrons of different spin and $V_{ij}$ for electrons of the same spin is

[TABLE]

Here, ${\mathop{\hat{n}}}_{i\sigma}={\mathop{\hat{a}}}_{i\sigma}^{\dagger}{\mathop{\hat{a}}}^{\phantom{\dagger}}_{i\sigma}$ is the particle number operator for orbital $i\in\{1,2\}$ and spin $\sigma\in\{\downarrow,\uparrow\}$ in the second quantization with creation (annihilation) operators ${\mathop{\hat{a}}}_{i\sigma}^{\dagger}$ (${\mathop{\hat{a}}}_{i\sigma}^{\phantom{\dagger}}$). In the AIM, the impurity is coupled to a bath of non-interacting fermions,

[TABLE]

${\mathop{\hat{c}}}_{ik\sigma}^{\dagger}$ (${\mathop{\hat{c}}}_{ik\sigma}^{\phantom{\dagger}}$) are the creation (annihilation) operators of the $k$th bath state of orbital $i$ with spin $\sigma$. For the on-site energy of the impurity, we use $\epsilon_{i\downarrow}=\epsilon-J$ and $\epsilon_{i\uparrow}=\epsilon$ with \epsilon=-0.25\leavevmode\nobreak\ eV and J=50\leavevmode\nobreak\meV, and

[TABLE]

with \tilde{U}=0.5\leavevmode\nobreak\eV. The bath parameters $\tilde{V}_{ik}$ and $\epsilon_{ik}$ are obtained from a flat bath hybridization function,

[TABLE]

defined by $-2\Im\left(\Delta_{i}(\omega)\right)=\Gamma\;\Theta(-1<\omega<1)$ with \Gamma=50\leavevmode\nobreak\meV. This set of parameters exhibits a generalized Kondo scenario with symmetry between SU(2) Haldane (1978) and SU(4) Filippone et al. (2014) with the corresponding Kondo temperatures of T_{\text{SU}(2)}=0.36\leavevmode\nobreak\K and T_{\text{SU}(4)}=20\leavevmode\nobreak\K. Due to the difference in the on-site energies of the impurity orbitals $J$, this AIM exhibits side peaks close to the Abrikosov-Suhl resonance. We present the spin-down and spin-up spectral functions of this AIM in Figs. 2 and 4.

V.2 Technical details of the methods

We solve the AIM in the imaginary-time domain using the CTQMC solver in hybridization expansion as implemented in the TRIQS library Werner et al. (2006); Parcollet et al. (2015); Seth et al. (2016). We performed 15 CTQMC runs at $\beta=400\leavevmode\nobreak\ \text{(eV)}^{-1}$ (T=29\leavevmode\nobreak\K), each on 20 node points and with $10^{6}$ measurements. The difference in the impurity on-site energies $J$ lifts the spin degeneracy but keeps the orbital degeneracy intact and, therefore, the 15 CTQMC runs give 30 independent samples. Based on this sample, we estimate reliable variances $\sigma_{i}$ for the QMC data without having to bother about possible autocorrelations.

In the following, we do not distinguish orbitals anymore and just discuss the spectral functions depending on the spin. The AC in the present paper is performed with MEM and BPAC. In both cases, we use $N_{\text{d}}=401$ data points on an equally spaced $\tau$ grid for $\tau\in[0,\beta]$ and the same amount of $\omega$ points equally spaced on the interval $\omega\in[-1,1]$. We applied the MEM of Ref. Jarrell and Gubernatis (1996) with an alternative evidence approximation von der Linden et al. (1999) and the preblur formalism Skilling (1991). BPAC is applied as explained in Sec. IV using a $N_{\text{w}}=1000$ walker and $\epsilon_{\lambda}=10^{-5}$.

Additionally, we compare the results with those obtained by the real-time FTPS impurity solver, which does not need any AC, since it calculates the Green’s function already on the real axis. In contrast to the CTQMC solver, the FTPS solver is a zero-temperature method, which has to be considered when comparing the results.

V.3 Comparison of the results

First, we employ the FTPS solver for the spin-down part of the AIM and show the corresponding spectrum in Fig. 2 (black line). As expected from the definition of the impurity model in Sec. V.1, the spectral function shows Hubbard satellites at approximately -0.3\leavevmode\nobreak\eV and 0.2\leavevmode\nobreak\eV and the Abrikosov-Suhl resonance at the chemical potential (0\leavevmode\nobreak\eV). Additionally, there is a peak at approximately -50\leavevmode\nobreak\meV contributed by the exchange coupling parameter $J$.

The CTQMC + MEM result (red dash-dotted line) shows a spurious peak at approximately 0.04\leavevmode\nobreak\eV, which does not appear in the FTPS solution. To find out whether this peak is supported by the CTQMC data or whether it is an artifact of the AC by MEM, we employ BPAC (blue line). First, we use a four-peak model $M_{4}=M({\bf C}=(1,1,1,1))$ where each peak consists of a single asymmetric Lorentzian, i.e., with peak complexity $C_{n}=1$. The NESA logarithmic data evidence yields $\ln(p({\bf G}_{\text{d}}|M_{4},\Sigma,{\cal I}))=2828.2\pm 0.3$. Tab. 1 in the App. IX.1 shows the prior ranges used and the parameters estimated. The four-peak BPAC solution (Fig. 2, blue line) does not show the additional peak slightly above $\omega=0$ by construction. The evaluation of the five-peak model $M_{5}=M({\bf C}=(1,1,1,1,1))$ where we introduce an additional peak at $\mu\in(0.03,0.07)$ produces a logarithmic data evidence $\ln(p({\bf G}_{\text{d}}|M_{5},\Sigma,{\cal I}))=2825.0\pm 0.3$. This yields a logarithmic Bayes factor of $3.2\pm 0.6$, which corresponds to a probability of $(93\%-98\%)$ that the four-peak model is preferred over the five-peak model. This is in agreement with the FTPS solution and demonstrates that the fifth peak at \sim 0.04\leavevmode\nobreak\eV is not supported by the CTQMC data and, therefore, an artifact of the MEM solution. In general, the spectral function obtained by BPAC depends on the choice of the model $M$. In the spirit of Bayesian probability theory, we can average over different models weighted by their corresponding model probability. Therefore, we actually should compute

[TABLE]

If one model is highly preferable, as $M_{4}$ in the present case, then $p({\bf A}|{\bf G}_{\text{d}},\Sigma,{\cal I})\approx p({\bf A}|{\bf G}_{\text{d}},\Sigma,M_{4},{\cal I})$. Hence, we plotted $p({\bf A}|{\bf G}_{\text{d}},\Sigma,M_{4},{\cal I})$ in Fig. 2 (blue line).

Starting point for the determination of the spectral function $A(\omega)$ are CTQMC data on the imaginary-time Green’s function $G(\tau)$, which we denote by $G_{\text{d}}(\tau)$. Inserting the reconstructed spectral function $A(\omega)$ into Eq. 3 yields the reconstructed $G(\tau)$, which allows to asses the misfit in data space. Likewise, we can apply Eq. 3 to the FTPS spectral function to obtain the corresponding $G(\tau)$. The reconstructed Green’s function for imaginary times $G(\tau)$ is compared with the CTQMC data $G_{\text{d}}(\tau)$ in Fig. 3 for the MEM, BPAC, and FTPS.

Even though the spectral functions of MEM and BPAC differ slightly, $G(\tau)$ of both solutions lies within the error of the CTQMC data. In the lower panel, the difference between $G(\tau)$ and $G_{\text{d}}(\tau)$ is shown on an enlarged scale, which reveals a systematic deviation between FTPS and CTQMC data. The reason is that the FTPS solver calculates the spectral function at T=0\leavevmode\nobreak\K, whereas $\beta=400\leavevmode\nobreak\ (\text{eV})^{-1}$ (T=29\leavevmode\nobreak\K) is used in the CTQMC simulation.

The spectrum of the spin-up part of the AIM obtained with the FTPS solver is presented in Fig. 4 (black line) and shows a two-peak structure as does the MEM (red dash-dotted line).

Hence, for parametrizing the spin-up spectral function in BPAC (blue line), we use a two-peak model $M_{2}=M({\bf C}=(1,1))$ with the prior ranges given in Tab. 2 in App. IX.1. Since we are not interested in specific questions about spurious peaks, we are satisfied with model $M_{2}$.

VI Application II: MEM + BPAC

In this section, we propose a combination of MEM and BPAC (MEM + BPAC) and apply the method to the impurity problem studied in Ref. Bauernfeind et al. (2017) for SrVO3. We present the details of the AIM in Sec. VI.1, give the technical details of MEM + BPAC in Sec. VI.2, and discuss the results in Sec. VI.3.

VI.1 The impurity problem

The multiorbital AIM discussed in Ref. Bauernfeind et al. (2017) arises from the ab initio treatment of SrVO3, which has become a test-bed material in DMFT. The solution of the AIM obtained with the FTPS solver shows a three-peak structure in the upper Hubbard band between 1.75\leavevmode\nobreak\eV and 4.25\leavevmode\nobreak\eV, see Fig. 5 (black line).

Ref. Bauernfeind et al. (2017) showed that CTQMC + MEM is not able to resolve these high-energy excitations. The question we want to address in this section is as follows: Is the absence of the three-peak structure a failure of MEM, or is it - due to the ill-posed inversion problem - generally impossible to recover certain high-energy details of the spectrum? To answer this question, we applied MEM + BPAC as explained in the following section.

VI.2 Technical details and MEM + BPAC

In order to obtain an answer to this question, we start from the FTPS real-frequency data, transform them to the imaginary-time axis, and add noise on the order of the CTQMC error ($\sigma=10^{-5}$). By this procedure, we ensure that we know precisely the error statistics of the data, and we know that the correct result has to have the three-peak structure.

We use an inverse temperature $\beta=200\leavevmode\nobreak\ \text{(eV)}^{-1}$ (T=58\leavevmode\nobreak\K), $N_{\text{d}}=501$ data points on an equally spaced $\tau$ grid in the interval $\tau\in[0,\beta]$ and the same amount of $\omega$ points equally spaced for $\omega\in[-4,6]$.

Instead of using BPAC as explained in the previous sections, here, we apply a combination of MEM and BPAC. MEM + BPAC takes the MEM solution for a given subinterval of the energy axis and applies BPAC only for the remaining interval. In that way, the number of parameters is small, which enables faster sampling in the calculation of the evidence integral with NESA. We take $\Omega\in(1,1.75)$ as an additional parameter and use the MEM solution for $\omega<\Omega$ and BPAC for $\omega\geq\Omega$. We use the prior ranges $1.75<\mu_{n}<4.25$ and $0<a_{n},\gamma_{n}^{\text{l}},\gamma_{n}^{\text{r}}<1$, and $C_{n}=1$ for each peak $n$. Furthermore, the remaining parameters describing the right tail are constrained by $4.25<\Omega^{r}<6$, $-2<\mu^{r}<4.25$, and $1<\nu^{r}<10$. In NESA, we use $N_{\text{w}}=2000$ walkers and $\epsilon_{\lambda}=10^{-5}$.

VI.3 Comparison of the results

We applied the MEM of Ref. Jarrell and Gubernatis (1996) with an alternative evidence approximation von der Linden et al. (1999) and the preblur formalism Skilling (1991) and were able to qualitatively reproduce the CTQMC + MEM solution in Fig. 5 of Ref. Bauernfeind et al. (2017), see Fig. 5 (red dash-dotted line, mostly covered by the blue line). The MEM spectral function does not show the three-peak structure in the upper Hubbard band.

We applied MEM + BPAC using one-, two-, and three-peak models $M_{1}$, $M_{2}$, and $M_{3}$. The obtained logarithmic data evidences $\ln(p({\bf G}_{\text{d}}|M_{i},\Sigma,{\cal I}))$ for $i\in\{1,2,3\}$ are $\{5036.7\pm 0.2,5036.0\pm 0.2,5035.3\pm 0.2\}$ and correspond to probabilities of 57 % for $M_{1}$, 29 % for $M_{2}$, and 14 % for $M_{3}$. It is interesting to note that the correct three-peak model actually has the lowest probability. Still, let us take a look at the result of $M_{3}$ shown in Fig. 5. Surprisingly, the three-peak model looks very similar to the MEM result, i.e., it is not even able to resolve the three-peak structure. Instead, it just shows one large peak in the energy region of the upper Hubbard band.

To elucidate this behavior, we consider the singular value representation of the kernel Pav (2018),

[TABLE]

Given the vector ${\bf A}$ of the discretized spectral function and the corresponding vector ${\bf G}$ of the Green’s function for discrete imaginary times as defined in Sec. II, we get

[TABLE]

The misfit defined in Eq. 9 can then be expressed in a very suggestive way. For simplicity, we assume constant noise $\sigma_{l}=\sigma$, $\forall l$. Then, the misfit to the data vector ${\bf G}_{\text{d}}$ is

[TABLE]

In Fig. 6, the singular values $E_{i}$ of the kernel matrix are plotted with decreasing magnitude on a logarithmic scale. We find that the singular values decrease exponentially and that, above $i=70$, the singular values are smaller than machine precision. We select the three singular modes $i=\{27,48,66\}$ corresponding to singular values of approximately $\{10^{-5},10^{-10},10^{-15}\}$, respectively. The corresponding modes ${\bf v}_{i}$, which are depicted in the upper inset, show an increasing number of nodes with increasing index $i$. Minimization of Eq. 24 with respect to the spectral function ${\bf A}$ yields the maximum likelihood solution,

[TABLE]

Since the singular values decay exponentially, their inverse increases exponentially and small noise in the coefficients \big{(}{\bf G}_{\text{d}}\cdot{\bf u}_{i}\big{)} becomes amplified. Fig. 7 demonstrates nicely that \big{(}{\bf G}_{\text{d}}\cdot{\bf u}_{i}\big{)} becomes dominated by noise above the mode index $i$ where the singular value $E_{i}$ reaches the order of the noise $\sigma$.

In order to see which part of the spectrum can, therefore, be reconstructed, we expand ${\bf A}$ into the mode vectors ${\bf v}_{i}$,

[TABLE]

In the lower inset of Fig. 6, we present the projected FTPS spectrum for $N=27,48$, and $66$. We observe that the three-peak structure is not resolvable at all with $N=27$ and the resulting spectrum (green dashed line in the lower inset) looks similar to the MEM and MEM + BPAC solutions in Fig. 5. $N=48$ allows to resolve two of the three peaks (blue dash-dotted line in the lower inset), whereas only $N=66$ resolves the full three-peak structure (red line in the lower inset). This demonstrates that, for a CTQMC error of $\sigma\geq 10^{-10}$, the AC kernel definitely does not allow to resolve the three-peak-structure. CTQMC errors of this magnitude, however, would imply enormous data-acquiring times, and even then only two of the three peaks would be visible.

VII Conclusions

We proposed a Bayesian parametric approach for the analytic continuation of noisy imaginary-time Green’s function data as, e.g. obtained by CTQMC. The commonly used Bayesian form-free reconstruction of QMC data is the MEM which is based on the entropic prior that uses a minimum amount of prior information, merely positivity and additivity. Due to the nature of the form-free reconstruction, there are typically as many unknown parameters as noisy data points. This, in combination with the ill-conditioned kernel, can lead to spurious features in the reconstructed spectrum. In many applications, however, we have additional prior knowledge, e.g., we know that there will be a small number of peaklike structures of a specific shape. This prior knowledge can be encoded by representing the spectrum in terms of suitably parametrized basis functions, encoding the spectrum with only a few parameters, much less than the number of data points. Our approach, which we denote BPAC, employs Bayesian parameter estimation to obtain the parametrized spectral function. In the present paper, we used asymmetric Lorentzians and additional tails. Of course, in other applications a different basis might be favorable.

Moreover, we employed Bayesian model comparison to validate different numbers of Lorentzian peaks. For Bayesian model comparison, the evaluation of high-dimensional evidence integrals is necessary. To this end, we employ NESA, which is particularly efficient for such integration problems.

We demonstrated the capability of BPAC in terms of the CTQMC data for an AIM that has a generalized Kondo scenario and compared the BPAC spectra to the MEM result as well as the spectral function obtained with the real-time FTPS impurity solver. It was shown that BPAC is able to tell true peaks from artifacts which are present in the MEM’s solution close to Abrikosov-Suhl resonance.

In a second application, we studied the AIM arising from the ab initio treatment of SrVO3. The spectral function obtained with the real-time FTPS solver shows a three-peak structure in the upper Hubbard band, which is not present in the MEM reconstruction of the CTQMC data. Although the MEM cannot resolve the three-peak structure, the rest of the spectrum is captured well. To start with a data set that definitely contains the three-peak structure, we generated imaginary-time data from the real-frequency FTPS spectrum. Adding noise to simulate the CTQMC error, we studied the MEM and BPAC reconstructions of this data set. To keep the number of parameters and, therefore, the numerical effort small, we employed BPAC focused on the structure in the upper Hubbard band. Therefore, we only described the upper part of the spectrum by Lorentzians whereas keeping the MEM reconstruction for the rest of the spectrum. Bayesian model comparison then allows to infer which details of the upper Hubbard band can reliably be inferred from the data. Considering the singular value decomposition of the kernel allows us to find rigorous arguments how numerical noise is propagated by the kernel.

Remarkably, we found that the information of the three-peak structure present in the real-frequency spectrum is attenuated by ten orders of magnitude during the transformation to imaginary-time space. It is, therefore, buried in the noise and impossible to be retrieved from the QMC data. This means that independent of the model chosen, we obtained a single large peak resembling the MEM solution. Although BPAC was not able to reconstruct the true shape of the upper Hubbard band, its advantage is that it reliably detects how many details of the spectrum are actually above the noise threshold in the data.

Therefore, we conclude that BPAC is a valuable addition to nonparametric reconstruction methods, such as MEM. The reconstruction can be performed either only with BPAC or with a MEM reconstruction that can be used first, and BPAC is employed to assess whether the data support specific features found in the MEM spectral function.

VIII Acknowledgments

The authors like to thank G. Kraberger and R. Triebl for discussions concerning the maximum entropy method. This work was partially supported by the Austrian Science Fund (FWF) through Grant No. SFB-ViCoM F41P04 and through the START Program No. Y746 as well as by NAWI Graz.

IX Appendix

IX.1 Prior ranges and estimated parameters

Tab. 1 shows the prior ranges and the estimated parameters for the four-peak model $M_{4}$ of the spin-down spectral function of the AIM discussed in Sec. V.

Tab. 2 shows the prior ranges and the estimated parameters for the two-peak model $M_{2}$ of the spin-up spectral function of the AIM discussed in Sec. V.