Spectral sum of current correlators from lattice QCD

We propose a method to use lattice QCD to compute the Borel transform of the vacuum polarization function appearing in the Shifman-Vainshtein-Zakharov QCD sum rule. We construct the spectral sum corresponding to the Borel transform from two-point functions computed on the Euclidean lattice. As a proof of principle, we compute the $s \bar{s}$ correlators at three lattice spacings and take the continuum limit. We confirm that the method yields results that are consistent with the operator product expansion in the large Borel mass region. The method provides a ground on which the OPE analyses can be directly compared with nonperturbative lattice computations.


I. INTRODUCTION
The spectral sum of hadronic correlation functions, such as the vacuum polarization function Π(q 2 ), of the form ds e −s/M 2 Im Π(s) (1) has often been introduced since the seminal work of Shifman, Vainshtein, and Zakharov [1,2]. The integral over invariant mass squared s smears out contributions of individual resonances so that one can use perturbative treatment of quantum chromodynamics (QCD) with quarks and gluons as fundamental degrees of freedom, as far as the Borel mass M , a parameter to control the typical energy scale, is sufficiently large. The integral (1) is a quantity effectively defined in the spacelike momentum region, and there would be no issue of the violation of the quark-hadron duality [3].
The integral (1) suppresses the contributions from the energy region above M and thus, is more sensitive to low-lying hadronic states. If one can find a window where M 2 is large enough to use perturbative expansion of QCD with nonperturbative corrections included by operator product expansion (OPE) and at the same time sufficiently small to be sensitive to lowest-lying hadronic states, the spectral sum (1) may be used to obtain constraints on the parameters of low-lying hadrons, such as their masses and decay constants. This method, called the QCD sum rule, has been widely applied to estimate masses, decay constants, and other properties of hadronic states in various channels [1,2]. However, an important question of how well the perturbative QCD with some nonperturbative corrections included through OPE can represent the spectral sum is yet to be addressed, especially when the correlation function is not always fully available from the experimental data, e.g. due to a limitation of accessible kinematical region.
In principle, the test of perturbative expansion and OPE can be performed using nonperturbatively calculated correlation functions using lattice QCD. Comparison of the lattice correlators at short distances with perturbative QCD may be found, e.g., in [4][5][6][7] for lighthadron current-current correlators and in [8,9] for charmonium correlators. The energy scale where the comparison is made has to be sufficiently low to avoid discretization effects in the lattice calculations, while the OPE analysis is more reliable at high energy scales. It has been pointed out that the convergence of OPE is a crucial problem in the energy region for which lattice QCD can provide reliable calculations by now [7,10].
In this work, we perform another test of perturbative QCD and OPE against nonperturbative lattice computation using the spectral sum of the form (1). It has an advantage that the OPE converges more rapidly compared to that applied for the correlator itself either in the coordinate space or in the momentum space. And, this is exactly the quantity that has been used in many QCD sum rule analyses; hence, it serves as a test of those sum rule calculations as well.
On the lattice, computation of the spectral sum (1) is highly nontrivial because it requires the spectral function ρ(q 2 ) ∝ Im Π(q 2 ) for all values of timelike q 2 above the threshold where a cut begins. Extraction of the spectral function from the lattice correlators is a notoriously difficult problem that requires solving an ill-posed inverse problem. Namely, one has to extract ρ(q 2 ) by solving with a lattice correlator C(t) of a current operator J calculated at a discrete set of time separations. There have been several methods developed to perform this inverse-Laplace transform, including the maximum entropy method (MEM) [11][12][13], Bayesian approach [14], Backus-Gilbert approach [15][16][17][18], the sparse modeling method [19], but none of them succeeded to achieve sufficiently precise extraction of ρ(ω 2 ) that can be used for the purpose of this work.
In this work, instead, we apply the method developed in [20]. It is based on a representation of the weight function e −ω 2 /M 2 in (1) as a polynomial of e −aω , which is then related to the transfer matrix e −aĤ defined on the lattice. (Here, a stands for the lattice spacing.) The method relates the spectral sum directly to the lattice correlators without explicitly solving the spectral function ρ(ω 2 ), so that the inverse-Laplace transformation can be avoided. The method has so far been applied to the B meson inclusive semileptonic decays [21] as well as the inelastic lepton-nucleon scatterings [22]. As we demonstrate in the next sections, the method allows us to construct the spectral sum with small and controlled systematic errors.
This paper is organized as follows. In Sec. II we introduce the spectral sum for the Borel transform in the continuum theory. We also introduce our lattice QCD setup for the evaluation in Sec. III. We discuss lattice calculations and their errors in Sec. IV. We show comparison with OPE and the ground state contribution in Sec. V. Section VI is devoted to our conclusion and outlook.

II. CURRENT CORRELATORS IN QCD AND THEIR SPECTRAL SUM
We briefly review the use of the spectral sum of QCD current correlators. More detailed reviews and discussions are found in the literature, e.g. [23,24].
We define the hadronic vacuum polarization (HVP) function as a Fourier transform of the current-current correlator, where J µ =qγ µ q is the quark vector current. Taking account of its analytical property, the HVP may be written in terms of a spectral function ρ(s), where Q 2 is the momentum squared, Q 2 = −q 2 . This integral diverges since the spectral function does not vanish in the limit s → ∞, and we can remove the divergence by subtracting, for instance, Π(q 2 0 ) at a certain point q 2 = q 2 0 , and define a subtracted HVP. In the QCD sum rule analyses, one introduces the Borel transform of HVP to enhance the contributions from low-lying hadronic states. The Borel transformation is defined as where M is the Borel mass that specifies a typical energy scale. The Borel transform of HVP may then be written as The exponential factor e −s/M 2 suppresses the contributions from high-energy states above M .
One can use OPE to evaluateΠ(M 2 ) including nonperturbative power corrections. We start from an expression of Π(−Q 2 ) as an expansion in 1/Q 2 , where α s (µ 2 ) is the strong coupling constant defined at a renormalization scale µ, m is the quark mass, and 0| αs π G 2 |0 and 0|qq|0 are the gluon and chiral condensates, respectively. Here, the four-quark condensate is represented by a vacuum saturation approximation (VSA) with a parameter κ 0 , which describes the violation of VSA when κ 0 = 1. By the Borel transformation, the logarithmic function and negative powers of Q 2 are transformed as where n is a positive integer. Therefore, the Borel transform of HVP can be expressed as follows: The perturbative coefficients of the leading order term, O(1/M 0 ), in the massless limit are known up to O(α 4 s ) [25], where the disconnected diagrams are neglected. The other corrections taken into account in this paper are summarized in Sec. V. Because of the factor 1/(n − 1)! in (10), the Borel transform is less affected by higher dimensional condensates, and the OPE is made more convergent than that for HVP itself (8).
Perturbative expansion ofΠ OPE (M 2 ) in the massless limit shows a good convergence. We set the renormalization scale µ 2 to M 2 e −γ E since the Borel transformation of the logarithmic function B M [log n (µ 2 /Q 2 )] appears as a polynomial of log(µ 2 /M 2 e −γ E ). (See Appendix A.) We showΠ pert 0 (M 2 ), which is the leading order of the 1/M 2 expansion, as a function of 1/M 2 in Fig. 1. We set Λ (n f =3) MS = 332 MeV for the coupling constant α s (µ 2 ). The running of α s (µ 2 ) is incorporated at five-loop level using RunDec [26,27]. Figure 1 indicates that the truncation error of the perturbative expansionΠ pert 0 (M 2 ) is not substantial for M > 1 GeV. Indeed, the O(α 4 s ) correction is at the level of 0.3% or smaller. For the next-to-leading order terms of OPE, i.e. the terms of m 2 /Q 2 , the perturbative coefficients are known to α 3 s [28],  where the renormalization scale µ is set at µ 2 = Q 2 and n f = 3. The numerical expressions for different n f 's are found, e.g., in [29]. We define the Borel transform of the correctioñ Applying the formula in (A2) and setting µ 2 = M 2 e −γ E , we found the expression, We plotΠ pert m 2 (M 2 ) in Fig. 2 (top). UnlikeΠ pert 0 (M 2 ), we observe significant dependence on the order of the perturbative expansion. To improve the convergence, we set the renormalization scale at µ 2 = 4M 2 e −γ E as shown in Fig. 2 (middle). The dependence on the scale µ is demonstrated in Fig. 2 (bottom), where the perturbative expansion truncated at should be independent of the renormalization scale up to truncation errors, we treat the variation due to the unphysical scale setting as the truncation error in the later sections.
In phenomenological studies, an ansatz for the spectral function of the form is often used. Here, m V and f V are a mass and a decay constant of the ground-state hadron, respectively. Excited states of hadrons are modeled by the continuum (or scattering) states calculated in perturbative QCD, and the spectral function of the continuum states ρ cont (s) is introduced above the threshold s th . This replacement amounts to assume the quark-hadron duality. The Borel transformation reduces the dependence on this assumption. The integral in (7) with ρ ph (s) corresponds to the OPE expression in (11). Namely, is used in the QCD sum rule analysis. Solving this equation for m V and f V , one can predict the mass and decay constant of this particular channel from the fundamental parameters of QCD, such as α s (µ 2 ), m as well as the condensates.
The QCD sum rule for the φ meson, which we mainly study in this work, is discussed in the literature, e.g., [2,30].

III. BOREL TRANSFORM OF THE SPECTRAL FUNCTION
We compute the Borel transformΠ(M 2 ) using lattice QCD. The weighted integral of the spectral function of the form (7) can be interpreted as a smeared spectral function. To compute the smeared spectrum in lattice QCD, we use the method proposed in [20], which is based on the expansion of the smearing kernel in terms of the transfer matrix on the lattice.
The method relates the smeared spectrum to the correlators computed on the lattice via the spectral representation. Applications to the inclusiveB s decay [21] and the inelastic lN scattering [22] have been discussed. We briefly review the key ideas of this method in the following. In this section, all parameters are in the unit of the lattice spacing a, unless otherwise stated.
We consider a current-current correlator with zero spatial momentum where J z stands for the z component of the vector current. Computation of such correlators as a function of the time separation t is straightforward in lattice QCD. The relation between the correlator and the spectral function is given by [31], We recall that ρ(ω 2 ) is defined in (5). Here, we make a change of variable ω = √ s. Estimation of the spectral function ρ(ω 2 ) from (17) is an ill-posed inverse problem because the functions e −ωt with different ω's are hard to distinguish numerically when ω's are close to each other. To avoid this problem, the method of [20] relates the correlator to the smeared spectral function such as (7), instead of the spectral function ρ(ω 2 ) itself.
We define the spectral density for a state |ψ , whereĤ is the Hamiltonian. The spectral densityρ(ω) evaluates the number of states having an energy ω. Setting |ψ = e −Ĥt 0 x J z (0, x)|0 , the Laplace transform of the spectral density may be written in terms of the correlators, Here, we introduce a small-time separation t 0 > 0 to avoid the contact term that potentially diverges at t 0 = 0. In this paper, we set t 0 = 1 not to lose high energy state contributions too much. The correlatorC(t) is normalized asC(0) = 1.
Let us now consider a smeared spectral function, with a smearing kernel S(ω), which will be specified later. One may approximate the smearing kernel in terms of the shifted Chebyshev polynomials T * j of e −ω , where N t stands for the truncation order of the approximation. The explicit form of the poly- . Note that the Chebyshev approximation is an orthogonal expansion and we do not impose any condition such as the one that e −ω being small for its convergence. We substitute this expression to (20). Then the smeared spectral function is written in terms of the transfer matrix e −Ĥ as where Here we replaced ω byĤ when sandwiched by the states ψ| and |ψ , and performed the integral over ω in (20). We can write the expectation value T * j (e −Ĥ ) using the correlators as where we use (16) and (18).
In practice, the lattice correlators contain statistical errors. Since (25) involves cancellations ofC(t) with different t's, the resulting expectation values T * j (e −Ĥ ) may induce large statistical errors. In particular, since we have to take an additional constraint [20], the statistical error causes a significant problem. We therefore determine T * j (j = 1, · · · , N t ) through a fit of correlators. Since T * Nt (x) includes terms up to x Nt , the data ofC(t) at t = 0-N t are used in the fit. Now we turn to the discussion of the Borel transform. The relation betweenρ(ω) and ρ(ω 2 ) is found as [see (17) and (19)] We therefore set S(ω) to be S(M, ω) as a function of the Borel mass M as to obtain the Borel transform as a smeared spectral function, where we change the variable as s = ω 2 . The smearing kernel (27) has an apparent problem of divergence at ω = 0, which induces divergences of the coefficients c * j (22). We therefore introduce a cutoff to regularize the integral (22). Since the spectrum ρ(s) vanishes below the energy of the lowest-lying state, any modification of the kernel below the lowest energy does not affect the final result. We therefore modify the smearing kernel, where ω 0 is set smaller than the mass of the ground state. The form of S cut (M, ω) is shown in Fig. 3. With ω 0 not much smaller than the lowest hadronic state, the modified smearing underestimates the smeared spectrum. In this work, we consider the ss states, for which the lowest energy state is the φ meson, whose mass is ∼ 1 GeV. We will discuss how the error due to the modified smearing can be corrected.
To summarize, we obtain the approximation between the smeared spectral function and where c * j (M ) is evaluated as (22) with S(ω) = S cut (M, ω).

IV. LATTICE CALCULATION
We compute two-point correlators of the vector current J µ =sγ µ s using lattice QCD.
In this work, we neglect the disconnected diagrams. We use ensembles with N f = 2 + 1 dynamical Möbius domain-wall fermions [32], where the gauge action is tree-level Symanzik improved. Parameters of the ensembles are listed in Table I Dirac eigenvalues [33], charmonium moments [9], short distance current-current correlators [5], topological susceptibility [34], and η meson mass [35]. Other details of the ensembles are available in [36,37].
We compute the Borel transform of the HVP using the technique outlined in the previous section. The estimate for the Chebyshev matrix elements T * j in (31) is obtained by a fit of lattice correlators. The fit is implemented using lsqfit [38], which is based on Bayesian statistics [39]. Following [20], we write the correlator at each temporal separation by the Chebyshev matrix elements as using the reverse formula of the shifted Chebyshev polynomials, The Chebyshev matrix elements T * j are determined such that they best reproduceC(t) under the given statistical error while satisfying the necessary condition | T * j | ≤ 1. Combining them with the coefficients c * j (M ), we obtainΠ cut (M 2 ) through (31). In order to match the lattice results with the counterpart in the MS scheme, the renormalization factor has to be multiplied. We use the renormalization constants of the vector current Z V = 0.955(9), 0.964 (6) In the following subsections, we discuss potential systematic effects due to the truncation of the Chebyshev expansion, the effect of the low-energy cut introduced in the smearing function, and the continuum extrapolation. to the long-distance correlator. We expect that higher-order polynomials are needed when the lattice spacing is small.
The truncation error can also be estimated through the coefficients c * j (M ) in (22) because T * j is bounded as | T * j | ≤ 1. In Fig. 7, we show the absolute values of the coefficients at various M 2 s at each lattice spacing. The plots demonstrate that the coefficients decrease exponentially for large j. When the scale M is large, the coefficient c * j (M ) drops more rapidly for high orders (larger j's). It implies that j c * j (M ) T * j is dominated by the lowerorder terms, which corresponds to shorter-distance correlators. At 1/M 2 ∼ 2 GeV −2 which corresponds to the lowest scale treated in this work, the coefficient c * j (M 2 ) is sufficiently small [∼ O(10 −4 )] already at j = 18 . We therefore set N t = 18 in the following.
In order to have another insight into the possible error due to the Chebyshev approximation, let us consider a simple model that has a single pole, with massm and decay constantf . The corresponding Euclidean correlator is and the normalized correlator (19) is given bȳ In this test, we ignore statistical errors and replace the expectation values T * j by the shifted Chebyshev polynomials T * j (e −am ) without introducing the fit. Combining the polynomials and the coefficients c * j (M ) determined by (22) with the smearing kernel S cut (M, ω), we obtain the Borel transformΠ pole (M 2 ). We can also analytically calculate the Borel transform of The difference between the original function and that with the cutoff remains when the pole mass is small,m = 1 GeV. We correct them as discussed in the following.

B. Correction for the low-energy cut of smearing function
The low-energy cut tanh(ω/ω 0 ) introduced to avoid the artificial divergence of the Chebyshev coefficients modifies the shape of the smearing kernel below ω ω 0 . If we set ω 0 sufficiently small, only the contribution from the ground state, i.e. the φ meson in our example, is significantly affected. We therefore correct for the error by using the information available for the ground state.
The contribution of the ground state ρ φ (s) for the spectral function is where f φ and m φ are the decay constant and the mass of the φ meson, respectively. The φ meson's contribution to the Borel transform is theñ Taking the limit ω 0 → 0, it recovers the physical result The difference between the Borel transform with and without the modification is then which we add back to the result ofΠ cut (M 2 ) as The deficit δΠ cut φ (M 2 ) can be computed using the values of f φ and m φ obtained for each lattice ensemble.
We show a typical threshold ω 0 dependence ofΠ lat (M 2 ) at certain values of M 2 in Fig. 9. Squares and circles denote theΠ cut (M 2 ) andΠ lat (M 2 ), respectively. As ω 0 increases,Π cut (M 2 ) decreases, as we expected. After the correction,Π lat (M 2 ) is insensitive to ω 0 . On the fine lattice, the small value of ω 0 enhances the statistical errors. To avoid large errors, we set ω 0 = 0.6 GeV for all lattice spacings in the following results. The error due to the low-energy modification is negligible after correcting for the ground state contribution.

C. Continuum limit
We take 50 points of 1/M 2 in the range 1/M 2 = 0.05-2.05 GeV −2 and computeΠ lat (M 2 ) for each lattice spacing. The results are shown in Fig. 10. We find that the results obtained at two coarser lattice spacing agree well except in the region 1/M 2 0.2 GeV −2 , where discretization effects are visible. The data at finest lattice spacing show a slight deviation from those at two coarser lattices, but we note that the strange quark mass is slightly mistuned on this ensemble and we have to correct that effect (see below.) We take the continuum limit ofΠ lat (M 2 ) using the data at three lattice spacings. Since both the statistical and systematic errors correlate highly among different values of 1/M 2 , we introduce an ansatz, with coefficients b 0 and b 1 to parametrize the discretization effect independent of 1/M 2 . We introduce a correction δΠ m to incorporate the mistuning of the valence quark mass m s . At tree level, the correction δΠ m is expressed as where m phys (µ 2 ) and m siml (µ 2 ) are the strange quark masses at the scale µ. We take   [34,41,[49][50][51][52][53], and the gluon condensate 0| αs π G 2 |0 = 0.0120(36) GeV 4 [1,2] (adding ±30% error).
In the calculation of the perturbative expansion and OPE, we set the renormalization scale µ 2 = 4M 2 e −γ E . The running of α s , m s , and 0|qq|0 are taken into account using RunDec [26,27] at five-loop level.
In OPE we include corrections up to mass-dimension six operators: where c 0 and c 2 stand for the perturbative expansion in the massless limit and the leading mass correction, respectively. The coefficient c 4 includes the gluon and quark condensates.
The coefficients c 0 and c 2 are already discussed in Sec. II. The coefficients c 4 and c 6 can 1 We use the chiral condensate evaluated in the massless quark limit, rather than the "strange quark condensate," which has also be evaluated using lattice QCD [47] as ss (2 GeV) = −(296(11) MeV) 3 .
The reason is that the difference from the massless limit involves a quadratic divergence and a renormalon ambiguity of order of m s Λ 2 QCD , which is the same order of the correction itself, is induced when the divergence is subtracted. In [47], the subtraction scheme is not explicitly shown, and in [48]  be computed by applying Eq. (A2) to the Wilson coefficients (see also [54]). Letting L M ≡ log(µ 2 e γ E /M 2 ), we can express the coefficients as where the gluon condensate 0| αs π G 2 |0 is defined in the MS scheme. The coefficient κ 0 in (47) parametrizes corrections to the VSA for the four-quark condensate. When the condensate is assumed to be fully factorized in the vacuum, κ 0 is equal to 1. There are studies that suggest the violation of VSA as large as κ 0 ∼ 6 [55]. We set κ 0 = 1 for the solid curve and incorporate the variation of κ 0 from 0 to 6 to estimate the error in Fig. 12. The higher dimensional condensates are neglected in this paper. We also include the renormalization scale dependence to estimate the truncation error as discussed in Sec. II, whereΠ pert 0 and Π pert m 2 correspond to c 0 and c 2 . We introduce the renormalization scales µ 0 and µ 2 for c 0 and c 2 , respectively, vary them in the range 2M 2 e −γ E ≤ µ 2 0 , µ 2 2 ≤ 8M 2 e −γ E separately, and take the maximal (minimum) value of c 0 + c 2 /M 2 as the upper (lower) limit of the band.  Fig. 12.

B. Extraction of the gluon condensate
As an application of the lattice calculation ofΠ(M 2 ), we try to determine the coefficient c 4 from the lattice data. Since the perturbative expansion and OPE converges reasonably well forΠ(M 2 ) ( although some uncertainty remains if κ 0 ∼ 5-6), the determination is less affected by the truncation error than that for the HVP function Π(q 2 ), and the systematic error of c 4 may be reduced. We consider corrections up to mass dimension six, since the To evaluate the systematic uncertainties, we use three sets of the renormalization scales (µ 2 0 , µ 2 2 ) = (4M 2 e −γ E , 4M 2 e −γ E ), (µ 2 0 , µ 2 2 ) = (2M 2 e −γ E , 8M 2 e −γ E ), and (µ 2 0 , µ 2 2 ) = (8M 2 e −γ E , 2M 2 e −γ E ), and take the maximum variants of the results as their systematic errors. We obtainc 4 = −0.34(7) +26 −19 . The first parenthesis gives the statistical error. The superscript (subscript) represents the upper (lower) systematic error.c 6 is not well constrained.
We subtract the contributions of the chiral condensate and the finite mass correction from c 4 , which are relatively well determined, and obtain 0| αs π G 2 |0 = 0.011(7) +22 −16 GeV 4 in the MS scheme at the scale µ = 2 GeV, which corresponds to 0| αs π G 2 |0 = 0.013(8) +27 −20 GeV 4 in the renormalization group invariant (RGI) scheme. They are related by (see also [56]) 0| α s π G 2 |0 RGI = 1 + 16 9 The first error includes the statistical errors of lattice calculations and inputs Λ n f =3 MS , m s , and 0|qq|0 . The second one corresponds to the systematic uncertainty associated with the perturbative expansion. It is known that the gluon condensate suffers from the renormalon ambiguity. (See, for instance, [57].) More precise determination of the gluon condensate will require more statistics and an improvement of the perturbative calculation.

VI. CONCLUSION AND OUTLOOK
The Borel transform has often been used in the QCD sum rule analyses in order to improve the convergence of OPE and to enhance the contribution of the ground state, which is of the main interest. A crucial question is then whether the theoretical uncertainty in the perturbative expansion and OPE is well under control. The uncertainty due to the modeling of the excited state and continuum contributions is another important issue in the QCD sum rule. In this work, we provide a method to compute the Borel transform utilizing the lattice QCD data for current correlators. Since the computation is fully nonperturbative in the entire range of the Borel mass M , one can use the result to verify the theoretical methods so far used in the QCD sum rule.
We find a good agreement between the lattice data and OPE in the region of M > 1.0 GeV. The OPE is truncated at the order 1/M 6 . Since the OPE involves unknown condensates, this comparison can be used to determine these parameters, provided that the lattice data are sufficiently precise. As the first example, we attempt to extract the gluon condensate, which appears in OPE at the order 1/M 4 . The size of the error is comparable to those of previous phenomenological estimates. With more precise lattice data in various channels, one would be able to determine the condensates of higher dimensions, which have not been determined well solely from phenomenological inputs.
Using baryonic current correlators, one may also study another side of the QCD phenomenology. Since there are no experimental inputs, the lattice data may play a unique role in the QCD sum rule analysis. For instance, the Ioffe formula for the nucleon mass m N [−2(2π) 2 0|qq|0 ] 1/3 [60] indicates a relation between the nucleon mass and chiral symmetry breaking, and it is interesting to study the baryonic correlator on the lattice to see if this relation comes out.
Another interesting application of the lattice calculation of the Borel transform is the determination of α s . A similar analysis has been performed directly on the current correlators [7], but it turned out that OPE does not converge sufficiently quickly to allow precise determination of α s from the perturbative expansion at the leading order of OPE. With the Borel transform, one expects that OPE converges more rapidly, and it may provide another way to extract α s , especially because the perturbative expansion is known to O(α 4 s ), i.e. among the best quantities for which high order perturbative expansion is available.
Our work provides a technique to relate two major tools to study nonperturbative aspects of QCD, i.e., the QCD sum rule and the lattice QCD. As outlined above, there are a number of applications, for which new insights into the QCD phenomenology are expected.