# ACFlow: An open source toolkit for analytical continuation of quantum Monte
Carlo data

Li Huang  Science and Technology on Surface Physics and Chemistry Laboratory,
P.O. Box 9-35, Jiangyou 621908, China

###### Abstract

The purpose of analytical continuation is to establish a real frequency
spectral representation of single-particle or two-particle correlation
function (such as Green’s function, self-energy function, and dynamical
susceptibilities) from noisy data generated in finite temperature quantum
Monte Carlo simulations. It requires numerical solutions of a family of
Fredholm integral equations of the first kind, which is indeed a challenging
task. In this paper, an open source toolkit (dubbed ACFlow) for analytical
continuation of quantum Monte Carlo data is presented. We at first give a
short introduction to the analytical continuation problem. Next, three primary
analytical continuation algorithms, including maximum entropy method,
stochastic analytical continuation, and stochastic optimization method, as
implemented in this toolkit are reviewed. And then we elaborate major
features, implementation details, and basic usage of this toolkit. Finally,
four representative examples are shown to demonstrate usefulness and
flexibility of the ACFlow toolkit.

###### keywords:

Quantum Monte Carlo simulation , Analytical continuation problem , Maximum
entropy method , Stochastic analytical continuation , Stochastic optimization
method

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: ACFlow
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: https://github.com/huangli712/ACFlow
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions (please choose one): GPLv3
Programming language: Julia
Supplementary material:
Journal reference of previous version:*
Does the new version supersede the previous version?:*
Reasons for the new version:*
Summary of revisions:*
Nature of problem (approx. 50-250 words): Most of the quantum Monte Carlo
methods work on imaginary axis. In order to extract physical observables and
compare them with the experimental results, analytical continuation must be
done in the post-processing stage to convert the quantum Monte Carlo simulated
data from imaginary axis to real axis.
Solution method (approx. 50-250 words): Three established analytical
continuation methods, including maximum entropy method, stochastic analytical
continuation, and stochastic optimization method, have been implemented in the
ACFlow toolkit.
Additional comments including restrictions and unusual features (approx.
50-250 words): The ACFlow toolkit is written by pure Julia language. It is
highly optimized and parallelized. It can be executed interactively in a
Jupyter notebook environment.

##  1 Introduction

It is well-known that quantum Monte Carlo (QMC) method is a powerful and exact
numerical approach, and has been widely used in many research fields, such as
nuclear physics [1], condense matter physics [2], and many-body physics [3].
In this paper, we just focus on the finite temperature QMC algorithms, which
are used to solve the interacting lattice models or quantum impurity models
[4]. Generally speaking, the simulated results of QMC methods are some sorts
of single-particle or two-particle correlation functions, which are usually
defined on imaginary time axis ($\tau\equiv-it$) or Matsubara frequency axis
($i\omega_{n}$). Therefore, they can’t be compared directly with the
correspondingly experimental results, including but not limited to the
electronic density of states $A(\omega)$, optical conductivity
$\sigma(\omega)$, dynamical structure factor $S(\mathbf{q},\omega)$, and so
on. It is necessary to convert the QMC simulated results from imaginary time
axis or Matsubara frequency axis to real axis (i.e. $\tau\to\omega$ or
$i\omega_{n}\to\omega$), which is the origin of the analytical continuation
problem.

Let’s concentrate on the following Fredholm integral equation of the first
kind:

$g(y)=\int K(y,x)f(x)~{}dx.$ (1)

Here, $K(y,x)$ is the known kernel function, $f(x)$ is the model function, and
$g(y)$ denotes the raw data. Given $f(x)$, it is quite easy to get $g(y)$ via
numerical integration. However, given $g(y)$, solving the Fredholm integral
equation reversely to get $f(x)$ is not as easy as expected. There is no
universal solution. Notice that the so-called analytical continuation problem
can be reformulated in terms of the Fredholm integral equation. Thus, its
objective is to seek a reasonable $f(x)$ to satisfy the above equation. The
QMC simulated data $g(y)$ are noisy and the kernel function $K(y,x)$ is ill
conditioned, which make analytical continuation of QMC simulated data a huge
challenge. In order to solve this problem, peoples have developed numerous
methods in the past decades. These methods include the least square fitting
method, singular value decomposition [5, 6], Padé approximation [7, 8, 9, 10],
Tikhonov-Philips regularization method, maximum entropy method [11, 12],
stochastic analytical continuation [13, 14], stochastic optimization method
[15, 16], sparse modelling method [17], and machine learning method [18, 19,
20], etc. However, each method has its pros and cons. None of these methods
can override the others. The analytical continuation problem is still far away
from being completely solved.

In recent years, quite a few analytical continuation codes have been released,
including maxent (by Mark Jarrell) [12], maxent (in ALPSCore) [21],
$\Omega$Maxent [22], ana_cont [23], SOM (in TRIQS) [24, 25], Stoch (in ALF)
[26], just to name a few. We note that the maximum entropy method has
dominated this field for quite a long time. Thus most of these codes only
implement the maximum entropy method [12, 21, 22, 23]. It is rather difficult
to crosscheck the simulated results obtained by various analytical
continuation methods. In addition, the features of the available codes are
quite limited and hard to examine new algorithms. In order to fill in this
gap, we would like to present a new open source toolkit, called ACFlow, for
analytical continuation. This toolkit implements three primary analytical
continuation methods, including the maximum entropy method, stochastic
analytical continuation, and stochastic optimization method, within an united
framework. It provides an easy-to-used library and application interface. Some
diagnostic and analytical tools are also available. With ACFlow, the users can
easily setup and execute analytical continuation calculations, and validate
the obtained results. We believe that this toolkit will play a vital role in
solving analytical continuation problems.

The rest of this paper is organized as follows. In section 2, background of
the analytical continuation problem is introduced. In section 3, basic
principles and key ingredients of the three analytical continuation methods as
implemented in the ACFlow toolkit are summarized. Section 4 gives a brief
overview about ACFlow’s main features and structures. Section 5 is the major
part of this paper, it explains basic usage, input and output files of ACFlow.
In order to demonstrate usefulness of this toolkit, four typical examples are
illustrated in section 6. Finally, section 7 serves as a short conclusion.

##  2 Problem

###  2.1 Finite temperature Green’s functions

Under the Wick’s rotation $t\to i\tau$, the time evolution operator in the
Heisenberg picture $e^{itH}$ will be replaced by $e^{-\tau H}$. Such a
transformation will increase efficiency of QMC random walking and suppress
numerical oscillation (when $t$ is large, the periodic oscillation of
$e^{itH}$ is quite obvious). This is an important reason why most of the
finite temperature QMC algorithms are formulated in imaginary time axis. The
outputs of finite temperature QMC simulations are usually single-particle or
two-particle correlation functions. For example, the single-particle Green’s
function $G(\tau)$ is defined as follows:

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

where $\tau$ denotes imaginary time, $\mathcal{T}_{\tau}$ denotes time-ordered
operator, and $d$ and $d^{\dagger}$ are annihilation and creation operators,
respectively. The Matsubara Green’s function $G(i\omega_{n})$ can be measured
by QMC simulations or constructed from $G(\tau)$ via direct Fourier
transformation:

$G(i\omega_{n})=\int^{\beta}_{0}d\tau~{}e^{-i\omega_{n}\tau}G(\tau),$ (3)
$G(\tau)=\frac{1}{\beta}\sum_{n}e^{i\omega_{n}\tau}G(i\omega_{n}).$ (4)

Here, $\beta$ means the inverse temperature ($\beta\equiv 1/T$) and
$\omega_{n}$ is the Matsubara frequency. Note that $\omega_{n}$ is equal to
$(2n+1)\pi/\beta$ for fermions and $2n\pi/\beta$ for bosons ($n$ is an
integer).

###  2.2 Spectral density

Clearly, neither $G(\tau)$ nor $G(i\omega_{n})$ can be observed
experimentally. We have to extract dynamical response function, i.e., the
spectral density $A(\omega)$, from them. $A(\omega)$ is indeed an observable
quantity. It is related to $G(\tau)$ via the following Laplace transformation:

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

where +(-) in the denominator is for fermionic (bosonic) system.
$G(i\omega_{n})$ and $A(\omega)$ manifest similar relation:

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

It is obvious that Eq. (5) and Eq. (6) are indeed two special forms of the
Fredholm integral equation of the first kind [see Eq. (1)]. So, the central
problem of analytical continuation is to search optimal $A(\omega)$ for given
$G(\tau)$ or $G(i\omega_{n})$.

Sometimes the spectral density $A(\omega)$ is called spectral function in the
references. It is tied to the imaginary part of real frequency Green’s
function $G(\omega)$:

$A(\omega)=-\frac{1}{\pi}\rm{Im}G(\omega).$ (7)

From Im$G(\omega)$, Re$G(\omega)$ could be calculated via the Kramers-Kronig
transformation:

$\mathrm{Re}G(\omega)=\frac{1}{\pi}\mathcal{P}\int_{-\infty}^{\infty}d\omega^{\prime}~{}\frac{\mathrm{Im}G(\omega^{\prime})}{\omega^{\prime}-\omega},$
(8)

where $\mathcal{P}$ means Cauchy principal value. Besides Eq. (5) and Eq. (6),
$A(\omega)$ has to obey some additional constraints or sum-rules. For
fermionic systems, the spectral functions must be positive:

$A(\omega)\geq 0.$ (9)

While for bosonic systems, the constraint becomes:

$\text{sign}(\omega)A(\omega)\geq 0.$ (10)

In addition, the spectral function $A(\omega)$ is always bounded,

$\int^{+\infty}_{-\infty}d\omega~{}A(\omega)<\infty.$ (11)

It can be utilized to normalize the resulting spectral function.

###  2.3 Kernel functions

Eq. (5) and Eq. (6) can be reformulated as follows:

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

and

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

where $K(\tau,\omega)$ and $K(\omega_{n},\omega)$ are the so-called kernel
functions. Their definitions are as follows:

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

and

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

where +(-) in the denominator of Eq. (14) stands for fermions (bosons).

As mentioned above, the kernel function is quite strange. The values of
$K(\tau,\omega)$ could change by tens of orders of magnitude. Especially, at
large positive and negative frequencies, $K(\tau,\omega)$ is exponentially
small. It implies that at large $|\omega|$ the features of $A(\omega)$ are
sensitive to the fine structures of $G(\tau)$. However, the data of $G(\tau)$
provided by QMC simulations are always fluctuant and noisy [27]. Tiny
deviations in $G(\tau)$ from its expected values can lead to enormous changes
in $A(\omega)$. Thus, analytical continuation is often characterized as an
ill-posed problem [12].

In principle, for incomplete and noise $G(\tau)$ or $G(i\omega_{n})$, the
number of spectral functions $A(\omega)$ that satisfy Eq. (12) and Eq. (13) is
infinite. So the question becomes which $A(\omega)$ should be chosen. Now
there are two different strategies to solve this problem. The first one is to
choose the most likely $A(\omega)$. The second one is to evaluate the average
of all the candidate spectral functions. In next section, we will introduce
three primary analytical continuation methods that follow the two strategies
and have been implemented in the ACFlow toolkit. For the sake of simplicity,
we will concentrate on analytical continuation of imaginary time Green’s
functions in main text.

##  3 Methods

###  3.1 Maximum entropy method

Perhaps the maximum entropy method is the most frequently used approach for
analytical continuation problems in the last decades [11, 12] because of its
high computational efficiency. Next, we will discuss the basic principle and
several variants of it.

####  3.1.1 Bayesian inference

Bayes’s theorem is the cornerstone of the maximum entropy method. Given two
events $a$ and $b$, Bayes’s theorem says:

$P[a|b]P[b]=P[b|a]P[a],$ (16)

where $P[a]$ is the probability of event $a$, $P[a|b]$ is the conditional
probability of event $a$ with given event $b$. In the scenario of analytical
continuation problem, $\bar{G}(\tau)$ and $A(\omega)$ are treated as two
events, where $\bar{G}(\tau)$ denotes the measured value of $G(\tau)$. So the
best solution for $A(\omega)$ is of course the one that maximizes
$P[A|\bar{G}]$, which is called the posterior probability. According to the
Bayes’s theorem, we get

$P[A|\bar{G}]=\frac{P[\bar{G}|A]P[A]}{P[\bar{G}]},$ (17)

where $P[\bar{G}|A]$ is the likelihood function, $P[A]$ is the prior
probability, and $P[\bar{G}]$ is the evidence. Since the evidence is a
normalization constant depending on the prior probability and the likelihood
function only, it is ignored in the following discussions. Thus,

$P[A|\bar{G}]\propto P[\bar{G}|A]P[A].$ (18)

####  3.1.2 Posterior probability

In the maximum entropy method, the likelihood function $P[\bar{G}|A]$ is
assumed to be in direct proportion to $e^{-\chi^{2}/2}$. Here, $\chi^{2}$ is
named as goodness-of-fit function. It measures the distance between
$\bar{G}(\tau)$ and reconstructed imaginary time Green’s function
$\tilde{G}(\tau)$:

$\chi^{2}=\sum^{L}_{i=1}\left[\frac{\bar{G}_{i}(\tau)-\tilde{G}_{i}(\tau)}{\sigma_{i}}\right]^{2},$
(19) $\tilde{G}_{i}=\sum_{j}K_{ij}A_{j}.$ (20)

Here, $L$ is number of imaginary time points, $\sigma$ denotes the error bar
(standard deviation) of $\bar{G}(\tau)$. $K_{ij}$ and $A_{j}$ are discrete
kernel and spectral functions, respectively. On the other hand, the prior
probability $P[A]$ is supposed to be in direct proportion to $e^{\alpha S}$,
where $\alpha$ is a regulation parameter and $S$ means entropy. Sometimes $S$
is also known as the Kullback-Leibler distance. Its formula is as follows:

$S=\int
d\omega\left(A(\omega)-m(\omega)-A(\omega)\log\left[\frac{A(\omega)}{m(\omega)}\right]\right),$
(21)

where $m(\omega)$ is the default model function.

According to the Bayes’s theorem, the posterior probability
$P[A|\bar{G}]\propto e^{Q}$ and

$Q=\alpha S-\frac{\chi^{2}}{2}.$ (22)

####  3.1.3 Algorithms of maximum entropy method

Now the original analytical continuation problem becomes how to figure out the
optimal $A(\omega)$ that maximizes $Q$. In other words, we have to solve the
following equation:

$\frac{\partial Q}{\partial A}\bigg{|}_{A=\hat{A}}=0,$ (23)

where $\hat{A}(\omega)$ is the optimal $A(\omega)$. Eq. (23) can be easily
solved by using standard Newton method. However, the obtained
$\hat{A}(\omega)$ is $\alpha$-dependent. That is to say, for a given $\alpha$,
there is always a $\hat{A}(\omega)$ that satisfies Eq. (23). So, new problem
arises because we have to figure out a way to construct the final spectral
function from these $\alpha$-resolved $\hat{A}(\omega)$. Now there exist four
algorithms, namely “historic”, “classic”, “bryan”, and “$\chi^{2}$kink”. Next
we will introduce them one by one.

_Historic algorithm_. The historic algorithm is quite simple. The $\alpha$
parameter will be adjusted iteratively to meet the following criterion:

$\chi^{2}=N,$ (24)

where $N$ is the number of mesh points for spectral density $A(\omega)$.

_Classic algorithm_. The basic equation for the classic algorithm reads:

$-2\alpha S(A_{\alpha})=\text{Tr}\left[\frac{\Lambda(A_{\alpha})}{\alpha
I+\Lambda(A_{\alpha})}\right],$ (25)

where $I$ is an identity matrix. The elements of $\Lambda$ matrix are
calculated as follows:

$\Lambda_{ij}=\sqrt{A_{i}}\left(\sum_{kl}K_{ki}[C^{-1}]_{kl}K_{lj}\right)\sqrt{A_{j}},$
(26)

where $C$ is the covariance matrix. Eq. (25) will be iteratively solved until
the optimal $\alpha$ and $\hat{A}(\omega)$ are determined.

_Bryan algorithm_. In both historic and classic algorithms, the spectral
function $\hat{A}(\omega)$ is always related to an optimal $\alpha$ parameter.
However, the spirit of the bryan algorithm [28] is completely different. It
tries to generate a series of $\alpha$ parameters and yield the corresponding
$A_{\alpha}(\omega)$. Then the final spectral function $A(\omega)$ is obtained
by evaluating the following integration:

$\overline{A(\omega)}=\int d\alpha~{}A_{\alpha}(\omega)P[\alpha|\bar{G}].$
(27)

_$\chi^{2}$ kink algorithm_. This algorithm was proposed by Bergeron and
Tremblay [22] recently. The first step is to generate a series of $\alpha$
parameters, and evaluate the corresponding spectral functions
$A_{\alpha}(\omega)$ and the goodness-of-fit functions $\chi^{2}[A_{\alpha}]$.
Then we plot $\log_{10}(\chi^{2})$ as a function of $\log_{10}(\alpha)$.
Usually this plot is split into three different regions: (1) Default model
region. In the limit of $\alpha\to\infty$, $\chi^{2}$ goes to a constant high
value. It means that the likelihood function $e^{-\chi^{2}/2}$ has negligible
weight, such that the prior probability $e^{\alpha S}$ becomes dominant and
minimizes $Q[A]$. At that time, the calculated $A(\omega)$ resembles the
default model function $m(\omega)$. (2) Noise-fitting region. In the limit of
$\alpha\to 0$, $\chi^{2}$ is relatively flat and approaches its global
minimum. In this region, the minimization algorithm tends to fit the noise in
$G(\tau)$. (3) Information-fitting region. $\alpha S$ is comparable with
$\chi^{2}/2$, so that $\chi^{2}$ is strongly dependent on $\alpha$. Bergeron
_et al._ suggested that the optimal $\alpha$ parameter situates in the
crossover between noise-fitting region and information-fitting region [22]. So
the second derivative of $\chi^{2}$ with respect to $\alpha$ is calculated,
and the maximum value in the resulting curve indicates the optimal value of
$\alpha$. Quite recently, Kaufmann and Held proposed a more numerically stable
and flexible approach to compute the optimal $\alpha$ [23]. They use the
following empirical function to fit dataset
$\\{\log_{10}(\alpha),\log_{10}(\chi^{2})\\}$:

$\phi(x;a,b,c,d)=a+\frac{b}{1+e^{-d(x-c)}},$ (28)

where $a$, $b$, $c$, and $d$ are fitting parameters. Then the optimal $\alpha$
is approximated by $10^{c-f/d}$, where $f$ is a numerical constant (Its
favorite value lies in $[2,2.5]$).

###  3.2 Stochastic analytical continuation

In principle, for given Green’s function $G$, there exists infinitely many
spectral densities $A(\omega)$ that can be used to reconstruct $G$ via Eq.
(12) and Eq. (13). The maximum entropy method tries to pick up the most likely
spectral function which maximizes $P[A|\bar{G}]$ (It actually maximizes $Q$)
[11, 12]. Here, we would like to introduce an alternative approach, namely the
stochastic analytical continuation [13, 14, 29, 30, 31, 32, 33, 34]. It is
argued that the weights for all the possible spectral densities are the same
if they can give rise to the same $\chi^{2}$. At first, a sequence of spectral
densities will be generated by stochastic method. Then an unbiased thermal
average of all possible spectra, Boltzmann weighted according to goodness-of-
fit function $\chi^{2}$, produces an average spectrum. Thus sometimes the
method was named as average spectrum method or stochastic sampling method in
the references [35, 36, 37, 38]. There are several variants for the stochastic
analytical continuation. Next we will introduce two representative algorithms
as proposed by A. W. Sandvik [13] and K. S. D. Beach [29], respectively.

####  3.2.1 Sandvik’s algorithm

Figure 1: Typical Monte Carlo field configurations for the stochastic
analytical continuation (A. W. Sandvik’s algorithm) [13]. Here, the $\delta$
functions reside at unrestricted frequencies $\\{\omega_{i}\\}$, but their
amplitudes $\\{\gamma_{i}\\}$ are equal and fixed. Note that different
parameterizations are also possible [34].

It was early on realized that a different way to achieve a smooth spectrum is
to average over many solutions with reasonable $\chi^{2}$ values [33]. Several
years later, A. W. Sandvik introduced the stochastic analytical continuation
in a slightly different form [13]. He suggested that the spectral function
$A(\omega)$ can be parameterized using $N$ $\delta$ functions (Please see
Figure 1 for a schematic diagram):

$A(\omega)=\sum^{N}_{i=1}\gamma_{i}\delta(\omega-\omega_{i}),$ (29)

where $\gamma_{i}$ and $\omega_{i}$ denote the amplitude and position of the
$i$-th $\delta$ function, respectively. Next, the Metropolis important
sampling algorithm is employed to sample the configuration space
$\mathcal{C}=\\{\omega_{i},\gamma_{i}\\}$. In practice, there are two
elementary Monte Carlo updates. One is to change the amplitudes of a pair of
$\delta$ functions under the constraint $\sum_{i}\gamma_{i}=1$. Another one is
to shift position of a randomly chosen $\delta$ function. Of course, block or
global updates can be implemented to improve ergodicity and sampling
efficiency [35, 36].

The transition probability of Monte Carlo updates reads:

$p(\mathcal{C}\to\mathcal{C}^{\prime})=\exp\left(-\frac{\Delta\chi^{2}}{2\Theta}\right),$
(30)

where the goodness-of-fit function $\chi^{2}$ can be evaluated by Eq. (19),
$\Theta$ is a regulation parameter which is similar to the $\alpha$ parameter
appeared in the maximum entropy method. Well, the remaining problem is how to
fix $\Theta$. Sandvik suggested to measure the following entropic term for a
series of $\Theta$:

$S(\Theta)=-\sum^{N}_{i=1}\gamma_{i}\log(\gamma_{i})K(0,\omega_{i}),$ (31)

where $K$ is the kernel function as defined above [13]. Then make a plot of
$S$ with respect to $\log(\Theta^{-1})$. Overall, when $\Theta$ is large, $S$
exhibits large fluctuations. When $\Theta$ is small, $S$ will approach its
global minimum steadily. A sharp drop in $S$ before the approach to a constant
value has been observed, and there is a local maximum at some
$\Theta=\hat{\Theta}$ preceding the drop. Thus, Sandvik postulated that
$\hat{\Theta}$ was the optimum value at which to accumulate and average the
spectral function. Syljuåsen _et al._ [38] suggested that let $\Theta=1$.
Fuchs _et al._ tried to fix $\Theta$ by using Bayesian inference. Such that
their approach was named as stochastic analytical inference [14]. Very
recently, Shao and Sandvik _et al._ proposed a smart method to determine the
optimal value of $\Theta$ [32, 30]. $\Theta$ is adjusted so that

$\langle\chi^{2}(\Theta)\rangle\approx\chi^{2}_{\text{min}}+c\sqrt{2\chi^{2}_{\text{min}}},$
(32)

where $c$ is a constant of order 1, $\chi^{2}_{\text{min}}$ is the minimum
value of $\chi^{2}$ at given $\Theta$. Note that $\chi^{2}_{\text{min}}$ can
be obtained in a simulated annealing process [39] to very low $\Theta$.

####  3.2.2 Beach’s algorithm

Figure 2: Typical Monte Carlo field configurations for the stochastic
analytical continuation (K. S. D. Beach’s algorithm) [29]. Note that the
amplitudes $\\{\gamma_{i}\\}$ of all the $\delta$ functions are not identical.
Both amplitudes $\\{\gamma_{i}\\}$ and positions $\\{r_{i}\\}$
($0.0<r_{i}<1.0$) can be sampled by using Monte Carlo method.

K. S. D. Beach proposed another variant of stochastic analytical continuation
in 2004 [29]. In his proposal, the analytical continuation problem is mapped
into a system of interacting classic fields at first. Then the classic field
is sampled using Monte Carlo method to obtain the final solution. He concluded
that the maximum entropy method is simply the mean field limit of the
stochastic analytical continuation. Next, this algorithm will be explained
concisely.

_Classic fields_. Recalled that the goodness-of-fit functional $\chi^{2}[A]$
measures how closely the Green’s function generated from $A(\omega)$ matches
the raw input data. Its expression is rewritten as follows:

$\chi^{2}[A]=\int^{\beta}_{0}\frac{1}{\sigma(\tau)^{2}}\left|\int
d\omega~{}K(\tau,\omega)A(\omega)-\bar{G}(\tau)\right|^{2}d\tau.$ (33)

At first, a new variable $x$ is introduced. The relation between $x$ and
$\omega$ is:

$x=\phi(\omega)=\int^{\omega}_{-\infty}d\omega^{\prime}~{}m(\omega^{\prime}),$
(34)

where $m(\omega)$ denotes the default model function. Clearly, the
$\phi(\omega)$ function defines a smooth mapping from $\mathbf{R}\to[0,1]$.
Since $\omega=\phi^{-1}(x)$, a dimensionless classic field $n(x)$ is created:

$n(x)=\frac{A(\phi^{-1}(x))}{m(\phi^{-1}(x))}.$ (35)

It is easy to prove that both $n(x)$ and $A(\omega)$ obey similar
normalization condition:

$\int d\omega~{}A(\omega)=\int^{1}_{0}dx~{}n(x)=1.$ (36)

Next, in analogy with the goodness-of-fit functional $\chi^{2}[A]$, the
Hamiltonian for the system of classic field $\\{n(x)\\}$ can be defined as
follows:

$H[n(x)]=\int^{\beta}_{0}\frac{d\tau}{\sigma(\tau)^{2}}\left|\int^{1}_{0}dx~{}K(\tau,x)n(x)-\bar{G}(\tau)\right|.$
(37)

Supposing $\alpha$ is an inverse temperature of the system, then the partition
function $Z$ is:

$Z=\int\mathcal{D}n~{}e^{-\alpha H[n]},$ (38)

where

$\int\mathcal{D}n=\int^{\infty}_{0}\left[\prod_{x}dn(x)\right]\delta\left(\int^{1}_{0}dx~{}n(x)-1\right).$
(39)

The thermally averaged value of the classic field is:

$\langle n(x)\rangle=\frac{1}{Z}\int\mathcal{D}n~{}n(x)e^{-\alpha H[n]}.$ (40)

Finally, according to the definition of the classic field, the averaged
spectral density $\langle A(\omega)\rangle$ can be expressed as:

$\langle A(\omega)\rangle=\langle n(\phi(\omega))\rangle m(\omega).$ (41)

So, by introducing the classic field $\\{n(x)\\}$, the analytical continuation
problem is converted into a statistical sampling of the classic field, which
is easily solved by using Monte Carlo method.

_Monte Carlo sampling_. Next we clarify how to sample the classic field.
Similar to Sandvik’s algorithm [13, 34], $n(x)$ is parameterized as a
superposition of many $\delta$ functions (see Figure 2 for a schematic
diagram):

$n_{\mathcal{C}}(x)=\sum_{i}\gamma_{i}\delta(x-r_{i}),$ (42)

where $\gamma_{i}$ and $r_{i}$ denote amplitude (weight) and position of the
$i$-th $\delta$ function, respectively. And $\mathcal{C}$ means a
configuration space formed by a set of $r_{i}$ and $\gamma_{i}$,

$\mathcal{C}=\\{r_{i},\gamma_{i}\\}.$ (43)

Note that $\gamma_{i}$ and $r_{i}$ satisfy the following constraints:

$\forall i,~{}\gamma_{i}>0,~{}\sum_{i}\gamma_{i}=1,~{}0\leq r_{i}\leq 1.$ (44)

Supposed that there is a transition from $\mathcal{C}$ to
$\mathcal{C}^{\prime}$
($\\{r_{i},\gamma_{i}\\}\to\\{r^{\prime}_{i},\gamma^{\prime}_{i}\\}$):

$r_{i}\to r^{\prime}_{i}=r_{i}+\sum_{\lambda\in\Lambda}\delta_{i\lambda}\Delta
r_{\lambda},$ (45)
$\gamma_{i}\to\gamma^{\prime}_{i}=\gamma_{i}+\sum_{\lambda\in\Lambda}\delta_{i\lambda}\Delta\gamma_{\lambda},$
(46)

where $\Lambda$ means a subset of the $\delta$ functions, then the Hamiltonian
of the system is changed from $H_{\mathcal{C}}$ to $H_{\mathcal{C}^{\prime}}$.
According to Eq. (37), $H_{\mathcal{C}}$, $H_{\mathcal{C}^{\prime}}$, and
their difference $\Delta H$ can be calculated by:

$H_{\mathcal{C}}=\int^{\beta}_{0}d\tau~{}h_{\mathcal{C}}(\tau)^{2},$ (47)
$H_{\mathcal{C}^{\prime}}=\int^{\beta}_{0}d\tau\left[h_{\mathcal{C}}(\tau)+\Delta
h(\tau)\right]^{2},$ (48) $\Delta
H=H_{\mathcal{C}^{\prime}}-H_{\mathcal{C}}=\int^{\beta}_{0}d\tau~{}\Delta
h(\tau)[2h_{\mathcal{C}}(\tau)+\Delta h(\tau)].$ (49)

Here,

$h(\tau)=\frac{1}{\sigma(\tau)}\left[\int^{1}_{0}dx~{}K(\tau,x)n(x)-\bar{G}(\tau)\right],$
(50)

and

$\Delta
h(\tau)=\frac{1}{\sigma(\tau)}\sum_{\lambda\in\Lambda}\left[\gamma^{\prime}_{\lambda}K(\tau,r^{\prime}_{\lambda})-\gamma_{\lambda}K(\tau,r_{\lambda})\right].$
(51)

Finally, the transition probability from $\mathcal{C}$ to
$\mathcal{C}^{\prime}$ reads

$p(C\to C^{\prime})=\exp(-\alpha\Delta H).$ (52)

_Parallel tempering_. The parallel tempering trick [40] is adopted to improve
the Monte Carlo algorithm as described above. It is possible to proceed
multiple simulations simultaneously for a sequence of inverse temperature
parameters $\\{\alpha_{1},\alpha_{2},\cdots,\alpha_{N}\\}$. The ratio for two
adjacent $\alpha$ parameters is a constant: $\alpha_{p+1}/\alpha_{p}=R$. Note
that the field configurations in all simulations evolve in parallel but not
independently. We can swap the field configurations between two adjacent
layers. Of course, the detailed balance is always preserved, and each
simulation will eventually settle into thermal equilibrium at given $\alpha$.
The transition probability of such a global Monte Carlo update is:

$p(\mathcal{C}\to\mathcal{C}^{\prime})=\exp[(\alpha_{p}-\alpha_{q})(H_{p}-H_{q})],$
(53)

where $p$ and $q$ are layer indices, and $p=q\pm 1$. Parallel tempering
eliminates the need for an initial annealing stage. Another advantage of
parallel tempering is that it yields a complete temperature profile of all the
important thermodynamic variables (such as specific heat and internal energy),
which can be used to estimate the critical $\alpha$ and the final spectral
function $\langle A(\omega)\rangle$.

_Critical inverse temperature_. Clearly, $\langle n(x)\rangle$ strongly
depends on the inverse temperature $\alpha$. How to use these
$\alpha$-dependent $\langle n(x)\rangle$ to construct the final spectral
function? Beach suggested a novel method [29]. During parallel tempering
process, the internal energy of the system is also measured in addition to
$\langle n(x)\rangle$:

$U(\alpha_{p})=\langle H[n]\rangle_{\alpha_{p}}.$ (54)

Let us plot $\log_{10}[U(\alpha)]$ as a function of $\log_{10}(\alpha)$. We
find that $\log_{10}[U(\alpha)]$ drops quickly at first when
$\log_{10}(\alpha)$ increases, and then it approaches to a constant value
slowly. The knee in $\log_{10}[U(\alpha)]$ function, occurring in the vicinity
of $\alpha=\alpha^{*}$ (the corresponding layer index $p=p^{*}$), signals a
jump in specific heat (a thermodynamic phase transition). Then the averaged
spectral function is constructed by:

$\langle\langle
n(x)\rangle\rangle=\frac{\sum^{N-1}_{p=p*}[U(\alpha_{p})-U(\alpha_{p+1})]\langle
n(x)\rangle_{\alpha_{p}}}{U(\alpha_{p*})-U(\alpha_{N})},$ (55)

where $N$ is the total number of $\alpha$, and $\alpha_{p*}$
($\equiv\alpha^{*}$) is the critical inverse temperature.

_Likelihood function_. Neither of the Sandvik’s and Beach’s algorithms needs
extra entropic term to regulate the spectral densities [13, 29]. All the
stochastically generated spectra are treated on the same footing. Thus, the
calculated spectral function retains more subtle structures than that obtained
by the maximum entropy method. Actually, in the stochastic analytical
continuation,

$\langle A\rangle=\int\mathcal{D}A~{}P[A|\bar{G}]A.$ (56)

The weight of the candidate spectral function $A$ is given by the likelihood
function $P[A|\bar{G}]$. Eq. (30) and Eq. (52) can be viewed as likelihood
functions in the stochastic analytical continuation.

###  3.3 Stochastic optimization method

Figure 3: Typical Monte Carlo field configurations for the stochastic
optimization method [15]. The spectral function is parameterized by multiple
rectangle functions. Here, $c_{i}$, $w_{i}$, and $h_{i}$ denote the center,
width, and height of the $i$-th rectangle, respectively.

A. O. Mishchenko _et al._ [15] proposed the stochastic optimization method.
Though it looks like the stochastic analytical continuation [13, 29], their
differences are quite apparent. The stochastic optimization method does not
need any likelihood function or Boltzmann distribution to weight the candidate
spectral functions. It generates a lot of spectral functions through Monte
Carlo samplings. For each candidate spectral function, the deviation $D$
between the reconstructed Green’s function $\tilde{G}$ and original Green’s
function $\bar{G}$ is measured. Those spectral functions with small deviations
$D$ are selected and averaged. Such that the desired spectral function is
obtained.

_Deviation function_. In the stochastic optimization method, the deviation
between reconstructed data $\tilde{G}$ and input data $\bar{G}$ is described
by:

$D[A]=\sum^{M}_{m=1}|\Delta(m)|,$ (57)

where $M$ is the number of input data, and $\Delta(m)$ is the deviation
function,

$\Delta(m)=\frac{\bar{G}(m)-\tilde{G}(m)}{S(m)}.$ (58)

Here, $S(m)=|G(m)|^{d}$ (where $0\leq d\leq 1$). Recently, Krivenko _et al._
suggested that it would be better to use the goodness-of-fit functional
$\chi^{2}[A]$ to replace $D[A]$ [24, 25].

_Spectral density_. The stochastic optimization method will try to accumulate
the candidate spectral functions that manifest small $D[A]$. Supposed the
Monte Carlo simulations are repeated for $L$ times. For the $i$-th Monte Carlo
simulation, the spectral density $A_{i}(\omega)$ and deviation $D[A_{i}]$ are
recorded. The minimum value of deviation is $\min\\{D[A_{i}]\\}$. Thus, the
final spectral density reads:

$A(\omega)=\frac{1}{L_{\text{good}}}\sum^{L}_{i=1}\theta(\alpha_{\text{good}}\min\\{D[A_{i}]\\}-D[A_{i}])A_{i}(\omega).$
(59)

Here, $\theta(x)$ is the Heaviside step function, and $\alpha_{\text{good}}$
is a adjustable parameter. $L_{\text{good}}$ denotes the number of “good”
spectral functions:

$L_{\text{good}}=\sum^{L}_{i=1}\theta(\alpha_{\text{good}}\min\\{D[A_{i}]\\}-D[A_{i}]).$
(60)

That is to say, only those spectral functions who satisfy the following
condition will be selected:

$D[A_{i}]\leq\alpha_{\text{good}}\min\\{D[A_{i}]\\}.$ (61)

Clearly, the larger $\alpha_{\text{good}}$ is, the more spectral functions are
included. It is usually set to 2.

_Rectangle representation_. Similar to the stochastic analytical continuation
[13, 29], the stochastic optimization method usually employs a few rectangle
functions to parameterize the spectral function:

$A(\omega)=\sum_{i}R_{\\{c_{i},w_{i},h_{i}\\}}(\omega),$ (62)

where $i$ is the index of rectangle function. The definition of rectangle
function $R_{\\{c_{i},w_{i},h_{i}\\}}(\omega)$ reads:

$R_{\\{c_{i},w_{i},h_{i}\\}}(\omega)=h_{i}\theta[\omega-(c_{i}-w_{i}/2)]\theta[(c_{i}+w_{i}/2)-\omega],$
(63)

where $c_{i}$, $w_{i}$, $h_{i}$ denote the center, width, and height of the
$i$-th rectangle, respectively. Pay attention to that the area of all
rectangles must be normalized to 1:

$\sum_{i}h_{i}w_{i}=1.$ (64)

_Monte Carlo sampling_. The parameters of all rectangle functions create a
configuration space:

$\mathcal{C}=\\{c_{i},w_{i},h_{i}\\}.$ (65)

Then the Metropolis algorithm is utilized to sample this configuration space.
Mishchenko _et al._ introduces seven Monte Carlo updates [15, 24], including:
(a) Insert a new rectangle, change width and height of another rectangle; (b)
Remove an existing rectangle, change width and height of another rectangle;
(c) Shift position of any rectangles; (d) Change widths of any two rectangles;
(e) Change heights of any two rectangles; (f) Split a rectangle into two new
rectangles; (g) Merge two adjacent rectangles into a new rectangle. The
transition probability of these Monte Carlo updates is:

$p(\mathcal{C}\to\mathcal{C}^{\prime})=\left(\frac{D[A_{\mathcal{C}}]}{D[A_{\mathcal{C}^{\prime}}]}\right)^{1+d}$
(66)

As compared to the maximum entropy method [11, 12], the likelihood function,
entropic term, and model function are absent in the stochastic optimization
method. As compared to the stochastic analytical continuation [13, 29], there
are no adjustable parameters, such as $\Theta$ in Sandvik’s algorithm and
$\alpha$ in Beach’s algorithm. Thus, the simulated results of the stochastic
optimization method are less affected by artificial parameters.

##  4 Overview

###  4.1 Major features

Now the ACFlow toolkit supports three analytical continuation methods as
introduced above. It includes four different analytical continuation solvers,
namely MaxEnt, StochAC, StochSK, and StochOM. Just as their names suggested,
the MaxEnt solver implements the maximum entropy method [12]. The StochAC and
StochSK solvers implement the K. S. D. Beach’s algorithm [29] and A. W.
Sandvik’s algorithm [13] of the stochastic analytical continuation,
respectively. The StochOM solver implements the stochastic optimization method
[15]. The ACFlow toolkit also provides a convenient library, which can be used
to prepare and carry out analytical continuation calculations flexibly. The
major features of the ACFlow toolkit are summarized in Table 1.

Features | MaxEnt | StochAC | StochSK | StochOM
---|---|---|---|---
Matrix-valued Green’s function | Y | N | N | N
Imaginary time grid | Y | Y | Y | Y
Matsubara frequency grid | Y | Y | Y | Y
Linear mesh | Y | Y | Y | Y
Nonlinear mesh | Y | Y | Y | Y
Fermionic kernel | Y | Y | Y | Y
Bosonic kernel | Y | Y | Y | Y
Self-defined model function | Y | N | N | N
Constrained analytical continuation | N | Y | Y | Y
Regeneration of input data | Y | Y | Y | Y
Kramers-Kronig transformation | Y | Y | Y | Y
Parallel computing | N | Y | Y | Y
Parallel tempering | N | Y | N | N
Interactive mode | Y | Y | Y | Y
Script mode | Y | Y | Y | Y
Standard mode | Y | Y | Y | Y
Table 1: Major features of the ACFlow toolkit. MaxEnt, StochAC, StochSK, and
StochOM are the four analytical continuation solvers as implemented in this
toolkit.

In Table 1, “Y” means yes while “N” means no. “Interactive mode”, “Script
mode”, and “Standard model” are the three running modes supported by the
ACFlow toolkit. We will introduce them in next section. The MaxEnt solver
supports the “historic”, “classic”, “bryan”, and “chi2kink” algorithms to
determine the $\alpha$ parameter. The StochAC solver is only compatible with a
flat model function, while the StochSK and StochOM solvers don’t rely on any
default model functions. The StochOM solver does not support analytical
continuation of fermionic imaginary time Green’s function for the moment.

###  4.2 Implementations

The ACFlow toolkit is developed with pure Julia language. Thanks to powerful
type system and multiple dispatch paradigm of the Julia language, the four
different analytical continuation solvers are integrated into an united
software architecture. Redundant codes are greatly reduced. It is quite easy
to implement new analytical continuation solver or add new features to the
existing solvers in the future. Distributed computing is a built-in feature of
Julia. So, it is straightforward to realize parallel calculations in the
ACFlow toolkit. Now except for the MaxEnt solver, all the other solvers are
parallelized.

Filename | Description
---|---
ACFlow.jl | Entry of the ACFlow module.
maxent.jl | Maximum entropy method.
sac.jl | Stochastic analytical continuation (K. S. D. Beach’s algorithm).
san.jl | Stochastic analytical continuation (A. W. Sandvik’s algorithm).
som.jl | Stochastic optimization method.
global.jl | Numerical and physical constants.
types.jl | Basic data structures and computational parameters.
base.jl | Driver for analytical continuation simulation.
inout.jl | Read input data and write calculated results.
config.jl | Parse configuration file and extract computational parameters.
math.jl | Root finding, numerical integration, interpolation, Einstein summation, and curve fitting.
util.jl | Some utility functions.
mesh.jl | Meshes for spectral density.
grid.jl | Grids for input data.
model.jl | Default model functions.
kernel.jl | Kernel functions.
Table 2: List of source codes of the ACFlow toolkit.

The source codes of the ACFlow toolkit are placed in the acflow/src folder.
Their functions are summarized in Table 2. The documentation of the ACFlow
toolkit is written by using the Markdown language and Documenter.jl package.
The source codes are placed in the acflow/docs folder. The users can build
documentation by themselves. Please see section 5 for how to do that. Or they
can read the latest documentation in the following website:

        https://huangli712.github.io/projects/acflow/index.html

Ten tests and four tutorials are also shipped with the ACFlow toolkit. Their
source codes are placed in the acflow/test folder. See acflow/test/test.md and
acflow/test/tutor.md for more details. The code repository of the ACFlow
toolkit is:

        https://github.com/huangli712/ACFlow

##  5 Getting started

In this section, we will discuss how to install and use the ACFlow toolkit.

###  5.1 Installation

It is an easy task to install the ACFlow toolkit. First, since it is written
in pure Julia language, it is necessary to install the Julia runtime
environment at first. The newest version of Julia is always preferred (version
$>$ 1.60). Since the core codes only rely on Julia’s built-in standard
library, no the third-party packages are needed. Second, just download source
codes of the ACFlow toolkit from its github repository. It should be a
compressed file, such as acflow.zip or acflow.tar.gz. Please uncompress it in
your favorite directory by using the following commands:

        $ unzip acflow.zip

or

        $ tar xvfz acflow.tar.gz

Third, the users have to declare a new environment variable ACFLOW_HOME.
Supposed that the root directory of the ACFLow toolkit is
/home/your_home/acflow, then ACFLOW_HOME should be setup as follows:

        $ export ACFLOW_HOME=/home/your_home/acflow/src

Finally, in order to generate the documentation, the users should type the
following commands in the terminal:

        $ pwd
        /home/your_home/acflow
        $ cd docs
        $ julia make.jl

After a few seconds, the documentation is built and saved in the
acflow/docs/build directory if everything is OK. The home page of the
documentation is acflow/docs/build/index.html. We can open it with any web
browsers.

###  5.2 Run

The ACFlow toolkit is designed to be flexible and easy-to-use. It provides
three running modes to facilitate analytical continuation calculations, namely
the interactive, script, and standard modes.

_Interactive mode_. With the ACFlow toolkit, the users can setup and carry out
analytical continuation simulations interactively in Julia’s REPL (Read-Eval-
Print Loop) environment. For example,

        julia> push!(LOAD_PATH, ENV["ACFLOW_HOME"])
        julia> using ACFlow
        julia> setup_args("ac.toml")
        julia> read_param()
        julia> mesh, Aout, Gout = solve(read_data())

Here, ac.toml is a configuration file, which contains essential computational
parameters. The return values of the solve() function (i.e., mesh, Aout, and
Gout) are mesh at real axis $\omega$, spectral density $A(\omega)$, and
reproduced Green’s function $\tilde{G}$, respectively. They can be further
analyzed or visualized by the users.

_Script mode_. The core functionalities of the ACFlow toolkit are exposed to
the users via a simple application programming interface. So, the users can
write Julia scripts easily by themselves to perform analytical continuation
simulations. A minimal Julia script (acrun.jl) is listed as follows:

        #!/usr/bin/env julia
        push!(LOAD_PATH, ENV["ACFLOW_HOME"])
        using ACFlow
        setup_args("ac.toml")
        read_param()
        mesh, Aout, Gout = solve(read_data())

Of course, this script can be extended to finish complex tasks. In section
6.1, a realistic example is provided to show how to complete an analytical
continuation of Matsubara self-energy function via the script mode.

_Standard mode_. In the standard mode, the users have to prepare the input
data manually. In addition, a configuration file must be provided. Supposed
that the configuration file is ac.toml, then the analytical continuation
calculation is launched as follows:

        $ /home/your_home/acflow/util/acrun.jl ac.toml

or

        $ /home/your_home/acflow/util/Pacrun.jl ac.toml

Noted that the acrun.jl script runs sequentially, while the Pacrun.jl script
supports parallel and distributed computing. As we can conclude from the
filename extension of configuration file (ac.toml), it adopts the TOML
specification. The users may edit it with any text-based editors. Next we will
introduce syntax and format of the input data files and configuration files.

###  5.3 Input files

The input files for the ACFlow toolkit can be divided into two groups: data
files and configuration files.

_Data files_. The input data should be store in CSV-like text files. For
imaginary time Green’s function, the data file should contain three columns.
They represent $\tau$, $\bar{G}(\tau)$, and standard deviation of
$\bar{G}(\tau)$. For fermionic Matsubara Green’s function, the data file
should contain five columns. They represent $\omega_{n}$, Re$G(i\omega_{n})$,
Im$G(i\omega_{n})$, standard deviation of Re$G(i\omega_{n})$, and standard
deviation of Im$G(i\omega_{n})$. For bosonic correlation function
$\chi(i\omega_{n})$, the data file should contain four columns. They represent
$\omega_{n}$, Re$\chi(i\omega_{n})$, and standard deviation of
Re$\chi(i\omega_{n})$.

_Configuration files_. The configuration file adopts the TOML format. It is
used to setup the computational parameters. It consists of one or more blocks.
Possible blocks (or sections) of the configuration file include [BASE],
[MaxEnt], [StochAC], [StochSK], and [StochOM]. The [BASE] block is mandatory,
while the other blocks are optional. A schematic configuration file (ac.toml)
is listed as follows:

⬇

[BASE]

finput = "giw.data"

solver = "StochOM"

...

[MaxEnt]

method = "chi2kink"

...

[StochAC]

nfine = 10000

...

[StochSK]

method = "chi2min"

...

[StochOM]

ntry = 100000

...

In the [BASE] block, the analytical continuation problem is defined. The
solver used to solve the problem must be assigned. The types of mesh, grid,
default model function, and kernel function are also determined. The [MaxEnt],
[StochAC], [StochSK], and [StochOM] blocks are used to customize the
corresponding analytical continuation solvers further. In Table 3-Table 6, all
the possible input parameters for these blocks are collected and summarized.
As for detailed explanations of these parameters, please refer to the user
guide of the ACFlow toolkit. The uses can find it in the acflow/docs
directory.

[BASE] block
---
Parameter | Type | Default | Description
finput | string | “green.data” | Filename for input data.
solver | string | “MaxEnt” | Solver for the analytical continuation problem.
ktype | string | “fermi” | Type of kernel function.
mtype | string | “flat” | Type of default model function.
grid | string | “ffreq” | Grid for input data (imaginary axis).
mesh | string | “linear” | Mesh for output data (real axis).
ngrid | integer | 10 | Number of grid points.
nmesh | integer | 501 | Number of mesh points.
wmax | float | 5.0 | Right boundary (maximum value) of mesh.
wmin | float | -5.0 | Left boundary (minimum value) of mesh.
beta | float | 10.0 | Inverse temperature.
offdiag | bool | false | Treat the off-diagonal part of matrix-valued function?
pmodel | array | N/A | Additional parameters for customizing the default model.
pmesh | array | N/A | Additional parameters for customizing the mesh.
exclude | array | N/A | Restriction of energy range of the spectrum.
Table 3: Possible parameters for the [BASE] block. [MaxEnt] block
---
Parameter | Type | Default | Description
method | string | “chi2kink” | How to determine the optimized $\alpha$ parameter?
nalph | integer | 12 | Total number of the chosen $\alpha$ parameters.
alpha | float | 1e9 | Starting value for the $\alpha$ parameter.
ratio | float | 10.0 | Scaling factor for the $\alpha$ parameter.
blur | float | -1.0 | Shall we preblur the kernel and spectrum?
Table 4: Possible input parameters for the [MaxEnt] block, which are used to
setup the solver based on the maximum entropy method [12, 11].  [StochAC]
block
---
Parameter | Type | Default | Description
nfine | integer | 10000 | Number of points of a very fine linear mesh.
ngamm | integer | 512 | Number of $\delta$ functions.
nwarm | integer | 4000 | Number of Monte Carlo thermalization steps.
nstep | integer | 4000000 | Number of Monte Carlo sweeping steps.
ndump | integer | 40000 | Intervals for monitoring Monte Carlo sweeps.
nalph | integer | 20 | Total number of the chosen $\alpha$ parameters.
alpha | float | 1.0 | Starting value for the $\alpha$ parameter.
ratio | float | 1.2 | Scaling factor for the $\alpha$ parameter.
[StochSK] block
Parameter | Type | Default | Description
method | string | “chi2min” | How to determine the optimized $\Theta$ parameter?
nfine | integer | 100000 | Number of points of a very fine linear mesh.
ngamm | integer | 1000 | Number of $\delta$ functions.
nwarm | integer | 1000 | Number of Monte Carlo thermalization steps.
nstep | integer | 20000 | Number of Monte Carlo sweeping steps.
ndump | integer | 200 | Intervals for monitoring Monte Carlo sweeps.
retry | integer | 10 | How often to recalculate the goodness-of-fit function.
theta | float | 1e6 | Starting value for the $\Theta$ parameter.
ratio | float | 0.9 | Scaling factor for the $\Theta$ parameter.
Table 5: Possible input parameters for the [StochAC] and [StochSK] blocks,
which are used to setup the two solvers based on the stochastic analytical
continuation (Beach’s and Sandvik’s algorithms) [13, 29].  [StochOM] block
---
Parameter | Type | Default | Description
ntry | integer | 2000 | Number of attempts to figure out the solution.
nstep | integer | 1000 | Number of Monte Carlo steps per try.
nbox | integer | 100 | Number of boxes to construct the spectrum.
sbox | float | 0.005 | Minimum area of the randomly generated rectangles.
wbox | float | 0.02 | Minimum width of the randomly generated rectangles.
norm | float | -1.0 | Is the norm calculated?
Table 6: Possible input parameters for the [StochOM] block, which are used to
setup the solver based on the stochastic optimization method [15].

###  5.4 Output files

Once the analytical continuation simulation is finished, the final spectral
function $A(\omega)$ is outputted to Aout.data. As is shown in Eq. (7),
$A(\omega)$ is equivalent to the imaginary part of real frequency Green’s
function Im$G(\omega)$. Then the ACFlow toolkit will automatically calculate
the corresponding real part Re$G(\omega)$ via the Kramers-Kronig
transformation [see Eq. (8)]. The full Green’s function at real axis
$G(\omega)$ is stored in Gout.data. The spectral function is also used to
reconstruct the imaginary time or Matsubara Green’s functions
[$\tilde{G}(\tau)$ or $\tilde{G}(i\omega_{n})$], which is stored in repr.data.
Besides the three output files, the ACFlow toolkit will generate quite a few
output files, which can be used to analyze and diagnose the calculated
results. All of the possible output files of the ACFlow toolkit are collected
and explained in Table 7.

Filename | Description
---|---
Aout.data | Final spectral function $A(\omega)$.
Gout.data | Full Green’s function at real axis $G(\omega)$.
repr.data | Reproduced Green’s function $\tilde{G}$ at imaginary time or frequency axis.
model.data | Default model function $m(\omega)$.
chi2.data |  $\log_{10}(\chi^{2})$ vs $\log_{10}(\alpha)$.
prob.data |  $P[\alpha|\bar{G}]$ vs $\alpha$ for the MaxEnt solver (bryan algorithm).
Aout.data.alpha_$i$ |  $\alpha$-resolved spectral function $A_{\alpha}(\omega)$ for the StochAC solver.
hamil.data |  $U(\alpha)$ vs $\alpha$ for the StochAC solver.
goodness.dat |  $\log_{10}(\chi^{2})$ vs $\log_{10}(\Theta)$ for the StochSK solver.
stat.data | Monte Carlo statistical information for stochastic sampling method.
Table 7: Possible output files of the ACFlow toolkit.

##  6 Examples

In order to demonstrate usefulness of the ACFlow toolkit, four examples are
illustrated in this section. These examples cover typical application
scenarios of the ACFlow toolkit, including analytical continuations of
Matsubara self-energy function, Matsubara Green’s function, imaginary time
Green’s function, and current-current correlation function within the script
mode or standard mode. All of the necessary source codes and data files, which
can be used to reproduce the results as shown in this section, are placed in
the /home/your_home/acflow/test/T* folders.

###  6.1 Matsubara self-energy function

Now let us consider the following single-band Hubbard model on a Bethe lattice
at first:

$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},$
(67)

where $t$ is the hopping parameter, $\mu$ is the chemical potential, $U$ is
the Coulomb interaction, $n$ is the occupation number, $\sigma$ denotes the
spin, $i$ and $j$ are site indices. This model is solved by using the
dynamical mean-field theory (dubbed DMFT) [41] with the hybridization
expansion continuous-time quantum Monte Carlo solver (dubbed CT-HYB) [3] as
implemented in the $i$QIST package [42, 43]. The parameters used in the DMFT +
CT-HYB calculation are $t=0.5$, $U=2.0$, $\mu=1.0$, and $\beta=10.0$. Once the
DMFT self-consistent calculation is finished, the Matsubara self-energy
function $\Sigma(i\omega_{n})$ is obtained. We are going to convert it to real
frequency self-energy function $\Sigma(\omega)$. The data of Matsubara self-
energy function $\Sigma(i\omega_{n})$ have been preprocessed and stored in
siw.data. This file contains five columns, which are used to record the
Matsubara frequency $\omega_{n}$, Re$\Sigma(i\omega_{n})$,
Im$\Sigma(i\omega_{n})$, error bar of Re$\Sigma(i\omega_{n})$, error bar of
Im$\Sigma(i\omega_{n})$, respectively. Only the first twenty Matsubara
frequency points are kept, because the high-frequency data are somewhat noisy.

The purpose of this example is to demonstrate usage of the MaxEnt solver and
the script mode of the ACFlow toolkit. Next we will explain the key steps in
detail. As for the complete Julia script, please refer to sigma.jl and
gendata.jl in the /home/your_home/acflow/test/T01/ folder.

First, we have to load the essential Julia packages. Both the DelimitedFiles
and Printf packages belong to Julia’s standard library. They are used to read
input data and write calculated results, respectively.

⬇

#!/usr/bin/env julia

push!(LOAD_PATH, ENV["ACFLOW_HOME"])

using DelimitedFiles

using Printf

using ACFlow

welcome() # Print welcome message only

Next, the data of Matsubara self-energy function are read from siw.data. The
Hartree term $\Sigma_{H}$ should be subtracted from its real part:

$\Sigma(i\omega_{n})\to\Sigma(i\omega_{n})-\Sigma_{H}.$ (68)

Note that $\Sigma_{H}$ is approximately equal to the asymptotic value of real
part of $\Sigma(i\omega_{n})$ when $n$ goes to infinite.

⬇

# Deal with self-energy function

#

# Read self-energy function

dlm = readdlm("siw.data")

#

# Get grid

grid = dlm[:,1]

#

# Get self-energy function

Sinp = dlm[:,2] + im * dlm[:,3] # Value

Serr = dlm[:,4] + im * dlm[:,5] # Error bar

#

# Subtract hartree term

Sh = 1.0

@. Sinp = Sinp - Sh

Next, the computational parameters are encapsulated into two dictionaries. The
dictionary B is for the [BASE] block, while the dictionary S is for the MaxEnt
solver. Then the setup_param() function is called, so that these parameters
take effect. Here, the MatEnt solver [12, 11] is employed to tackle the
analytical continuation problem. But the other stochastic sampling solvers are
also applicable. The default model function is gaussian. The mesh for spectral
density is non-uniform (A tangent mesh). The number of used $\alpha$
parameters is 15, and the optimal $\alpha$ parameter is determined by the
$\chi^{2}$kink algorithm [22].

⬇

# Setup parameters

#

# For [BASE] block

# See types.jl/_PBASE for default setup

B = Dict{String,Any}(

"solver" => "MaxEnt", # Choose MaxEnt solver

"mtype" => "gauss", # Default model function

"mesh" => "tangent", # Mesh for spectral density

"ngrid" => 20, # Number of input points

"nmesh" => 801, # Number of output points

"wmax" => 8.0, # Right boundary of mesh

"wmin" => -8.0, # Left boundary of mesh

"beta" => 10.0, # Inverse temperature

)

#

# For [MaxEnt] block

# See types.jl/_PMaxEnt for default setup

S = Dict{String,Any}(

"nalph" => 15, # Number of alpha

"alpha" => 1e12, # Starting value of alpha

"blur" => -1.0, # Enable preblur or not

)

#

# Let the parameters take effect

setup_param(B, S)

It is quite easy to start the analytical continuation calculation. Just call
the solve() function and pass the grid, input data, and error bar data to it.
The return values of this function call are real frequency mesh, spectral
density, and reconstructed Matsubara self-energy function.

⬇

# Call the solver

mesh, Aout, Sout = solve(grid, Sinp, Serr)

Finally, the real frequency self-energy function must be supplemented with the
Hartree term, and then the final results are written into sigma.data.

⬇

# Calculate final self-energy function on real axis

#

# Add hartree term

@. Sout = Sout + Sh

#

# Write self-energy function to sigma.data

open("sigma.data", "w") do fout

for i in eachindex(mesh)

z = Sout[i]

@printf(fout, "%20.16f␣%20.16f␣%20.16f\n",

mesh[i], real(z), imag(z))

end

end

Figure 4: Analytical continuation of Matsubara self-energy function by using
the maximum entropy method. (a) Real part of real frequency self-energy
function. (b) Imaginary part of real frequency self-energy function. (c)
$\chi^{2}$ as a function of $\alpha$. The vertical bar indicates the optimal
$\alpha$ parameter chosen by the $\chi^{2}$kink algorithm. (d) Reproduced and
original data for imaginary part of the Matsubara self-energy functions.

The calculated results are displayed in Fig. 4. Fig. 4(a) and (b) show the
real and imaginary parts of the real frequency self-energy function,
respectively. Near the Fermi level, Re$\Sigma(\omega)$ exhibits quasi-linear
behavior, with which the quasiparticle weight $Z$ and effective mass of
electron $m^{*}$ can be easily evaluated [41]. As for the imaginary part,
Im$\Sigma(0)$ is finite, which indicates that the electron-electron scattering
is not trivial. Fig. 4(c) shows the $\alpha$-dependent $\chi^{2}$. The
vertical bar in this figure indicates the optimal $\alpha$ is around 102.154.
In Fig. 4(d), the reproduced and raw Matsubara self-energy functions are
compared. It is apparent that they are consistent with each other.

###  6.2 Matsubara Green’s function

The purpose of the second example is to treat the Matsubara Green’s function
by using the StochOM solver.

At first, please consider the following spectral density with two gaussian
peaks:

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

with $A_{1}=1.0$, $A_{2}=0.3$, $\epsilon_{1}=0.5$, $\epsilon_{2}=-2.5$,
$\Gamma_{1}=0.2$, and $\Gamma_{2}=0.8$. Then the Matsubara Green’s function
$G(i\omega_{n})$ is evaluated by using Eq. (13) with $\beta=10.0$. Random
noises, generated by the formula $0.0001r_{1}\exp(i2\pi r_{2})$ where $r_{1}$
and $r_{2}$ are pseudo random numbers in ($0.0<r_{1}$, $r_{2}<1.0$), are added
to $G(i\omega_{n})$. The error bar of $G(i\omega_{n})$ is fixed to 1e-4. The
generated data for $G(i\omega_{n})$ are written in giw.data.

Next, we are going to use the standard mode, such that a configure file
(ac.toml) must be prepared. It is listed as follows. Since the StochOM solver
is chosen, the [BASE] and [StochOM] blocks must be present.

⬇

[BASE]

finput = "giw.data"

solver = "StochOM"

ktype = "fermi"

mtype = "flat"

grid = "ffreq"

mesh = "linear"

ngrid = 10

nmesh = 501

wmax = 5.0

wmin = -5.0

beta = 10.0

offdiag = false

[StochOM]

ntry = 100000

nstep = 1000

nbox = 100

sbox = 0.005

wbox = 0.02

norm = -1.0

Figure 5: Analytical continuation of Matsubara Green’s function by using the
stochastic optimization method. (a) Simulated and exact spectral functions.
(b) Reconstructed and synthetic Matsubara Green’s functions. Only the
imaginary parts are presented in this figure.

Then we use the acrun.jl or Pacrun.jl script to perform analytical
continuation simulation. The calculated results are shown in Fig. 5. As is
seen in Fig. 5(a), both the sharp peak around 0.5 eV and the broad peak around
-2.5 eV are correctly reproduced by the StochOM solver. In Fig. 5(b), the
reconstructed Matsubara Green’s function agrees quite well with the raw input
data.

###  6.3 Imaginary time Green’s function

In this example, analytical continuation of imaginary time Green’s function
will be tested. Note that this example is borrowed from Reference [29]
directly.

The exact spectral function reads:

$A(\omega)=\begin{cases}\frac{1}{W}\frac{|\omega|}{\sqrt{\omega^{2}-\Delta^{2}}},~{}\quad&\text{if}~{}\Delta<|\omega|<W/2.\\\
0,&\text{otherwise}.\end{cases}$ (70)

Here, $W$ denotes bandwidth, and $\Delta$ is used to control size of the
energy gap. Let $W=6$ and $2\Delta=1$. This spectrum should exhibit flat
shoulders, steep peaks, and sharp gap edges at the same time. Actually, it is
a typical spectrum of a BCS superconductor.

First, the imaginary time Green’s function $G(\tau)$ is generated using Eq.
(12). Then a normally-distributed random noise is add to $G(\tau)$. Maximum
amplitude of the noise is 1e-4. The error bar of $G(\tau)$ is fixed to 1e-3.
The data are written in gtau.data.

Next, we try to prepare the configure file (ac.toml). In this case, we would
like to benchmark the StochAC solver, so the solver parameter is set to
“StochAC” and the grid parameter is set to “ftime”. Furthermore, the exclude
parameter is enabled to impose some _a priori_ constraints to the spectrum.
The full ac.toml is listed as follows:

⬇

[BASE]

finput = "giw.data"

solver = "MaxEnt"

ktype = "fermi"

mtype = "flat"

grid = "ffreq"

mesh = "linear"

ngrid = 10

nmesh = 501

wmax = 5.0

wmin = -5.0

beta = 10.0

offdiag = false

exclude = [[-5.0,-3.0], [-0.5,0.5], [3.0,5.0]]

[StochAC]

nfine = 10000

ngamm = 512

nwarm = 4000

nstep = 10000000

ndump = 40000

nalph = 40

alpha = 1.00

ratio = 1.20

We perform analytical continuation simulation by using the acrun.jl or
Pacrun.jl script. In order to obtain smooth spectral density, it is useful to
increase number of $\delta$ functions (See ngamm parameter) and number of
Monte Carlo sampling steps (See nstep parameter).

Figure 6: Analytical continuation of imaginary time Green’s function by using
the stochastic analytical continuation (Beach’s algorithm). (a) Simulated and
exact spectral functions. (b) $\alpha$-dependent spectral functions. (c)
Internal energy $U$ as a function of $\alpha$. The vertical bar indicates the
optimal $\alpha$ parameter. (d) Simulated and exact imaginary time Green’s
functions.

Figure 6 shows the calculated results. In Fig. 6(a), the exact spectral
function is compared with the simulated spectrum. Note that besides the
StochAC solver, the other three solvers are also tested. Their results are
also plotted in this figure for a direct comparison. It is remarkable that the
StochAC and StochSK solvers do a superior job of modelling the spectrum. The
major characteristics of the spectrum, including flat regions, steep peaks,
and sharp gap edges, are well captured by the two solvers. Especially, we have
finished more tests without any constraints on the spectral density. The gap
in the spectrum can be reproduced as well. On the other hand, the spectra
obtained by the MaxEnt and StochOM solvers are much too smooth, and show extra
shoulder peaks around $\pm$2.0. Figure 6(b) shows $\alpha$-resolved spectral
functions $A_{\alpha}(\omega)$ for selected $\alpha$ parameters. Fluctuation
in the flat regions of the calculated spectral density grows when $\alpha$
increases. Figure 6(c) shows internal energy $U$ as a function of $\alpha$.
From this figure, the critical $\alpha$ is estimated, which is indicated by
the vertical bar. Finally, the reproduced Green’s function $\tilde{G}(\tau)$
agrees quite well with the raw input data, as is shown in Fig. 6(d).

###  6.4 Current-current correlation function

The previous three examples only involve fermionic correlators. How about
bosonic correlation functions? In this example, we will demonstrate how to
perform analytical continuation simulation for a typical bosonic correlator,
the current-current correlation function $\Pi(\tau)$, to obtain the optical
conductivity $\sigma(\omega)$. Note that this example is taken from Reference
[5] directly.

The exact optical conductivity $\sigma(\omega)$ reads:

$\sigma(\omega)=\left\\{\frac{W_{1}}{1+(\omega/\Gamma_{1})^{2}}+\frac{W_{2}}{1+[(\omega-\epsilon)/\Gamma_{2}]^{2}}+\frac{W_{2}}{1+[(\omega+\epsilon)/\Gamma_{2}]^{2}}\right\\}\frac{1}{1+(\omega/\Gamma_{3})^{6}},$
(71)

where $W_{1}=0.3$, $W_{2}=0.2$, $\Gamma_{1}=0.3$, $\Gamma_{2}=1.2$,
$\Gamma_{3}=4.0$, and $\epsilon=3.0$. The current-current correlation function
$\Pi(\tau)$ can be evaluated from $\sigma(\omega)$ by using the following
equation:

$\Pi(\tau)=\int^{\infty}_{-\infty}K(\tau,\omega)\sigma(\omega)~{}d\omega,$
(72)

where the kernel function $K(\tau,\omega)$ is different from Eq. (14). It
reads:

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

In this case, $\beta$ is fixed to be 20.0.

At first, we use Eq. (71) $\sim$ Eq. (73) to prepare $\Pi(\tau)$. A normally-
distributed random noise is add to $\Pi(\tau)$. Maximum amplitude of the noise
is 1e-4. The error bar of $\Pi(\tau)$ is fixed to 1e-4. The data of
$\Pi(\tau)$ are written in chit.data.

Next, we conduct analytical continuation simulation as usual. The used
configuration file is attached as follows. Here, the StochSK solver is
adopted, so the solver parameter is “StochSK” and the grid parameter is
“btime”. And the Shao-Sandvik algorithm [30] is applied to seek optimal
$\Theta$, so the method parameter is “chi2min”. The users can further increase
the values of nfine, ngamm, and nstep parameters to improve computational
accuracy.

⬇

[BASE]

finput = "chit.data"

solver = "StochSK"

ktype = "bsymm"

mtype = "flat"

grid = "btime"

mesh = "linear"

ngrid = 501

nmesh = 801

wmax = 8.0

wmin = 0.0

beta = 20.0

offdiag = false

[StochSK]

method = "chi2min"

nfine = 40000

ngamm = 1000

nwarm = 1000

nstep = 1000000

ndump = 200

retry = 10

theta = 1e+6

ratio = 0.90

Figure 7: Analytical continuation of current-current correlation function by
using the stochastic analytical continuation (Sandvik’s algorithm). (a)
Simulated and exact optical conductivities $\sigma(\omega)$. (b) Simulated and
exact current-current correlation functions $\Pi(\tau)$.

The calculated results are illustrated in Fig. 7. From Fig. 7(a), it is clear
that the main features of optical conductivity are successfully captured by
the StochSK solver. Both the sharp Drude peak at $\omega=0$ and a broad
satellite peak around $\omega=3.0$ are well reproduced. As is seen in Fig.
7(b), the reconstructed $\tilde{\Pi}(\tau)$ coincides with the original
$\Pi(\tau)$.

##  7 Concluding remarks

In this paper, a full-fledged analytical continuation toolkit named ACFlow is
presented. It implements several primary analytical continuation methods,
including the maximum entropy method, stochastic analytical continuation, and
stochastic optimization method. It provides quite a few validation and
diagnostic tools. It can be used with great flexibility for the analytical
continuations of arbitrary fermionic and bosonic correlation functions
generated by finite-temperature quantum Monte Carlo simulations.

Note that analytical continuation problem is a hotspot in computational
physics and many-body physics all the time. Many efforts have been devoted to
solve it in recent years. Noticeable achievements include maximum quantum
entropy method [44], Nevanlinna analytical continuation [45, 46], blocked-mode
sampling and grid point sampling in stochastic analytical continuation [35,
36], constrained stochastic analytical continuation [31, 34], machine learning
assisted analytical continuation [18, 19], and so on. We would like to
incorporate these new progresses into the ACFlow toolkit in the near future.

Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

Data availability

The data that support the findings of this study will be made available upon
reasonable requests to the corresponding author.

Acknowledgement

This work is supported by the CAEP Foundation (under Grant No. CX100000) and
the National Natural Science Foundation of China (under Grants No. 11874329
and No. 11934020).

## References

  * [1] J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt, R. B. Wiringa, Quantum monte carlo methods for nuclear physics, Rev. Mod. Phys. 87 (2015) 1067–1118.  doi:10.1103/RevModPhys.87.1067.
  * [2] W. M. C. Foulkes, L. Mitas, R. J. Needs, G. Rajagopal, Quantum monte carlo simulations of solids, Rev. Mod. Phys. 73 (2001) 33–83.  doi:10.1103/RevModPhys.73.33.
  * [3] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, P. Werner, Continuous-time monte carlo methods for quantum impurity models, Rev. Mod. Phys. 83 (2011) 349–404.  doi:10.1103/RevModPhys.83.349.
  * [4] J. Gubernatis, N. Kawashima, P. Werner, Quantum Monte Carlo Methods: Algorithms for Lattice Models, Cambridge University Press, 2016.  doi:10.1017/CBO9780511902581.
  * [5] O. Gunnarsson, M. W. Haverkort, G. Sangiovanni, Analytical continuation of imaginary axis data for optical conductivity, Phys. Rev. B 82 (2010) 165125.  doi:10.1103/PhysRevB.82.165125.
  * [6] C. E. Creffield, E. G. Klepfish, E. R. Pike, S. Sarkar, Spectral weight function for the half-filled hubbard model: A singular value decomposition approach, Phys. Rev. Lett. 75 (1995) 517–520.  doi:10.1103/PhysRevLett.75.517.
  * [7] H. J. Vidberg, J. W. Serene, Solving the eliashberg equations by means ofn-point padéapproximants, J. Low Temp. Phys. 29 (3) (1977) 179–192.  doi:10.1007/BF00655090.
  * [8] K. S. D. Beach, R. J. Gooding, F. Marsiglio, Reliable padé analytical continuation method based on a high-accuracy symbolic computation algorithm, Phys. Rev. B 61 (2000) 5147–5157.  doi:10.1103/PhysRevB.61.5147.
  * [9] i. c. v. Osolin, R. Žitko, Padé approximant approach for obtaining finite-temperature spectral functions of quantum impurity models using the numerical renormalization group technique, Phys. Rev. B 87 (2013) 245135.  doi:10.1103/PhysRevB.87.245135.
  * [10] J. Schött, I. L. M. Locht, E. Lundin, O. Grånäs, O. Eriksson, I. Di Marco, Analytic continuation by averaging padé approximants, Phys. Rev. B 93 (2016) 075104.  doi:10.1103/PhysRevB.93.075104.
  * [11] J. E. Gubernatis, M. Jarrell, R. N. Silver, D. S. Sivia, Quantum monte carlo simulations and maximum entropy: Dynamics from imaginary-time data, Phys. Rev. B 44 (1991) 6011–6029.  doi:10.1103/PhysRevB.44.6011.
  * [12] M. Jarrell, J. Gubernatis, Bayesian inference and the analytic continuation of imaginary-time quantum monte carlo data, Phys. Rep. 269 (3) (1996) 133–195.  doi:https://doi.org/10.1016/0370-1573(95)00074-7.
  * [13] A. W. Sandvik, Stochastic method for analytic continuation of quantum monte carlo data, Phys. Rev. B 57 (1998) 10287–10290.  doi:10.1103/PhysRevB.57.10287.
  * [14] S. Fuchs, T. Pruschke, M. Jarrell, Analytic continuation of quantum monte carlo data by stochastic analytical inference, Phys. Rev. E 81 (2010) 056701.  doi:10.1103/PhysRevE.81.056701.
  * [15] A. S. Mishchenko, N. V. Prokof’ev, A. Sakamoto, B. V. Svistunov, Diagrammatic quantum monte carlo study of the fröhlich polaron, Phys. Rev. B 62 (2000) 6317–6336.  doi:10.1103/PhysRevB.62.6317.
  * [16] O. Goulko, A. S. Mishchenko, L. Pollet, N. Prokof’ev, B. Svistunov, Numerical analytic continuation: Answers to well-posed questions, Phys. Rev. B 95 (2017) 014102.  doi:10.1103/PhysRevB.95.014102.
  * [17] J. Otsuki, M. Ohzeki, H. Shinaoka, K. Yoshimi, Sparse modeling approach to analytical continuation of imaginary-time quantum monte carlo data, Phys. Rev. E 95 (2017) 061302.  doi:10.1103/PhysRevE.95.061302.
  * [18] R. Fournier, L. Wang, O. V. Yazyev, Q. Wu, Artificial neural network approach to the analytic continuation problem, Phys. Rev. Lett. 124 (2020) 056401.  doi:10.1103/PhysRevLett.124.056401.
  * [19] H. Yoon, J.-H. Sim, M. J. Han, Analytic continuation via domain knowledge free machine learning, Phys. Rev. B 98 (2018) 245101.  doi:10.1103/PhysRevB.98.245101.
  * [20] L.-F. Arsenault, R. Neuberg, L. A. Hannah, A. J. Millis, Projected regression method for solving fredholm integral equations arising in the analytic continuation problem of quantum physics, Inverse Problems 33 (11) (2017) 115007\.  doi:10.1088/1361-6420/aa8d93.
  * [21] R. Levy, J. LeBlanc, E. Gull, Implementation of the maximum entropy method for analytic continuation, Comput. Phys. Commun. 215 (2017) 149–155.  doi:https://doi.org/10.1016/j.cpc.2017.01.018.
  * [22] D. Bergeron, A.-M. S. Tremblay, Algorithms for optimized maximum entropy and diagnostic tools for analytic continuation, Phys. Rev. E 94 (2016) 023303.  doi:10.1103/PhysRevE.94.023303.
  * [23] J. Kaufmann, K. Held, ana_cont: Python package for analytic continuation, Comput. Phys. Commun. 282 (2023) 108519.  doi:https://doi.org/10.1016/j.cpc.2022.108519.
  * [24] I. Krivenko, M. Harland, Triqs/som: Implementation of the stochastic optimization method for analytic continuation, Comput. Phys. Commun. 239 (2019) 166–183.  doi:https://doi.org/10.1016/j.cpc.2019.01.021.
  * [25] I. Krivenko, A. S. Mishchenko, Triqs/som 2.0: Implementation of the stochastic optimization with consistent constraints for analytic continuation, Comput. Phys. Commun. 280 (2022) 108491.  doi:https://doi.org/10.1016/j.cpc.2022.108491.
  * [26] F. F. Assaad, M. Bercx, F. Goth, A. Götz, J. S. Hofmann, E. Huffman, Z. Liu, F. P. Toldin, J. S. E. Portela, J. Schwab, The ALF (Algorithms for Lattice Fermions) project release 2.0. Documentation for the auxiliary-field quantum Monte Carlo code, SciPost Phys. Codebases (2022) 1doi:10.21468/SciPostPhysCodeb.1.
  * [27] L. Boehnke, H. Hafermann, M. Ferrero, F. Lechermann, O. Parcollet, Orthogonal polynomial representation of imaginary-time green’s functions, Phys. Rev. B 84 (2011) 075145.  doi:10.1103/PhysRevB.84.075145.
  * [28] R. K. Bryan, Maximum entropy analysis of oversampled data problems, Eur. Biophys. J. 18 (3) (1990) 165–174.  doi:10.1007/BF02427376.
  * [29] K. S. D. Beach, Identifying the maximum entropy method as a special limit of stochastic analytic continuation (2004).  arXiv:0403055.
  * [30] H. Shao, Y. Q. Qin, S. Capponi, S. Chesi, Z. Y. Meng, A. W. Sandvik, Nearly deconfined spinon excitations in the square-lattice spin-$1/2$ heisenberg antiferromagnet, Phys. Rev. X 7 (2017) 041072.  doi:10.1103/PhysRevX.7.041072.
  * [31] A. W. Sandvik, Constrained sampling method for analytic continuation, Phys. Rev. E 94 (2016) 063308.  doi:10.1103/PhysRevE.94.063308.
  * [32] Y. Q. Qin, B. Normand, A. W. Sandvik, Z. Y. Meng, Amplitude mode in three-dimensional dimerized antiferromagnets, Phys. Rev. Lett. 118 (2017) 147207\.  doi:10.1103/PhysRevLett.118.147207.
  * [33] S. R. White, D. J. Scalapino, R. L. Sugar, N. E. Bickers, Monte carlo calculation of dynamical properties of the two-dimensional hubbard model, Phys. Rev. Lett. 63 (1989) 1523–1526.  doi:10.1103/PhysRevLett.63.1523.
  * [34] H. Shao, A. W. Sandvik, Progress on stochastic analytic continuation of quantum monte carlo data (2022).  arXiv:2202.09870.
  * [35] K. Ghanem, E. Koch, Average spectrum method for analytic continuation: Efficient blocked-mode sampling and dependence on the discretization grid, Phys. Rev. B 101 (2020) 085111.  doi:10.1103/PhysRevB.101.085111.
  * [36] K. Ghanem, E. Koch, Extending the average spectrum method: Grid point sampling and density averaging, Phys. Rev. B 102 (2020) 035114.  doi:10.1103/PhysRevB.102.035114.
  * [37] K. Vafayi, O. Gunnarsson, Analytical continuation of spectral data from imaginary time axis to real frequency axis using statistical sampling, Phys. Rev. B 76 (2007) 035115.  doi:10.1103/PhysRevB.76.035115.
  * [38] O. F. Syljuåsen, Using the average spectrum method to extract dynamics from quantum monte carlo simulations, Phys. Rev. B 78 (2008) 174429.  doi:10.1103/PhysRevB.78.174429.
  * [39] S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671–680.  doi:10.1126/science.220.4598.671.
  * [40] E. Marinari, Optimized monte carlo methods (1996).  arXiv:9612010.
  * [41] A. Georges, G. Kotliar, W. Krauth, M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68 (1996) 13–125.  doi:10.1103/RevModPhys.68.13.
  * [42] L. Huang, Y. Wang, Z. Y. Meng, L. Du, P. Werner, X. Dai, iqist: An open source continuous-time quantum monte carlo impurity solver toolkit, Comput. Phys. Commun. 195 (2015) 140–160.  doi:https://doi.org/10.1016/j.cpc.2015.04.020.
  * [43] L. Huang, iqist v0.7: An open source continuous-time quantum monte carlo impurity solver toolkit, Comput. Phys. Commun. 221 (2017) 423–424.  doi:https://doi.org/10.1016/j.cpc.2017.08.026.
  * [44] J.-H. Sim, M. J. Han, Maximum quantum entropy method, Phys. Rev. B 98 (2018) 205102\.  doi:10.1103/PhysRevB.98.205102.
  * [45] J. Fei, C.-N. Yeh, E. Gull, Nevanlinna analytical continuation, Phys. Rev. Lett. 126 (2021) 056402.  doi:10.1103/PhysRevLett.126.056402.
  * [46] J. Fei, C.-N. Yeh, D. Zgid, E. Gull, Analytical continuation of matrix-valued functions: Carathéodory formalism, Phys. Rev. B 104 (2021) 165111.  doi:10.1103/PhysRevB.104.165111.