Towards NNLO Accuracy in the QCD Sum Rule for the Kaon Distribution Amplitude

We calculate the $O(\alpha_s)$ and $O(\alpha_s^2)$ gluon radiative corrections to the QCD sum rule for the first Gegenbauer moment $a_1^K$ of the kaon light-cone distribution amplitude. The NNL0 accuracy is achieved for the perturbative term and quark-condensate contributions to the sum rule. A complete factorization is implemented, removing logarithms of $s$-quark mass from the coefficients in the operator-product expansion. The sum rule with radiative corrections yields $a_1^K(1 \GeV)=0.10\pm 0.04$.

1. Light-cone distribution amplitudes (DA's) of hadrons enter various factorization formulae used for description of exclusive processes in QCD. The concept of DA's allows to describe collinear partons in an energetic hadron, separating long-distance dynamics from the perturbatively calculable hard-scattering amplitudes.
The set of DA's with a growing twist is especially useful for the pion and kaon because their intrinsically small masses make the collinear description more efficient. The lowest twist-2 DA has a transparent physical interpretation, describing the longitudinal momentum distribution in the quark-antiquark Fock-state of a meson. Switching from the pion to kaon, one encounters the SU(3) f l -symmetry violation effects, which originate from the quark mass difference m s − m u,d . These effects have to be accounted as accurate as possible, in order to assess the SU(3) f l symmetry relations between the hadronic amplitudes with pions and kaons. Important examples are the relations between B → ππ and B → πK, KK charmless decay amplitudes employed in the studies of CP-violation and quark-flavour mixing.
The most essential SU(3) f l -violating effects in the kaon twist-2 DA include the ratio of the decay constants f K /f π and the difference between the longitudinal momenta of strange and nonstrange quark-partons. This difference is proportional to the first moment a K 1 in the decomposition of the kaon twist-2 DA in Gegenbauer polynomials, whereas a π 1 vanishes in the isospin (G-parity) symmetry limit. In addition, the ratio of the second Gegenbauer moments a K 2 /a π 2 can also deviate from unity; the effects related to a K n at n ≥ 3 are usually neglected.
In this paper we concentrate on the determination of the asymmetry parameter a K 1 (µ) for the kaon, at a low scale µ ∼ 1 GeV. The method originally suggested in [1] and based on QCD sum rules [2] is employed. The most recent sum rule estimates of a K 1 were obtained in [3] and [4], where, in addition to the known leading-order (LO) results, the next-to-leading (NLO), O(α s ) correction to the quark-condensate contribution are taken into account. These calculations, together with the estimates [5,6] based on the operator identities, yield the interval (quoted as a best estimate in [7]): a K 1 (1GeV) = 0.06 ± 0.03 . The positive sign of a K 1 corresponds, as expected, to a larger average momentum of the heavier valence s-quark in the kaon.
The aim of this work is to upgrade the precision of the QCD sum rule for a K 1 . We calculate the gluon radiative corrections to the perturbative and quark-condensate contributions in NNLO, including O(α s ) and O(α 2 s ) corrections. This task is technically feasible, due to the currently achieved state-of-the-art in the calculations of multiloop effects in the two-point correlation functions with strange and nonstrange quarks. For the correlation function with the scalar and pseudoscalar currents the O(α 4 s ), five-loop accuracy has recently been achieved [8], and used, e.g., for the QCD sum rule determination of the strange quark mass [9,10,11]. In this case, in the perturbative expansion the O(α 2 s ) terms are important numerically, which is one motivation to include these terms also in the sum rule for a K 1 . The correlation functions underlying the sum rules for Gegenbauer coefficients are however different, because the currents contain derivatives. Therefore, the calculation reported in the present paper, involves a certain technical novelty. In addition, we clarify and take into account the mixing of operators that is necessary for the complete factorization of small and large scales in the correlation function. Our result for the first Gegenbauer moment of the kaon DA is: In what follows, we introduce the correlation function, present the expressions for the new radiative corrections, derive the resulting QCD sum rule for a K 1 , including the new O(α s ) and O(α 2 s ) terms and perform the numerical analysis.
2. The twist-2 DA of the kaon enters the standard expression for the light-cone expansion of the vacuum-kaon bilocal matrix element (we take K − for definiteness): where the s-andū quarks carry the momentum fractions u andū , and µ is the normalization scale determined by the interval z 2 near the light-cone. We use the compact notation A ρ = g s A a ρ λ a /2 for the gluon field and the covariant derivative is defined as D ρ = ∂ ρ − iA ρ . In (2), the twist-2 DA ϕ K (u) is normalized to unity, so that in the local limit z → 0 one reproduces the definition of the kaon decay constant f K .
As usual, ϕ K (u) is expanded in the Gegenbauer polynomials with the coefficients a K n (µ) (Gegenbauer moments). The first Gegenbauer moment a K 1 is proportional to the average difference between the longitudinal momenta of the strange and nonstrange quarks in the two-parton state of the kaon. Expanding both parts of Eq. (2) around z = 0 in local operators and using the decomposition (3) with C 3/2 1 (x) = 3x one relates a K 1 to the vacuum-to-kaon matrix element of a local operator with one derivative: . The Gegenbauer moments a π,K n (µ) are known to be multiplicatively renormalizable only at the one-loop level. Generally, this property is lost at higher orders in α s , e.g., the two-loop renormalization of a π 2 calculated in [12] includes operator-mixing effects. Still the a K 1 case is special, in so far as the underlying operatorsγ ν γ 5 i ↔ Dλ u can only mix with ∂ λ (sγ ν γ 5 u), as there is no other local operator with the same dimension and flavour content. However, the above two operators have opposite G (s) -parities, where G (s) is the analog of the isospin G-parity for the SU(2) subgroups of SU(3) f l involving s quark (V or U-spins). Naturally, G (s) -conservation is only realized in the m s = m u,d limit. Note, however, that the ultraviolet renormalization in MS-scheme is a mass-independent procedure. Hence it is legitimate to consider the SU(3) f l limit, while performing the renormalization, so that the G (s) -conservation protects the operators from mixing with each other. As a result, a K 1 remains multiplicatively renormalizable at any order in perturbation theory. For completeness, we present the well-known expression for scale-dependence of a K 1 with the two-loop (NLO) accuracy, written in an unexpanded form: where γ 0 = 8/9, γ 1 = 590/243 are the anomalous dimensions [13] and β 0 = 9/4, β 1 = 4 are the coefficients of β-function for n f = 3.
As originally suggested in [1], the few first Gegenbauer moments of DA's can be calculated employing operator-product expansion (OPE) and QCD sum rules for two-point vacuum correlation functions. The method works well only for the first two coefficients a π,K 1,2 . In the sum rules for a π,K n≥3 the condensate contributions grow fast with n and the control over OPE is lost 2 .
To obtain a QCD sum rule for a K 1 , it is convenient to take the so called "diagonal" correlation function: (6) where, for brevity, only the relevant kinematical structure is shown. In fact, (6) is not quite diagonal, because one of the operators, the same as in (4), contains a derivative. A different choice is to correlate the operator in (4) with the pseudoscalar currentūiγ 5 s. This, so called "nondiagonal" correlation function was tried in [15], but produces an unstable sum rule (see the discussion in [3]).
3. In LO, the operator product expansion (OPE) for the correlation function (6) includes the quark-loop diagram at O(α 0 s ) and the contributions of the vacuum condensates, calculated in the deep spacelike region, The corresponding diagrams are collected in Fig. 1. One obtains a generic expansion for the invariant amplitude Π(q 2 ) defined in (6) in inverse powers of the variable Q 2 : where the coefficients A d>2 contain condensate densities with growing dimensions. In the above, µ is the ultraviolet renormalization scale in the loop diagrams. The OPE is applicable at sufficiently large Q 2 , provided that the coefficients A d are proportional to the powers of the light-quark masses and/or to the condensate densities, the latter being of O(Λ QCD ) in some power. Hereafter we neglect the u, d-quark masses with respect to m s . In the expansion (7), the terms proportional to 1/Q d>6 are also neglected, since already the contribution of the d = 6 term is quite small. Hence, only the vacuum condensates with dimension d ≤ 6 are taken into account. The gluon radiative corrections can be systematically included in each term of the OPE (7). Diagrammatically, the α s -(α 2 s -) corrections correspond to all possible one-gluon (two-gluon) insertions in the Fig. 1 diagrams. The second small parameter entering the OPE of the correlation function (6) is the ratio m 2 s /Q 2 . Hence, each of the coefficients in (7) can be cast into a form of a generic double expansion: where the coefficients a It is important to assess the numerical role of the small parameters in the combined expansion (8). We expect to use OPE at Q 2 ≃ 1 GeV 2 . Taking α s (1 GeV) = 0.47 [16] and a conservative upper limit for the strange quark mass, m s (1 GeV) < 150 MeV, one has m 2 s /Q 2 ≤ 0.02 ≪ α s /π ≃ 0.15. Hence, the perturbative O(α k s ) contributions to the OPE are expected to be more important than the O((m 2 s /Q 2 ) k ) terms with the same power. In particular, in the first line of (8) the second-order, α 2 s -correction is expected to be of order of the m 2 s /Q 2 -term. Moreover, the observed hierarchy allows one to neglect all "mixed" In what follows, we restrict ourselves to the previously known first-order corrections in m 2 s /Q 2 in the d = 2, 4 terms of OPE, neglecting them in the d = 6 term.
So far only the O(α s ) correction to the quark-condensate contribution A 4 was calculated [3,4], that is, the coefficient a and a (20) 4 , respectively, in (8). Hence, for the largest d = 2, 4 terms of the OPE (7) the NNLO accuracy in α s is achieved. For the subleading d = 6 term in the OPE we retain the known LO result.
The gluon radiative corrections are computed in the standard MS-scheme of renormalization. The well developed techniques of loop calculations are employed, in particular, the programs QGRAF [17], FORM [18], and MINCER [19,20]. The following results are given for n f = 3. Including the new O(α s ) and O(α 2 s ) contributions, we obtain the d = 2 term in (7) originating from the perturbative contribution: The d = 4 contribution in (7) generated by the quark-condensate term has the following expression to O(α 2 s ): where qq ≡ 0|qq|0 , (q = s, u) is the quark-condensate density. Finally, the d = 6 term in (7) contains the LO quark-gluon, gluon and four-quark condensate contributions [3,4]: where sGs ≡ 0|sσ µν g s G a µν (λ a /2)s|0 , G 2 ≡ 0| αs π G a µν G a µν |0 , and the four-quark condensate is factorized into the products of two quark condensate densities, assuming isospin symmetry d d = ūu .
Comparing with the previous calculations [3,4], we observe a difference in the O(α s ) s-quark condensate term. More specifically, the constant term 112/27 in the first line of (10) has to replace 124/27 in the corresponding expressions obtained in [3] and [4]. 3 This difference has, however, only a minor influence on the numerical results.
The form of the coefficient function multiplying m 4 s /Q 2 in (9) and m 2 s G 2 in (11) deserves a separate discussion. One has to emphasize that, in order to achieve a complete factorization of small and large scales in OPE, all logarithms of the small parameter m s have to be removed from the coefficient functions leaving only the powers of m 2 s /Q 2 . This procedure, understood long ago [21,22,23,24,25], generates terms proportional to ln(µ 2 /Q 2 ) instead of ln(m 2 s /Q 2 ), e.g., in the LO parts of A 2 and A 6 in (9) and (11), respectively. These logarithms were not properly treated in previous calculations. Note that a simple replacement ln(m 2 s /Q 2 ) → ln(µ 2 /Q 2 ) in the coefficient functions can miss a constant term which has to be added to the logarithm. Let us, for example, explain the calculation of the coefficient function in (11), taking into account the mixing of sGs and m s G 2 terms under renormalization in MS scheme.
We isolate the contributions of the quark-antiquark-gluon and gluon condensates to the correlation function (6) and write them in the following convenient form: where all other contributions indicated by ellipses are not important for this discussion. The two terms shown in (12) are generated by the renormalized d = 6 operators m ss Gs and m 2 s G 2 . To calculate the coefficient functions C 1,2 (Q 2 ) at D = 4 − 2ε it is sufficient to write the pattern of the mixing of these two operators in the form where the index "nr" indicates that the operators on r.h.s. are constructed from bare, nonrenormalized quark and gluon fields. The Z factor of the multiplicative renormalization of (sGs) nr can be put to unity in this approximation. To proceed, we need the expansion of the operator product in (6) before the vacuum average is taken: where we use (13). The coefficient functions C 1,2 (Q) can now be obtained by projecting both parts of (14) onto the suitable free quark-gluon states. Sandwiching (14) between the sg| and |s states, and computing the relevant diagrams we find: Note that one has to retain the O(ε) term for the finite tree-level diagrams. Furthermore, projecting (14) onto the two-gluon state, we obtain: After substituting C 1 from (15) to (16), the 1/ε poles cancel each other and the desired results for the renormalized coefficient functions are obtained at ε → 0: We again emphasize that one has to keep the O(ε) term in the expression for the tree-level coefficient function C 1 (Q 2 ) in order to get the non-logarithmic contribution to C 2 (Q 2 ). The absence of 1/ε poles in the above calculation can also be interpreted as a result of the cancellation between ultraviolet and infrared divergences. Indeed, the 1/ε pole in the mixing of operators emerges as an ultraviolet divergence, whereas 1/ε in the loop diagrams used to calculate the coefficient functions has an infrared origin. A similar derivation is applied to the mixing of O(m s ss ) and O(m 4 s ) terms in OPE, yielding the constant term 5/2 that accompanies ln(µ 2 /Q 2 ) in the coefficient function (9).

4.
After the correlation function is calculated, the derivation of the sum rule follows the standard procedure. The OPE (7) with the coefficients given in (9), (10) and (11) is equated to the hadronic dispersion relation In the above, possible subtractions are ignored in anticipation of the Borel transformation. The residue of the kaon pole is obtained by combining the matrix element (4) with the definition of the kaon decay constant. The spectral density ρ h (s) includes the contributions of hadronic continuum and resonances with J P = 0 − , 1 + and strangeness: Kππ, K * π, Kρ, K 1 (1270), K 1 (1400),... . Accordingly, the lower limit of integration is s h = (m K + 2m π ) 2 , the invariant mass squared of the lightest continuum state in this channel. To approximate ρ h (s), we employ the quark-hadron duality approximation: where s K 0 is the effective threshold, and the spectral density: with l s ≡ ln(µ 2 /s), is obtained by calculating the imaginary part of Π(Q 2 ) at positive s = −Q 2 . In the above, the running parameters (α s , m s ) are taken in the MS scheme and normalized at the scale µ. The next step is the Borel transformation of (18), which eliminates the subtraction terms and suppresses the integral over ρ h (s) , so that the resulting relation becomes less sensitive to the duality approximation. The transformed invariant amplitude Π(q 2 ) has the following form: where l M ≡ ln(µ 2 /M 2 ) and γ E is the Euler constant. Finally, the sum rule for the kaon Gegenbauer moment a K 1 obtained from (18) using (19) reads: where the functions Π(M 2 ) and ρ OP E (s) are given in (21) and (20), respectively. An equivalent form of the sum rule where the entire r.h.s. is represented as a duality integral over the spectral density, was used in [3,4].

5.
To perform the numerical analysis of the sum rule (22), we specify the relevant input parameters, starting from the kaon mass m ± K = 493.58 MeV and decay constant f K = 159.8 ± 1.4 ± 0.44 MeV [16]. For the strange quark mass we adopt: which covers the intervals m s (2 GeV) = 105 ± 6 ± 7 MeV [9] (with uncertainties added in quadrature) and m s (2 GeV) = 92 ± 9 MeV [10] from the recent QCD sum rule determinations with O(α 4 s ) accuracy. The running of the coupling α s (m Z ) = 0.1176 ± 0.002 [16] to the lower scale µ is taken with the 4-loop accuracy applying the program presented in [26]. The default value of the renormalization scale is µ = 1 GeV, where the expansion parameter is α s (1 GeV)/π = 0.15 ± 0.01 and m s (1 GeV) = 128 ± 21 MeV is obtained from (23).
The accuracy of the condensate densities with higher dimensions is less important. In particular, we take for the quark-gluon condensate density the standard parameterization sGs = m 2 0 ss (1 GeV) with m 2 0 = 0.8 ± 0.2 GeV 2 [30,29] (neglecting the running), and a very wide interval for the gluon condensate density G 2 = 0.012 ± 0.012 GeV 4 . We also allow for a varying factor 0.1 ÷ 1 multiplying the square of the quark condensate density in the factorization relation for the four-quark condensates.
The Borel parameter interval M 2 = 1.0 − 2.0 GeV 2 is adopted, so that the d = 6 contribution in Π(M 2 ) remains lower than 30%. The upper limit is taken to avoid too large contributions of the excited and continuum states estimated with duality approximation. The dependence of a K 1 on M 2 for s K 0 = 1.05 GeV 2 and at central values of all input parameters is plotted in Fig. 2, showing a remarkable stability. There we also plot separate contributions to the sum rule originating from the d = 2, 4, 6 terms of the expansion (7). The threshold parameter s K 0 = 1.05 GeV 2 was fixed in [3] from the sum rule for f K . In fact, the s K 0 -dependence of the sum rule result turns out to be rather weak, if one varies the threshold in rather wide limits s K 0 = 0.9 − 1.4 GeV 2 , the dependence of a K 1 on M 2 remaining flat. The reason is that the spectral density (20) is small in the region above the duality interval for the kaon. Therefore, although the pattern of hadronic states (resonances and continuum states) in this region is rather complicated, we expect that the integral over these states is also small. Hence the sum rule result will not noticeably change if one modifies the duality ansatz for the hadronic spectral density, e.g. by adding more resonances and correspondingly increasing the threshold parameter (as shown e.g., in [3]). The numerical prediction of the sum rule is represented in the form: where the first ("SR") error is the combined uncertainty of a K 1 due to the variation of M 2 and s K 0 , representing a sort of an "intrinsic " uncertainty of the sum rule. The subsequent errors correspond to the individual variations of α s (m Z ), m s (2 GeV) and m u,d (2 GeV) (that is, qq ), within the adopted intervals. The last error ("cond") shown in (24) is a combined uncertainty due to variation of all remaining condensate parameters. Finally, adding the individual uncertainties in quadrature we obtain the interval (1).
Let us now discuss the structure of the perturbative series as it follows from the numerical analysis. First, we observe that gluon radiative corrections substantially enhance the perturbative d = 2, O(m 2 s ) term and suppress the quark condensate term. Numerically, for the O(m 2 s ) term in (21) at µ = M = 1 GeV the following α s -expansion is obtained (in MS-scheme): revealing a poor convergence in α s at α s (1 GeV)/π ≃ 0.15. It is not surprising, if one recalls the higher-order perturbative corrections calculated for the m s determination from τ decays or from QCD sum rules based on scalar/pseudoscalar correlation functions, where the situation is in fact similar. For instance, in the O(m 2 s ) part of the Borel-transformed scalar/pseudoscalar correlation function [8,9] one has (see e.g., eq.(16) in [9]): The ratio of (26) (taken with O(α 2 s ) accuracy) and (25) has a much better convergence: where the corrections cancel each other to a great extent.
In the quark condensate contribution, the radiative corrections are less sizable: at µ = 1 GeV the numerical hierarchy of the terms multiplying the strange quark condensate density in (10) where we again find that the second-order correction in α s is important. We conclude that after taking into account the NNLO perturbative corrections, the numerical pattern of the sum rule for a K 1 drastically changes: the coefficient of the d = 2 term gets enhanced, whereas the d = 4 term decreases. Note also that including higherorder perturbative corrections in the sum rule, makes more consistent the use of the input parameters, such as m s , determined with a high accuracy, up to O(α 4 s ). Finally, the d = 6 subleading contributions to the sum rule play an important role in providing the Borel stability. The mixed quark-antiquark-gluon condensate dominates numerically in (11) yielding a negative contribution to the Borel-transformed correlation function and stabilizing the whole sum of OPE terms. Computation of the radiative corrections to this term is a difficult task beyond our scope. Having in mind the large uncertainty of the mixed condensate density, we expect that the α s -corrections to the coefficient function of the sGs term will hardly improve the overall accuracy of its contribution.
Finally, we checked that due to the enhanced precision in α s , the dependence of the sum rule prediction for a K 1 on the renormalization scale µ becomes small. Taking µ > 1 GeV we calculate a K 1 (1 GeV) by rescaling the sum rule result with the NLO scale dependence (5). The LO and NLO logarithmic dependences naturally cancel out in d = 2, 4 terms, leaving a very mild residual scale-dependence due to the unaccounted α s -correction to the sGs term.
6. Concluding, we have calculated the NNLO gluon radiative corrections to the QCD sum rule for the first Gegenbauer moment of the kaon distribution amplitude. The corrections turned out to be numerically important, they change the relative magnitude of the d = 2 (loop diagrams) and d = 4, 6 (condensate) terms in the OPE, improving also the Borel stability of the sum rule.
The uncertainty of a K 1 is still large and amounts up to 40%, due mainly to the limited precision of the light quark masses: m s directly entering the sum rule and m u,d determining the quark-condensate densities via Gell-Mann-Oakes-Renner relation. A better determination of the ratio of strange and nonstrange condensates and of the mixed quark-gluon condensate are another possibilities of reducing the theory error. The calculation of the radiative corrections to the quark-gluon condensate contribution can also improve the accuracy. This however requires a dedicated computational effort including an analysis of the whole basis of the dimension-six operators. The weak dependence of the sum rule on the threshold parameter s K 0 indicates that the one-resonance (kaon) duality ansatz is quite satisfactory. Still, an additional analysis, including the axial-vector and radially excited kaon resonances could provide a better understanding of the duality pattern in this channel.
Our result for a K 1 is somewhat larger than the previous estimates [3,4,7], and the uncertainty, we believe, is more realistic. For comparison we also quote the two recent lattice QCD determinations of this parameter: a K 1 (2 GeV) = 0.0453 ± 0.0009 ± 0.0029 [31] and a K 1 (2 GeV) = 0.048 ± 0.003 [32], which have achieved a rather small error. By evolving our result (1) to this scale with the help of the scale-dependence (5) we find a K 1 (2 GeV) = 0.08 ± 0.04, that is, within uncertainties, only a marginal agreement with the lattice results. Let us note that both lattice determinations use a linear extrapolation in m s (kaon mass squared) inspired by ChPT in the leading order (see also [33]). We would like to stress again that in our analysis the contribution to the sum rule proportional to m 2 s is enhanced while the term proportional to m s ss is suppressed, both enhancement and suppression being caused by the radiative corrections. In the language of ChPT, this observation could indicate an important role of the next-to-leading terms in the expansion of a K 1 in the kaon mass.