Transition form factors of the pion in light-cone QCD sum rules with next-to-next-to-leading order contributions

The transition pion-photon form factor is studied within the framework of Light-Cone QCD Sum Rules. The spectral density for the next-to-leading order corrections is calculated for any Gegenbauer harmonic. At the level of the next-to-next-to-leading (NNLO) radiative corrections, only that part of the hard-scattering amplitude is included that is proportional to the $\beta$-function, taking into account the leading zeroth-order harmonic. The relative size of the NNLO contribution in the prediction for the form factor $F^{\gamma^{*}\gamma\pi}(Q^2)$ has been analyzed, making use of the BLM scale-setting procedure. In addition, predictions for the form factor $F^{\gamma^{*}\rho\pi}$ are obtained that turn out to be sensitive to the endpoint behavior of the pion distribution amplitude, thus providing in connection with experimental data an additional adjudicator for the pion distribution amplitude. In a note added, we comment on the preliminary high-$Q^2$ BaBar data on $F^{\gamma^{*}\gamma\pi}$ arguing that the significant growth of the form factor between 10 and 40 GeV$^2$ cannot be explained in terms of higher-order perturbative corrections at the NNLO.


I. Introduction
Although higher-order calculations in QCD perturbation theory have already a long history, little is known about exclusive processes at the next-to-leading order (NLO) level [1,2,3,4,5,6], and beyond [7,8], because these are quite complex in detail. In view of more and more high-precision experimental data for a variety of hadronic processes becoming gradually available, the importance of such higher-order calculations exceeds the pure theoretical interest and acquires phenomenological relevance. In particular, processes with two photons in the initial state, one far off-shell and the other quasi real, provide a useful tool to access (after their fusion) the partonic structure of the produced hadronic states, e.g., pseudoscalar mesons.
Experimentally, the photon-to-pion transition form factor within this class of two-photon processes has been measured by the CLEO Collaboration [9] with high precision and extending the range of Q 2 up to 9 GeV 2 , as compared to the previous low-momentum CELLO data [10]. Theoretically, this high precision allows one to test models and fundamental quantities, like the pion distribution amplitude (DA), the applicability of QCD factorization, etc.-see [5,6,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27] and references cited therein. Moreover, one can determine [26] a compatibility region between the CLEO data and constraints derived from lattice simulations on the second moment of the pion DA [28,29]. This information can then be used to extract a range of values of the fourth moment of the pion DA that would simultaneously fulfil both constraints (CLEO and lattice). This prediction [26] can provide a guide for the determination of this moment on the lattice, a task that has not been accomplished yet.
For two highly virtual photons, perturbative QCD works well because factorization at some factorization scale µ 2 F applies, so that the process can be cast into the form of a convolution F γ * γ * π (Q 2 , q 2 ) = C Q 2 , q 2 , µ 2 F , x ⊗ ϕ π x, µ 2 F + O(Q −4 ) , (1.1) which contains a hard part C, calculable within perturbation theory, and a wave-function part ϕ π that is the (leading) twist-two pion distribution amplitude [30] and has to be modeled within some nonperturbative framework (or be extracted from experiment). Here, the omitted twist-four contribution represents subleading terms in the operator product expansion (OPE), which are suppressed by inverse powers of the photon virtualities.
If both virtualities, Q 2 and q 2 , are sufficiently large, the T -product of the currents can be expanded near the light cone (z 2 = 0) by virtue of the OPE to obtain the well-known where we have used the abbreviationsx ≡ 1 − x and (1. 4) In contrast, the kinematics probed in the CELLO [10] and the CLEO [9] experiments involves a quasi-real photon with q 2 → 0. At such a low virtuality, the hadronic content of the quasi-real photon, i.e., its long-distance structure [14] becomes important (see the diagram in the right panel of Fig. 1), thus preventing a straightforward QCD calculation [13] of the form factor F γ * γπ (Q 2 , q 2 → 0) ≡ F γ * γπ (Q 2 ) on the ground of factorization. The method of Light-Cone QCD Sum rules (LCSR for short) allows one to avoid this problem by providing the means of performing all QCD calculations at sufficiently large q 2 (γ * ) and then use a dispersion relation to "approach" the mass-shell photon (γ) with zero virtuality. This calculational scheme, which can accommodate the large-distance properties of the photon, i.e., its hadronic content, was proposed by Khodjamirian in [14] and the form factor F γ * γπ (Q 2 ) was calculated at the LO level of the LCSR including also twist-four contributions.
The core ingredient of the LCSR approach is the spectral density, which provides a powerful tool for a quantitative description of hadronic processes in QCD in terms of a dispersion relation: ρ phen (Q 2 , s) s + q 2 .
(1. 5) [Note that here the label "phen" abbreviates phenomenological]. Effects due to the longdistance dynamics of the γ * γ → π 0 process are partly contained in the form factor F γ * ρπ (Q 2 ) (which can be obtained by means of quark-hadron duality) and also in the π distribution amplitudes of different twists [31,32,33,34]. Within the LCSR approach the form factor F γ * ρπ (Q 2 ) appears inevitably because one assumes that the spectral density, entering the dispersion relation, can be approximated by the ground states of vector mesons [14,35], like the ρ and the ω. 1 At this point, two important remarks are in order. (i) Contrary to previous calculations, e.g., in [5,6], we will not adopt a zero-width approximation here, but use instead a more realistic Breit-Wigner ansatz for the effective resonances (see Sec. V). (ii) The scaled form factor Q 4 F γ * ρπ (Q 2 ), obtained in the framework of LCSRs, depends at large Q 2 mainly on the differential pion characteristic d dx ϕ π (x)| x=ǫ , with ǫ ∼ s 0 Q 2 ≪ 1 being a small neighborhood around the origin, where s 0 is the duality interval entering the model for ρ phen . This feature appears to be opposite to the case of the Q 2 F γ * γπ (Q 2 ) (scaled) form factor, that depends mainly on (though it is not directly proportional to) the inverse moment [22] x −1 π = 1 0 ϕ π (x; µ 2 )x −1 dx, cf. Eq. (1. 3) evaluated at the scale q 2 → 0, because the latter is an integral characteristic of the pion DA [18,19]. Hence, Q 4 F γ * γρπ (Q 2 ) can provide complementary information on the pion DA and help discriminate among various proposed pion DA models.
The structure of the paper is the following. In the next section, we recall the formalism of LCSRs for the form factors F γ * γπ (Q 2 ), F γ * ρπ (Q 2 ) and construct the spectral density in a systematic way. This calculation is extended beyond the LO in Sec. III and an explicit expression for the spectral density at the NLO for any index n-the latter indicating the order of the expansion in Gegenbauer harmonics-is derived. The further extension to the NNLO is also given in this section. Actually, we include only the β 0 -proportional contributions that can be obtained from the corresponding terms of the hard-scattering amplitudes and denote it by NNLO β . The effects of the NNLO contributions are discussed in Sec. IV in connection with the BLM prescription (and its modifications [36]). Our predictions for the form factors F γ * γπ (Q 2 ), F γ * ρπ (Q 2 ) are presented in Sec. V, where we also provide a comparison of F γ * γπ (Q 2 ) with the experimental data. Section VI contains our conclusions emphasizing our main results. Important technical details are provided in three dedicated appendices. In a Note Added, we point out that the new BaBar data on F γ * γπ , which show a significant growth of the form factor beyond 10 GeV 2 , cannot be described within the QCD convolution approach.
Here we present the theoretical description of the transition form factor F γ * γπ (Q 2 ) of the exclusive process γ * γ → π 0 , employing the framework of light-cone sum rules [14,37,38,39] beyond the next-to-leading-order (NLO) of perturbative QCD. 2 An integral part of this sort of calculation is the form factor F γ * ρπ (Q 2 ), describing the transition γ * ρ → π, which will, therefore, be computed in parallel. 1 For the sake of simplicity, one sets the masses of the two vector mesons equal and appeals to isospin symmetry to treat both particles in terms of a combined effective resonance. 2 For the sake of clarity, we use the following notation for the form factors: The associated reaction is denoted by superscripts, whereas the calculational context, e.g., QCD, is marked by subscripts. Note that our notation differs from abbreviated notations used in the cited works.

A. Factorization
The calculation of F γ * γπ (Q 2 ) proceeds through the following main steps: (i) First, the form factor F γ * γ * π QCD (Q 2 , q 2 ) at large Euclidean virtualities of the photons, Q 2 , q 2 ≥ 1 GeV 2 is calculated. (ii) Then, an appropriate realistic model for the spectral density at low s, based, for instance, on quark-hadron duality, is constructed. (iii) Finally, the dispersion relation for the form factor F γγ * π (Q 2 , q 2 ) is exploited.
Applying the factorization theorems [31,32,33,34], the form factor F γ * γ * π QCD (Q 2 , q 2 ) can be cast in the form The hard-scattering amplitude T for this process, written below in the square brackets, is calculable within perturbative QCD. The symbols µ R and µ F denote, respectively, the scale of the renormalization of the theory and the factorization scale of the process. The pion DA ϕ (2) π (x; µ 2 F ) of twist two, entering the convolution with the hard-scattering amplitude, is inaccessible to perturbative QCD and demands the application of nonperturbative methods (see Section V A). We quote the well-known result for F γ * γ * π QCD at LO in a s that also includes the twist-four contribution [14], viz., and where ϕ (4) π (x) is a naturally appearing combination of twist-four pion DAs (for more details, see [40] and Sec. V A). An explicit expression for T 1 has been obtained in [1,2,3]. More recently [7], the general structure of T 2 was investigated and its β-part, b 0 · T β , was calculated. For convenience, both amplitudes are included in Appendix A. In fact, the hard-scattering amplitudes T 0 , T 1 , T 2 also determine the spectral densities ρ (0) , ρ (1) , ρ (2) , as one can see from the following generic expression The core issues of the hard-scattering amplitudes are listed below, while characteristic Feynman graphs for each order of the perturbative expansion are depicted for illustration in Fig.  2: In this order of the expansion we have T 0 (Q 2 , q 2 ; x) ⊗ ϕ(x; µ 2 F ) and the factorization scale cannot be determined uniquely.
2. NLO perturbative QCD: T 1 . A typical diagram is shown in Fig. 2 Here the hard-scattering amplitude starts depending upon the factorization scale µ F . In explicit terms this reads and the µ F -dependence enters via the logarithm whereas the LO Efremov-Radyushkin-Brodsky-Lepage (ERBL) kernel [31,32,41] V (0) is given in the next section with more details provided in Appendix B. As regards the term T (1) (y, x), defined in (A.7), it contains contributions from the partial kernel V b (which enters V (0) and is generated by contracting the diagram in Fig. 2(b)) and also the kernel g stemming from V (1) and defined in (A.5).
3. NNLO perturbative QCD: T 2 with an example depicted graphically in Fig. 2(c). In this case, things are more complicated. First of all, we reiterate that only the β 0 -proportional contributions to T 2 (termed T β ) will be considered, making use of the calculations performed in [7]. This is, because that part can provide the right size of the whole NNLO contribution using the scheme-independent BLM scale-setting procedure to eliminate the appearance of the β-function in the perturbative series. To discuss the structure of T β , let us first present it explicitly (postponing details to the appropriate sections to follow): Second, as one sees from this expression, the hard-scattering amplitude depends explicitly also on the renormalization scale, µ 2 R , owing to the renormalization of the strong running coupling. As a result, additional logarithms of the form ln (µ 2 F /µ 2 R ) appear. These terms are controlled by the renormalization-group (RG) equation and can be resummed by applying the BLM procedure [42] which amounts to a rescaling of the argument of the running coupling to another value. This will be discussed later in more detail. Third, the logarithms Ln(y) and Ln 2 (y), which bear the µ 2 F -dependence [cf. Eq. (2.7)], are accompanied by elements of the kernel V (1) and the kernel V (0) governing the NLO and LO evolution, respectively.

B. Construction of the spectral density
We continue here with the systematic construction of the spectral density. In doing so, we will make use of another set of variables, (x, Q 2 ), instead of the usual set (s, Q 2 ). In the LO approximation, ρ (0) (x) follows from the definition (2.6) upon inserting in Eq. (2.5) the explicit expression for T 0 (Q 2 , q 2 ; x) [5,14]. Then, we obtain [cf. first item in Appendix A]. One notices that ρ (0) is directly proportional to the pion DA of leading twist two, ϕ π , and the derivatives with respect to x of the twist-four contribution, ϕ (4) . This observation allows one to simplify the subsequent analysis by introducing the normalized spectral densityρ(x, Q 2 ) via Eq. (2.9b).
In this representation, all scale dependence is contained in the coefficients a n (µ 2 ) and is controlled by the ERBL equation. For this reason, it is convenient to project the leadingtwist part of the spectral density ρ(Q 2 , s) on the same basis of eigenfunctions {ψ n } and introduce the partial density ρ n , which has the form where x = Q 2 /(Q 2 + s) and n (x) + . . . . (2.12) Therefore, from definition (2.11) and Eq. (2.9a) it follows The spectral density, discussed above, allows us now to construct the phenomenological spectral density (2.14) which consists of two parts. The first term, ρ h , (with "h" denoting hadronic) encodes the hadronic content of the spectral density, reexpressed in terms of the γ * ρ → π transition form factor F γ * ρπ : where we assumed that the ρ and ω resonances have the same mass and can be represented by a δ-function. Later on, we are going to show how this simple ansatz can be improved to include a finite width [14]. The second term in Eq. (2.14) represents its perturbative part 16) in which F γ * γ * π QCD (Q 2 , q 2 ) will be computed according to Eq. (2.2) (i.e., in convolution form) including contributions up to the NNLO β .
Then, substituting ρ phen in the dispersion relation for F γ * γπ (Q 2 , q 2 ) in (1.5) [14], we obtain where we have assumed that the hadronic spectral density, ρ h , in the dispersion integral can be approximated by the expression which determines the structure of F γ * ρπ (Q 2 ). Here, and below, s 0 = 1.5 GeV 2 denotes the duality interval in the ρ-meson channel. In order to derive the LCSR for F γ * ρπ (Q 2 , M 2 ) [14], we first perform a Borel transformation of Eq. (2.18) with respect to the virtuality q 2 and then insert the result into Eq. (2.17). Finally, taking the limit q 2 → 0, we arrive at where M 2 ≈ 0.7 GeV 2 is the typical value of the Borel mass parameter and m ρ is the ρ-meson mass. The first term in Eq. (2.19), which is proportional to F γ * ρπ , expresses the hadronic content of the quasi-real photon, whereas the second term encodes the pointlike subprocesses governed by QCD perturbation theory. The spectral density, given by Eq. (2.6), will allow us to obtain both parts of F γ * γπ LCSR (Q 2 ); viz., and where in correspondence with our remarks below Eq. (2.19). Employing this notation, F γ * ρπ (Q 2 ) from Eq.(2.18) can be recast in the form For the sake of completeness and for future use, we also display the spectral images of Eqs.
where x 0 = Q 2 /(Q 2 + s 0 ) and s = Q 2x /x. The spectral densityρ n beyond the LO will be constructed and discussed in the next section. Let us remark at this point that in the following we are going to consider also mixed forms of these sets of variables, i.e., a spectral density which depends on both variables s and x. This should not cause any confusion because it is understood that one has to replace each time the appropriate variable, i.e., either s = Q 2x /x or x = Q 2 /(Q 2 + s).

III. STRUCTURE OF THE SPECTRAL DENSITY BEYOND LO
This section extends our calculation of the spectral density beyond the leading order of perturbation theory. In the following two subsections we will take into account the radiative corrections at the NLO and also at the NNLO; in the latter case only the β 0 -proportional terms will be included.

A. Spectral density in NLO
Let us first write down the final result forρ (1) n for any index n and then proceed with a systematic discussion of its structure: The partial cases n = 0, 2, 4 coincide, after some algebraic manipulations, with the results obtained before by Schmedding and Yakovlev (SY) [5]. 3 The above expression-besides being valid for any n-also reveals how the radiative corrections manifest themselves in its various terms, as we will now explain. The quantities v b (n) and v(n) in (3.1a) and (3.1b) are the eigenvalues of the corresponding parts of the one-loop ERBL evolution kernel V (0) , notably, V b + , and V b + + V a + . These are given by where the eigenvalues of the partial kernels V a,b are obtained from The latter definition reflects the vector-current conservation, while further details pertaining to the above equations are relegated to Appendix B. Note that the term (3.1b) originates from the logarithm in T 1 [cf. Eq. (2.7)]. Therefore, its contribution to the discontinuity contains, in comparison with ρ (0) in (2.9a), a new element which is discussed in Appendix A, item 2, Eq. (A.12). The corresponding kernel in Eq. (3.1b) maps each monomial term again onto a monomial one, e.g., where a notation deviating from the usual + prescription has been used. Keeping in mind that, ultimately, we have to integrate the spectral density ρ(x) over x, it is particularly useful to recast Eq. (3.1b) in terms of an expansion over an orthogonal polynomial basis, e.g., over the eigenfunctions {ψ n }: The expansion coefficients b nl are given in Appendix A, item 2. As regards the term (3.1c), it originates from the kernel g(y, u), introduced in [43] and discussed in [7] (with some technical details being provided for the reader's convenience in Appendix A, item 2). This kernel, termed g in [43], is not diagonal with respect to the {ψ n }-basis and is responsible for the apparent breaking of the conformal symmetry in the MS-scheme [43]. Notice that also the kernel in Eq. (3.1c) maps each term ψ n onto a sum of ψ l (y), l ≤ n. Therefore, (3.1c) can be cast in the form of the following algebraic expansion Note that the various terms in the spectral densityρ (1) have a one-to-one correspondence to the kernel V (0) and its elements, and partly also to the kernel V (1) via its element g.
(i) For the special case n = 0, ψ 0 (x) = 6xx so that the dependence ofρ on the factorization scale µ 2 F disappears, owing to the fact that the asymptotic DA does not evolve in this approximation. Indeed, in this case, the following chain of evident simplifications is induced This expression agrees with the result obtained in [5] for ρ (1) 0 (Q 2 , s). (ii) For the general case of an arbitrary n, the "nondiagonal" (in ψ n ) part ofρ (1) in the second line of the expression (note that s = Q 2x /x) has been rewritten in terms of the known coefficients G nl and b nl , supplied in Appendix A in terms of Eqs. (A.9) and (A.13). This way, we have achieved that Eq. (3.12) is a purely algebraic expression for ρ. To return to the n = 0 case, one should set in Eq. (3.12) v a (0) = v(0) = 0 and G 00 = 1. Expressions (3.1a)-(3.1c), or, equivalently, (3.12), provide an effective tool for analyzing any model of the pion DA within the LCSR approach. At this point it is worth comparing the NLO contribution, originating from two different approaches, with respect to the case when one of the photons becomes real. First, we have the expression obtained directly from the factorization formula in Eq. (A.7) for q 2 = 0 [3,6,44], i.e., Second, starting from the dispersion relation given by Eq. (2.25) for H, one can return to the previous expression by setting s 0 = 0, x 0 = 1 to get where we have usedρ The crucial observation here is that the dispersion method yields an expression that contains the leading squared logarithm with a negative sign-in contrast to the result one finds with the factorization approach [cf. (3.13)]. In this second case, it is more involved to show [44] that the leading logarithm provides suppression in the relevant integration region and can therefore be associated with Sudakov effects.

B. Spectral density in NNLO. β 0 -proportional contributions
The β-dependent part of the partial spectral density,ρ (2β) n , can be obtained from the corresponding part of the whole amplitude T 2 , given by Eqs. (A.19), i.e., from T β ⊗ ψ n . The calculation of the discontinuity of the latter expression is a rather technical task and is, therefore, relegated to the two Appendices A and C. It is important to realize that the structure of the spectral density,ρ (2β) n , resembles the structure of the analogous term in the NNLO β -amplitude that is proportional to T β ⊗ ψ n , as one may appreciate by comparing the following two expressions: For convenience, we have made use of the abbreviation Ln [cf. (2.7)], omitting arguments and separating for emphasis such terms from other functions by a dot. Recalling the results of the previous subsection, one appreciates that the spectral densityρ (1) n (x, s/µ 2 F ) follows from the first term on the RHS of Eq. (3.15a). On the other hand, the new contribution R (2) n (x, s/µ 2 F ) in Eq. (3.16b) derives from Eq. (3.15b) and represents one of the main results of this investigation.
To continue, consider the important partial case n = 0 for which the corresponding expression forρ (2β) 0 becomes significantly simplified. Indeed, the leading logarithmic term contributing toR (2) 0 and stemming from Ln 2 (y) [cf. (3.15b)] cancels out because it is proportional to v(n). Actually, this is a general property of the leading logarithmic terms in all orders of the expansion that first reveals itself at the NLO level-see Eq. (3.1a). On the other hand, the subleading logarithmic term ∼ ln (s/µ 2 F ) survives, because it originates from the V where the individual ingredients of this equation are the following  [45] and g [43].
Note that just this term leads to the breaking of the (x ↔x)-symmetry of the spectral Fig. 3, we compare the contributions to the spectral density from the NLO and NNLO β at one single scale Q 2 = µ 2 F = µ 2 SY , the latter scale µ 2 SY = (2.4 GeV) 2 corresponding to the typical average momentum [5] measured by the CLEO Collaboration [9] in the Q 2 region [1.5 − 8] GeV 2 . This is done in terms of 0 (x,x/x) for the NNLO β spectral density (dashed green line) in comparison withρ (1) 0 (x) calculated at the NLO (solid red line).
0 (x,x/x) at the typical CLEO reference scale (labeled by the acronyms of Schmedding and Yakovlev [5] To render our presentation more transparent, we have relegated the discussion of the elements T (3.18)], to the Appendices A and B. The expressions for the partial densitiesR (2) n (x, s/µ 2 F ), appearing as convolutions of these elements with ψ n , are supplied in Appendix C.

IV. SIZE OF THE NNLO CONTRIBUTION TO THE FORM FACTORS
In this section we discuss the effects of the NNLO β corrections to the transition form factors for the ρ and π mesons, relying upon the BLM prescription and its modifications. Employing the results obtained in the previous section, we are going to simplify the expressions for the spectral density by adopting the so-called 'default' scale setting µ 2 F = µ 2 R = Q 2 . Then, the expression for ρ (2β) n , given by Eq. (3.16b), reduces to Fig. 3 It turns out that the NNLO β contribution, calculated here, is negative (dashed line in Fig. 3).
Hence, taken together with the already known NLO contribution, which is also negative, the total effect of the radiative corrections at the considered level of the perturbative expansion is to decrease the magnitude of the form factor. Following the BLM procedure, the last term in Eq. (4.2) determines the "shift" of the scale in the argument of the running coupling from the value Q 2 to the BLM-scale [42], Q 2 BLM , according to . As a result, a s (Q 2 ) → a s (Q 2 BLM ) > a s (Q 2 ) at Q 2 BLM (Q 2 ) < Q 2 and, hence, the form factor given by Eq. (4.2) assumes the BLM-improved form The main contribution to F γ * γπ LCSR is provided by the asymptotic DA, Fig. 4, left panel). This scale may be considered as the borderline for applying perturbative QCD, defining this way some minimal scale for the BLM scheme-denoted Q 2 min . It is remarkable (though accidental) that, as mentioned above, this scale corresponds approximately to µ 2 SY = 5.76 GeV 2 . Obviously, below this particular scale, the BLM prescription, expressed through Eq. (4.3), would entail a renormalization scale that will be out of the region where perturbation theory can be safely applied.
The scale Q 2 ⋆ , determined above, still belongs to the perturbative regime and, therefore, the ratio can serve as a crude measure for the relative weight of the NNLO contribution. In Fig. 4 (left panel), we show Q 2 BLM (Q 2 ) vs. Q 2 (dashed red line in the middle), following from Eq. (4.3), in comparison with a simpler linear dependenceQ 2 BLM = Q 2 exp(−1.811) (solid green line) that results from the standard factorization formula of the perturbative approach employed in [7]. One observes that both results are rather close to each other in the important CLEO-data region. To be able to apply the BLM procedure below the minimal scale Q 2 min , we shall use a somewhat improved version of this procedure-termed BLM-introduced in [36]. Above the perturbation-theory borderline Q 2 ≥ Q 2 min , this modified BLM procedure coincides on the RHS of Eq. (4.4) with the standard one. But for Q 2 < Q 2 min , the BLM scale is frozen at Q 2 ⋆ and the expanded expression contains only a tail of the NNLO β correction provided by the third term in the equation below Substituting this expression for F γ * γπ BLM into Eq. (4.5), one obtains the quantity δ BLM , shown in the right panel of Fig. 4. The numerical value, estimated in Eq. (4.5), turns out to overestimate the size of the NNLO contribution, as one observes from this figure. In mathematical terms this becomes evident by glancing at the ratio and recalling that the magnitude of this contribution to the total form factor is a few times smaller than δ BLM at the moderate values of Q 2 characterizing the CLEO-data region. The size of the NNLO correction seems to be rather important and at the level of about 10% at low Q 2 , even by taking it into account only in the incomplete form of Eq. (4.2). On the other hand, at higher momenta, inspection of the δ 2 behavior in the right panel of Fig. 4 reveals that the size of the corrections rapidly decreases to the level of 5% around the scale µ 2 SY .

V. PREDICTIONS AND COMPARISON WITH EXPERIMENTAL DATA
This section contains a discussion of the implications of our theoretical findings on the transition form factors F γ * γπ and F γ * ρπ vis-a-vis the experimental data for the former. To understand the influence of the NNLO radiative correction on these hadronic observables, we have to analyze its relative weight with respect to the NLO contribution. Moreover, we have to discuss the interplay between perturbative corrections and nonperturbative ingredients, notably, the quark virtuality λ 2 q and the twist-four scale δ 2 . This discussion can be further substantiated by comparing the calculated photon-to-pion transition form factor with the available experimental data from measurements by the CLEO [9] and the CELLO [10] Collaborations. From this comparison, we can extract valuable information as to what extent our calculation can describe the data in the whole measured Q 2 region. From the theoretical point of view, we can use these data in order to estimate what is still missing on the theoretical side. We will focus below not on the exact phenomenological description of the mentioned data, but analyze instead the ramifications they imply on the theoretical approach and its various elements.

A. Nonperturbative input
We discuss first the nonperturbative ingredients of our analysis. The main one is the leading twist-two pion DA, ϕ π , which can be derived with the help of various methods. These include-among others-QCD sum rules and lattice simulations. In addition, one has to model the twist-four component of the pion DA, see, e.g., the discussion in [19,21,22]. Lacking a complete derivation of the full pion DA from first principles of QCD, we are actually forced to reverse-engineer its structure from calculations of its first few moments ξ n π .
To be more precise, one can calculate the moments of ϕ π with standard QCD SR [46] and also with those employing nonlocal condensates (NLC) [18,47,48,49,50]. On the other hand, one can use LCSR to analyze the high-precision CLEO data [9] on F γ * γπ and extract rather strict constraints on the first Gegenbauer coefficients a 2 and a 4 [5,6,19]. More recently, two independent Collaborations have published results for the first coefficient a 2 by measuring the first moment ξ 2 π of ϕ π on the lattice [28,29,51]. The lattice calculation of the second coefficient a 4 ( or, equivalently, the fourth moment ξ 4 π ) is still lacking, but a compatibility region between the CLEO data [9] and the a 2 -lattice constraints was worked out in [26] to predict a rather narrow interval for the ξ 4 π values.  [28] and solid ones [29] are also shown together with predictions obtained from nonlocal QCD sum rules (slanted green rectangle) [18]-the latter corresponding to the vacuum quark virtuality λ 2 q = 0.4 GeV 2 . Right: CLEO data in comparison with various theoretical models and lattice results listed in Table I. The dashed green 1σ ellipse (together with the shifted best-fit point &) describes the effect of including the twist-four contribution to the pion DA via renormalons [22].
We provide below a short overview of our present knowledge of ϕ (2) π from different sources, omitting specific details for which we refer the interested reader to the original literature. The models shown below in the figures are summarized in Table I. Using QCD sum rules with NLC, we derived a "bunch" of admissible pion DAs [18], taking into account in the expansion (2.10) only the first two terms with a 2 and a 4 (details can be found in [35]). Pion DA models, like the Chernyak-Zhitnitsky (CZ) model [46], or the Braun-Filyanov (BF) [38] one, also used these two harmonics for modeling the pion DA. Note, however, that in the NLC approach [18] this is not the result of an arbitrary truncation of the Gegenbauer expansion after n = 4, but follows from the fact that all calculated higher-order coefficients up to n = 10 turn out to be compatible with zero. Hence, from a pragmatic point of view, in order to capture the main characteristics of the pion DA, it is sufficient to restrict the analysis of the experimental data on the photon-to-pion transition to two-parameter models. Such a data analysis of the CLEO measurements was first performed by Schmedding and Yakovlev (SY) [5] within the framework of LCSR. In a subsequent series of papers [6,19,22,25,26,35] (nicknamed BMS), this type of data processing was further pursued with results confirming the previous SY findings while also improving the theoretical accuracy. These results are displayed in Fig. 5 with details being provided in Table I. In this graphics (left panel), the slanted shaded (green) rectangle represents the BMS DA "bunch" from the NLC approach together with the middle point ", which corresponds to the BMS model with the coefficients a BMS 2 (µ 2 SY ) = +0.14, a BMS 4 (µ 2 SY ) = −0.09. The CLEO-data constraints are shown in the form of error ellipses: 1σ (thick solid green line), 2σ (solid blue line), and 3σ (dashed-dotted red line). The range of values of a 2 , determined recently on the lattice by two independent Collaborations, is denoted by vertical dashed lines [28] and solid ones [29]. All constraints and predictions shown have been evolved to the scale µ 2 SY using the NLO ERBL evolution equation.
In the right panel of Fig. 5 we show the effect of the twist-four contribution to the pion DA, taken into account via the renormalon approach of [60], in comparison with the standard obtained by different methods and for several pion DA models. The designations correspond to those used in Fig. 5. Also included are the theoretical constraints derived from QCD SRs with NLC and estimates derived from an analysis [6] of the CLEO data [9] using LCSRs. The lattice measurements of [29], [28], and [51] are also shown. [Note that the uncertainties of the coefficients a 2 and a 4 are correlated. Here, the rectangular limits of the fiducial ellipse [6,19] are shown.] DA methods/models Symbols a 2 (µ 2 SY ) a 4 (µ 2 SY ) Lattice UKQCD/RBC [29] 0.215 ± 0.07 -QCDSF/UKQCD [28] 0.19 ± 0.11 - [51] 0.233 ± 0.145 - one, which is based on the asymptotic DAs of twist four [14], namely,

Data analysis
Here δ 2 (1GeV 2 ) = (0.19 ± 0.02) GeV 2 . This estimate has been obtained in [6] under the assumption that λ 2 q ≡ q (ig σ µν G µν ) q /(2 qq ) = 0.4 ± 0.05 GeV 2 (see Appendix A in [6]). A full-fledged analysis of the renormalon-model corrections has been given in [22]-see also [21]. The net effect is to shift the whole 1σ error ellipse (broken contour in Fig. 5) further away from the asymptotic pion DA and hence to higher values of a 2 (observe the shift of the best-fit point to its new position denoted by &). The original 1σ error ellipse (solid contour) is also shown in this figure and one appreciates that the variation of the twist-four term can have substantial influence on the transition form factor.
Let us conclude this subsection with the following remarks: (i) The most striking message of Fig. 5 (to be read in conjunction with Table I) is that the CLEO 1σ error ellipse, the two recent lattice calculations, and the fiducial region of pion DAs extracted from nonlocal QCD sum rules all have a common region of validity. Moreover, the BMS model DA " is entirely within this area. (ii) In contrast, the asymptotic pion DA x is outside the 3σ ellipse, while the CZ pion DA s, as well as the BF model w, are outside the 4σ error ellipse of the CLEO data. (iii) All models inside the "rhombus", determined in [58], are also more or less outside the 1σ error ellipse. This fact gains more weight in view of the reduced interval of a 2 computed on the lattice [29] (in the form of ξ 2 π ) which seems to exclude all these models as well. (iv) Analogous considerations for the moment ξ 4 π have been discussed in [26].

B. Phenomenological input of the light-cone sum rules
The next task concerns the model of the resonances entering the dispersion integral for the transition form factor. As we mentioned before, we are going to refine the simple δfunction ansatz by taking into account a finite width in terms of the Breit-Wigner model. To this end, we replace ρ h in ρ phen (s, with f ω ≃ f ρ /3, and F γ * ωπ ≃ 3F γ * ρπ . Then, we trade the simple δ-function resonance model for the Breit-Wigner ansatz [14], given by also taking into account the difference in the masses and widths of the ρ and the ω vector mesons, m ρ = 0.7693 GeV, m ω = 0.7826 GeV, and Γ ρ = 0.1502 GeV, Γ ω = 0.00844 GeV, respectively. To continue, in order to obtain F γ * ρπ , we appeal to the duality between the hadronic part of the spectral density and its perturbative counterpart that we express via The main effect of using Eq. (5.4) on the LHS of Eq. (5.5) is a slight suppression of the integral relative to the outcome of the simple δ-function ansatz.
After the Borel transformation of Eq. (5.5), one arrives at an expression for F γ * ρπ in terms of V [cf. (2.23)], notably, where k ≈ 0.932 weakly depends on the Borel parameter M 2 in a region containing the standard scale M 2 ≃ 0.7 GeV 2 . Therefore, the Breit-Wigner ansatz for the resonances (abbreviated in what follows by BW) supplies a factor k −1 > 1, entailing an increase of the value of the form factor as compared to the simple δ-function ansatz. A similar enhancement effect in analogy to Eq. (2.20) applies to the total form factor F γ * γπ LCSR , leading to where k 1 ≈ 0.984. Crudely speaking, the net enhancement amounts to 2-4% in that momentum region, where the resonance part prevails over the pointlike one.
C. NNLO effects on Q 2 F γ * γπ -integral characteristics Having set up the framework for calculating the pion-photon transition form factor, including also its nonperturbative ingredients, we now turn our attention to the specific effects due to the NNLO β radiative corrections according to Eq. (4.2).
Being interested mainly in the magnitude of the form factor, it is actually sufficient to include only the ψ 0 term. An analysis including higher harmonics will be presented in a future publication.
The final results for the photon-to-pion transition form factor, including the NNLO β corrections, within the improved LCSR approach are displayed in Figs. 6 and 7. The presented results have been obtained along the lines followed in our previous dedicated and detailed works in [6,19,22] and will not be repeated here. Recalling the features of the pion DAs discussed in Subsec. V A, we reduce our discussion to only three models: As (lower black line), CZ (upper red line), and BMS (green line in the middle). The solid lines represent the scaled form factor with the NNLO β and the BW-ansatz corrections included, whereas the FIG. 6: Theoretical predictions from LCSR for the form factor Q 2 F γ * γπ with the NNLO β corrections (together with the BW resonance model) included (solid lines) and without them (dashed lines). All predictions shown are evaluated with the twist-four parameter value δ 2 Tw−4 = 0.19 GeV 2 [6,19]. They correspond to selected pion DAs. These are: CZ model-upper (red) line [46], BMS-model-middle (green) line [18], and As DA-lower (black) line. For comparison, the corresponding predictions for each pion DA model without these corrections are displayed as dashed lines. The experimental data shown are from the CELLO (diamonds, [10]) and the CLEO (triangles, [9]) Collaborations.
dashed lines (with the corresponding color for each model) show the result without these corrections. These theoretical predictions are displayed in the background of the experimental data from the CLEO [9] (triangles) and the CELLO (diamonds) [10] Collaborations.
The main lesson from these figures is that the effect on Q 2 F γ * γπ , induced by the radiative corrections, amounts to about (−10 ÷ −5)% (see Sec. IV), whereas that caused by the use of the more realistic BW-resonance ansatz provides a small growth of approximately (+4÷+2)%. Combining both effects, results into a net reduction of the magnitude to within 7% at small Q 2 ≃ 2 GeV 2 . In the case of the As DA, the size of this reduction rapidly drops to 2.5% for Q 2 ≥ µ 2 SY , as one may appreciate by comparing the solid line (NNLO β ) with the dashed line (NLO) at the bottom of Fig. 6. These results have been obtained with the spectral density ρ On the other hand, the transition form factor is mainly determined via the inverse moment [35], which belongs to the integral characteristics of the pion DA, i.e., that can be easily understood from the expression for the H part in Eq. (2.25) setting Q 2 ≫ s 0 , x 0 → 1. This means that for pion DAs, like the BMS one, which have the particular property a 2 ∼ −a 4 (see the anti-diagonal in Fig. 5), the corrections associated with the higher harmonics ψ 2 and ψ 4 mutually cancel, leaving only the correction due to ψ 0 . The same conclusion can be drawn for all pion DAs belonging to the BMS "bunch"-shown as a green strip in Fig. 7 FIG. 7: The same designations as in the previous figure apply. The shaded (green) strip denotes the form-factor predictions derived with the help of the BMS "bunch" [18]. The dashed blue line represents a dipole-form interpolation formula [9] of the CLEO data.
We observe from Figs. 6 and 7 that the improved theoretical predictions which contain the discussed effects fail to describe both sets of the experimental data at low Q 2 , say, below 4 GeV 2 . Actually, the data can be grouped into three regions: (i) The low Q 2 region is mainly covered by the CELLO data and extends to momenta up to approximately 1 GeV 2 . This region, where hadronization is immanent, is virtually inaccessible to the methods used in our analysis-let us call it therefore the unknowable regime. (ii) There is some intermediate momentum region between 1 and 3 GeV 2 , where the existing data are underestimated by our theoretical predictions. In contrast to the previous case, here we can figure out what the origin of this drawback may be. Hence, let us call this intermediate domain the unknown regime, because the missing contributions can, in principle, be estimated. (iii) Finally, above about 3 GeV 2 , the agreement between the CLEO data and the shaded BMS strip is fairly good, as one also realizes by comparing these data points with the dipole fit (dashed blue line) used by the CLEO Collaboration. In fact, the dipole curve and the prediction due to the BMS model (solid curve within the BMS strip) almost coincide in this Q 2 region. Therefore, it is safe to claim that this high-momentum domain is well-reproduced by our techniques-hence the term known regime. 5 4 Note that the width of the BMS strip is somewhat narrower compared to our previous results in [35] because here we use a smaller range for the δ 2 uncertainties in the twist-four contribution. 5 However, one may face the challenge posed by the recently released BaBar data [61] exactly in this momentum-transfer region-see Note Added below.
Let us now discuss the reasons for the discrepancy in the unknown regime.
• One reason for the observed discrepancy may be traced to the fact that we used for the evolution of the twist-four contribution only the one-loop anomalous dimension. Considering the evolution effect at the two-loop level, could potentially reduce its size, thus rendering the reduction of the transition form factor less severe.
• The uncalculated remnant of the NNLO contribution could eventually enhance its total magnitude-should this part appear with the opposite sign with respect to the calculated β-function part.
• Still another source of uncertainty lies with the value of δ 2 which controls the size of the twist-four contribution. A smaller value of this parameter would entail less suppression of the form factor, given that this nonperturbative contribution has a negative sign, just like the calculated NNLO β contribution. Here there is a subtlety. If one assumes a smaller value of δ 2 , then this would automatically mean that also λ 2 q should assume a smaller value because δ 2 ≈ λ 2 q /2, i.e., these two nonperturbative parameters work synergistically. But a decrease of the latter parameter would yield to an enhancement of the leading-twist contribution, meaning that the whole BMS strip (made up on the basis of the BMS "bunch" [18,35]) would move somewhat upward as a whole. As a result, one would obtain an increase of the transition form factor exactly in that regime between 2 and 4 GeV 2 , providing a better agreement with the CLEO and CELLO data.
D. NNLO effects on Q 4 F γ * ρπ (Q 2 )-differential characteristics The scaled form factor Q 4 F γ * ρπ (Q 2 ), obtained in the framework of LCSRs, see Eq. (5.6), is determined by the hadronic part V , defined in Eq. (2.24). At moderate Q 2 around the scale µ 2 SY , V is mainly formed by the ϕ π -dependent leading-twist contribution and to a lesser extent by the twist-four one. The combined effect of the BW-ansatz and the NNLO β corrections is rather important and amounts to approx. −10% in the region limited from above by the scale µ 2 SY . For still higher momentum values, and up to Q 2 ≈ 15 GeV 2 , it reaches the level of +9%-see left panel of Fig. 8. More specifically, the NNLO β corrections reduce the value of the form factor by an amount of approx. 20%, starting at Q 2 = 2 GeV 2 , and become significantly smaller at the end of this region (see dashed line in Fig. 8, left  panel), whereas the use of the BW-ansatz leads to an overall increase of 8%. The ratio of the combined effect of the NNLO β and the BW-ansatz contributions, relative to the total transition form factor, is illustrated in Fig. 8 in terms of a solid red line, whereas the normalized contribution of the NNLO β term alone is depicted in the same figure by the dashed blue line.
The form factor Q 4 F γ * ρπ (Q 2 ) is presented in Fig. 8 (right panel) using different pion DAs. We observe from this figure that this process can actually be used to discriminate among different pion DAs of twist two, as first pointed out by Khodjamirian in [14] because it is sensitive to the particular shape of the pion DA. Especially the momentum region between 2 GeV 2 and 8 GeV 2 seems to be particularly convenient for this task because (i) the obtained predictions are clearly distinguishable and (ii) because it corresponds to the range probed already by the CLEO experiment for the pion-photon transition [9]. The measurement of both transition form factors in this momentum region would provide definitive clues for the underlying pion DA. Let us now make some detailed remarks on the structure of the BMS strip, obtained from the "bunch" of pion DAs following from NLC QCD sum rules. The first observation from Fig. 8 (right panel) is that the BMS strip (which includes the BMS model: long-dashed line) [18] crosses around 11 GeV 2 the prediction obtained with the As DA. We argue that this effect can be traced back to the differential pion characteristics of the pion DA, viz., The usefulness of the slope of the pion DA was first discussed in [39] in connection with the pion's electromagnetic form factor. In our case, it is instructive to recall the definition of V in Eq. (2.24), from which one appreciates that the lower limit of the dispersion integral x 0 tends to the upper one, 1, when Q 2 ≫ s 0 . Therefore, near the upper integration limit, this integral is determined by the endpoint behavior ofρ(x)| x∼1 , which in LO is proportional to ϕ ′ π (0) [cf. Eq. (2.9a)], entailing for the form factor a behavior like ∼ ϕ ′ π (0) This becomes visible in the vicinity of the upper limit 1, say, for a value 0.1 = s 0 /(s 0 + Q 2 * ), which for the BMS model corresponds to Q 2 * ≈ 14 GeV 2 , while for the upper part of the strip the value of Q 2 * is much larger. Thus, it becomes possible to investigate the endpoint behavior of different pion DAs-mentioned also in [14,38]-in terms of this transition form factor by probing the slope of the pion DA at the origin.
(m + 1)(m + 2)/2 a m (µ 2 ) . (5.9) Then, we get the following results: These simple estimates are qualitatively responsible for the associated form-factor values at Q 2 * (slopes and altitudes) in Fig. 8 (right panel). It turns out that for the CZ DA the value of the form factor is approximately 3 to 4 times larger than the one for the As DA, while for the BMS DA it is smaller or quite close to it [cf. Eq. (5.10)]. Note that this feature appears to be opposite to the Q 2 F γ * γπ (Q 2 ) form factor that mainly depends on the inverse moment [22] x −1 π , the latter being an integral pion characteristic [18,19]. Indeed, one sees from Fig. 7 that the results for Q 2 F γ * γπ (Q 2 ) (green-shaded strip) are always larger than the predictions obtained from the asymptotic DA. Hence, Q 4 F γ * ρπ (Q 2 ) can provide complementary information about the pion DA and help discriminate among various proposed pion DA models.
Going beyond the leading order ofρ, a new behavior near the endpoints emerges due to the hard-gluon exchange leading to a 1/Q 2 -behavior of the form factor in the far asymptotic domain. This is, because in NLO, the expressionρ (1) in Eqs. (3.11) and (3.12) contains a term proportional to − ln 2 (x/x) that accumulates this effect.
A second issue related to the BMS strip in the right panel of Fig. 8, which deserves to be considered as well, concerns the fine structure of its envelopes. Indeed, close inspection reveals that the upper envelope has a dip at Q 2 ≈ 9 GeV 2 , whereas the analogous irregularity of the lower envelope is much less pronounced and appears at a scale close to 10.5 GeV 2 . The origin of these irregularities might be related to the crossing of the predictions obtained from different DAs inside the BMS "bunch" that have different initial values of the Gegenbauer coefficients a 2 and a 4 (exemplified by the intersecting lines inside the BMS strip). When it happens that such a crossing point lies just near the upper or lower boundary of the BMS strip, the corresponding envelope is "bent" inwards. Physically, this means that the range of the calculated uncertainties is somewhat smaller there relative to the regions where no crossing points appear close to the boundaries. On the other hand, we cannot exclude that these irregularities are a spurious effect created by the algorithm we used to obtain our predictions. In any case, they have no influence on our analysis.

VI. CONCLUSIONS
Results have been presented for the photon-to-pion and the ρ-meson-to-pion transition form factors using light-cone sum rules and including those NNLO contributions which are proportional to the β function. Let us summarize our main findings.
1. The spectral density of the LCSR at the NLO level was systematically constructed for any Gegenbauer harmonic of order n and was presented in compact form in Eqs. (3.1) and (3.12).
2. Using the NNLO hard-scattering amplitude T (2) β , proportional to the β function and calculated before in [7], we derived all necessary ingredients of the NNLO spectral density ρ  (3.20). These quantities enter the dispersion integral which determines the structure of the form factors F γ * ρπ and F γ * γπ within the light-cone sum-rule approach.
3. Predictions for the form factor F γ * γπ were presented which include the NLO and the β-part of the NNLO radiative correction, an improved resonance contribution, and also estimates of the uncertainties stemming from the twist-four contributions.
The first contribution has a negative sign and amounts to a reduction of up to 10%, whereas the Breit-Wigner ansatz provides a small enhancement below 5%, and the twist-four contributions can yield to a small enhancement by adopting a smaller value of the parameter δ 2 which controls their size. This, however, would entail a parallel enhancement of the leading-order contribution to F γ * γπ via the parameter λ 2 q , because these two parameters are correlated (δ 2 ≈ λ 2 q /2). Comparison with the CLEO and the CELLO data shows that at 2 GeV 2 , and below, the calculated form factor falls slightly short compared to the data, leaving room for additional soft and higher perturbative contributions that are unknown at present. 4. We have also given predictions for the F γ * ρπ form factor using different pion distribution amplitudes including the asymptotic one, the whole set of pion DAs derived from nonlocal QCD sum rules, and the CZ model. In this case, the combined effect of the negative NNLO β radiative corrections and the use of the Breit-Wigner ansatz lead to a net reduction of the form-factor magnitude varying between −10%, below the scale µ 2 SY , and growing to the level of +9% at Q 2 ≈ 15 GeV 2 , whereas the twistfour contributions amount to −11% even up to the scale 15 GeV 2 . Once there will be experimental data for this reaction, it will provide an additional and useful tool to discriminate among various pion DAs-especially between those of the CZ type that receive strong endpoint enhancement on the one hand and such with their endpoints being suppressed-like the BMS one-on the other. Moreover, because this form factor is sensitive to the differential characteristics of the pion DA-expressed via ϕ ′ π (0)-one can use the endpoint behavior of the underlying pion DAs as an additional adjudicator in selecting the optimal pion DA.
To conclude, we have extended the analysis of the pion-to-photon transition process to the NNLO level in a systematic though partial way. Nevertheless, our analysis provides the possibility to estimate the size of the associated form factors in that region of momenta which is accessible to present measurements. It will be interesting to pursue and complete this sort of calculation in the future by including the whole NNLO contribution.

Note Added
After completion of this work we became aware (thanks to Dieter Müller and Maxim Polyakov) of new preliminary data of the BaBar Collaboration 6 on two-photon-induced processes in the momentum range 4 < Q 2 < 40 GeV 2 which show above 10 GeV 2 "a powerlaw growth behavior that contradicts most models for the pion DA". Meanwhile, these data have been officially released [61].
Taking the high-Q 2 data points of BaBar at face value, one might even come to the conclusion that they are incompatible with QCD per se because the indicated power-law enhancement at large momentum values cannot be explained by any known QCD effect, like higher-order radiative corrections or contributions due to higher-twist effects. Our presented analysis serves to prove that the inclusion of the main part of the NNLO radiative corrections provides suppression-not enhancement-at the 10% level, while the twist-four contribution gives a small enhancement of a few percent in the CLEO region. The expected size of the uncalculated NNLO remainder-even if it should have a positive sign-is not expected to exceed the few-percent level (< 10%). An enhancement of the twist-four contribution due to two-loop evolution is possible but it is expected to be of the order of a few percent as well. The size of all these QCD corrections is far less than indicated by the high-Q 2 BaBar data, so that the observed enhancement at high Q 2 cannot be explained by higher-order perturbative QCD and power corrections.
Note that endpoint-enhanced pion DAs, like the CZ model, are also in conflict with these data, despite opposite claims in [61], because the corresponding Q 2 F γ * γπ (Q 2 ) prediction (see Fig. 9) scales with Q 2 (analogously to all other pion DA models in the convolution scheme) above approximately 15 GeV 2 , in sharp contrast with the significant growth of the BaBar data in this momentum range. In this context, we find it remarkable that the BaBar data points in the momentum range already probed by the CLEO Collaboration are compatible with the CLEO dipole fit and the asymptotic QCD prediction √ 2f π , being also within the BMS strip. Even more significantly, two more data points-outliers-at about 14 GeV 2 and 27 GeV 2 turn out to be just on the upper boundary of the BMS strip and, hence, in compliance with the QCD expectations. (see Fig. 9). Hence, one may divide the BaBar data into two branches: one containing the data points in the CLEO region plus the two outliers, the other consisting of the remaining 10 high-Q 2 data points. The first branch supports the QCD predictions with NNLO radiative corrections and twist-four contributions. The other branch is in clear conflict with the convolution scheme of QCD. In view of this data structure, odds are that the BaBar data may bear some intrinsic inconsistency.
We end this discussion with the following key observations: (i) The main NNLO radiative corrections and the twist-four contributions do not provide enhancement to Q 2 F γ * γπ (Q 2 ) in the range of momentum transfer 10-40 GeV 2 , exclusively covered at present by the BaBar experiment. Hence, the observed behavior of Q 2 F γ * γπ (Q 2 ) growing with Q 2 in this momentum region with an almost constant slope cannot be explained within the convolution scheme of QCD. (ii) Within the QCD convolution scheme, all pion DA models, which have a convergent projection onto the eigenfunctions (Gegenbauer polynomials) of the meson evolution equation and hence vanish at the endpoints 0 and 1, are conflicting with the BaBar data for II: Deviation in terms ofχ 2 ≡ χ 2 /ndf (ndf = number of degrees of freedom) of Q 2 F γ * γπ (Q 2 ) predictions for the asymptotic (Asy), the BMS DA, and the CZ one. For a direct comparison with the BaBar analysis, we employ the NLO approximation of perturbative QCD, including also twistfour contributions. The first column shows the results for the combined sets of the CLEO [9] and the BaBar [61] data. The second column refers only to the BaBar data, while the third column takes into account only the last 10 high-Q 2 BaBar data, starting with the data point at 10 FIG. 9: Predictions for Q 2 F γ * γπ (Q 2 ) derived from three different pion DA models: Asymptotic, BMS, and CZ. The same designations as in Fig. 7 are used. The momentum-transfer range is extended to 40 GeV 2 in order to include the new BaBar data [61] (shown as thick bullets with error bars). The displayed theoretical results include the NNLO β radiative corrections and the twist-four contributions. The horizontal dashed line represents the asymptotic QCD prediction √ 2f π .
The bottom line: The anomalous high-Q 2 -behavior of the pion-photon transition form factor, found by BaBar, could potentially be of great importance, if confirmed by independent measurements, e.g., by the BELLE Collaboration, and would demand a new framework of analysis within QCD. Some scenarios to explain the BaBar effect have already been proposed [62,63,64].
APPENDIX A: DETERMINING THE SPECTRAL DENSITY FROM THE DIS-CONTINUITY OF THE HARD-SCATTERING AMPLITUDE Let us first define the task ahead in a bottom-to-top approach, starting with the LO case and finishing with the NNLO one. First, we analyze the structure of the amplitudes themselves, then, we construct the discontinuity for their elements in order to pair each T k with the spectral density, ρ (k) , associated with this order of the perturbative expansion. At the LO level we encounter just one single contribution to the discontinuity, i.e., that one arising from T 0 ⊗ f (where f (x) is some appropriate case-dependent test function). In NLO, T 1 ⊗ f , a new element appears: T 0 ⊗ Ln(x) ⊗ f (x). Finally, at the NNLO level of the expansion, T 2 , we are faced with the discontinuity of still one more quantity, notably, [Recall that Ln(x) denotes the logarithm of the photon momenta over the factorization scale-cf. (2.7).] Before we continue, we mention parenthetically in this context that the contributions to ρ (k) , just mentioned, can be also obtained from a generating function R(s; ε), which we display below: The total sum of these contributions,R(a s , x; s, Q 2 ), (following from Eq. (A.1)) enters the spectral density in the form of a convolution with V (0) + . Therefore, the resummed expression does not contribute to the ψ 0 part by virtue of the current conservation v(0) = 0.
1. From the discontinuity of T 0 in LO we have [5,14]

2.
In NLO we obtain [6,7] where the notations of Ref. [7] and the following abbreviation have been used: HereQ 2 = − [(q 1 − q 2 )/2] 2 = (Q 2 + q 2 )/2, and µ 2 R and µ 2 F denote, respectively, the scale of the renormalization of the theory and the factorization scale of the process. The kernels V a and V b are diagonal with respect to ψ n and are defined in Appendix B, whereas the kernel g reads and is not diagonal with respect to the Gegenbauer-polynomials. This kernel is responsible for the apparent breaking of conformal symmetry in the MS-scheme [7]. Substituting (A.5a) and (A.5c) into Eq. (A.4), we arrive at One observes that in this expression a simpler logarithm Ln(y) ≡ ln [(Q 2 y + q 2ȳ ) /µ 2 F ] appears in comparison to LN in Eq. (A.5c). From the former expression and Eq. (A.5b), we get the following result for the convolution T 1 ⊗ ψ n = C F T 0 (Q 2 , q 2 ; y)⊗ − ψ n (y) · 3 1 + v b (n) + [g + ⊗ ψ n ] (y) + ψ n (y) · 2v(n) Ln(y) , where the second term in Eq. (A.8) reads The RHS of Eq. (A.9) can be further evaluated using the following relation The discontinuity of T 1 ⊗ ψ n in Eq. (A.10) consists of one part, determined by the discontinuity of T 0 and following directly from Eq. (A.3), and a second nontrivial part entailed by the logarithm, as one can also verify from the analysis in [5]. Then one has  3. In NNLO we extract the b 0 -proportional contribution by collecting all terms in the general structure of T 2 computed in [7]. The result is where T (2) β (ω, x) = T 0 (Q 2 , q 2 ; y) ⊗ T (2) β (y, x) + LN(ω, y) (V We use here for the elements of T 2 the notation T (2) β , T β , T (1) and V (1) β , introduced in [7], that differs from the previous one by a factor of 2. We have T (2) β (x, y)=

12
V a +V a − 209 36 V − 7 3V − 1 4V + 19 6 g +ġ + (x, y) − 6δ(x − y) , (A. 16) while T (1) is defined in (A.5b). One observes that these T -kernels are calculated in terms of the kernels V and g, which enter the evolution kernel V (1) , and the derivativesV ,V that will be discussed in Appendix B. The origin ofġ has been clarified in Sec. 3 of Ref. [7]. Under the assumption that T 2 → b 0 · T β , it turns out that the entire LN(ω, y) dependence appears only inside the term Ln(y) [cf. Eq.(2.7)], in analogy to the NLO case, see (A.7): (A.17a) T β = C F T 0 (Q 2 , q 2 ; y) ⊗ − ln This important property of the T β structure has already been mentioned in [7]. While all terms in (A.17b) contribute to the discontinuity of T β , only the last term, which contains the square of a logarithmic expression, contributes a new type of discontinuity. On the other hand, the term proportional to ln(µ 2 F /µ 2 R ) equals T 1 , as it can be seen from Eq. (A.7). Thus, the final structure of T β assumes the form From this expression, one sees that the first term, which contains no logarithm at all, has no influence on the discontinuity of T 0 . On the other hand, the remaining terms contain logarithms, which do affect the discontinuity, with the last one being the new contribution first appearing at the NNLO level. Then, the partial amplitude T β ⊗ ψ n reads Substituting identity [5] 1 π Im ln 2 (y − y 0 ) y − y 0 = δ(y − y 0 ) ln 2 (y 0 ) − π 2 /3 − 2 θ(y 0 > y) ln(y − y 0 ) y − y 0 + (A. 21) into Eq. (A.20), we obtain for its RHS 1 2(Q 2 + s) This Appendix compiles the crucial properties of the evolution kernels borrowing results from [65,66] and [7]. Following [67,68], we introduce the auxiliary kernels V a,b (x, y; λ), viz., and their sum V (x, y; λ) = V a + (x, y; λ) + V b + (x, y; λ) that includes the main of the logarithmic contributions-generated by the one-loop renormalization of the running couplingaccompanied by the factor a s ln(x/y). Just the effects of this renormalization lead to those contributions that are proportional to b 0 and pertain to the derivatives of the auxiliary kernel V (x, y; λ), i.e., 2V ′ (x, y; λ) | λ=0 ≡V (x, y) = 2 θ(y > x) x y 1 + 1 y − x ln x y + x →x y →ȳ + , (B. These contributions can be obtained from the generating kernel [66] V β (x, y|λ) = (1 + λ)V a (x, y; λ) + V b (x, y; λ) + C(λ) (B.6) for any order of b 0 . Here C(λ) is an analytic function in the variable λ with C(0) = 1.
To obtain the b 0 -contribution at any desired fixed order of the parameter a s b 0 , one has to expand the kernel V β (x, y|a s b 0 ) in a Taylor series with respect to a s up to this order. Hence V β (x, y|0) = V + (x, y), while the first differentiation of V β with respect to a s , β (x, y) , (B.7) leads to V (1) β . The generalized kernel V β (x, y|λ) has been derived from the diagrams for the ordinary one-loop kernels, generated by replacing single gluon lines by a sum of renormalonchain insertions. The coefficient C(λ) in Eq. (B.7) accumulates non-logarithmic parts of these renormalon-chain contributions and was determined in [65,66].

APPENDIX C: MAIN ELEMENTS OF THE NNLO PARTIAL AMPLITUDE
We present now the main elements of the partial amplitudes T β ⊗ ψ n and (2V (1) β ) + − T (1) ⊗ ψ n entering Eq. (A. 19) in Appendix A. We split each of the expressions below in two parts: the singular part with x → 0 is extracted in an explicit form that turns out to be proportional to ψ n , whereas the other part contains an integration over longitudinal momentum fractions.