Stochastic pole expansion method

The single-particle imaginary time Green’s function $G(\tau)$ reads:

$$G(\tau)=\langle\mathcal{T}_{\tau}c(\tau)c^{\dagger}(0)\rangle,$$

(2)

where $\mathcal{T}_{\tau}$ is the imaginary time ordering operator, $c^{\dagger}$ and $c$ are creation and annihilation operators in the Heisenberg representation, respectively. For fermions, $G(\tau)$ must fulfil the anti-periodicity condition, i.e., $G(\tau)=-G(\tau+\beta)$. While for bosons, $G(\tau)$ must be $\beta$-periodic, i.e., $G(\tau)=G(\tau+\beta)$. Here $\beta$ denotes the inverse temperature of the system ($\beta=1/T$). The Matsubara Green’s function $G(i\omega_{n})$ can be derived from $G(\tau)$ via Fourier transformation,

$$G(i\omega_{n})=\int^{\beta}_{0}d\tau~{}e^{-i\omega_{n}\tau}G(\tau),$$

(3)

and vice versa,

$$G(\tau)=\frac{1}{\beta}\sum_{n}e^{i\omega_{n}\tau}G(i\omega_{n}).$$

(4)

Note that $\omega_{n}=(2n+1)\pi/\beta$ and $2n\pi/\beta$ for fermions and bosons, respectively, where $n\in\mathbb{Z}$.

Supposed that the spectral density of the single-particle Green’s function is $A(\omega)$, then we have:

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

(5)

with the positive (negative) sign for fermionic (bosonic) operators. Similarly,

$$G(i\omega_{n})=\int^{+\infty}_{-\infty}d\omega\frac{1}{i\omega_{n}-\omega}A(\omega).$$

(6)

These equations denote the spectral representation of single-particle Green’s function. We notice that the SPX method, as well as the other analytic continuation methods that are classified as ASM, are closely related to the spectral representation. Next, we would like to make further discussions about this representation for the fermionic and bosonic correlators.

Fermionic correlators. The spectral density $A(\omega)$ is defined on $(-\infty,\infty)$. It is positive definite, i.e., $A(\omega)\geq 0$. Eq. (5) and Eq. (6) can be reformulated as:

$$G(\tau)=\int^{+\infty}_{-\infty}d\omega~{}K(\tau,\omega)A(\omega),$$

(7)

and

$$G(i\omega_{n})=\int^{+\infty}_{-\infty}d\omega~{}K(\omega_{n},\omega)A(\omega),$$

(8)

respectively. The kernel functions $K(\tau,\omega)$ and $K(\omega_{n},\omega)$ are defined as follows:

$$K(\tau,\omega)=\frac{e^{-\tau\omega}}{1+e^{-\beta\omega}},$$

(9)

and

$$K(\omega_{n},\omega)=\frac{1}{i\omega_{n}-\omega}.$$

(10)

Bosonic correlators. The spectral density $A(\omega)$ obeys the following constraint: $\text{sign}(\omega)A(\omega)\geq 0$. Thus, it is more convenient to define a new function $\tilde{A}(\omega)$ where $\tilde{A}(\omega)=A(\omega)/\omega$. Clearly, $\tilde{A}(\omega)$ is always positive definite. As a result Eq. (5) and Eq. (6) can be rewritten as:

$$G(\tau)=\int^{+\infty}_{-\infty}d\omega~{}K(\tau,\omega)\tilde{A}(\omega),$$

(11)

and

$$G(i\omega_{n})=\int^{+\infty}_{-\infty}d\omega~{}K(\omega_{n},\omega)\tilde{A}(\omega),$$

(12)

respectively. Now the bosonic kernel $K(\tau,\omega)$ becomes:

$$K(\tau,\omega)=\frac{\omega e^{-\tau\omega}}{1-e^{-\beta\omega}}.$$

(13)

Especially, $K(\tau,0)=1/\beta$. As for $K(\omega_{n},\omega)$, its expression is:

$$K(\omega_{n},\omega)=\frac{\omega}{i\omega_{n}-\omega}.$$

(14)

Especially, $K(0,0)=-1$. Besides the bosonic Green’s function, typical correlator of this kind includes the transverse spin susceptibility $\chi_{+-}(\tau)=\langle S_{+}(\tau)S_{-}(0)\rangle$, where $S_{+}=S_{x}+iS_{y}$ and $S_{-}=S_{x}-iS_{y}$.

Bosonic correlators of Hermitian operators. There is a special case of the previous observable kind with $c=c^{\dagger}$. Here, $A(\omega)$ becomes an odd function, and equivalently, $\tilde{A}(\omega)$ is an even function [i.e., $\tilde{A}(\omega)=\tilde{A}(-\omega)$]. Therefore the limits of integrations in Eq. (5) and Eq. (6) are reduced from $(-\infty,\infty)$ to $(0,\infty)$. So the two equations can be transformed into:

$$G(\tau)=\int^{+\infty}_{0}d\omega~{}K(\tau,\omega)\tilde{A}(\omega),$$

(15)

and

$$G(i\omega_{n})=\int^{+\infty}_{0}d\omega~{}K(\omega_{n},\omega)\tilde{A}(\omega),$$

(16)

respectively. The corresponding $K(\tau,\omega)$ reads:

$$K(\tau,\omega)=\frac{\omega\left[e^{-\tau\omega}+e^{-(\beta-\tau)\omega}\right]}{1-e^{-\beta\omega}}.$$

(17)

Especially, $K(\tau,0)=2/\beta$. And $K(\omega_{n},\omega)$ becomes:

$$K(\omega_{n},\omega)=\frac{-2\omega^{2}}{\omega_{n}^{2}+\omega^{2}}.$$

(18)

Especially, $K(0,0)=-2$. Perhaps the longitudinal spin susceptibility $\chi_{zz}(\tau)=\langle S_{z}(\tau)S_{z}(0)\rangle$ and the charge susceptibility $\chi_{ch}(\tau)=\langle N(\tau)N(0)\rangle$ are the most widely used observables of this kind.

It is well known that the finite temperature many-body Green’s functions can be expressed within the Lehmann representation Negele and Orland (1998):

$$G_{ab}(z)=\frac{1}{Z}\sum_{m,n}\frac{\langle n|d_{a}|m\rangle\langle m|d_{b}^{\dagger}|n\rangle}{z+E_{n}-E_{m}}\left(e^{-\beta E_{n}}\pm e^{-\beta E_{m}}\right),$$

(19)

where $a$ and $b$ are the band indices, $d$ ($d^{\dagger}$) denote the annihilation (creation) operators, $|n\rangle$ and $|m\rangle$ are eigenstates of the Hamiltonian $\hat{H}$, and $E_{n}$ and $E_{m}$ are the corresponding eigenvalues, $Z$ is the partition function ($Z=\sum_{n}e^{-\beta E_{n}}$). The positive sign corresponds to fermions, while the negative sign corresponds to bosons. The domain of this function is on the complex plane, but the real axis is excluded ($z\in\{0\}\bigcup\mathbb{C}~{}\backslash~{}\mathbb{R}$). If $z=i\omega_{n}\in i\mathbb{R}$, $G_{ab}(i\omega_{n})$ is the Matsubara Green’s function. If $z=\omega+i0^{+}$, $G_{ab}(\omega+i0^{+})=G_{ab}^{R}(\omega)$ is called the retarded Green’s function.

At first, we focus on the diagonal cases ($a=b$). For the sake of simplicity, the band indices are ignored in the following discussions. Let $A_{mn}=\langle n|d|m\rangle\langle m|d^{\dagger}|n\rangle\left(e^{-\beta E_{n}}+e^{-\beta E_{m}}\right)/Z$ and $P_{mn}=E_{m}-E_{n}$, then $G(z)=\sum_{m,n}A_{mn}/(z-P_{mn})$ Huang et al. (2023). Clearly, only those nonzero elements of $A_{mn}$ contribute to the Green’s function. If the indices $m$ and $n$ are further compressed into $\gamma$ (i.e, $\gamma=\{m,n\}$), then Eq. (19) is simplified to:

$$G(z)=\sum^{N_{p}}_{\gamma=1}\frac{A_{\gamma}}{z-P_{\gamma}}.$$

(20)

Here, $A_{\gamma}$ and $P_{\gamma}$ mean the amplitude and location of the $\gamma$-th pole, respectively. $N_{p}$ means the number of poles, which is equal to the total number of nonzero $A_{mn}$. Such an analytic expression of Green’s function is called the pole expansion. It is valid for both fermionic and bosonic correlators.

Fermionic correlators. For fermionic systems, the pole representation for Matsubara Green’s function can be recast as:

$$G(i\omega_{n})=\sum^{N_{p}}_{\gamma=1}\Xi(\omega_{n},P_{\gamma})A_{\gamma}.$$

(21)

Here, $\Xi$ is called the kernel matrix. It is evaluated by:

$$\Xi(\omega_{n},\omega)=\frac{1}{i\omega_{n}-\omega}.$$

(22)

Note that $A_{\gamma}$ and $P_{\gamma}$ should satisfy the following constraints:

$$\forall\gamma,~{}0\leq A_{\gamma}\leq 1,~{}\sum_{\gamma}A_{\gamma}=1,~{}P_{\gamma}\in\mathbb{R}.$$

(23)

Bosonic correlators. For bosonic systems, the pole representation for Matsubara Green’s function can be defined as follows:

$$G(i\omega_{n})=\sum^{N_{p}}_{\gamma=1}\Xi(\omega_{n},P_{\gamma})\tilde{A}_{\gamma}.$$

(24)

Here, $\Xi$ is evaluated by:

$$\Xi(\omega_{n},\omega)=\frac{G_{0}\omega}{i\omega_{n}-\omega},$$

(25)

where $G_{0}=-G(i\omega_{n}=0)$, which should be a positive real number. Be careful, $\Xi(0,\omega)=-G_{0}$. $\tilde{A}_{\gamma}$ is the renormalized amplitude of the $\gamma$-th pole:

$$\tilde{A}_{\gamma}=\frac{A_{\gamma}}{G_{0}P_{\gamma}}.$$

(26)

Note that $\tilde{A}_{\gamma}$ and $P_{\gamma}$ should satisfy the following constraints:

$$\forall\gamma,~{}0\leq\tilde{A}_{\gamma}\leq 1,~{}\sum_{\gamma}\tilde{A}_{\gamma}=1,~{}P_{\gamma}\in\mathbb{R}.$$

(27)

Bosonic correlators of Hermitian operators. Its pole representation can be defined as follows ($\forall\gamma$, $A_{\gamma}>0$ and $P_{\gamma}>0$):

	$\displaystyle G(i\omega_{n})$	$\displaystyle=$	$\displaystyle\sum^{N_{p}}_{\gamma=1}\left(\frac{A_{\gamma}}{i\omega_{n}-P_{\gamma}}-\frac{A_{\gamma}}{i\omega_{n}+P_{\gamma}}\right)$		(28)
		$\displaystyle=$	$\displaystyle\sum^{N_{p}}_{\gamma=1}\Xi(\omega_{n},P_{\gamma})\tilde{A}_{\gamma}.$

Thus, the kernel matrix $\Xi$ reads:

$$\Xi(\omega_{n},\omega)=\frac{-G_{0}\omega^{2}}{\omega^{2}_{n}+\omega^{2}}.$$

(29)

Especially, $\Xi(0,0)=-G_{0}$. The renormalized weight $\tilde{A}_{\gamma}$ reads:

$$\tilde{A}_{\gamma}=\frac{2A_{\gamma}}{G_{0}P_{\gamma}}.$$

(30)

The constraints for $\tilde{A}_{\gamma}$ and $P_{\gamma}$ are the same with Eq. (27).

As for the off-diagonal cases ($a\neq b$), it is lightly to prove that $\sum_{\gamma}A_{\gamma}=0$. It implies that there exist poles with negative weights. Hence we can split the poles into two groups according to the signs of their amplitudes. The Matsubara Green’s function can be expressed as follows:

	$\displaystyle G(i\omega_{n})$	$\displaystyle=$	$\displaystyle\sum^{N^{+}_{p}}_{\gamma=1}\frac{A^{+}_{\gamma}}{i\omega_{n}-P^{+}_{\gamma}}-\sum^{N^{-}_{p}}_{\gamma=1}\frac{A^{-}_{\gamma}}{i\omega_{n}-P^{-}_{\gamma}}$		(31)
		$\displaystyle=$	$\displaystyle\sum^{N^{+}_{p}}_{\gamma=1}\Xi(\omega_{n},P^{+}_{\gamma})A^{+}_{\gamma}-\sum^{N^{-}_{p}}_{\gamma=1}\Xi(\omega_{n},P^{-}_{\gamma})A^{-}_{\gamma}.$

Here, $\Xi(\omega_{n},\omega)$ is already defined in Eq. (22). The $A^{\pm}_{\gamma}$ and $P^{\pm}_{\gamma}$ are restricted by Eq. (23). In addition,

$$N_{p}=N^{+}_{p}+N^{-}_{p},$$

(32)

and

$$\sum^{N^{+}_{p}}_{\gamma=1}A^{+}_{\gamma}-\sum^{N^{-}_{p}}_{\gamma=1}A^{-}_{\gamma}=0.$$

(33)

Supposed that the input Matsubara Green’s function is $\mathcal{G}(i\omega_{n})$, where $n=$ 1, 2, $\cdots$, $N_{\omega}$, the objective of analytic continuation is to fit the (possibly noisy and incomplete) Matsubara data into Eq. (20) under some constraints. In mathematical language, we should solve the following multivariate optimization problem:

$$\mathop{\arg\min}\limits_{\left\{A_{\gamma},P_{\gamma}\right\}^{N_{p}}_{\gamma=1}}\chi^{2}\left[\left\{A_{\gamma},P_{\gamma}\right\}^{N_{p}}_{\gamma=1}\right].$$

(34)

Here, $\chi^{2}\left[\left\{A_{\gamma},P_{\gamma}\right\}^{N_{p}}_{\gamma=1}\right]$ is the so-called goodness-of-fit function or loss function. Its definition is as follows:

$$\chi^{2}\left[\left\{A_{\gamma},P_{\gamma}\right\}^{N_{p}}_{\gamma=1}\right]=\frac{1}{N_{\omega}}\sum^{N_{\omega}}_{n=1}\left|\left|\mathcal{G}(i\omega_{n})-\sum^{N_{p}}_{\gamma=1}\frac{A_{\gamma}}{i\omega_{n}-P_{\gamma}}\right|\right|^{2}_{F},$$

(35)

where $||\cdot||_{F}$ denotes the Frobenius norm. The minimization of Eq. (34) is highly non-convex. Traditional gradient-based optimization methods, such as the non-negative least squares method, conjugate gradient method, Newton and quasi-Newton methods, are frequently trapped in local minima Nocedal and Wright (2006). Their optimized results strongly depend on the initial guess. The semi-definite relaxation (SDR) fitting method Mejuto-Zaera et al. (2020); Huang et al. (2023), adaptive Antoulas-Anderson (AAA) algorithm Nakatsukasa and Trefethen (2020); Nakatsukasa et al. (2018), and conformal mapping plus Prony’s method Ying (2022a, b), which have been employed to search the locations of poles in previous works Huang et al. (2023), are also tested. But these methods usually fail when $N_{p}$ is huge [$N_{p}\sim O(10^{3})$] or the Matsubara data are noisy.

In order to overcome the above obstacles, we employ the simulated annealing method Kirkpatrick et al. (1983) to locate the global minimum of $\chi^{2}$. The core idea is as follows: First of all, a set of $\{A_{\gamma},P_{\gamma}\}$ parameters are generated randomly. These parameters form a configuration space $\mathcal{C}=\{A_{\gamma},P_{\gamma}\}$. Second, this configuration space is sampled by using the Metropolis Monte Carlo algorithm. In the present SPX method, four Monte Carlo updates are supported (see Figure 1). They include: (i) Select one pole randomly and shift its location. (ii) Select two poles randomly and shift their locations. (iii) Select two poles randomly and change their amplitudes. The sum-rules, such as Eq. (23) and Eq. (27), should be respected in this update. (iv) Select two poles randomly and exchange their amplitudes. Assumed that the current Monte Carlo configuration is $\mathcal{C}=\{A_{\gamma},P_{\gamma}\}$, the new one is $\mathcal{C}^{\prime}=\{A^{\prime}_{\gamma},P^{\prime}_{\gamma}\}$, and $\Delta\chi^{2}=\chi^{2}(\mathcal{C}^{\prime})-\chi^{2}(\mathcal{C})$, then the transition probability reads:

$$p(\mathcal{C}\to\mathcal{C}^{\prime})=\left\{\begin{array}[]{lr}\exp\left(-\frac{\Delta\chi^{2}}{2\Theta}\right),&\text{if}~{}\Delta\chi^{2}>0,\\ 1.0,&\text{if}~{}\Delta\chi^{2}\leq 0,\end{array}\right.$$

(36)

where $\Theta$ is an artificial system temperature and $\chi^{2}$ is interpreted as energy of the system. Third, the above two steps should be restarted periodically to avoid being trapped by local minima. Fourth, once all the Monte Carlo sampling tasks are finished, we should pick up the “best” solution which exhibits the smallest $\chi^{2}$, or select some “good” solutions with small $\chi^{2}$ and evaluate their arithmetic average. Finally, with the optimized $N_{p}$, and $A_{\gamma}$, and $P_{\gamma}$ parameters, the retarded Green’s function $G^{R}(\omega)$ can be easily evaluated by replacing $i\omega_{n}$ with $\omega+i\eta$ in Eqs. (21), (24), (28), and (31), where $\eta$ is a positive infinitesimal number. And the spectral density $A(\omega)$ is calculated by:

$$A(\omega)=-\frac{1}{\pi}\text{Im}G^{R}(\omega).$$

(37)

In the SPX method, the amplitudes and locations of the poles should be optimized by the Monte Carlo algorithm under some constraints [see Eqs. (23), (27), and (33)]. We note that these constraints are from the canonical relations for the fermionic and bosonic operators Negele and Orland (1998); Coleman (2016). They must be satisfied, or else the causality of the spectrum can not be guaranteed. But beyond that, more constraints are allowable. Further restrictions on the amplitudes and locations of the poles can greatly reduce the configuration space that needs to be sampled and enhance the possibility to reach the global minimum of the optimization problem. The possible strategies include: (1) Restrict $\{A_{\gamma}\}$ only; (2) Restrict $\{P_{\gamma}\}$ only; and (3) Restrict both $\{A_{\gamma}\}$ and $\{P_{\gamma}\}$ at the same time. These extra constraints can be deduced from a priori knowledge about the Matsubara Green’s function $G(i\omega_{n})$ and the spectral density $A(\omega)$. For example, for a molecule system, the amplitudes of the poles are likely close. On the other hand, if we know nothing about the input data and the spectra, we can always try some constraints. The universal trend is that the more reasonable the constraints, the smaller the $\chi^{2}$ function. This is the so-called constrained sampling algorithm. By combining it with the SPX method (dubbed C-SPX), the ability to resolve fine features in the spectra will be greatly enhanced. To the best of our knowledge, the constrained sampling algorithm was first proposed by A. W. Sandvik Sandvik (2016). And then it is broadly used in analytic continuations for spin susceptibilities of quantum many-body systems Shao et al. (2017); Ma et al. (2018). Quite recently, Shao and Sandvik summarized various approaches to mount the constraints and benchmark their performances in a comprehensive review concerning the SAC method Shao and Sandvik (2023). Due to the similarities of the SPX and SAC methods, it is believed that all the constraint tactics as suggested in Reference [Shao and Sandvik, 2023] should be applicable for the SPX method.

In analogy to the SAC method, the poles in the SPX method are distributed randomly in a real frequency grid Beach (2004); Sandvik (1998). This grid must be extremely dense and is usually linear. But in principle a nonuniform grid is possible. For example, Shao and Sandvik Shao and Sandvik (2023) have suggested a nonlinear grid with monotonically increasing spacing for the $\delta$ functions which are used to parameterize a spectrum that exhibits a sharp band edge. Since a spectral density can be viewed as a probability distribution Jarrell and Gubernatis (1996) and we notice that the distribution of the poles looks quite similar to the spectrum. So, it is natural to adjust the frequency grid dynamically to make sure that the grid density has an approximate distribution with the spectral density as obtained in the previous run. We adopt the following algorithm to manipulate the desirable frequency grid: (1) Calculate the integrated spectral function $\phi(\epsilon)$:

$$\phi(\epsilon)=\int^{\epsilon}_{\omega_{\text{min}}}A(\omega)d\omega,~{}\epsilon\in[\omega_{\text{min}},\omega_{\text{max}}].$$

(38)

(2) The new frequency grid $f_{i}$ is evaluated by:

$$f_{i}=\phi^{-1}(\lambda_{i}),~{}i=1,\cdots,N_{f},$$

(39)

where $\lambda_{i}$ is a linear mesh in $[\phi(\omega_{\text{min}}),\phi(\omega_{\text{max}})$], and $N_{f}$ denotes the number of grid points. Next, we should perform the analytic continuation simulation again to yield a new spectrum, then repeat steps (1) and (2). We find that the $\chi^{2}$ drops quickly in the first few iterations, and then approaches a constant value slowly. The spectrum is refined simultaneously. During the iterations, the frequency grid for the poles is adaptively refreshed via Eqs. (38) and (39), thus we call it the self-adaptive sampling algorithm. It is actually a new variation of the constrained sampling algorithm Sandvik (2016). More importantly, it is quite effective. According to our experiences, 3 $\sim$ 5 iterations are enough to obtain a convergent solution. In practice, we often adopt the spectrum generated by the MaxEnt method to initialize the frequency grid, and then employ the SPX method (dubbed SA-SPX) to refine this spectrum further.

The SPX method, together with the MaxEnt method Jarrell and Gubernatis (1996); Gubernatis et al. (1991), have been implemented in an open source software package, namely ACFlow Huang (2022), which is a full-fledged analytic continuation toolkit. This package supports analytic continuation for both fermionic and bosonic correlators.

The flowchart of the SPX method is illustrated in Figure 2. Next, we would like to explain some essential steps. (1) Initialize the kernel matrix $\Xi$. It is defined in Eqs. (22), (25), and (29). We should evaluate and save them during the initialization stage to improve the computational efficiency. (2) Initialize or reset $\{A_{\gamma},P_{\gamma}\}$. For each Monte Carlo run, $\{A_{\gamma},P_{\gamma}\}$ must be reinitialized to avoid getting trapped in local minima. (3) Calculate $\chi^{2}$ and $G(i\omega_{n})$. Here, $G(i\omega_{n})$ means the reproduced Matsubara data. It can be derived by using the present $\{A_{\gamma},P_{\gamma}\}$ via Eqs. (21), (24), (28), and (31). The loss function $\chi^{2}$ measures the distance between the input Matsubara data $\mathcal{G}(i\omega_{n})$ and the reproduced Matsubara data $G(i\omega_{n})$. It is evaluated by Eq. (34). (4) Try to update $A_{\gamma}$ or $P_{\gamma}$ randomly. Four Monte Carlo updates (Move 1 $\sim$ Move 4) are introduced to change the randomly chosen $A_{\gamma}$ or $P_{\gamma}$. (5) Accept new $A_{\gamma}$ or $P_{\gamma}$. When $\Delta\chi^{2}<0$, the Monte Carlo proposal is always accepted. (6) Accept new $A_{\gamma}$ or $P_{\gamma}$ probably. When $\Delta\chi^{2}>0$, it is still possible to accept the Monte Carlo proposal. The transition probability is determined by Eq. (36). In a standard simulated annealing algorithm, $\Theta$ should decline gradually from high temperature to low temperature Marinari (1996). However, in the present implementation, we just set $\Theta$ to a large value (about $10^{6}\sim 10^{9}$). We find that such a large $\Theta$ is essential to enable the Monte Carlo walker to escape the local minima and visit as many configurations as possible. (7) Backup $\chi^{2}$ and $\{A_{\gamma},P_{\gamma}\}$. For each Monte Carlo run, we always record the best solution, i.e., the smallest $\chi^{2}$ and the corresponding $\{A_{\gamma},P_{\gamma}\}$. (8) Enough runs? The Monte Carlo procedure is repeated many times in order to find out the “true” global minimum of $\chi^{2}$. (9) Choose the best $\{A_{\gamma},P_{\gamma}\}$. Once the outer iteration is finished, we obtain a collection of $\chi^{2}$ and $\{A_{\gamma},P_{\gamma}\}$. For the molecule cases, we should pick up the smallest $\chi^{2}$ and the corresponding $\{A_{\gamma},P_{\gamma}\}$. However, for the condensed matter cases, it would be better to select some “good” solutions and evaluate their averaged value. (10) Output $G(i\omega_{n})$ and $G^{R}(\omega)$. Finally, the Matsubara Green’s function $G(i\omega_{n})$ and retarded Green’s function $G^{R}(\omega)$ are calculated via Eq. (20) using the optimized $\{A_{\gamma},P_{\gamma}\}$.

In the present implementation, several numerical issues need to be emphasized.

• •

  The kernel matrix $\Xi(\omega_{n},\omega)$ should be computed at advance and kept in memory. This allows us to calculate the transition probability $p$ and the Green’s function $G(i\omega_{n})$ efficiently. The Matsubara frequencies $\omega_{n}$ are taken from the input data directly. And we have to create a frequency grid on the real axis for $\omega$ ($\omega\in[\omega_{\text{min}},\omega_{\text{max}}]$). In principle, the poles are distributed in a continuous frequency space. So in order to improve the computational accuracy, the number of grid points ($N_{f}$) should be as large as possible.

• •

  According to Eq. (36), even when $\chi^{2}(\mathcal{C}^{\prime})>\chi^{2}(\mathcal{C})$, the transition from $\mathcal{C}^{\prime}$ to $\mathcal{C}$ is not always rejected. Though the transition probability $p$ could be quite small, we still have the chance to accept a less optimal solution than what we currently have and escape the local minima. This transition probability is largely controlled by the $\Theta$ parameter, which behaves as the system’s annealing temperature Kirkpatrick et al. (1983); Marinari (1996). In accordance with the spirit of the simulated annealing algorithm, $\Theta$ should be decreased gradually. However, in the present implementation, $\Theta$ is fixed and considered as a user-supplied parameter. According to our experiences, the preferred value of $\Theta$ is between $10^{6}$ and $10^{9}$.

• •

  In the SPX method, the Monte Carlo sampling procedure should be repeated from scratch many times. In each run, the $\chi^{2}$ and the corresponding Monte Carlo configuration $\mathcal{C}$ will be tracked. We find that when the final spectral density exhibits multiple $\delta$-like peaks, it is wise to pick up the “best” solution whose $\chi^{2}$ is the smallest. However, when the spectral density is expected to be broad and smooth, it is better to calculate an arithmetic average of some selected “good” solutions. We just adopted the following strategy to filter the solutions Krivenko and Harland (2019); Krivenko and Mishchenko (2022). At first, we try to calculate the median or mean value of the collected $\chi^{2}$ data, i.e., $\langle\chi^{2}\rangle$. Then, only the solutions, whose $\chi^{2}$ are smaller than $\langle\chi^{2}\rangle/\alpha_{\text{good}}$, are retained. Here, $\alpha_{\text{good}}$ is an adjustable parameter. Its optimal value is around $1.0\sim 1.2$.

• •

  The SPX method has been parallelized to improve computational efficiency. It is straightforward to launch multiple Monte Carlo processes simultaneously to accelerate the calculation.

We at first benchmark the SPX method on the condensed matter cases, in which the spectral functions are usually broad and smooth. The spectral functions are assumed to be the superposition of a few Gaussian peaks:

$$A(\omega)=\sum^{N_{g}}_{i=1}w_{i}\exp\left[\frac{-(\omega-\epsilon_{i})^{2}}{2\Gamma^{2}_{i}}\right],$$

(41)

where $N_{g}$ is the number of Gaussian peaks, $w_{i}$, $\epsilon_{i}$, and $\Gamma_{i}$ denote the weight, position, and broadening of the $i$-th Gaussian peak, respectively. The spectral functions are then back-continued to the Matsubara frequency axis by using Eq. (8) with $\beta=10.0$. The Matsubara data, supplemented with random noises, are used as inputs for the SPX method.

Figure 3 shows the analytic continuation results by the SPX method for three typical scenarios: (i) two broad Gaussian peaks separated by a sizable gap ($\Delta_{\text{gap}}\approx 2.0$, $N_{g}=2$, $w_{1}=w_{2}=0.5$, $\epsilon_{1}=-\epsilon_{2}=2.5$, $\Gamma_{1}=\Gamma_{2}=0.5$), (ii) two Gaussian peaks with a “pseudogap” ($N_{g}=2$, $w_{1}=1.0$, $w_{2}=0.3$, $\epsilon_{1}=0.5$, $\epsilon_{2}=-2.5$, $\Gamma_{1}=0.2$, $\Gamma_{2}=0.8$), and (iii) two Gaussian peaks plus a sharp quasiparticle resonance peak ($N_{g}=3$, $w_{1}=1.0$, $w_{2}=0.3$, $w_{3}=0.4$, $\epsilon_{1}=0.0$, $\epsilon_{2}=-2.5$, $\epsilon_{3}=1.5$, $\Gamma_{1}=0.02$, $\Gamma_{2}=0.8$, $\Gamma_{3}=0.5$). As is evident in Fig. 3(a)-(c), the major features of the synthetic spectral functions, including the energy gap, positions and weights of the peaks, are well recovered by the SPX method. The only exception is that the full width at half maximum and height of the quasiparticle resonance peak for scenario (iii) [see Fig. 3(c)] are somewhat overestimated. The MaxEnt method Jarrell and Gubernatis (1996); Gubernatis et al. (1991) leads to slightly better results for scenarios (i) and (ii). But it is also unable to recover the sharp quasiparticle resonance peak for scenario (iii). Fig. 3(d)-(f) illustrate the distributions of poles for the three scenarios. Not surprisingly, they manifest approximate characteristics to the correspondingly spectral functions. Overall, for the condensed matter cases, the performance of the SPX method is comparable with the other fitting-based analytic continuation methods.

Now let us turn to the molecule cases, in which the spectral functions usually exhibit multiple discrete $\delta$-like peaks Fei et al. (2021b); Ying (2022b). To construct the Matsubara data, Eq. (21) is utilized ($\beta=20.0$). Random noises are added to the synthetic Matsubara data as described above.

Three typical scenarios are prepared to examine the SPX method: (i) an off-center $\delta$-like peak ($N_{p}=1$, $A_{1}=1.0$, $P_{1}=1.0$), (ii) four low-frequency $\delta$-like peaks with equal amplitudes ($N_{p}=4$, $A_{1}=A_{2}=A_{3}=A_{4}=0.25$, $P_{1}=-P_{2}=4.0$, $P_{3}=-P_{4}=1.0$), and (iii) eight $\delta$-like peaks distributed over a wide frequency range ($N_{p}=8$, $A_{1}=A_{2}=A_{3}=A_{4}=A_{5}=A_{6}=A_{7}=A_{8}=0.125$, $P_{1}=8.0$, $P_{2}=4.0$, $P_{3}=-P_{5}=0.5$, $P_{4}=0.0$, $P_{6}=-1.0$, $P_{7}=-5.0$, $P_{8}=-9.0$). In the calculations, the number of poles is considered as a priori knowledge. Besides Eq. (23), there are no additional limitations for the amplitudes and locations of the poles. Especially, to obtain reasonable solutions for scenarios (ii) and (iii), the numbers of individual Monte Carlo runs are increased to $2\times 10^{5}$. Figure 4 shows the analytic continuation results. Clearly, both the amplitudes and locations of the peaks (poles) are resolved accurately by the SPX method. Not only the low-frequency multiplets but also the high-frequency sharp peaks are well reproduced. For the three scenarios, the traditional analytic continuation methods, such as MaxEnt Jarrell and Gubernatis (1996); Gubernatis et al. (1991), fail to distinguish the multiple peaks. They can only give rise to a envelop curve for the $\delta$-like peaks.

Though the three pole models are well solved by the SPX method, quite a lot of computational resources are consumed, especially for scenarios (ii) and (iii). Just as mentioned before, extra constraints could be imposed within the parameterization of the poles, which is a widely used trick for the SAC method and its variants Sandvik (1998); Beach (2004); Mishchenko et al. (2000); Sandvik (2016); Vafayi and Gunnarsson (2007); Fuchs et al. (2010); Shao et al. (2017); Ghanem and Koch (2020a); Shao and Sandvik (2023); Ghanem and Koch (2020b); Syljuåsen (2008). Such constraints, which represent some kind of innate knowledge [e.g., existence of the sharp band edge, and properties of the poles (including number, amplitudes, and approximate locations)], could significantly improve fidelity of the calculated spectrum and computational efficiency of the SPX method.

The C-SPX method is at first examined by using the four-pole model [i.e. scenario (ii)]. Because the analytic continuation calculation without extra constraints has been done before (see Figure 4), we try to limit the locations of the poles with $P_{\gamma}\in[-4.5,-3.5]\bigcup[-1.5,-0.5]\bigcup[0.5,1.5]\bigcup[3.5,4.5]$ and rerun the calculations (We will discuss the tricks about how to speculate the restricted zones for the poles in Section VI). The number of individual Monte Carlo is reduced to 1000. The calculated results are displayed in Figure 5. Figure 5(a) shows the distribution of the poles. The horizontal and vertical dashed lines denote the exact amplitudes and locations of the poles, respectively. Figure 5(b) shows the collected $\chi^{2}$ for individual Monte Carlo runs. We find that once the constraints are imposed, the computational accuracy of not only the locations but also the amplitudes of the poles is greatly improved, and it becomes easier to figure out the global minimization.

And then the eight-pole model [i.e. scenario (iii)] is used to test the C-SPX method. The amplitudes of the eight poles are the same, but they reside in a wide frequency range. It is rather difficult to get the correct solution by using the traditional analytic continuation methods. As for the SPX method, it is tough to arrive at the global minimum unless the number of individual Monte Carlo runs is increased up to $10^{6}$, which is very time-consuming. However, if we assume that the amplitudes of the poles are known and try to optimize their locations by the C-SPX method, it is easy to reproduce the correct spectrum. Now the number of individual Monte Carlo runs can be reduced to 1000. Figure 6 shows the benchmark results. For the SPX method, though the peaks at the low-frequency region are roughly resolved, it fails to reproduce the high-frequency peaks. However, by using the C-SPX method, both the high-frequency and low-frequency peaks are accurately resolved. Furthermore, we find that the values of the goodness-of-fit function $\chi^{2}$ are more scattered when the SPX method is used [See Fig. 6(b)]. It implies that the SPX method is easier trapped by the local minima than the C-SPX method.

The third example concerns the analytic continuation of Matsubara Green’s function of a BCS superconductor Beach (2004). Its spectral function reads:

$$A(\omega)=\left\{\begin{array}[]{lr}\frac{1}{W}\frac{|\omega|}{\sqrt{\omega^{2}-\Delta^{2}}},{}&\text{if}~{}\Delta<|\omega|<W/2.\\ 0,&\text{otherwise}.\end{array}\right.$$

(42)

Here, $W$ denotes the bandwidth ($W=6.0$), and $\Delta$ is used to control the gap’s size ($\Delta=0.5$). This spectrum is comprised of flat shoulders, steep peaks, and sharp gap edges. These unusual features pose severe challenges to the existing analytic continuation methods. Figure 7(a)-(b) shows the analytic continuation results obtained by the SPX method without any constraints. Clearly, the calculated spectrum exhibits extra shoulder peaks around $\pm 2.0$ and long tails for $|\omega|>3$. The energy gap and the sharp band edges are not well captured. Then constraints are applied to the locations of the poles. They are allowed to appear in a restricted frequency range ($[-3.0,0.5]\bigcup[0.5,3.0]$). The calculated results are shown in Fig. 7(c)-(d). It is clear that the major characteristics of the exact spectrum are well reproduced by the C-SPX method. The excrescent peaks around $\pm 2.0$ are mostly suppressed, leaving small ripples.

In this subsection, let us concentrate on a realistic case, Matsubara Green’s function from quantum many-body simulation. We just consider the following single-band half-filling Hubbard model:

$$H=-t\sum_{\langle ij\rangle,~{}\sigma}c^{\dagger}_{i\sigma}c_{j\sigma}-\mu\sum_{i}n_{i}+U\sum_{i}n_{i\uparrow}n_{i\downarrow},$$

(43)

where $t$ is the hopping parameter ($t=0.5$), $\mu$ is the chemical potential ($\mu=1.0$), $U$ is the Coulomb interaction ($U=2.0$), $n$ is the occupation number ($n=1.0$), $\sigma$ denotes the spin, $i$ and $j$ are site indices. The inverse system temperature is $\beta=10.0$. This model is defined in a Bethe lattice. It was solved by the single-site dynamical mean-field theory (dubbed DMFT) Georges et al. (1996) in combination with the hybridization expansion continuous-time quantum Monte Carlo impurity solver (dubbed CT-HYB) Gull et al. (2011). All the calculations were done by using the iQIST software package Huang et al. (2015); Huang (2017). The major outputs of the DMFT + CT-HYB calculations are the Green’s function and the self-energy function at the imaginary axis. In this example, the Matsubara Green’s function is analytically continued to extract the spectral function by the MaxEnt method and the SPX method. The Matsubara self-energy function will be treated in the next subsection.

The analytic continuation results are shown in Figure 8. Since the Coulomb interaction is comparable with the bandwidth ($W=4t=2.0$), the ground state of this model should be metallic, but close to the Mott metal-insulator transition Imada et al. (1998). Thus, it is not strange that the spectral function exhibits a three-peak structure (i.e. the quasiparticle resonance peak in the vicinity of the Fermi level + lower and upper Hubbard bands), which is a hallmark of the strongly correlated metallic systems Georges et al. (1996). We can observe that the spectral functions given by the MaxEnt method and the SPX method [see Fig. 8(a)] match with each other, so the results should be reliable. Figure 8(b)-(c) show the real and imaginary parts of the reproduced Green’s functions, which are compared with the input data. Apparently, the original Matsubara data are well reproduced, which implies the derived pole expansion formula [see Eq. (21)] is reasonable.

Besides one-particle Green’s functions, it often is necessary to analytically continue the self-energy functions after the DMFT simulations Georges et al. (1996). Now we would like to demonstrate how to convert the self-energy function from Matsubara to real frequencies by using the SPX method. The high-frequency asymptotic behaviors of self-energy are different from those of one-particle Green’s function. When approaching the high-frequency limit, the real part is a non-zero constant, and the imaginary part decays like $\Sigma_{1}/(i\omega_{n})$ where $\Sigma_{1}$ is the first moment of the self-energy Kaufmann and Held (2023). So the input data of self-energy should be preprocessed. The original data are taken from the DMFT solution of Eq. (43). At first, the Hartree term $\Sigma_{H}$ is subtracted from $\Sigma(i\omega_{n})$:

$$\tilde{\Sigma}(i\omega_{n})=\Sigma(i\omega_{n})-\Sigma_{H}.$$

(44)

Note that $\Sigma_{H}$ is approximately equal to the asymptotic value of the real part of $\Sigma(i\omega_{n})$ when $n$ goes to infinite. It is also the zeroth moment of the self-energy. Sometimes $\tilde{\Sigma}(i\omega_{n})$ could be further normalized by division through its first moment $\Sigma_{1}$ Kaufmann and Held (2023); Kraberger et al. (2017),

$$\tilde{\Sigma}(i\omega_{n})=\frac{\Sigma(i\omega_{n})-\Sigma_{H}}{\Sigma_{1}},$$

(45)

but it is not necessary. Then, we perform analytic continuation by using the SPX method as usual and take $\tilde{\Sigma}(i\omega_{n})$ as the input data. The analytic continuation procedure will convert $\tilde{\Sigma}(i\omega_{n})$ into $\tilde{\Sigma}(\omega)$. Finally, $\tilde{\Sigma}(\omega)$ is supplemented with the Hartree term $\Sigma_{H}$ to get $\Sigma(\omega)$:

$$\Sigma(\omega)=\tilde{\Sigma}(\omega)+\Sigma_{H}.$$

(46)

The benchmark results are illustrated in Fig. 9. The results obtained by the MaxEnt method Jarrell and Gubernatis (1996); Gubernatis et al. (1991) are also displayed as a comparison. When $\omega>0$, the results obtained by both methods are consistent with each other. If $\omega<0$, the results obtained by the SPX method exhibit an additional peak near $\omega=-2.0$ [see Fig. 9(a)-(b)]. On the contrary, the MaxEnt method only yields a smooth spectrum. Anyway, the reproduced Matsubara data generated by both methods are quite close to the original Matsubara data [see Fig. 9(c)-(d)].

The first example we address here is a Gaussian model. The analytic structure of the exact spectral function is:

$$\frac{A(\omega)}{\omega}=\alpha_{1}\exp\left[-\frac{(\omega-\epsilon_{1})^{2}}{2\gamma^{2}_{1}}\right]+\alpha_{2}\exp\left[-\frac{(\omega-\epsilon_{2})^{2}}{2\gamma^{2}_{2}}\right],$$

(47)

where $\alpha_{1}=\alpha_{2}=0.5$, $\gamma_{1}=\gamma_{2}=0.5$, $\epsilon_{1}=-\epsilon_{2}=2.5$. Clearly, $A(\omega)$ is an odd function. It exhibits two distinct peaks at energies $\epsilon_{1}$ and $\epsilon_{2}$. And the two Gaussian peaks are antisymmetric about $\omega=0.0$. There is a huge gap between the two peaks ($\Delta_{\text{gap}}\approx 3.0$). The Matsubara data are generated by using Eqs. (12) and (14). The inverse system temperature $\beta$ is 10.0.

In Fig. 10(a), the exact spectrum is drawn and compared to the calculated spectra as obtained by the SPX and MaxEnt methods. Both methods could capture the major characteristics, including the two-peaks structure and the big gap, of the spectrum. The lower panel of Fig. 10 shows a typical distribution of the poles. The overall profile of this distribution is also Gaussian-like. There are a small fraction of poles in the band gap region, but their contributions are trivial.

For the second example, we consider a four-pole model. Its analytic form is as follows:

$$G(z)=\sum^{4}_{i=1}\frac{\alpha_{i}}{z-\epsilon_{i}},$$

(48)

where $\alpha_{1}=\alpha_{3}=0.5$, $\alpha_{2}=\alpha_{4}=-0.5$, $\epsilon_{1}=-\epsilon_{2}=4.0$, and $\epsilon_{3}=-\epsilon_{4}=1.0$. The exact spectrum exhibits four off-centered $\delta$-like peaks at $\epsilon_{1}$, $\epsilon_{2}$, $\epsilon_{3}$, and $\epsilon_{4}$. These peaks are also antisymmetric with respect to $\omega=0.0$. Note that similar spectra are usually seen in molecule systems (such as the Hubbard dimer) and the momentum-resolved Kohn-Sham eigenvalues (i.e. the “band structure”) of condensed matter systems Fei et al. (2021a); Huang et al. (2023). We just use this equation or Eq. (24) to generate the Matsubara data. The inverse system temperature $\beta$ is 20.0.

In this example, we respect the symmetries of the Green’s function and spectral function. In other words, $G(z)$ is treated as a bosonic correlator of a Hermitian operator, and the corresponding spectrum can be defined in the positive half-axis only. The exact and calculated spectral functions are shown in Fig. 11(a). As for the MaxEnt method, it resolves the low-frequency peak at $\omega=1.0$ quite well. However, it fails to resolve the high-frequency peak at $\omega=4.0$. The resulting high-frequency peak is much broader and smoother than the exact one. On the other hand, the SPX method does an excellent job. Not only the positions but also the heights of the two $\delta$-like peaks are precisely reproduced. In Fig. 11(b) the poles belonged to the “best” solution are visualized. Clearly, the weights and locations of the poles agree quite well with the exact values.

For the third example, we consider a resonance model. Its spectral function is defined as follows:

$$A(\omega)=\left\{\begin{array}[]{lr}\frac{1}{W}\frac{\omega}{\sqrt{\omega^{2}-\Delta^{2}}},{}&\text{if}~{}\Delta<\omega<W/2.\\ 0,&\text{otherwise}.\end{array}\right.$$

(49)

It is obviously a variant of Eq. (42). The $W$ and $\Delta$ parameters are the same with those used in M₀₃ (see Section V.3). This spectrum has finite weights only at the positive half-axis. It should exhibit a sharp resonance peak at $\omega=\Delta$ and a broad platform [see Fig. 12(a)]. We try to generate the bosonic Green’s function via Eqs. (16) and (18). The inverse system temperature is $\beta=10.0$.

As usual, both the MaxEnt and SPX methods are employed to solve the analytic continuation problem. The results are visualized in Fig. 12(a). We observe that the simulated spectra exhibit two broad humps at $\omega\approx$ 0.8 and 2.1, respectively. However, none of the important features of the ideal spectrum is captured. Therefore we resort to the C-SPX method again. Since there are two band edges in the true spectrum, we introduce two cutoff parameters, namely $\Lambda^{l}$ and $\Lambda^{r}$, where $\omega_{\text{min}}<\Lambda^{l}<\Lambda^{r}<\omega_{\text{max}}$. Here, the superscripts “$l$” and “$r$” mean the left and right band edges, respectively. Thus, the forbidden region for the poles becomes $[\omega_{\text{min}},\Lambda^{l}]\bigcup[\Lambda^{r},\omega_{\text{max}}]$. Next, we try to calibrate the two parameters and evaluate the corresponding $\chi^{2}(\Lambda^{l})$ and $\chi^{2}(\Lambda^{r})$. The calculated results are presented in Fig. 12(b). We find that neither $\chi^{2}(\Lambda^{l})$ nor $\chi^{2}(\Lambda^{r})$ is monotonic function. They both exhibit minima:

	$\displaystyle\Lambda^{l}_{c}$	$\displaystyle=$	$\displaystyle\mathop{\arg\min}\limits_{\Lambda^{l}}\chi^{2}(\Lambda^{l}),$		(50)
	$\displaystyle\Lambda^{r}_{c}$	$\displaystyle=$	$\displaystyle\mathop{\arg\min}\limits_{\Lambda^{r}}\chi^{2}(\Lambda^{r}).$		(51)

We can conclude that $\Lambda^{l}_{c}=0.5$ and $\Lambda^{r}_{c}=3.0$. They are the optimal cutoff parameters. Such that the restrictions are fixed. We conduct the analytic continuation simulation with the SPX method again. The simulated spectrum is shown in Fig. 12(a) as well. It is apparent that the C-SPX method outperforms the SPX method and the MaxEnt method in this example. Both the sharp resonance peak at $\omega=0.5$ and the broad platform at $1.0<\omega<3.0$ are well reproduced. The only deviation is the small ripples in the platform region as seen in the simulated spectrum. But these ripples can be further suppressed by using more poles and collecting more “good” solutions.

For the fourth example, we consider a spin model in the XY chain. Its Hamiltonian reads:

$$H=J_{xy}\sum_{\langle ij\rangle}\left(S^{i}_{x}S^{j}_{x}+S^{i}_{y}S^{j}_{y}\right),$$

(52)

where $S_{\alpha}=\frac{\sigma_{\alpha}}{2}$ are the spin-$\frac{1}{2}$ operators ($\alpha=x,~{}y,~{}z$), $\langle ij\rangle$ denotes the nearest neighbors, and $J_{xy}=1$. This model can be exactly solved by the Jordan-Wigner transformation. The energy spectrum is given by:

$$\epsilon_{k}=J_{xy}\cos(ka),$$

(53)

where $a$ is the lattice spacing. The local spin-spin correlation function $\chi_{zz}(\tau)=\langle{S_{z}(\tau)S_{z}(0)}\rangle$ is basically a density-density correlator in the sense of spinless fermions. Its exact expression reads Katsura et al. (1970):

$$\chi_{zz}(\tau)=\int_{-\pi}^{\pi}\frac{\mathrm{d}k\mathrm{d}k^{\prime}}{(2\pi)^{2}}\frac{e^{\tau\epsilon_{k}}e^{(\beta-\tau)\epsilon_{k^{\prime}}}}{(1+e^{\beta\epsilon_{k}})(1+e^{\beta\epsilon_{k^{\prime}}})}.$$

(54)

Therefore the correspondingly spectral function is:

$$A(\omega)=\int_{-\pi}^{\pi}\frac{\mathrm{d}k\mathrm{d}k^{\prime}}{2\pi}\frac{(1-e^{-\beta\omega})e^{\beta\epsilon_{k^{\prime}}}}{(1+e^{\beta\epsilon_{k}})(1+e^{\beta\epsilon_{k^{\prime}}})}\delta(\omega-\epsilon_{k^{\prime}}+\epsilon_{k}).$$

(55)

It is obvious that there are thresholds at $\omega=\pm 2J_{xy}$ because $\max|\epsilon_{k}-\epsilon_{k^{\prime}}|$ is $2|J_{xy}|$.