Pion-photon transition form factor in LCSR and tests of asymptotics

We study the pion-photon transition form factor (TFF) $F^{\gamma*\gamma\pi^0}(Q^2)$ using a state-of-the art implementation of light cone sum rules (LCSRs) within fixed-order QCD perturbation theory. The spectral density in the dispersion relation includes all currently known radiative corrections up to the next-to-next-to-leading-order (NNLO) and all twist contributions up to order six. Predictions for the TFF are obtained for various pion distribution amplitudes (DAs) of twist two, including two-loop evolution which accounts for heavy-quark mass thresholds. The influence of the main theoretical uncertainties is quantified in order to enable a more realistic comparison with the data. The characteristics of various pion DAs are analyzed in terms of the conformal coefficients $a_2$ and $a_4$ in comparison with the $1\sigma$ and $2\sigma$ error regions of the data and the most recent lattice constraints on $a_2$ with NLO and NNLO accuracy. Our results provide more stringent bounds on the variation of the pion DA and illuminate the corresponding asymptotic behavior of the calculated TFF.


I. INTRODUCTION
In this work we consider the pion-photon transition form factor F γ * γ * π 0 (q 2 1 , q 2 2 ) for the process γ * (q 2 1 )γ * (q 2 2 ) → π 0 with q 2 1 = −Q 2 and q 2 2 = −q 2 assuming Q 2 ≫ q 2 and adopting a single-tagged experimental set-up. In that case, one measures the differential cross section dσ(Q 2 , q 2 = 0)/dQ 2 for the above exclusive process and selects events in which the π 0 and one final-state electron (or positron)-the "tag"-are registered, while the other lepton remains undetected.
A self-consistent calculation of the TFF within QCD encompasses various regimes of dynamics from low Q 2 1 GeV 2 , where perturbation theory is unreliable and nonperturbative effects are eventually more important but poorly known, up to high Q 2 values where one would expect that the perturbative contributions in terms of a power-series expansion in the strong coupling prevail and provide an accurate dynamical picture within perturbative QCD (see Fig. 1). There are mainly three different sources of nonperturbative effects related to confinement that pertain to the TFF: (i) mass generation due to Dynamical Chiral Symmetry Breaking (DCSB), (ii) the bound-state dynamics of the pion encoded in a light-cone parton distribution amplitude (DA), and (iii) the hadronic content of the quasireal photon that is emitted from the untagged electron (or positron) at large distances and interacts nonperturbatively with the pion. We do not address DCSB in this work, but we refer to other approaches which account for this and use their results in * Electronic address: stefanis@tp2.ruhr-uni-bochum.de  FIG. 1: Illustration of the single-tag π 0 production in a twophoton process with one highly virtual photon γ * (Q 2 ) and a quasireal photon γ(q 2 ∼ 0) emitted from the untagged electron (or positron). The TFF γ * γ → qq → π 0 is shown as the convolution of the hard quark-gluon subprocesses within fixed-order perturbative QCD (see the text) with the pion light-cone distribution amplitude for the pion (shaded oval). the analysis. A reliable theoretical scheme able to include the other two nonperturbative ingredients, together with perturbative radiative corrections and nonperturbative higher-twist contributions, is the method of light-cone sum rules (LCSRs) [1,2] in combination with fixed-order perturbation theory (FOPT) within QCD. This scheme provides computational techniques which can be used in connection with various pion DAs and is particularly useful for the analysis of the experimental data [3,4] that are eventually indicating discrepant observations applying to the same phenomenon; see [5,6] for a detailed comparison of various theoretical approaches and a classification scheme of the predictions.
In this work we present a LCSR-based calculation of the π − γ TFF which contains several new elements relative to previous approaches: (i) The twist-two spectral density includes all presently known radiative corrections up to the next-to-nextto-leading order (NNLO), i.e., up to the order of O α 2 s β 0 [7]. The LCSR also contains the twist-four term and the twist-six contribution [8]. Recently, the method of LCSRs was combined with the solution of the renormalization-group (RG) equation [9] to perform a summation over the radiative corrections and extend its application to momenta Q 2 < 1 GeV 2 . This momentum regime is outside the scope of the present investigation. (ii) The hadronic content of the quasireal photon γ(q 2 ) is included in the LCSR by employing a physical spectral density which models the vector-meson properties of the quasireal photon in terms of ρ/ω resonances by means of a Breit-Wigner form.
(iii) Several pion DAs are used as nonperturbative input and their characteristics are investigated in the (a 2 , a 4 ) plane, where a 2 and a 4 are the first nontrivial coefficients in the conformal expansion of the pion DA. Comparison is given with the 1σ and 2σ error regions created within the LCSR from the combined sets of the CELLO [10], CLEO [11], Belle [4], and BABAR( 9 GeV 2 ) [3] data. The graphical representation of the TFF predictions also includes the recently released preliminary BESIII data [12,13], see also [14]. (iv) All pion DAs used in the TFF predictions are evolved from their normalization scale to the measured momenta using a NLO (two-loop) evolution scheme which takes into account heavy-quark mass thresholds.
(v) A crucial attenuation effect in the conformal expansion of the TFF within the LCSR approach is worked out, which marks a crucial difference to perturbative QCD and is of particular importance with respect to the behavior of the TFF at large Q 2 . (vi) To facilitate the discussion of the asymptotic characteristics of the TFF, a new quantity is introduced which measures the scaling rate of the scaled TFF with Q 2 .
The paper is organized as follows. In Sec. II we present the theoretical formalism to carry out the TFF calculations. We specify the applied LCSR and discuss the important attenuation effect related to the hadronic structure of the quasireal photon. Our results for the pion DAs and the TFF predictions are presented in Sec. III. This section includes a dedicated discussion of the TFF asymptotics. Our conclusions are summarized in Sec. IV. The involved evolution scheme to handle the scale dependence of the pion DA at the two-loop order by including heavy-quark thresholds is explained in Appendix A. Appendix B completes the paper by providing a compilation of the experimental data together with the corresponding TFF values and uncertainties for the BMS [15] DAs. For the first time, the analogous values for the platykurtic (pk) DA [16] are also given.

A. QCD Factorization
The amplitude T µν describing the process γ * (q 1 )γ * (q 2 ) → π 0 (P ) can be defined by the correlation function where j µ = 2 3ū γ µ u − 1 3d γ µ d is the quark electromagnetic current. Expanding the T-product of the composite (local) current operators in terms of Q 2 and q 2 (assuming that they are both sufficiently large), one gets by virtue of the factorization theorem, the LO term denoting the pion DA of twist two. For vanishing q 2 this expression reduces to where we have recast the inverse moment 1/x π in terms of the projection coefficients a n on the set {ψ n } of the eigenfunctions of the one-loop Efremov-Radyushkin-Brodsky-Lepage (ERBL) evolution equation: Here ψ 0 (x) = 6x(1 − x) ≡ 6xx is the asymptotic pion DA ϕ asy π and the higher eigenfunctions are given in terms of the Gegenbauer polynomials ψ n (x) = 6xxC The pion DA parameterizes the matrix element where the path-ordered exponential (the lightlike gauge link) [z, 0] = P exp ig z 0 t a A µ a (y)dy µ ensures gauge invariance. It is set equal to unity by virtue of the lightcone gauge z · A = 0 adopted in this work. Higher-twist DAs in the light cone operator product expansion of the correlation function in (1) give contributions to the TFF that are suppressed by inverse powers of Q 2 . Physically, ϕ (tw-2) π (x, Q 2 ) describes the partition of the pion's longitudinal momentum between its two valence partons, i.e., the quark and the antiquark, with longitudinalmomentum fractions x q = x = (k 0 + k 3 )/(P 0 + P 3 ) = k + /P + and xq = 1−x ≡x, respectively. It is normalized to unity, 1 0 dxϕ (tw-2) π (x) = 1, so that a 0 = 1. The expansion coefficients a n (µ 2 ) are hadronic parameters and have to be determined nonperturbatively at the initial scale of evolution µ 2 , but have a logarithmic Q 2 development via α s (Q 2 ) governed by the ERBL evolution equation, see, for instance, [17] for a technical review. The one-loop anomalous dimensions γ (0) n are the eigenvalues of ψ n (x) and are known in closed form [18]. The ERBL evolution of the pion DA at the two-loop order is more complicated because the matrix of the anomalous dimensions is triangular in the {ψ n (x)} basis and contains off-diagonal mixing coefficients [8,[19][20][21][22][23][24][25]. To obtain the TFF predictions in the present work, we employ a two-loop evolution scheme (App. A), which updates the procedure given in Appendix D of [24] by including the effects of crossing heavy-quark mass thresholds in the NLO anomalous dimensions γ (1) n and also in the evolution of the strong coupling, see, e.g., [26][27][28].
The TFF can be expressed in more general form to read [18,29] where µ F is the factorization scale between short-distance and large-distance dynamics and h.tw. denotes highertwist contributions. The hard-scattering amplitude T has a power-series expansion in terms of the strong coupling a s ≡ α s (µ 2 R )/4π, where µ R is the renormalization scale. In order to avoid scheme-dependent numerical coefficients, we set µ F = µ R ≡ µ (default choice) relegating the discussion of the scheme dependence and the factorization/renormalization scale setting of the TFF to [30,31].
Then we have where the short-distance coefficients on the right-hand side can be computed within FOPT in terms of Feynman diagrams as those depicted in Fig. 1. In our present calculation we include the following contributions, cast in convolution form via (6) with ⊗ ≡ 1 0 dx, where the abbreviation L ≡ ln Q 2 y + q 2ȳ /µ 2 has been used [7,32].
The dominant term is [32,33] T β = T where β 0 = 11 3 C A − 4 3 T R N f is the first coefficient of the QCD β function with T R = 1/2, C F = 4/3, C A = 3 for SU (3) c and N f is the number of active flavors.
Recently, two more contributions to the NNLO radiative corrections have been calculated in [7] to which we refer for their explicit expressions and further explanations. These are while the term T (2) c in (8c) has not been computed yet and is considered in this work as the main source of theoretical uncertainties. Finally, suffices to say that V

B. Light cone sum rules
Let us now turn to the description of the TFF using a dispersion relation within the LCSR approach.
The TFF for one highly virtual photon with the hard virtuality Q 2 and one photon with a small virtuality q 2 ≪ Q 2 can be expressed in the form of a dispersion integral in the variable q 2 → −s, while Q 2 is kept fixed, to obtain where ρ(Q 2 , s) is the spectral density The first term ρ h (Q 2 , s) models the hadronic (h) content of the spectral density, while ρ pert (Q 2 , s) denotes the QCD part in terms of quarks and gluons, calculable within perturbative QCD, Each of these terms can be computed from the convolution of the associated hard part with the corresponding DA of the same twist [2]. Below some effective hadronic threshold in the vector-meson channel, the photon emitted at large distances is replaced in F γ * V π 0 by a vector meson V = ρ, ω, etc., using for the corresponding spectral density a phenomenological ansatz, for instance, a δ-function model.
Thus, after performing the Borel transformation 1/(s+ q 2 ) → exp −s/M 2 , with M 2 being the Borel parame-ter, one obtains the following LCSR (see [7,8,32] for more detailed expositions) where the spectral density is given bȳ For simplicity, we have shown above the LCSR expression for the simple δ-function model to include the ρ-meson resonance into the spectral density. However, the actual calculation of the TFF predictions to be presented below, employs a more realistic Breit-Wigner form, as suggested in [2] and used in [32]. This reads where the masses and widths of the ρ and ω vector mesons are given by m ρ = 0.770 GeV, m ω = 0.7826 GeV, Γ ρ = 0.1502 GeV, and Γ ω = 0.00844 GeV, respectively. The other parameters entering (15) are s =xQ 2 /x with x ≡ 1 − x, x 0 = Q 2 / Q 2 + s 0 , and the effective threshold in the vector channel is s 0 ≃ 1.  (15) includes in an effective way the nonperturbative long-distance properties of the real photon in terms of the duality interval s 0 and the masses of the vector mesons that are absent in the pQCD formulation of the TFF, but play an important role in the kinematic region Q 2 s 0 and x 0 0.5 (cf. the first term in Eq. (15)). The real-photon limit q 2 → 0 can be taken in (15) by simple substitution because there are no massless resonances in the vector-meson channel. Thus, this equation correctly reproduces the behavior of the TFF for a highly virtual and a quasireal photon from the asymptotic limit Q 2 → ∞ down to the hadronic normalization scale of Q 2 ∼ 1 GeV 2 , as measured in singletag experiments. For still lower momenta, outside the validity range of the standard LCSR scheme, other dispersive approaches may be more preferable [9,[36][37][38][39][40]. The pertinent role of subleading power corrections to the TFF has been investigated in [41,42].
Using the conformal expansion for ρ tw-2 , the spectral density can be expressed in the form wherē with the elementsρ (i) n being given in Appendix B of Ref. [7].
The dispersive analysis here includes the twist-four and twist-six spectral densities in explicit form. Theρ tw-4 spectral density is given bȳ where the twist-four coupling parameter takes values in the range δ 2 tw-4 (µ 2 = 1 GeV 2 ) ≈ λ 2 q /2 = 0.19 ± 0.04 GeV 2 and is closely related to the average virtuality λ 2 q of vacuum quarks [43][44][45][46][47], defined by λ 2 q ≡ q(igσ µν G µν )q /(2 qq ) = 0.4±0.05 GeV 2 . Details on its estimation and evolution can be found in [24], whereas the sensitivity of the TFF to its variation was examined in [48]. In the present analysis the evolution of δ 2 tw-4 is also included. Expression (20) is evaluated with the asymptotic form of the twist-four pion DA [2] while more complicated forms were considered in [49,50]. The twist-six part of the spectral density, i.e., ρ tw-6 (Q 2 , x) = (Q 2 + s)ρ tw-6 (Q 2 , s), was first derived in [8]. An independent term-by-term calculation in [7] confirmed this result. We quote it here in the form where the plus prescription [f (x, . To obtain detailed numerical results for the TFF F (Q 2 ) using (15), we employ several DAs from different approaches with various shapes encoded in their conformal coefficients a n . The latter are determined at their native normalization scale (as quoted in the referenced approaches) by means of the moments of the pion DA where ξ = x −x and N = 2, 4, . . .. The expansion coefficients a n can be expressed in terms of the moments ξ N π as follows

C. Attenuation effect
We now turn our attention to an effect that marks a crucial difference between the perturbative approach and the use of a dispersion relation and has important consequences for the scaling behavior of the TFF.
As shown in [7] (see Fig. 1 and Eq. (20) there), and discussed here further, the leading-twist expression for a given harmonic n of the TFF within the LCSR scheme is not the same as in pQCD. The TFF considered in this work is where the terms in the square brackets represent the twist-two contribution. While the result based on factorization is given by the inverse moment ofx with respect to ψ n (x) the analogous LCSR expression deviates from that because of the hadronic structure of the quasireal photon taken into account via the vector meson dominance. It gives instead where the involved parameters are defined below Eq. (17) and the spectral density is given by Eq. (18).
As long as Q 2 ≫ s 0 and x 0 = 1 + s 0 /Q 2 −1 → 1, the hadronic part of the quasireal photon in the spectral densityρ(Q 2 , x) is suppressed at large Q 2 . As a result, all harmonics contribute at once like in perturbative QCD and Q 2 F LCSR n (Q 2 ) → 3 on account of 6 both terms in Eq. (27) contribute with comparable magnitudes. Thus, nonperturbative higher-twist contributions controlled by s 0 in Eq. (27) are no more suppressed. At the same time, the vector-meson generated factor entails an attenuation effect of the conformal expansion so that the ψ n harmonics contribute successively in pace with Q 2 , see Fig. 2. The upshot of this nonperturbative attenuation effect is that the LCSR-based TFF predictions can deviate significantly from those obtained in perturbative QCD.
The core observations from the graphics in Fig. 2 are the following: (i) The strongest contribution to the form factor stems from the zeroth order term F 0 (Q 2 ) which grows uniformly. (ii) The higher partial terms oscillate with zero crossings clustering at Q 2 1 GeV 2 . These oscillations cause an attenuation effect in the sense that harmonics of higher order start to contribute at larger and larger Q 2 values one following the other as n grows. For instance, ψ 12 starts to grow uniformly only beyond the zero crossing around 15 GeV 2 (right panel). (iii) Therefore, the more harmonics with positive coefficients are included in the conformal expansion of the pion DA, the stronger the attenuated enhancement of the form factor becomes as Q 2 grows. (iv) The decrease of the conformal coefficients due to ERBL evolution is only log- arithmic and is thus insufficient to compensate for this enhancement, though at asymptotically large Q 2 values it finally prevails. (v) These considerations apply not only to DAs with a large number of expansion coefficients, they are also valid to a less degree for DAs with a few number of coefficients but having an inverse hierarchy [35]. (vi) On the other hand, at Q 2 < 1 GeV 2 only the term F 0 (Q 2 ) contributes so that, irrespective of how many conformal coefficients are included in the pion DA representation, the total form factor will be dominated by the ψ 0 harmonic.

III. ANALYSIS OF THE RESULTS
In this section we present our results for the pion DA and the TFF predictions calculated within the LCSR approach described in the previous section.

A. Pion DAs
To obtain predictions for the TFF, we evaluate (18) in the LCSR (15) using for the physical spectral density in (13) the Breit-Wigner form (17) and employing various conformal coefficients a n at their normalization scale. We consider two such scales µ 1 = 1 GeV and µ 2 = 2 GeV, depending on the particular pion DA. If µ 2 is not the native normalization scale, the ERBL evolution scheme discussed in App. A is applied. This scheme, presented here for the first time, works for any polynomial order of the conformal expansion and accounts for the crossing of heavy-quark flavors. The coefficients a 2 , a 4 , a 6 of various models for the pion DA are given in Table I at both scales µ 1 and µ 2 . Using these values, one can readily compute the corresponding moments (23) using Eq. (24). This table also includes the values of the inverse moment at the scale µ 2 . Because broad, concave distributions cannot be adequately represented in terms of only the lowest three coefficients, the corresponding inverse moments of the DSE-DB, DSE-RL, and the holographic AdS/QCD DAs are calculated within the α − representation given by Eq. (38) [52,53]. Here the abbreviation DSE means Dyson-Schwinger equations with the label DB referring to the use of the most advanced Bethe-Salpeter kernel while RL denotes the rainbow ladder approximation.
The graphical representation of the pion DAs is displayed in Fig. 3 at the scale µ 2 = 2 GeV in terms of the DA projections on the plane (a 2 , a 4 ) using the symbols given in Table I. The experimental constraints are expressed in the form of 1σ (solid line) and 2σ (dashed line) error regions generated from the combined analysis of the CELLO [10], CLEO [11], Belle [4], and BABAR( 9 GeV 2 ) [3] data within LCSRs. The two slanted rectangles represent the constraints imposed by the QCD sum rules with nonlocal condensates used in [15] in connection with the determination of the Bakulev-Mikhailov-Stefanis (BMS) DAs. The larger one corresponds to the average vacuum quark virtuality λ 2 q (µ 2 ≈ 1 GeV 2 ) = 0.4 GeV 2 , whereas the smaller rectangle was determined in [16,60] using the slightly larger but still admissible value λ 2 q (µ 2 ≈ 1 GeV 2 ) = 0.45 GeV 2 (see [61] and references cited therein). It contains pion DAs with a characteristic platykurtic profile [16] (see Sec. III B).
This figure also contains the lattice constraints on a 2 at the scale µ 2 from [56] (red vertical lines further to the left) as well as those from [57] (blue vertical lines). The presented intervals in both cases are calculated by combining errors in quadrature. The results at the scale µ 2 are a 2 = 0.101 +0.024 −0.024 (NNLO) and a 2 = 0.078 +0.031 −0.029 (NLO) [56], whereas a 2 = 0.136 ± 0.021 [57]. A linear combination of errors would slightly overestimate the combined uncertainties yielding somewhat larger intervals of a 2 values. The results from [56], quoted in Table  I, were obtained from a combined chiral and continuum limit extrapolation at the NNLO and NLO level. This treatment differs from that applied in [57], where no extrapolation to the continuum limit was carried out. This  I: Conformal coefficients a2, a4, a6 for various pion DAs discussed in the text at two typical normalization momentum scales µ1 = 1 GeV and µ2 = 2 GeV. If µ2 is not the initial scale, NLO ERBL evolution in the global scheme is employed, see App. A. The range of the BMS and platykurtic DAs is related to the determination of a2 and a4 from QCD sum rules with nonlocal condensates using λ 2 q = 0.40 GeV 2 and λ 2 q = 0.45 GeV 2 , respectively. They cause the variation of the TFF predictions shown in the form of a green shaded band in Fig. 4. The coefficient a2 of the CZ DA was originally given at the scale µ = 0.5 GeV: a CZ 2 = 2/3 [54]. For the extrapolation to higher scales see [24]. Higher conformal coefficients up to and including a12 for the DSE-DB and DSE-RL DAs at the scale µ2 can be found in [53]. The coefficients up to and including a20 at the scale µ1 of the holographic AdS/QCD DA ϕ hol π (x) = (8/π) √ xx are tabulated in [55]. They were calculated here by means of the expression ξ 2n AdS/QCD [B(x, y) being the Euler Beta function] in combination with Eq. (24). The lattice results of [56] with NNLO (two loops) and NLO (one loop) matching to the MS scheme are quoted separately, where the subscript r denotes the systematic uncertainty due to the nonperturbative renormalization. They were obtained from a combined extrapolation to the chiral and continuum limit with associated uncertainties labeled by the subscripts m and a, respectively. The statistical errors of the data after extrapolation are given in sub-and super-scrip form. The question mark (?) in the lattice result of [57] indicates that it was not extrapolated to the continuum limit. - Table I by the question mark (?). The chiral extrapolation was included in the first parenthesis together with the statistical error, while the renormalization error is given in the second parenthesis. The general tendency of the new lattice estimates seems to favor DAs with a smaller value of a 2 . Recall that the second moment ξ 2 π (or equivalently a 2 ) gives information only on the variance statistic σ 2 [ϕ] = 1 4 ξ 2 π of the pion DA and contains no information about its shape in the central region. To this end, one needs the kurtosis statistic β 2 [ϕ] = ξ 4 π ( ξ 2 π ) 2 which measures the peakedness or flatness of a distribution in terms of the fourth moment, see [60] for a quantitative discussion.

B. Platykurtic pion DA
The particularity of the platykurtic pion DA [16] derives from the fact that it contextually and mathematically encapsulates in its profile the two consequences of confinement: (i) the appearance of nonlocal vacuum expectation values whose expansion in terms of local operators involves the virtuality λ 2 q of the vacuum quarks [43] and (ii) dynamical chiral symmetry breaking (DCHB) and mass generation [62]. The first feature entails quark correlations at a finite distance 1/λ q ∼ 0.3 fermi, while the second one entails the mass dressing of the confined quark propagator, see, for instance, [63]. As argued in [16,60], these effects induce distinctive geometrical characteristics of the pion DA. While the first one leads to the suppression of the endpoint regions x = 0, 1, the second one enhances the central region around x = 1/2.
The net result of this competition is a unimodal distribution with a unique short-tailed platykurtic profile. Take away the quark correlations and the endpoints of the DA get too strong resulting into a broad concave distribution in the whole x range similar to that obtained with DSE [53]. Leave aside the mass dressing and the central region is more or less depleted giving rise to a bimodal distribution whose bimodality strength is controlled by the nonlocality parameter λ 2 q . When λ 2 q = 0, one gets an infinite correlation length corresponding to the use of local condensates in the QCD sum rules. Such a situation gives rise to the Chernyak-Zhitnitsky (CZ) pion DA [54]. In contrast to bimodal DAs, like CZ and BMS, the platykurtic DA yields in the middle point ϕ (tw-2) π/pk (x = 1/2, µ 1 ) = 1.264 in agreement with the LCSR calculation [64] ϕ (tw-2) π (x = 1/2, µ 1 ) = 1.2 ± 0.3. For the derivation of the platykurtic DA and its range (Table I and small (green) strip in Fig. 3), we refer to [16,60].  (Table I) shown in comparison with the 1σ (solid line) and 2σ (dashed line) error regions created from the combined analysis of the CELLO [10], CLEO [11], Belle [4], and BABAR( 9 GeV 2 ) [3] data within LCSRs. The larger rectangle shows the range of values for BMS-like DAs from [15] and the smaller one the analogous region for the platykurtic DAs derived in [16,60]. The vertical lines mark the constraints on a2 from two lattice determinations: [56] NLO (dashed red lines) and NNLO (solid red lines); [57] (solid blue lines). The other designations are given in Table I with further explanations in the text.
From this figure we make the following striking observations: (i) The platykurtic strip shows a positive correlation between the coefficients a 2 and a 4 , while the BMS-type DAs (larger green rectangle) have coefficients with a negative correlation between them. (ii) Also the arrangement of the 1(2)σ error regions exhibits an anticorrelation pattern between a 2 and a 4 . (iii) Nevertheless, the platykurtic strip overlaps with the data regions at its upper right corner where the coefficients are given by a 2 (µ 2 ) ≈ 0.08, a 4 (µ 2 ) ≈ −0.009 (28) corresponding to the moments ξ 2 π ≈ 0.229, ξ 4 π ≈ 0.106 .
(iv) Remarkably, just there it also enters the range of the NNLO lattice constraints on a 2 from [56], while it mostly overlaps with the analogous NLO region. The observed agreement extends to the values of the second moment ξ 2 π . One has from [56] ξ 2 NNLO π (µ 2 ) = 0.234 +6 −6 (4) r (4) a (2) m , (30a) which gives after adding the errors in quadrature the values 0.234±0.0085 and 0.227±0.0095, respectively. These values conform with the platykurtic range for the second moment [60] ξ 2 pk π (µ 2 ) = 0.220 +0.009 −0.006 , ξ 4 pk π (µ 2 ) = 0.098 +0.008 −0.005 , (31) while the fourth moment is also given for the sake of comparison with other models. (v) All positively correlated DAs are unimodal but have enhanced tails, except the platykurtic one which shares tail suppression with the anticorrelated BMS-like DAs. (vi) Moreover, as we will see shortly, the platykurtic DA yields a TFF in good agreement with all data compatible with strict scaling at large Q 2 without crossing the pQCD asymptotic limit, cf. (33).

C. TFF predictions
We now outline the calculational procedure to obtain predictions for the scaled TFF Q 2 F γπ (Q 2 ) using the pion DAs given in Table I. This discussion relies upon Table  II in correspondence with the formalism exposed in Sec.
II. An extended discussion can be found in [67]. The results are shown in Fig. 4 in comparison with all existing data collected in App. B. The recently released preliminary data of the BESIII Collaboration [12,13] (see also [14]) are also included, keeping in mind that the probed momentum range extends below 1 GeV 2 , where our predictions are expected to be less reliable.
The calculated twist-two form factor is given in explicit form by with indications showing the sign of these contributions. The label (?) marks the only uncalculated NNLO term. To facilitate the use of Table II, we briefly describe the TFF calculation using as a reference model the set of the BMS DAs determined in [15]. These DAs are sufficiently well parameterized by means of the two lowest coefficients a 2 and a 4 , whereas higher coefficients a n>4 , (n = 6, 8, 10) can be ignored because they were found to be negligible albeit bearing large uncertainties [15]: a 6 ≈ a 2 /3; a 8 ≈ a 2 /4; a 10 ≈ a 2 /5. The variation of a 2 , a 4 allowed by the employed QCD sum rules gives rise to the narrower (green) strip of predictions in Fig. 4. It includes at the twist-two level the LO, NLO, and NNLO-T β0 contributions to the short-distance coefficients, cf. (8), (9), and involves the eigenfunctions {ψ 0 , ψ 2 , ψ 4 }. The ψ 0 eigenfunction yields the largest (negative) NNLO-T β0 contribution as we have seen in Fig. 2. Therefore, also the term T ∆V ≪ T β0 , cf. (10a), is taken into account only via the zero harmonic ψ 0 , whereas the term NNLO-T L vanishes for ψ 0 , cf. (10b) [7]. The remaining NNLO term T c is unknown and this unknownness induces the dominant theoretical uncertainty in the TFF prediction, shown in terms of the broader (blue) band enveloping the narrower (green) one. To gauge it, we assume that this term may be comparable in magnitude to the leading NNLO term T β0 (likely overestimating its significance), and obtain the uncertainties shown in the last column of Table III. Estimates of further theoretical errors-not considered here-can be found in [6,7].
The total TFF also comprises in the spectral density (18) the contributions (20) (twist four) and (22) (twist six) and includes NLO evolution with heavy-quark crossings, see App. A. The calculation of the TFF with the platykurtic DA is similar (black solid line in Fig. 4). The analogous computations with the other considered pion DAs include more conformal coefficients as indicated in Table II. The last row in this table provides information on the consistency of each TFF prediction with the asymptotic limit from pQCD [18,68], where we used the convenient notation The calculated TFF predictions are collected in Fig. 4 (left panel) and are shown in comparison with all data currently available: CELLO [10], CLEO [11], BABAR [3], Belle [4] (see Table III). The BESIII data are given in [12,13]. Apart from the results obtained with the BMS/pk DAs, already mentioned, this figure contains the following curves. The dashed red line farthest to the top shows F (Q 2 ) for the DSE-RL DA [52], whereas the solid red line below it gives the result for the DSE-DB DA [52] using {a 2 , a 4 , . . . , a 12 }. To exhibit the influence of the attenuation effect on the TFF, we also show the DSE-DB DA-based prediction using only a 2 and a 4 (longdashed-dotted-dotted red line). Obviously, this prediction agrees better with the asymptotic limit (horizontal black line). The reason is that it does not receive uninhibited contributions from the higher harmonics with a n>4 at higher Q 2 , see Fig. 2. The computation of the TFF for the holographic DA [55] yields the prediction represented by the dashed blue line running close to the asymptotic limit. Figure 4 also includes the form-factor predictions derived with the DA from the light-front quark model [58] (solid pink line) and the nonlocal chiral quark (NLχQM) (model 3 in Table II [59])-dashed-dotted red line. Both lines run below all data above 8 GeV 2 .

D. TFF Asymptotics
The behavior of the form factor is known theoretically in two limits. At Q 2 → 0 and in the chiral limit of quark masses, one obtains from the axial anomaly [70,71] lim where f π = 132 MeV is the leptonic decay constant of the pion. On the other hand, the asymptotic behavior of the form factor is given by Eq. (33) which is an exact  Table II and using various pion DAs defined in Table  I. The innermost (green) shaded strip shows the range of predictions obtained with the bimodal BMS DAs from [15]. The thick black line inside it denotes the result for the platykurtic DA [16]. expression from pQCD. The TFF in the Q 2 range between the aforementioned limits can be phenomenologically described by the interpolation formula of Brodsky and Lepage [68], In order to study the scaling behavior of the calculated TFF at large Q 2 more quantitatively, it is convenient to define the quantity which provides a normalized measure of the deviation of the scaled form factor from the asymptotic value-the "baseline". The graphical representation of the theoretical results for this quantity versus Q 2 is shown in comparison with the data in the right panel of Fig. 4. Inspection of the graphics in Fig. 4 leads to the following observations: (1) Most of the Belle data-except at Q 2 = 27.33 GeV 2come within errors close to the baseline clearly indicating that the TFF approaches an asymptotic value. Indeed, using the dipole formula F (Q 2 ) = BQ 2 /(C + Q 2 ) the Belle Collaboration determined B = 0.209 ± 0.016 GeV, which is slightly larger than the exact result √ 2f π , but still compatible.
(2) Several BABAR data points above 10 GeV 2 do not indicate TFF saturation as they move away from F ∞ .
Because also their error bars do not reach the F ∞ baseline, one may think that this data deviation is systematic and self-generated. It was shown in [6] that a dipole best-fit to the BABAR data yields B = 0.23 GeV and C = 2.6 GeV 2 with χ 2 = 1.7, being unable to reproduce the Belle data with acceptable accuracy. The origin of these discrepancies is the subject of several theoretical investigations (see [5] for a detailed discussion and references).
(3) The Q 2 intervals above 20 GeV 2 , probed by BABAR and Belle, are only scarcely populated and have a rather poor statistics. (4) The recently released preliminary BESIII data [12,13] cover the momentum range [0.3 − 3.1] GeV 2 and exceed the statistical accuracy of the CELLO data at Q 2 1 GeV 2 considerably, though their error bars become larger in the range 1.5 3.1 GeV 2 . These results are important for the hadronic light-by-light scattering calculations [36,37,39,72,73]. As regards the TFF below 1 GeV 2 , other dispersion approaches may be more adequate [9,[74][75][76] than the LCSR scheme applied in this work. (5) The broader (blue) band of predictions obtained with the BMS DAs (including their main uncertainties) approaches F ∞ gradually from lower values without reaching it (similarly also the platykurtic prediction). In the momentum interval [10 − 11] GeV 2 , the BMS TFF takes the value F BMS (Q 2 ) = 0.1604 −0.124 +0.0128 GeV, while the pk TFF reaches this value at [13.5 − 15.0] GeV 2 (see Table III). Both TFFs grow very slowly towards higher Q 2 values indicating that they have already entered the pre-asymptotic regime due to saturation. (6) The TFF predictions associated with the LFQMbased DA [58] (solid line) and the NLχ QM DA [59] (dashed-dotted red line), respectively, follow the trend of the platykurtic-generated prediction but have smaller magnitudes so that the onset of pre-asymptotic behavior in the TFF is shifted to much higher momenta. This is because both models have a negative a 6 coefficient (see Table I) that slows the TFF saturation. (7) The prediction based on the holographic AdS/QCD DA [55] provides a rather good agreement with the Belle data, while it disagrees with the BABAR data above 10 GeV 2 , where these start to grow. However, it crosses the baseline at still higher Q 2 around 57 GeV 2 . Variants of the holographic DA in [77] either yield similar results or tend to cross the baseline quite fast. (8) Still broader concave pion DAs, like DSE-DB and DSE-RL, lead to predictions that reach the preasymptotic regime in a way sensitive to the power α − in the "Gegenbauer-α" representation. This employs Gegenbauer polynomials of variable dimensionality α = α − + 1/2 [52,78] and gives TFF predictions with magnitudes growing in inverse proportion to α − . Employing the set {a 2 , a 4 , . . . , a 12 }, both TFF predictions cross the F ∞ line already at Q 2 ≈ 4 GeV 2 (DSE-RL) and Q 2 ≈ 10 GeV 2 (DSE-DB) and continue to grow. The reduced DSE-DB, which uses only a 2 , a 4 , leads to a TFF with a better asymptotic behavior (dashed-dotted-dotted red line) entering the pre-asymptotic regime around 10 GeV 2 and remaining then close to F ∞ but above it. (9) The theoretical predictions depend crucially on the DA models involved in the calculation. An attempt to extract the asymptotic behavior of the TFF directly from the data was given in [67].

IV. SUMMARY AND OUTLOOK
In this work we carried out a comprehensive analysis of the pion-photon transition form factor in QCD using the method of LCSRs within FOPT to NNLO and twist-six accuracy. The presented predictions for this exclusive observable are of considerable interest for two different reasons: (i) they provide a handle on the involved pion distribution amplitude and (ii) they represent a powerful tool to study the onset of scaling at high Q 2 in presentdays experiments.
To analyze in detail the Q 2 behavior of the TFF within the LCSR approach in comparison to pQCD, we studied the attenuation effect of the partial components of the TFF related to the conformal expansion of the pion DA and worked out how they contribute as Q 2 grows. We showed that broad, concave DAs relying on many positive conformal coefficients will tend to exceed the asymptotic limit F ∞ because higher components F n (Q 2 ) will start to contribute far beyond 10 GeV 2 .
A good agreement with F ∞ presumes saturation of the scaled TFF at large Q 2 and entails the onset of scaling. The big unknown is at which momentum scale this becomes obvious [79,80]. The over-all agreement of the TFF predictions, obtained in this work with the set of the BMS DAs [15] (including the platykurtic one [16]), with the CLEO [11], BABAR(< 9 GeV 2 ) [3], and Belle [4] data is good and any discrepancies are within the corresponding experimental errors, see Fig. 4. However, they disagree with the BABAR data above 10 GeV 2 . These predictions are shown as shaded bands and include the principal theoretical uncertainties of the conformal coefficients (narrower band) and the incomplete knowledge of the NNLO radiative corrections (wider band) in these figures. Saturation is observed in the interval [10 − 14] GeV 2 , where the TFF reaches the value 0.16 GeV, indicating the onset of the pre-asymptotic regime.
More ambitiously, the combination of the 1σ and 2σ error regions of all data compatible with asymptotic scaling with the most recent lattice constraints from [56] at the NNLO and NLO level supports a platykurtic pion DA with a 2 ≈ 0.08, a 4 ≈ −0.009 [16,60], though the sign and magnitude of a 4 require further consideration in the lattice context or otherwise. From the experimental side, it would be very helpful to have more data in regular steps of 1 GeV 2 above 10 GeV 2 . This would enable a reliable data-driven analysis based on the statespace reconstruction method, proposed in [67], and we hope that the Belle-II Collaboration will perform such measurements.
Acknowledgments I would like to thank Sergey Mikhailov and Alex Pimikov for a fruitful collaboration during the last years. I am indebted to Alex Pimikov for help with the numerical computations and their graphical representation. I am grateful to Gunnar Bali for useful remarks.
Appendix A: NLO evolution of the pion DA including heavy-quark thresholds In this Appendix (done partly in collaboration with S. V. Mikhailov and A. V. Pimikov) we discuss the NLO (i.e., the two-loop) ERBL evolution of the pion DA with an arbitrary number of Gegenbauer coefficients taking into account heavy-quark flavors (also known as global QCD scheme, see, e.g., [27] and references therein). This scheme employs the global coupling α glob s (Q 2 , Λ 2 N f ) that depends on the number of flavors N f through the QCD scale parameter Λ N f . This procedure was used in this work to derive the results given in Tables I, II, and III and obtain the predic-tions shown in all figures. It takes into account the heavyquark mass thresholds and thus requires the matching of the strong coupling in the Euclidean region of Q 2 at the corresponding heavy-quark masses when one goes from N f → N f + 1. Note that the dependence on N f in Appendix D of [24], which provided the basis for the NLO evolution of the pion DA in our earlier works, was ignored assuming a fixed number of flavors. The new scheme has already been used in our more recent investigations [5][6][7]81], but without exposing the underlying formalism in final form. This task will be accomplished here including further refinements. The NLO evolution of the pion DAs with two conformal coefficients a 2 and a 4 at the initial scale µ 2 ≃ 1 GeV 2 and a varying number of heavy flavors has also been applied in [82] (see Appendix D there). Our technical exposition below extends this treatment to any number of conformal coefficients and more heavy-flavor thresholds. For some specific details and references, we refer to [27,83].
Let us start with a fixed number of flavors and supply some basic formulas from [24]. The ERBL evolution equation for the pion DA is given by and is driven by the kernel with a s = α s /(4π).
To perform the pion DA evolution, while ignoring quark-mass thresholds, we make use of the evolution matrix E with the components E nk . Expanded over the basis {ψ n } of the Gegenbauer harmonics, this matrix assumes the following triangular form [84] E nk (N f ; Q 2 , µ 2 ) = P (n, where the coefficients d nk (Q 2 , µ 2 ) will be defined shortly, and where µ 2 and Q 2 refer to the initial and observation scale, respectively. The factor P (n, Q 2 , µ 2 ) in Eq. (A4) denotes the diagonal part of the evolution matrix that dominates the renormalization-group (RG) controlled evolution of the ψ n -harmonics in the conformal expansion ϕ RG π (x, Q 2 ) = n a n (µ 2 ) P (n, Q 2 , µ 2 ) ψ n (x) + a s (Q 2 ) Then, the diagonal part of the evolution exponential at the two-loop level can be given explicitly, where a s (µ 2 ) = α glob;(2) s (µ 2 , Λ 2 3 )/(4π) and c 1 = b 1 /b 0 , with b i being the expansion coefficients of the QCD βfunction. The evolution exponent of the coupling is de- The second term in the brackets in Eq. (A4) represents the nondiagonal part of the evolution equation to the order O(a 2 s ) induced by renormalization and encodes the mixing of the higher Gegenbauer harmonics for indices k > n related to the conformal-symmetry breaking at NLO [23]. Notice that all components on the right-hand side of Eqs. (A4) and (A7) depend on N f , which changes to N f + 1, when the next quark-mass threshold is crossed. The explicit form of the mixing coefficients is given by [24] where the values of the first few elements of the matrix M nk (k = 2, 4 ≥ n = 0, 2) read Analytic expressions for M nk have been obtained in [22]. The values in Eq. (A9) reproduce the exact results with a deviation less than about 1%. To make our further exposition more compact, we make use of the parameter vectors A(µ 2 ) and Ψ(x) defined at the reference momentum scale µ 2 as follows A = (1, a 2 , a 4 , · · · , a 2(N −1) ) , (A10a) Ψ = (ψ 0 , ψ 2 , · · · , ψ 2(N −1) ) , where their dimension and the dimension of the matrix E depends on the parameter N . Then, the evolution of the pion DA can be carried out in terms of the Gegenbauer coefficients a i with i = 2, 4, . . . , 2(N − 1). For a fixed number of flavors, one gets where E T is the transposed matrix of E, while A(µ 2 0 ) is the vector of the Gegenbauer coefficients defined at some initial scale µ 2 0 . In the global QCD scheme, the evolution of the pion DA defined at the initial scale µ 2 0 , is implemented by means of the threshold interval factors E i in the following step-by-step procedure, where the matrices E i and E i (µ 2 ) are given by and the thresholds are defined [27] by the heavy-quark masses m c ∼ M 4 = 1.65 GeV, m b ∼ M 5 = 4.75 GeV, and m t ∼ M 6 = 172.5 GeV, while M 2 3 ≡ µ 2 0 sets the initial scale taken to be either µ 0 = µ 1 = 1 GeV or µ 0 = µ 2 = 2 GeV, see Table I. Note that the global evolution matrix, Eq. (A12), is presented for µ 0 < M 4 and µ > µ 0 . No matching at the mass thresholds is needed in the case of equal initial and final momentum scales, i.e., E(N f ; Q 2 , Q 2 ) = 1 because of the independence of the evolution matrix on the number of flavors. For example, at the threshold M 4 , we have E 4 (M 2 4 ) = 1 ensuring the continuity of the global evolution matrix E glob . It is worth noting that our NLO evolution scheme in terms of Eq. (A12), has the following improvements relative to that used in [82] (see Appendix D there): (a) It is applicable to DAs with any number of Gegenbauer harmonics. (b) The number of heavy-quark thresholds is extended to four flavors. (c) When the interval of evolution contains two or more mass thresholds, our method can still incorporate contributions from the non-diagonal part of the evolution matrix, removing the restriction of using only the first two Gegenbauer coefficients a 2 and a 4 as in [82].
We reiterate that the matching of the coupling constants at the quark-mass thresholds requires the readjust-ment of the value of the QCD scale parameter Λ to Λ (N f ) . A detailed description of the matching procedure of the running coupling in the global scheme can be found, for instance, in [27]. For definiteness, we quote here the twoloop Λ (2) (N f ) values used in our code: Λ whereas the global evolution of the pion DA assumes the form ϕ glob π (x, µ 2 ) = A glob (µ 2 )Ψ(x) = A(µ 2 0 )Ψ(x; µ 2 ) = N −1 n=0 a glob 2n (µ 2 )ψ 2n (x) .