Reconstruction of smeared spectral function from Euclidean correlation functions

We propose a method to reconstruct smeared spectral functions from two-point correlation functions measured on the Euclidean lattice. Arbitrary smearing function can be considered as far as it is smooth enough to allow an approximation using Chebyshev polynomials. We test the method with numerical lattice data of Charmonium correlators. The method provides a framework to compare lattice calculation with experimental data including excited state contributions without assuming quark-hadron duality.

typeset using P T P T E X.cls 1 Introduction Reconstruction of hadron spectral function from Euclidean correlation functions is a notoriously difficult problem. In lattice quantum chromodynamics (LQCD) computations, which have so far been the only practical method to calculate non-perturbative quantities with errors under control, physical quantities are extracted from n-point correlation functions obtained on an Euclidean lattice. It means that all momenta inserted are space-like, so that physical amplitudes, especially those for on-shell particles or resonances, have to be read off from this unphysical setup. The ground state contribution can be obtained relatively easily by measuring an exponential fall-off of the correlator at long distances, while excited states are much harder to identify because the exponential function of the form exp(−Et) with an energy E and a time separation t is numerically very similar for different E's especially when many energy levels are close to each other as in the experimental situation. Typically, then, only one or even none of the excited states can be identified.
Still, because the spectral function is of phenomenological interest for various applications, several groups developed methods to extract it with some extra pieces of information.
The maximum entropy method [1,2] is one of such attempts, where one assumes that unknown functions (or its parameters) are statistically equally distributed within the model space and tries to determine the most "likely" function. A slightly different statistical approach based on Bayesian statistics was also proposed [3]. Unfortunately, they are not free from uncertainties since the statistical distribution of the spectral function does not really have theoretical basis. Similar problem would also remain in a recent attempt to use machine learning to reconstruct the spectral function [4].
Another attempt to approach the problem is the use of the Backus-Gilbert method, which is a deterministic method to obtain a smeared spectral function [5]. The smearing kernel is automatically determined by the data, requiring that the width of the smearing is minimized. In practice, one has to relax the minimization by adding an extra term to the function to be minimized in order to avoid numerical instability. The smearing kernel is therefore unknown until one actually performs such analysis. There is also a proposal to arrange the Backus-Gilbert method such that the smearing kernel obtained becomes close to what one inputs [6]. Several such methods have been tested with mock data led from a known spectral function [7], and there is no obvious best solution found so far.
We propose an alternative method to reconstruct the spectral function with smearing specified by arbitrary kernels. The Chebyshev polynomials are introduced to approximate the kernel function, and a smeared spectral function under this approximation is obtained from lattice data of temporal correlator. The procedure is deterministic and the systematic 2 error due to a truncation of the Chebyshev polynomials can be estimated. Limitation of the method comes from statistical error of the lattice data, which prevents one from using high order polynomials.
The smeared spectral function can offer an intermediate quantity that can be used to compare experimental data with non-perturbative theoretical calculation provided by LQCD.
For instance, let us consider the R-ratio R(s) defined for the e + e − → qq cross section as R(s) = σ e + e − →qq (s)/σ e + e − →µ + µ − (s). It cannot be directly compared with perturbative calculations in the resonance region where perturbative expansion does not converge. The LQCD calculation as it stands is not useful either, because it can only calculate low-lying hadron spectrum and scattering phase shift for specified final states, such as π + π − or KK, but fully inclusive rate is unavailable. Our method concerns how to extract the information for inclusive processes, such as qq, from lattice results of hadron correlators without specifying any final states. The R-ratio as a function of invariant mass of the final states is not directly accessible in our method, but the function that is smeared with some kernel can be related to quantities calculable on the lattice.
The smeared spectral function was considered in the early days of perturbative QCD by Poggio, Quinn and Weinberg [8]. They considered a smearing of the form whose kernel approaches a delta function δ(s − s ) in the limit of the width ∆ s → 0. The R-ratio (or the spectral function) can be related to the vacuum polarization function through the optical theorem R(s) = (1/π)ImΠ(s) and the dispersion relation 1 The smeared spectrum can then be written using the vacuum polarization function at complex values of q 2 : The main observation was that, unlike the imaginary part ImΠ(s) evaluated on the cut, one can avoid non-perturbative kinematical region and calculate Π(s + i∆ s ) in perturbation theory as long as the smearing range ∆ s is large enough. This is the argument behind the quark-hadron duality. The width parameter ∆ s is typically chosen larger than the QCD scale Λ QCD , but there is no a priori criteria of how large ∆ s should be for perturbation theory to 1 In practice, one needs to use the subtracted version to avoid ultraviolet divergences.
3 work to a desired accuracy. One can also consider a Gaussian smearing instead of (16), and arrive at the same conclusion [9].
In many applications of perturbative QCD, such as deep inelastic scattering or inclusive hadron decays, the smearing is not as transparent as in this example. Some smearing over kinematical variables is involved depending on the setup of the problems, and the question of how much smearing is introduced is more obscure. Yet, one usually assumes that perturbation theory works; systematic error due to this duality assumption is unknown.
Another commonly used form of smearing is the Laplace transform which is related to the Borel transform involved in QCD sum rule calculations [10]. Sincẽ Π(M 2 ) can be written using the vacuum polarization function Π(Q 2 ) in the space-like domain, the non-perturbative kinematical region is avoided. The effective range of smearing is controlled by the parameter M 2 .
Of course, one can view the dispersion relation (2) for a space-like value of q 2 = −Q 2 : or its subtracted version as a sort of smearing. The range of smearing is effectively infinity as the weight function decreases only by a power of s. This is another way of keeping away from the resonance region, and perturbation theory is expected to be applicable. One still needs to include power corrections using the operator product expansion, which involves unknown parameters (or condensates); LQCD calculation is desirable to eliminate such uncertainties.
We develop a formalism of LQCD calculation that allows us to compute the quantities mentioned above, i.e. those obtained by applying some smearing on the spectral function.
The smearings are designed to escape from the resonance region to some extent, so that perturbation theory is applicable, but some uncertainty still remains as mentioned above.
By using LQCD, on the other hand, fully non-perturbative calculation can be achieved and no remnant uncertainty due to the singularities nor the power corrections remains. In other words, one can entirely avoid the assumption of quark-hadron duality.
Through this method, the comparison between experimental data and lattice calculation would provide a theoretically clean test of QCD. It can also provide a testing ground for 4 the perturbative QCD analysis including the operator product expansion against the fully non-perturbative lattice calculation.
This paper is organized as follows. In Section 2 we introduce the method to reconstruct the smeared spectral function. It includes an approximation using the Chebyshev polynomials, whose performance is demonstrated for several cases with toy examples in Section 3.
Then, the method is tested with actual lattice data of Charmonium correlator in Section 4.
The discussion section (Section 5) lists possible applications of the methods including phenomenological ones as well as those of theoretical interests related to the QCD sum rule. Our conclusions are given in Section 6.

Reconstruction of smeared spectral function
We are interested in the spectral function for some state |ψ . It does not have to be an eigenstate of HamiltonianĤ, but can be created by applying some operator on the vacuum, e.g. x J µ (x)|0 for the case of e + e − → qq. Here J µ (x) is the electromagnetic current and the sum over x gives the (spatial) zero-momentum projection. An extension to the case of different initial and final states should be possible.
The spectral functionρ(ω) in (7) is normalized such that it becomes unity when integrated over all possible energy ω.
On the lattice, we calculate the temporal correlation function which is normalized to one at zero time separation. It can be rewritten using the spectral function asC thus the conventional spectral decomposition of a correlator.
In practice, we can take the state |ψ as with some (small, but non-zero) time separation t 0 in order to avoid any potential divergence due to a contact term when evaluating ψ|ψ . For instance, by taking t 0 = 1 in the lattice unit, we can identify (8) as C(t + 2)/C(2) for a vector correlator where µ stands for a spatial direction and not summed over. We note that the HamiltonianĤ in (8) is not explicitly written in lattice QCD simulations, but we assume that it exists so that the time evolution is written by a transfer matrixẑ = e −Ĥ . We also assume that the eigenvalues ofĤ are non-negative and equivalently that the eigenvalues ofẑ are constrained to lie between 0 and 1. Strictly speaking, many lattice actions currently used in numerical simulations do not satisfy the reflection positivity, which is a necessary condition for a Hermitian Hamiltonian to exist. Any problem due to the violation of the reflection positivity is expected to disappear in the continuum limit. Our assumption is, therefore, that our lattice calculations are sufficiently close to the continuum limit. Using the transfer matrixẑ, the correlator in (8) is simply written as For the vacuum polarization function due to electromagnetic currents, the spectral function is often defined as ρ(s) = (1/π)ImΠ(s). The Euclidean correlator is then expressed as (See [11], for instance.) This is slightly different from (8), but can be related by redefining the spectral function. Namely, we defineC(t) asC(t) ≡ C(t + 2t 0 )/C(2t 0 ), so that the spectral functionρ(ω) in (9) isρ(ω) = (1/C(2t 0 ))ω 2 ρ(ω 2 )e −2ωt 0 . We note that the normalization factor C(2t 0 ) is explicitly calculable on the lattice. Now we define a smeared spectral functionρ ∆ (ω) for a smearing kernel S ∆ (ω, ω ) as Then, the matrix element to be evaluated is ψ|S ∆ (ω,Ĥ)|ψ . The form of the smearing kernel S ∆ (ω, ω ) is arbitrary; we can consider the choices discussed in the previous section.
In the following, to be explicit, let us assume a specific form for S ∆ (ω, ω ) as where ∆ represents a range of smearing.
We consider a polynomial approximation of S ∆ (ω,Ĥ) of the form where T j (x) stands for the Chebyshev polynomial and the sum is up to its maximal order N . The first few terms of the Chebyshev polynomials are T 0 (x) = 1, to construct the following ones. According to the general formula of the Chebyshev approximation, the coefficients c j (ω) in (17) can be obtained as where the function f ω (x) is written as Here, x corresponds to eigenvalues ofẑ = e −Ĥ , so that the function f ω (x) represents the smearing function (16) with ω replaced byĤ. This standard formula for the Chebyshev approximation is written for a function f ω (x) defined in [−1, 1]. Here we use it only between [0, 1] and assume that f ω (x) vanishes for x ≤ 0. Technically, the numerical integral (18) becomes unstable for large j's due to a divergence of the integrand as x → 1. Instead, one may use an alternative formula evaluation of which is more stable for large j's. For the range of x between [0, 1], this integral is up to π/2.
To optimize the Chebyshev approximation, one can use a modified form written in terms The corresponding formula for the coefficients appearing in the Chebyshev approximation is The approximation formula (17) is unchanged other than replacing c j (ω)T j (ẑ) by c * j (ω)T * j (ẑ). Since the range of x is narrower, this series gives a better approximation of the original function for a given order N . 7 Finally, remember that the matrix elements of the transfer matrixẑ and its powerẑ t can be written asC(t) = ψ|ẑ t |ψ / ψ|ψ . Then, we arrive at an expression forρ ∆ (ω): where the last term T * j (ẑ) may be constructed from the correlatorC(t) by replacing the power of the transfer matrixẑ t appearing in T * j (ẑ) byC(t) = C(t + 2t 0 )/C(2t 0 ). Therefore, the first few terms are obtained as which is a deterministic procedure.
The general expression (22) for an approximation ofρ ∆ (ω) is valid for any smearing kernel and for any value of ω, as long as the coefficients c * j (ω) are calculated appropriately. As it is well known, the Chebyshev approximation provides the best approximation of any function defined in 0 ≤ x ≤ 1. It is the best among any polynomials at a given order N in the sense that the minmax error, the maximum deviation from the true function in the same range, is minimum; in order to achieve a better approximation, one needs a higher polynomial order N . Since the (shifted) Chebyshev polynomials T * j (x) are oscillating functions between 0 and 1, it is necessary to use larger N in order to better approximate detailed shape of the original spectral functionρ(ω) by narrowing the width ∆ of the smearing kernel. The approximation is demonstrated in the next section by taking a few examples.
The shifted Chebyshev approximation works only when the argument x is in 0 ≤ x ≤ 1.
In our case, it corresponds to the condition that the eigenvalues ofẑ are in [0, 1], which should be satisfied becauseẑ is the transfer matrixẑ = e −Ĥ . For a given eigenvalue z i ofẑ, each polynomial T * j (z i ) takes a value between −1 and 1, and if the state |ψ is decomposed as |ψ = i a i |i with a normalization i |a i | 2 = ψ|ψ , the individual polynomial becomes This provides a non-trivial constraint that must be satisfied by the correlator C(t).

3 Chebyshev polynomial approximation: examples
First, we demonstrate how well the smearing kernel S ∆ (ω, ω ), Eq. (16), is approximated by the Chebyshev polynomials. Setting ω = ω 0 = 1, in some unit, say the lattice unit, we draw a curve of S ∆ (ω, ω 0 ) in Fig. 1 (left). The Chebyshev approximation of the form (17), replacingẑ by e −ω 0 in this equation, is also plotted for N = 10 (dotted), 15 (dot-dashed) and 20 (dashed curve). On the right panels, we plot the error of the approximation, namely a difference from the true function S From the plots one can confirm that the smearing function with a larger width ∆ = 0.3 is well approximated by a limited order of the polynomials. Namely, the polynomials up to order N = 20 or even 15 give a nearly perfect approximation; the deviation is at a few per cent level. Apparently, the approximation becomes poorer when the function is sharper, ∆ = 0.2 or 0.1. One needs higher order polynomials to achieve better approximation. We limit ourselves to N = 10-20, because these are the orders that can be practically used for the analysis of lattice data, as we discuss in the next section.
When the target energy ω 0 is lower, ω 0 = 0.5, we observe a very similar pattern as shown in Fig. 2. An important difference is, however, that the approximation is better than those for ω 0 = 1.0, as one can see by comparing the size of the error S In order to see the rate of convergence of the Chebyshev approximation, we plot the coefficients c * j (ω) as a function of j in Fig. 3. As an example, the point ω = ω 0 is taken because this is where the error is largest. (It is not always the case especially when the approximation is already good. See Figs. 1 and 2.) One can see that the magnitude of the coefficient |c * j (ω)| decreases roughly exponentially as j. When the approximation is better (larger ∆), the decrease of |c * j (ω)| is faster. Since the Chebyshev polynomial | T * j (ẑ) | is bounded from above by 1, this shows (the upper limit of) the rate of convergence. The Calculation of c * j (ω) is numerically inexpensive, and one can easily estimate the error of the approximation due to a truncation at the order N by ±|c * N +1 (ω)|. As another test, we consider the Laplace transform (4), which is achieved by a smearing kernel   used for an extraction of the charm quark mass from the Charmonium temporal moments [12]. (The same lattice ensembles are also used for calculations of the Dirac spectrum [13,14], short-distance current correlators [15], topological susceptibility [16], η meson mass [17].) Among 15 ensembles generated at various lattice spacings and lattice sizes, we use the one at a lattice spacing a = 0.080 fm, a lattice size 32 3 × 64, and bare quark masses am ud = 0.007 and am s = 0.040. The corresponding pion mass is 309(1) MeV. We take 100 gauge configurations and calculate the Charmonium correlator with a tuned charm quark mass am c = 0.44037, and compute charm quark propagators from Z 2 noises distributed over a time slice to construct Charmonium correlators. We then construct the Charmonium correlators, which correspond to those of local currents in the pseudo-scalar (PP) and vector (VV) channels.
This calculation has been repeated for 8 source time slices to improve statistical signal, so that the total number of measurements is 800. Three spatial polarizations are averaged for the VV channel.   We construct the Chebyshev matrix elements T * j (ẑ) defined in (22) from the current correlators. The calculation is straightforward. Namely, we replace the term ofẑ t in the polynomials by C(t + 2t 0 )/C(2t 0 ) and calculate the linear combinations of them with the coefficients of the shifted Chebyshev polynomials. We take t 0 = 1 in the lattice unit. The results are shown in Fig. 6; the error is calculated using the jackknife method. It turned out that the Chebyshev matrix elements are precisely determined up to j = 11, which corresponds to t = 13. Beyond that point, the statistical error grows rapidly, and the results eventually get out of the range of ±1, which must be satisfied for the Chebyshev polynomials. Statistical error of T * j (ẑ) as a function of j. Circles and squares represent the data for PP and VV channels, respectively.
In fact, the statistical error grows exponentially for higher-order polynomials as shown in Fig. 7. Its growth rate is about a factor of three to proceed by another order in j, which means that 10 times larger statistical samples would be needed to include yet another order to improve the Chebyshev approximation. This is not unreasonable because we try to construct the quantities of O(1) as a linear combination of terms of exponentially different orders. For the Charmonium correlator at the lattice spacing chosen in this work, the terms ofẑ t are suppressed roughly by e −1.2t , which is however compensated by the Chebyshev coefficients growing even faster. In the end, a strong cancellation among different powers ofẑ gives a number between −1 and +1, and the noise is relatively enhanced. For instance, at the order j = 12, a cancellation of four orders of magnitudes takes place and it becomes even harder for higher orders.
Since the Chebyshev approximation drastically fails outside the domain 0 ≤ x ≤ 1, we are not able to use it beyond the order where | T * j (ẑ) | exceeds 1. Instead, we introduce a fit to determine these matrix elementsT j ≡ T * j (ẑ) in such a way that they are consistent withC(t) while satisfying a constraint |T j | ≤ 1. To do so, we can use the reverse formula of the shifted Chebyshev polynomials [18]:  where the prime on the sum indicates that the term of r = 0 is to be halved. The relation to be satisfied is thenC We take T * r 's as free parameters to be determined and fit the lattice dataC(t) with a constraint |T j | ≤ 1. Statistical correlations of C(t) among different t's are taken into account in the fit using the least-squares fit package lsqfit by Lepage [19]. The resulting values ofT j are listed in Table 2. Since the numbers of inputs,C(t)'s, and unknowns,T j 's, are the same, the condition (26) may be solved as a system of linear equations unless the constraints are introduced. In fact, the results are unchanged from the direct determination through (23) within the statistical error for small j's up to j 10. Beyond that, they are affected by the constraints. The error becomes large, of order 1, for large j's, where they are essentially undetermined by the fit but still kept within ±1.
Once a set of estimates for T * j (ẑ) , i.e.T j , is obtained, the remaining task is to use (22) to estimateρ ∆ (ω). The results are shown in Fig. 8. Three bands corresponding to the polynomial order N = 12 (red), 14 (blue), 16 (orange) are overlaid. We find that the overall shape is unchanged by adding more terms, i.e. from 12 to 14 or to 16, while the size and shape of the statistical error is affected. When the polynomial order is lower, some wiggle structure is observed, while the necks, the positions of small statistical error, are fatten by adding more terms and eventually the error becomes nearly uniform over ω. Beyond N = 16, the results are essentially unchanged, since the higher order coefficients c * j (ω) are exponentially suppressed.
In Fig. 9 we show the results for the smeared spectral functionρ ∆ (ω) obtained with N = 16. The smearing kernel is S ∆ (ω, ω ) with ∆ = 0.1 (top), 0.2 (middle) and 0.3 (bottom). The results are compared with the expected contributions from the ground state and first excited state (red and blue curves). They are drawn assuming δ-function distributions,ρ(ω ) = i A i e −ω t 0 δ(ω − E i ), with the fitted values of energy levels E i and their amplitudes A i given in Table 1. The factor e −ω t 0 is introduced to take account of the time evolution from 0 to t 0 , which is included in the definition of the state |ψ = e −Ht 0 J µ |0 .
When the smearing width is large, ∆ = 0.3 (bottom panels of Fig. 9), we observe that the reconstructed smeared spectral function follows the expected form from the low-lying states in the lower region of ω. As ω increases, the spectral function indicates more contributions from higher excited states. This is exactly what we expected. The lattice data in the short time separations contain the information of such states, which is properly extracted with our method. From perturbation theory, one expects a constant proportional to the number of color degrees of freedom N c = 3 for the (unsmeared) spectral function ρ(ω). This constant is slightly distorted because of the difference between ρ(ω) andρ(ω). For t 0 = 1 (in the lattice unit), this should give an exponentially decreasing spectral functionρ ∆ (ω) as ∼ ω 2 e −2ω at large ω, which is indeed observed in the results. On the other hand, the resonance structure This error due to the truncation is not included in the band shown in Fig. 9, but taking account of this marginal size of error the reconstructed ρ ∆ (ω) looks reasonable for ∆ = 0.2 (middle panels). Namely, it follows the expected curve of the ground and the first excited states up to around the peak of the latter and then drops slowly due to higher excited state contributions.
The truncation error increases to 50-100% for ∆ = 0.1 (upper panels of Figures 1 and  2), and we should not take the results (top panels of Fig. 9) too seriously. If it were precisely calculated, we would be able to resolve the resonance structures as the curve of the lowlying state contributions suggests. To do so, we need to include higher order terms of the Chebyshev approximation, which requires much better precision of the simulation data. This reflects the fact that the reconstruction of the full spectral function from Euclidean lattice data is an ill-posed problem. One needs ridiculously high precision in order to achieve a full reconstruction as emphasized in [6].

Discussions
As already mentioned earlier, our proposal to calculate the smeared spectral function is not limited to the case just discussed. Any sort of weighted integral of the spectral function can be considered. A well-known example is the contribution of quark vacuum polarization to  Results are shown by orange band, while the ground-state contribution (dot-dashed) and the ground-state and plus excited-state contribution (dashed) are plotted assuming they have δ-function structures. the muon anomalous magnetic moment g − 2. Phenomenologically, one employs the optical theorem and dispersion relation to relate the vacuum polarization function in the Euclidean domain Π(Q 2 ) to a weighted integral of the experimentally observed R-ratio, or the spectral function. Then, an integral of Π(Q 2 ) with an appropriate weight gives the contribution to g − 2. In this case, however, a direct expression in terms of the time correlator C(t) is known [11], and we do not really need the approximation method developed in this work.
Another phenomenologically interesting example is the hadronic τ decay. Using the finiteenergy sum rule [20] the hadronic width can be written as where S EW is a short-distance electroweak correction and |V ud | is a CKM matrix element.
(Here, only the ud contribution is considered. An extension to the us contribution is straightforward.) The spectral function ρ V +A (s) denotes a sum of V V and AA channels. The integral (27) reflects a particular kinematics of τ decay and has a complicated form, but our method can be applied for such a case in principle. A practical question would be, however, whether a good enough approximation can be achieved with a limited number of terms.
The Laplace transform (24) is often considered in QCD sum rule analyses [10] because the corresponding Borel transform of perturbative series makes it more convergent. Our method allows to calculate two-point function after the Borel transform directly using lattice QCD.
It may be useful to test the perturbative expansion and the operator product expansion involved in the QCD sum rule calculations. Conversely, it can also be used to validate lattice QCD calculations especially in the short-distance region. Such a test has been performed using short-distance correlators in the coordinate space [15], and it is interesting to do the test for the Borel transformed quantities.
The dispersion integral of the form (2) is of course another type of applications of our method. Although the vacuum polarization function Π(Q 2 ) in the space-like momenta Q 2 can be directly obtained by a Fourier transform of the lattice correlators, the Chebyshev approximation offers a method to extract it at arbitrary values of Q 2 corresponding to the momenta of non-integer multiples of 2π/L. The Laplace transform as introduced in [21] can do this too, but requires the information of large time separations in the integral of the form ∞ 0 dt e ωt C(t). The method developed in this work accesses only to relatively short time separations, but we need to introduce an approximation. Systematic error in each case has to be carefully examined.
Extension to the cases of more complicated quantities, such as the nucleon structure function as measured in deep inelastic scattering or the inclusive hadron decays, can also 20 be considered. One of the authors has proposed an analysis to use the dispersion integral to relate the inclusive decay rate to an amplitude in space-like momenta [22]. This method contains a difficulty of requiring the information at unphysical momentum regions, which may be avoided by more flexible integral transformation proposed in this work.
For this class of applications, there are two independent kinematical variables, q 2 and p · q, with p the momentum of a decaying particle (or the initial nucleon) and q the momentum transfer. In order to apply the method outlined in this work, we need to fix one of these kinematical variables and introduce a smearing on the other variable. More complicated integral might be useful for b → u ν decay analysis for which one introduces elaborate kinematical cuts in order to avoid backgrounds from b → c ν. The flexibility of our method would allow such analyses.
Going to high tempertature QCD, our method would not work as it is, because the correlator can not be simply written as ψ|ẑ t |ψ due to the contribution from the opposite time direction. The operator is then no longer a power of the transfer matrix, direct estimate of the Chebyshev matrix elements is not available.

Conclusions
Precise reconstruction of the spectral function from lattice data remains a difficult problem. Instead, we calculate a smeared counterpart, which contains some information of the spectral function after smearing out its detailed structures. For the Charmonium spectral function we obtain a reasonably precise result for a smearing width ∆ = 0.3, which is about 700 MeV in the physical unit. Since the mass splittings experimentally observed are narrower than this, we are not able to resolve the details of the spectrum. Taking the limit of small smearing width requires exponentially better statistical precision, and it comes back to the original problem. Still, our proposal has advantages compared to previously available methods.
In contrast to the Bayesian approach [3] or the maximum entropy method [1,2], our method allows a reliable estimate of the systematic errors, since it does not assume any statistical distribution of unknown function. In principle, the method is deterministic once the input lattice data are given. The least-squared fit involved in order to enforce the constraint that eigenvalues of Hamiltonian is positive plays only a minor role that becomes irrelevant when the lattice data are made precise.
Compared to the Backus-Gilbert method [5], the method proposed in this paper is more flexible as it allows any predefined smearing function, while it is automatically determined in the Backus-Gilbert method and therefore is uncontrollable. The variant of the Backus-Gilbert method [6] also has this flexibility. Our method also allows systematic improvements since the approximation is achieved by a series of exponentially decreasing coefficients. As the statistical precision of the input correlator is improved, one can include higher order terms and thus improve the approximation.
The method can be used in the analysis of inclusive processes to define intermediate quantities for which fully non-perturbative lattice calculation is possible. Its potential application is not limited to the spectral function for two-point correlators, but other processes such as deep inelastic scattering and inclusive B meson decays can be considered. Since the method does not rely on perturbation theory, the processes with small momentum transfer can be calculated on a solid theoretical ground, which was not available so far. More importantly, we do not have to rely on the assumption of quark-hadron duality.