Reconstructing lattice QCD spectral functions with stochastic pole expansion and Nevanlinna analytic continuation

The reconstruction of spectral functions from Euclidean correlation functions is a well-known, yet ill-posed inverse problem in the fields of many-body and high-energy physics. In this paper, we present a comprehensive investigation of two recently developed analytic continuation methods, namely stochastic pole expansion and Nevanlinna analytic continuation, for extracting spectral functions from mock lattice QCD data. We examine a range of Euclidean correlation functions generated by representative models, including the Breit-Wigner model, the Gaussian mixture model, the resonance-continuum model, and the bottomonium model. Our findings demonstrate that the stochastic pole expansion method, when combined with the constrained sampling algorithm and the self-adaptive sampling algorithm, successfully recovers the essential features of the spectral functions and exhibits excellent resilience to noise of input data. In contrast, the Nevanlinna analytic continuation method suffers from numerical instability, often resulting in the emergence of spurious peaks and significant oscillations in the high-energy regions of the spectral functions, even with the application of the Hardy basis function optimization algorithm.


I. INTRODUCTION
Lattice QCD (LQCD) is a well-established first-principles and non-perturbative approach for studying strong interactions [1][2][3].It serves as a valuable tool in understanding the genesis and evolution of the quark-gluon plasma (QGP) [4] and mapping out the phase diagram of strong-interaction matter [5][6][7].In LQCD, spectral functions play a vital role in scrutinizing and elucidating high-energy physical phenomena that involve quarks and gluons, such as the melting of heavy quarkonium [8][9][10][11][12][13][14][15] and the transport properties [16][17][18] of the QGP formed through relativistic heavy-ion collisions.However, accessing the spectral functions and other dynamical properties of the QCD medium from lattice simulations remains challenging due to LQCD's typical formulation on a discrete Euclidean space-time grid [1][2][3].Therefore, researchers must reconstruct the spectral functions from numerically computed Euclidean correlation functions on the lattice to understand the relevant physics and compare the theoretical results with corresponding experimental data obtained from the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) [19].
Mathematically, the Euclidean correlators G(t) and the spectral functions ρ(s) are connected through a Fredholm equation of the first kind: G(t) = K(t, s) ⊛ ρ(s), where K(t, s) represents a continuous kernel function and ⊛ signifies convolution.Mapping the spectral functions to the Euclidean correlators is a straightforward process that can be easily accomplished using numerical integration.However, extracting spectral functions from Euclidean correlators through analytic continuation poses a formidable challenge [19].We observe that similar inverse problems are quite common in many-body and high-energy physics [1,20].They are considered ill-posed.There are two main reasons for this statement.
Firstly, the Euclidean correlators are evaluated at a finite number of points due to the space-time discretization in LQCD [1][2][3].Secondly, as LQCD simulations rely on stochastic Monte Carlo sampling, the resulting Euclidean correlators are inherently noisy [21].Even small deviations or fluctuations in the Euclidean correlators result in significant uncertainties in the spectral functions.As a result, the majority of resulting spectral functions exhibit high oscillations and lack physical significance, thereby inhibiting a reliable comparison between the theoretical spectra and experimental data.To address these issues, researchers have developed a plethora of analytic continuation approaches in the past decades.Next, we will briefly introduce some of the most commonly employed methods in LQCD simulations.
Maximum entropy method.The maximum entropy method (MaxEnt) is perhaps the most popular analytic continuation tool, and it has dominated this field for a long time.In this method, the spectral density is interpreted as a probability distribution [22][23][24][25].The primary objective is to extract the most probable spectral density ρ from the correlation function G and maximize the posterior probability Pr [ρ|G].According to Bayes's theorem, Pr[ρ|G] ∝ Pr[G|ρ]Pr [ρ], where Pr [G|ρ] is the likelihood function and Pr[ρ] is the prior probability.It is important to incorporate analytical knowledge related to spectral properties in LQCD, such as positive definiteness or even the presence of pole structures within the spectra, into the probability distribution.A significant portion of the prior information could be encoded within the prior probability Pr [ρ], which is proportional to exp (αS ), where α is a regulation parameter and S denotes entropy.It is worth emphasizing that the entropic term αS is not unique and typically takes the form of the generalized Shannon-Jaynes entropy [22,23].Another popular alternative is the Bayesian reconstruction entropy [25,26].While the MaxEnt method is widely recognized for its efficiency and noise-tolerance, it sometimes struggles to faithfully recover sharp, subtle, and high-frequency features within the spectral functions.
Stochastic analytic continuation.In the past decade, the arXiv:2309.11114v1[hep-lat] 20 Sep 2023 stochastic analytic continuation method (SAC) and its variants have emerged as formidable contenders to surpass the Max-Ent method [27][28][29][30][31][32][33][34][35][36][37].Unlike the MaxEnt method, the SAC method treats all spectral functions equally instead of selecting the most probable one.Initially, the spectral functions are parameterized with hundreds or thousands of δ-like functions.And then these parameters, such as the amplitudes and locations of the δ functions, are stochastically sampled at a fictitious temperature Θ using a Boltzmann-like weight function, which essentially serves as a likelihood function [27,28].Finally, the gathered spectral functions are filtered and averaged.We note that Shao and Sandvik et al. have proven the equivalence in a generalized thermodynamic limit (large number of degrees of freedom) of the average spectrum and the maximum entropy solution [38].Although the SAC method supplements the MaxEnt method by enabling the resolution of subtle structures in spectra, it requires significant computational resources [32].
Machine learning approaches.In recent years, several machine learning aided methods have been developed to address the analytic continuation challenges in LQCD simulations.These methods include the deep neural networks (DNN) [39,40], radial basis functions network (RBFN) [41], entropy variational autoencoder (SVAE) [42], kernel ridge regression (KRR) [43], automatic differentiation (AD) [44], Gaussian processes regression (GPR) [45], and many others.Although these methods may offer improved performance in certain situations, their universality is not guaranteed.Furthermore, some studies have adopted supervised approaches to train machine learning network models, which incorporate prior knowledge from specific physics insights into the training sets.However, caution must be exercised as there is a risk of introducing biases in the training data.
Recently, two new analytic continuation methods, namely the stochastic pole expansion (SPX) [46] and the Nevanlinna analytic continuation (NAC) [47,48], have been proposed.The SPX method inherits the spirit of the SAC method, where the Matsubara Green's function is initially parameterized with hundreds or thousands of poles.Subsequently, the amplitudes and positions of these poles are optimized using a stochastic algorithm that based on stimulated annealing [49].The SPX method is applicable to fermionic and bosonic systems.It has been extended to support analytic continuation of matrixvalued Green's functions [46].On the other hand, similar to the Padé approximation (PA) [50][51][52][53], the NAC method aims to interpolate the Matsubara data in the complex plane using some form of continued fraction expansion [47,48].It takes the "Nevanlinna" analytic structure of the Matsubara Green's function into consideration, ensuring that the calculated spectral functions are inherently positive and normalized.However, this method is highly sensitive to the noise level in the raw Matsubara data.With noiseless data as input, it can successfully resolve complex spectral functions across a wide energy range with unprecedented accuracy.Unfortunately, when the Matsubara data contains noise, the Nevanlinna interpolants may not exist, and the resulting spectral functions are not guaranteed to be causal.
Both the SPX and NAC methods have not yet been em-ployed in addressing the issue of analytic continuation in the LQCD simulations, and it remains uncertain whether they are applicable for the analytic continuation of Euclidean correlation functions [19].Therefore, the purpose of this study is to fill this knowledge gap and expand the potential applications of these methods.We first generate noisy Euclidean data using four representative models: the Breit-Wigner model, the Gaussian mixture model, the resonance-continuum model, and the bottomonium model.These synthetic data sets are then processed using the SPX, NAC, and MaxEnt methods.Finally, we conduct a comprehensive comparison between the calculated spectra and the exact solutions if available.The results suggest that the SPX method manifests comparable or even superior performance in comparison to the commonly used MaxEnt method.On the other hand, the NAC method tends to suffer from numerical instability, even in the absence of noise in the input data.
The structure of the remaining sections of this paper is as follows: Section II provides an introduction to the basic formulations of Euclidean correlation functions and offers a brief overview of the SPX and NAC methods.In Section III, we present the computational setups, demonstrate, and discuss four representative examples.Section IV explores the robustness of the SPX method and numerical instability of the NAC method in the presence of noisy Euclidean data.Additionally, we analyze the effects of the constrained sampling algorithm and the self-adaptive sampling algorithm on the SPX method, along with the impact of the Hardy basis function optimization algorithm on the NAC method.Finally, our findings are summarized in Section V.

II. METHOD A. Euclidean correlation function
At finite temperature, the Euclidean correlation function G(τ) is related to the spectral function ρ(ω) through [19] Here β = 1/T represents the inverse system temperature, and τ represents the Euclidean (imaginary) time interval (τ ∈ [0, β]).In momentum space, the expression for the Euclidean correlation function is [54] where p represents the Euclidean (Matsubara) frequency, the value of G(p) can be derived by performing a discrete Fourier transform on G(τ).In the references, Eq. ( 2) is also known as the Källén-Lehmann (KL) spectral representation [55].By applying the process of analytic continuation, the retarded propagator G R (ω) can be attained, enabling the extraction of the spectral function through the following expression: where ω = −ip.It is important to note that ρ(ω) must be an odd function for bosonic system, i.e., ρ(−ω) = −ρ(ω).

B. Stochastic pole expansion
According to the textbooks of many-body physics, the Lehmann representation of the finite temperature many-body Green's functions is given by the following formula [20,56]:  4) can be simplified as: It is evident that only terms where A mn 0 can be present.The indices m and n can also be compressed as γ, resulting in the following expression: Eq. ( 6) is referred to as the pole representation of the manybody Green's functions [20].In this representation, N p denotes the number of poles, and A γ and P γ mean the amplitude and position of the γ-th pole.For the Euclidean correlation in momentum space, its pole representation can be reformulated as: Here, Ξ represents the kernel matrix, which is calculated using the following equation: Ãγ is the renormalized amplitude of the γ-th pole, given by: It can be easily proven that Ãγ and P γ must satisfy the following constraints: We assume that the input Euclidean correlation function is denoted as G(p n ), and the input data consists of N frequency points.We then utilize Eq. ( 7) to approximate the Euclidean data.To assess the discrepancy between G(p n ) and G(p n ), we introduce the so-called goodness-of-fit function χ 2 .Its definition is as follows: where || • || F represents the Frobenius norm.Hence, the objective of the analytic continuation is to solve the subsequent multivariate optimization problem: Once the optimized parameters N p , Ãγ , and P γ are determined, evaluating the retarded Green's function is straightforward by substituting p with ω + i0 + in Eq. ( 7).Additionally, the spectral function ρ(ω) is computed using Eq. ( 3).It is important to note that this optimization problem [i.e., Eq. ( 12)] is highly non-convex.Traditional gradient-based optimization methods typically fail to identify the global minimum unless the initial solution is of high quality [57].Therefore, in the SPX method, we employ the simulated annealing algorithm [49] to optimize the Ãγ and P γ parameters subject to the constraints defined by Eq. ( 10).For technical details regarding the possible Monte Carlo random walking rules in the configuration space C = { Ãγ , P γ }, please refer to Ref. [46].The advantages of the SPX method include its ability to derive approximate expressions for correlation functions and its ease of extension to support the analytic continuation of bosonic systems, two-particle Green's functions, matrix-valued Green's functions, and so on.Application to noisy Matsubara data suggests that the SPX method can accurately resolve both continuum spectra for condensed matter cases and multiple δ-like peaks for molecule cases.Notably, it performs well in reproducing sharp high-frequency features.

C. Nevanlinna analytic continuation
It is well known that the retarded Green's function, denoted as G R (ω + i0 + ), and the Matsubara Green's function, denoted as G(iω n ), can both be consistently represented as G(z), where z ∈ C\R.The NAC method utilizes the fact that the negative fermionic Green's function, denoted as f (z) = −G(z), belongs to the class of Nevanlinna functions.By applying the invertible Möbius transform h(z) = (z − i)/(z + i) to the function value of f (z), the Nevanlinna function is mapped in a oneto-one fashion to a contractive function θ(z) = h[ f (z)].This contractive function θ(z) can be expressed in the form of a continued fraction expansion, and an iterative algorithm can be constructed accordingly [47].The recursion relation between two steps θ j (z) and θ j+1 (z) is given by: In this equation, h j (z) = (z − Y j )/(z + Y j ), Y j = iω j represents the j-th Matsubara frequency used, and γ j = θ j (Y j ) represents the function value of the j-th contractive function at the point Y j .The final expression of the recursive function θ(z) can be written as [58]: where with j increasing from left to right.Here N s is the overall iteration step, which is equivalent to the number of data points.
After obtaining θ(z), one can immediately get the Green's function by an inverse Möbius transform as Note that the Pick criterion [59] should be fulfilled for the existence of the Nevanlinna interpolation.Additionally, it is worth noting that there is flexibility in choosing θ N s +1 (z), which can be used to select the most desirable spectral function.In Reference [47], θ N s +1 (z) is expanded in the Hardy basis and chosen in such a way that it achieves the smoothest possible spectral function [60].The loss function employed in this selection process is given by: This loss function consists of two terms.The first term enforces the proper sum rule, while the second term incorporates the smoothness condition.λ is an adjustable parameter.By preserving the "Nevanlinna" analytic structure of Green's functions, the NAC method automatically generates positive and normalized spectral functions [47].However, it is important to emphasize that the method is sensitive to noise, and either a large number of data points N or a high Hardy order H can potentially lead to numerical instabilities.Although the NAC method has been extended to support the analytic continuation of matrix-valued Green's functions [48], it cannot be directly applied to bosonic systems in its original formalism.Quite recently, Nogaki et al. suggest an ingenious trick to work around this limitation.Their basic idea is to introduce an auxiliary fermionic function [61].Let us start with a bosonic Green's function G(τ) that satisfies the periodic condition G(τ+β) = G(τ).One can construct an artificial antiperiodic fermionic Green's function G(τ) as follows: Clearly, this auxiliary fermionic Green's function exhibits the same value as the bosonic Green's function in the range 0 < τ < β.It is easy to prove that the relation between the bosonic spectral function ρ(ω) and the auxiliary fermionic spectral function ρ(ω) is as follows: Furthermore, the sum rule for ρ(ω) is given by: Given G(τ), it is easy to construct G(iν n ) via direct Fourier transformation, where ν n = (2n+1)π/β are the fermionic Matsubara frequencies.Since one can perform analytic continuation for G(iν n ) via the standard NAC method to get ρ(ω).And then the bosonic spectral function ρ(ω) can be derived according to Eq. ( 18).This procedure has been outlined in Figure 1 in Ref. [61].

III. BENCHMARKS A. Computational setups
To benchmark the SPX and NAC methods, we consider four typical models, namely the Breit-Wigner model, the Gaussian mixture model, the resonance-continuum model, and the bottomonium model in the present investigation.The first three models provide analytic formulas for generating the exact spectral functions, denoted as ρ(ω).For the bottomonium model, we take the "exact" spectral function from Ref. [62] as input.Using these spectral functions, one can create clean Euclidean data, denoted as G clean , using Eq. ( 2) [63].To mimic the noise present in LQCD simulations [21], we manually add multiplicative Gaussian noises to the clean Euclidean data.The formula we use is as follows: where δ measures the noise level of the input data and N C (0, 1) represents complex-valued normal Gaussian noise [64].In our subsequent analysis, unless explicitly stated, we set δ = 10 −4 for the SPX method and δ = 0.0 for the NAC method.We then supply the noisy Euclidean data, denoted as G noisy , into the SPX and NAC codes to extract the spectral functions.Finally, we compare the calculated spectral functions with the corresponding exact solutions.
The SPX method has been implemented within the ACFlow package [65].In this study, the number of poles (N p ) is fixed to 2000.For each test, we perform a total of 2 × 10 3 individual SPX runs.Each SPX run consists of 2 × 10 5 Monte Carlo sampling steps.The spectral functions generated in all SPX runs are gathered and the corresponding χ 2 values are recorded.Assuming the mean value of the collected χ 2 is denoted as ⟨χ 2 ⟩, we only retain the solutions whose χ 2 values are smaller than ⟨χ 2 ⟩/α good , and apply them to calculate the averaged spectrum.Here, α good is an adjustable parameter (α good ≥ 1.0).Its optimal value is about 1.2.
Regarding the NAC method, we utilize another open-source toolkit, namely the Nevanlinna.jlpackage [66].In the present calculations, the lowest 100 Matsubara frequencies are kept as input.In order to avoid breaking the Pick criterion [59], the optimal number of data points is automatically determined, which is denoted as N opt , through a "Pick Selection" procedure in the algorithm.Usually a typical value for N opt is 10.During the simulated process, the Hardy basis function optimization algorithm is always enabled.The highest Hardy order, denoted as H max , is set to be 50.To ensure numerical stability, the cutoff value of Hardy order, denoted as H cut , should be determined automatically on a case-by-case basis.The λ parameter, as seen in Eq. ( 16), is set to be 10 −6 .
In addition to the SPX and NAC methods, the classic Max-Ent method [24] is also employed to yield analytic continuation results for comparison.We again utilize the ACFlow package, which provides a state-of-the-art implementation of the MaxEnt method [65].For the MaxEnt simulations, we usually choose a flat default model and calibrate the regularization parameter α by using the χ 2 -kink algorithm [67].The starting and ending values of α are 10 16 and 10 1 , respectively.The ratio between two successive α parameters, i.e. α i /α i+1 , is 10.For the prior probability (the entropic term), we adopt the generalized Shannon-Jaynes entropy [22,23], but the Bayesian reconstruction entropy [25,26] is also examined.The simulated results are in good agreement with each other.Thus, we only present results obtained with the generalized Shannon-Jaynes entropy in the following discussions.

B. Breit-Wigner model
The Breit-Wigner spectral function, obtained from a parameterization derived directly from one-loop perturbative quantum field theory [19,39], is expressed as follows: where M represents the mass of the corresponding state, Γ is the width, and A is a positive constant.We start with a superposed collection of Breit-Wigner peaks.Specially, two typical scenarios are investigated in the present work: The system temperature T is fixed to be 0.02 GeV.The synthetic Euclidean data (G noisy ) consists of 50 frequency points for both the SPX method and the MaxEnt method, and 100 frequency points for the NAC method.The analytical continuation results obtained by the SPX, NAC, and MaxEnt methods are illustrated in Figure 1.
Overall, the SPX method demonstrates better performance for the 1BW model.It captures precisely not only the height, but also the position of the Breit-Wigner peak.In comparison, the NAC and MaxEnt methods tend to overestimate the width and underestimate the height of the Breit-Wigner peak.Furthermore, the NAC method leads to an obvious oscillation phenomenon around 0.5 GeV.For the 2BW model, all three methods are able to reproduce the low-energy peak around 1.0 GeV, but encounter difficulty in resolving the high-energy peak near 3.0 GeV.Although the SPX method successfully identifies the position of the high-energy peak, its weight is not accurately resolved.The spectrum obtained by the NAC method is quite sensitive to the λ parameter in this case.If λ = 10 −4 (it is the default choice of the code), the highenergy peak is shifted towards a higher energy (∼ 4.5 GeV) and broadened significantly.Only when λ = 10 −6 , the highenergy peak is well reproduced.The spectrum obtained by the MaxEnt method is also not ideal.The high-energy peak is smeared out and replaced with a shoulder-like feature.

C. Gaussian mixture model
Just as its name implies, the spectral function of the Gaussian mixture model [24] is a superposition of some Gaussian peaks.It can be expressed by the following equation: where A i , M i , and Γ i represent the amplitude, position, and broadening of the i-th Gaussian peak.In this example, we consider a three Gaussian peaks model.The specific values for the model parameters are as follows: GeV, M 3 = 6.5 GeV, Γ 1 = 0.01 GeV, Γ 2 = 0.2 GeV, and Γ 3 = 1.5 GeV.The mock Euclidean data consist of 100 points, and T = 0.02 GeV.In the simulations, we have enhanced the SPX method by utilizing the self-adaptive sampling algorithm, which we refer to as SA-SPX [46].The results of the analytical continuation are presented in Figure 2.
It is expected that the true spectrum will exhibit three welldefined peaks.We find that all three methods are successful in recovering the low-energy sharp peak at M 1 .However, resolving the two high-energy peaks presents some challenges.Specifically, for the SA-SPX method, it is able to roughly resolve the peak at M 2 , but it tends to overestimate the width of the peak at M 3 .In order to save the computational resources, the iteration number of the self-adaptive sampling algorithm is fixed to 10.Nevertheless, it should be emphasized that increasing the iteration number could further reduce the peak's width at M 3 .In Section IV B, we will delve into the combination of the self-adaptive sampling algorithm with the SPX method and discuss the usefulness of the SA-SPX method in resolving complex LQCD spectral functions.As for the NAC method, it produces a sharp peak around 2.1 GeV and a broad peak around 4.5 GeV, which are both in incorrect positions.More specifically, this method underestimates the energies of the two high-energy peaks.We also adjust the λ parameter, but it doesn't help.Regarding the MaxEnt method, it can recover the peak at M 2 , although with a larger width.However, it fails to resolve the peak at M 3 , instead exhibiting a broad hump around 7.0 ± 3.0 GeV.As for the SPX method, the self-adaptive sampling algorithm is enable to obtain a more reasonable spectrum.Once the self-adaptive sampling algorithm is turned off, the obtained spectrum by the SPX method will resemble the one by the MaxEnt method.The spectrum should become smoother and almost featureless in the high energy region (ω > 4.0 GeV).The NAC method is enhanced by the Hardy basis function optimization algorithm and N opt = 11.Please see the main text for more details.

D. Resonance-continuum model
The resonance-continuum model is a physics-motivated model borrowed from References [37,44].The spectral function of the resonance-continuum model can be viewed as a nonlinear combination of the resonance part (ρ r ) and the con-tinuum part (ρ c ): Here, ξ 1 and ξ 2 are the mixing coefficients.Their definitions are as follows: and where ξ is a cutoff function: It is used to smooth out the constructed spectral function.The resonance part of the spectral function is given by: which follows a relativistic Breit-Wigner form.The continuum part of the spectral function is expressed as: In this example, the model parameters are  For the resonance-continuum model [see Fig. 3(a)], it is evident that the resonance peak at M r is approximately reproduced.However, in the continuum part (ω > 0.4 GeV), both the SPX and MaxEnt methods exhibit moderate oscillations.These oscillations will decay when ω is increased.The NAC method results in huge oscillations, especially when ω is large.This unphysical feature can not be eliminated or suppressed by the Hardy basis function optimization algorithm.This fact suggests that the three analytic continuation methods do not accurately describe the continuum part.In the case of the resonance model [see Fig. 3(b)], all three methods successfully recover the location, width, and height of the resonance peak.As for the continuum model [see Fig. 3(c)], all three methods produce oscillating spectra.In particular, the NAC method leads to more pronounced oscillations.Even worse, these oscillations are enhanced with the increment of ω.The only useful information that can be extracted from the calculated spectra is the location of the band edge.

E. Bottomonium spectrum
The "exact" bottomonium spectrum utilized in this study is taken directly from References [62,70].It is generated by a N f = 2 + 1 LQCD calculation.The temperature employed in the LQCD simulation is 201 MeV, which exceeds the deconfinement crossover temperature (T c ).The spectrum is specifically for the Υ channel.Initially, we synthesize the Euclidean data by Eq. ( 2) for the first 100 Matsubara frequencies.Then, random Gaussian noises are added by Eq. (21).In this example, we adopt the combination of the SPX method with the constrained sampling algorithm (dubbed C-SPX) [46] to improve the performance.Specifically, locations for the randomly generated poles (P γ ) are restricted to the energy range: ω ∈ [9.5 GeV,16.0GeV] [68].For the NAC method, the Hardy optimization trick is applied [47].Regarding the MaxEnt method, its default model is a shifted Gaussian function [69].Here, the terminology "C-SPX" implies that the positions of the poles are restricted in the SPX simulations [46,68].The Hardy basis function optimization algorithm is enabled for the NAC method and N opt = 8.A shifted Gaussian function, instead of a constant, is used as the default model for the MaxEnt method [69].See the main text for more details.
The analytic continuation solutions, together with the "exact" bottomonium spectrum, are illustrated in Fig. 4.
As is seen in Fig. 4, the ideal bottomonium spectrum consists of a single resonance peak at approximately 9.6 GeV and a "rise-and-decay" feature with two sizable bumps around 10.8 GeV and 12.0 GeV.By employing the constrained algorithm, the SPX method successfully resolves the left boundary of the resonance peak and captures the long tail of the "rise-and-decay" feature.However, it falls short in resolving the resonance peak and the two bumps.In Section IV B, we will demonstrate the application of the self-adaptive sampling algorithm to cure this issue partly [46].For the NAC method, it accurately reproduces the resonance peak.But it fails to recover the desired "rise-and-decay" feature, which is instead replaced by two distinct peaks located at approximately 11.3 GeV and 14.0 GeV.According to our experience, though the spectrum obtained by the NAC method is not accurate enough, it can provide some hints about the energy range of the true spectrum (such as the position of the resonance peak in this case).We can use this information to refine subsequent simulations by imposing more reasonable constraints for the C-SPX method and more appropriate default models for the MaxEnt method.By employing the MaxEnt method, the resonance peak and the "rise-and-decay" pattern are smoothed out.Such that the calculated spectrum only features a Gaussian-like peak centered around 11.5 GeV.In addition to the shifted Gaussian model, we also benchmark alternative default models, such as the flat model, the shifted Lorentzian model, and the two Lorentzians model, etc.However, they do not contribute to improving the results.

IV. DISCUSSIONS A. Robustness with respect to noisy data
In the previous work, the robustness of the SPX method in the presence of noisy data from quantum Monte Carlo simulations has been demonstrated [46].This study aims to reexamine the noise resilience of the SPX method when applied to noisy LQCD data.Let us take the 1BW model as an example (please refer to Section III B for the model parameters) to address this issue.The noise level is varied from δ = 10 −8 to δ = 10 −2 .The analytic continuation results are displayed in Fig. 5. Just as expected, the SPX method is highly robust to variations in noise levels of the input Euclidean data.For low noise levels (10 −4 ≤ δ ≤ 10 −8 ), the calculated spectra closely approximate the exact solution and manifest minimal deviation.At a moderate noise level (δ = 10 −3 ), the calculated spectrum shows a light fluctuation around ω = 0.5 GeV, and the main peak is shifted slightly (approximately 0.1 GeV) towards the lower energy region.For a high noise level (δ = 10 −2 ), three sharp peaks emerge in the calculated spectrum.Apart from the peak at 1.8 GeV, the other peaks at 0.5 GeV and 3.0 GeV are unphysical.In Fig. 5(h), a plot of log(⟨χ 2 ⟩) against log(δ −1 ) is shown.Initially, log(⟨χ 2 ⟩) decreases linearly as log(δ −1 ) increases from 2.0 to 5.0.Subsequently, it approaches to a constant value (approximately -8.1) when log(δ −1 ) ≥ 6.0.This benchmark suggests that the SPX method remains robust when applied to noisy LQCD data, even with a moderately elevated noise level.Nonetheless, minimizing the noise level can enhance the performance of the SPX method.
Fei and Gull et al. have pointed out that the NAC method requires high-precision input data to ensure the Pick criterion is not violated and the existence of the Nevanlinna interpolants [47,64].Thus, in the previous calculations, we just assume that the input Euclidean data are noiseless for the NAC method (δ = 0.0).Now let us examine the noise resilience of the NAC method for synthetic LQCD data.For the sake of simplicity, the 2BW model is taken as a test-bed.The model parameters are presented in Section III B. The noise level δ is fixed to be 10 −8 .For each NAC run, the noisy part of the mock Euclidean data is always refreshed.We repeat the NAC simulations for 100 times with and without Hardy basis function optimization [60].Then we collect the calculated spectra and evaluate their arithmetic average.The analytic continuation results are shown in Fig. 6.We confirm again that the NAC method is extremely sensitive to noise, irrespective of the Hardy basis function optimization algorithm.Small fluctuations in the input Euclidean data can lead to huge variations in the resulting spectra.Even though the noise level is quite small, the performance of the NAC method is not good.It fails to capture the major characteristics of the 2BW model, and produces some spurious peaks.If the noise level is further increased, its performance should deteriorate ulteriorly (not shown in Figure 6).

B. Self-adaptive sampling algorithm for the SPX method
In the SPX method, the poles should be placed in a very dense frequency grid.In general, such a frequency grid can be either uniform or non-uniform.So, some prior knowledge about the spectrum and the physical system could be encapsulated in the form of the grid to improve the performance and usefulness of the SPX method.This has led to the development of the C-SPX method [46].Previous works have suggested that by modifying boundaries and grid interval distribution of the grid, the SPX method is capable of capturing complicated features in the spectra.However, it can be observed in Fig. 4 that the C-SPX method, as well as the NAC and Max-Ent methods, fail to resolve the major characteristics of the bottomonium spectrum [62,70].It implies that simple constraints on the spectral boundaries (or limitations on scopes of the poles) are not enough.We have to figure out a systematic way to refine the probability distribution of the poles to approximate the true spectrum.Next, we will demonstrate how to achieve this goal by a combination of the self-adaptive sampling algorithm and the SPX method (dubbed SA-SPX) [46].
The main principle behind the SA-SPX method is to iteratively adjust the grid interval distribution.In consequence, the probability distribution of the poles is coordinated to approximate the true spectral density and the corresponding goodness-of-fit functional [see Eq. ( 11)] is automatically minimized.This can be achieved by using the spectral density obtained from a previous SPX run or from other analytic continuation methods to update the grid.Now let us concentrate on the bottomonium model again [62,70].Its parameters can be found in Section III E. Initially, we generate the first frequency grid for the poles using the spectral functions obtained by the NAC and MaxEnt methods.From the calculated spectra, one can conclude that there is likely a sharp resonance peak around ω = 9.6 GeV (with a band edge at approximately 9.5 GeV) and a broad feature ranging from 10.0 GeV to 16.0 GeV (see Fig. 4).Keeping these hints in mind, we try to de- sign a pseudo-spectrum consisting of two Gaussian peaks: where A 1 = 5.0,A 2 = 1.80,M 1 = 9.60, M 2 = 11.5, Γ 1 = 0.01, and Γ 2 = 5.0.This pseudo-spectrum is illustrated in Fig. 7(a).
The first peak at M 1 originates from the resonance peak identified by the NAC method, while the second peak at M 2 is inspired by the spectrum obtained by the MaxEnt method.It is worth noting that these spectral parameters can be further adjusted to mimic more accurately the results obtained by the NAC and MaxEnt methods.This pseudo-spectrum serves as a reference model.To generate the frequency grid for the poles, we execute the following steps: (1) Calculate the integrated spectral function ϕ(ϵ) via the equation: Here, ω min and ω max are the left and right boundaries of the spectrum, respectively.And ϵ represents a point within the interval [ω min , ω max ]. (2) Evaluate the new frequency grid f i by using the equation: where λ i is a linear mesh in the interval [ϕ(ω min ), ϕ(ω max )], and N f is the number of grid points.Now the boundaries for the grid are set as ω min = 9.5 GeV and ω max = 16.0GeV.The resulting new grid, as displayed in Fig. 7(b), is compared with the standard linear grid.Next, the newly generated grid is utilized to perform a SPX simulation from scratch.The calculated spectrum is then used to generate a newer frequency grid [at this time, the ρ pseudo (ω) in Eq. ( 31) should be replaced with the calculated spectrum ρ(ω)], and the SPX simulation is repeated.This iterative procedure is carried out until the obtained spectrum and frequency grid are converged.In our experience, 5 ∼ 10 iterations are typically sufficient for achieving convergence.Figure 7(c) shows the results obtained by using the C-SPX and SA-SPX methods, as well as the exact spectrum.It is evident that the spectrum obtained with the SA-SPX method comes closer to the exact spectrum than that obtained with the C-SPX method.The sharp resonance peak, the small bump near 12.0 GeV, and the long tail of the rise-and-decay feature are well reproduced by using the SA-SPX method.The only missing characteristic is the bump near 10.8 GeV.Additionally, an error analysis about the reconstructed Euclidean data is presented in Fig. 7(d).The goodness-of-fit function of the NAC method is the largest (χ 2 ≈ 0.01 to 0.1), while those of the MaxEnt, C-SPX, and SA-SPX methods are comparable (χ 2 ≈ 10 −6 to 10 −4 ).This indicates that the spectrum obtained by the NAC method for this particular case is not reliable.

C. Hardy basis function optimization for the NAC method
As mentioned before, the iterative interpolation algorithm for the NAC method allows for the selection of an arbitrary contractive function, denoted as θ N s +1 , at the final step.In the literature, Fei and Gull et al. proposed an optimized approach [47], in which θ N s +1 is expanded using the Hardy basis and its conjugate generate functions [60], and the coefficients for the expansion are determined by minimizing a smoothness norm [see Eq. ( 16)].They demonstrated that using a constant value for θ N s +1 is apt to yield spectral functions with oscillations, whereas the optimized algorithm is useful for eliminating these oscillations and generating smoother spectral functions.Until now we just employ the optimized NAC approach for analytic continuations.However, we wonder whether the Hardy basis optimization is always better than the standard option.In order to answer this question, we perform additional tests for the Breit-Wigner model by using the NAC method.We compare the analytic continuation results obtained with a constant θ N s +1 and the optimized θ N s +1 (see Fig. 8).As anticipated, the optimized θ N s +1 suppresses oscillations in a large measure and yields smoother spectra, especially for the 2BW model.However, we also notice that the performance of the Hardy basis function optimization algorithm strongly depends on the value of the λ parameter.The optimized NAC method tends to make a wrong estimation about the location of the high-energy peak if the λ parameter is not reasonable.Therefore, if we know nothing about the basic features of the spectra, perhaps a constant θ N s +1 is a much safer choice.

V. CONCLUSION
In the present work, we conduct a systematic investigation of two newly developed methods, namely the SPX method and the NAC method, for analytically continuing for the mock LQCD data.We treat four exact spectral functions, which are derived from physically motivated models or realistic LQCD simulations, including the Breit-Wigner model, the Gaussian mixture model, the resonance-continuum model, and the bottomonium model.We use the exact spectral functions to build clean Euclidean data by numerical integration.And later the statistical noises are added.The synthetic Euclidean data are used as input and then transformed back to the real axis using different analytic continuation methods.By comparing the results with the exact spectra, we are able to assess the accuracy of these methods.
The SPX method is generally capable of resolving the major features of the spectral functions involved in this study.However, it encounters difficulties when dealing with spectra that exhibit a wide platform, such as the continuum model (see Sec. III D), or when two features are too close together, such as the bottomonium spectrum (see Sec. IV B).We believe that these difficulties can be partially overcome with the help of the constrained sampling algorithm and the selfadaptive sampling algorithm.The SPX method demonstrates good noise tolerance and exhibits robustness with respect to moderate noise levels.Overall, the performance of the SPX method is comparable to that of the commonly used MaxEnt method.In cases where the spectral function is complicated, the SPX method could outperform the MaxEnt method due to its ability to incorporate prior information about the spectrum into the frequency grid for the poles.This grid could be iteratively refined to obtain better spectrum.
As for the NAC method, it is found to be numerically unstable even for input Euclidean data with extremely low noise level (δ = 10 −8 ).This drawback greatly limits the application of the NAC method in the LQCD simulations.Additionally, we observe that the Hardy basis optimization for θ N s +1 sometimes produces worse results when compared to those obtained with constant θ N s +1 .Although the Hardy basis optimization can suppress possible oscillations in the spectrum, it tends to yield incorrect estimations for the positions of the high-energy peaks if the λ parameter is not optimal.Therefore, better basis functions for expanding θ N s +1 are highly desirable.Or else we need a smart algorithm to determine the optimal λ.Nonetheless, the NAC method still proves its usefulness in analytic continuations of LQCD simulation data, as it allows for quick yet accurate estimations of the positions of the low-energy band edges and the resonance peaks.These important clues can then be used to construct a reference model for the probability distribution of the poles, which is subsequently utilized by the constrained SPX method.

FIG. 1 .FIG. 2 .
FIG. 1. Analytic continuations of the Breit-Wigner models.For the NAC method, the Hardy basis function optimization algorithm is adopted.(a) Single Breit-Wigner peak.N opt = 13.(b) Two Breit-Wigner peaks.N opt = 14.In panel (b) the spectra are scaled by a factor of 0.5 for a better view.
M r = 0.10 GeV, M c = 0.05 GeV, C r = 2.0 GeV, C c = 2.10 GeV, and Γ = 0.06 GeV.The synthetic Euclidean data consist of 100 frequency points, and T = 0.02 GeV.In this study, we consider three different cases: (i) the resonance-continuum model, (ii) the resonance model, and (iii) the continuum model.The analytic continuation results are shown in Figure 3.

FIG. 3 .
FIG. 3. Analytic continuations of the resonance-continuum models.The NAC method is enhanced with the Hardy basis function optimization algorithm.(a) Spectra of the resonance-continuum model.N opt = 12.(b) Spectra of the resonance model.N opt = 10.(c) Spectra of the continuum model.N opt = 13.

FIG. 4 .
FIG.4.Analytic continuations of the bottomonium correlation function.Here, the terminology "C-SPX" implies that the positions of the poles are restricted in the SPX simulations[46,68].The Hardy basis function optimization algorithm is enabled for the NAC method and N opt = 8.A shifted Gaussian function, instead of a constant, is used as the default model for the MaxEnt method[69].See the main text for more details.

FIG. 5 .
FIG. 5. Robustness of the SPX method with respect to the noisy LQCD data.Here we just consider the Breit-Wigner model (1BW model).The noise level δ is varied from 10 −8 to 10 −2 .The other model parameters can be found in Sec.III B. (a)-(g) Dependence on noise level δ of calculated spectral functions.(h) The goodness-of-fit function χ 2 as a function of the noise level δ.The horizontal bar indicates the asymptotic value of log(⟨χ 2 ⟩).

FIG. 6 .
FIG. 6. Robustness of the NAC method with respect to the noisy LQCD data.Here we just consider the Breit-Wigner model (2BW model).The noise level δ = 10 −8 .The other model parameters can be found in Sec.III B. The darker shaded region denotes the window for the reconstructed spectral functions.The blue and yellow solid lines mean the exact spectrum and the averaged spectrum, respectively.(a) Without Hardy basis function optimization.(b) With Hardy basis function optimization.

FIG. 7 .
FIG. 7. Analytic continuations of the bottomonium correlation function.(a) Spectra obtained by the MaxEnt and optimized NAC methods.The information extracted from the two spectra is used to construct a reference model [the pseudo-spectrum ρ pseudo (ω), see the solid blue line].(b) Standard linear grid and non-uniform grid used in the C-SPX and SA-SPX simulations, respectively.Note that the non-uniform grid is constructed from the reference model via Eqs.(31) and (32).(c) Spectra obtained by the C-SPX and SA-SPX methods.(d) Error analysis for the reconstructed Euclidean data from the MaxEnt, optimized NAC, C-SPX, and SA-SPX methods.
In this expression, d and d † represent the annihilation and creation operators, respectively.|n⟩ and |m⟩ are the eigenstates of the Hamiltonian Ĥ, and E n and E m are the corresponding eigenvalues.Z = n exp(−βE n ) is the partition function.z ∈ C\R.The positive sign corresponds to fermions, while the negative sign corresponds to bosons.By introducing A mn = ⟨n|d|m⟩⟨m|d † |n⟩ e −βE n ± e −βE m /Z and P mn = E m − E n , Eq. (