Subleading power corrections to the pion-photon transition form factor in QCD

We reconsider QCD factorization for the leading power contribution to the $\gamma^{\ast} \gamma \to \pi^0$ form factor $F_{\gamma^{\ast} \gamma \to \pi^0} (Q^2)$ at one loop using the evanescent operator approach, and demonstrate the equivalence of the resulting factorization formulae derived with distinct prescriptions of $\gamma_5$ in dimensional regularization. Applying the light-cone QCD sum rules (LCSRs) with photon distribution amplitudes (DAs) we further compute the subleading power contribution to the pion-photon form factor induced by the"hadronic"component of the real photon at the next-to-leading-order in ${\cal O}(\alpha_s)$, with both naive dimensional regularization and 't Hooft-Veltman schemes of $\gamma_5$. Confronting our theoretical predictions of $F_{\gamma^{\ast} \gamma \to \pi^0} (Q^2)$ with the experimental measurements from the BaBar and the Belle Collaborations implies that a reasonable agreement can be achieved without introducing an"exotic"end-point behaviour for the twist-2 pion DA.


Introduction
Hard exclusive processes play a prominent role in exploring the strong interaction dynamics of hadronic reactions in the framework of QCD. The pion-photon transition form factor γ * γ → π 0 at large momentum transfers (Q 2 ) serves as one of the simplest exclusive processes for testing the theoretical predictions based upon perturbative QCD factorization. The hard-collinear factorization theorem for the pion-photon form factor F γ * γ→π 0 (Q 2 ) can be demonstrated at leading power in 1/Q 2 utilizing both diagrammatic approaches [1][2][3] and effective field theory techniques [4]. The hard coefficient function entering the leading-twist factorization formula has been computed at one loop [5][6][7], and at two loops [8] in the large β 0 approximation. In virtue of the fact that the twist-2 pion distribution amplitude (DA) is defined by an axialvector light-ray operator, a subtle issue in evaluating QCD corrections to the hard function in dimensional regularization lies in the definition the Dirac matrix γ 5 in the complex Ddimensional space demanding a new set of algebraic identities and various prescriptions for the treatment of γ 5 have been proposed to meet the demand of precision QCD calculations in different contexts (see [9,10] for an overview and [11][12][13][14][15] for more discussions). Employing the trace technique, the γ 5 ambiguity of dimensional regularization was resolved by adjusting the way of manipulating γ 5 in each diagram to preserve the axial-vector Ward identity [6], which is less straightforward (systematic) for the higher-order QCD calculations of hadronic reactions. One of our major objectives of this paper is to demonstrate the equivalence of factorization formulae for the pion-photon transition form factor constructed with naive dimensional regularization (NDR) and 't Hooft-Veltman (HV) schemes of the γ 5 matrix, using the spinor decomposition technique [16][17][18] and the evanescent operator approach [19,20].
Confronting the theoretical predictions with the precision experimental measurements of the π 0 γ * γ form factor at accessible Q 2 evidently necessities a better understanding of the subleading power terms in the large momentum expansion, due in particular to the scaling violation implied by the BaBar data [21]. The significance of the power suppressed contributions to F γ * γ→π 0 (Q 2 ) was highlighted by evaluating the soft correction to the leading twist effect with the dispersion approach [22,23] and turned out to be crucial to suppress the contributions from higher Gegenbauer moments of the twist-2 pion DA (see also [24,25]). An attractive advantage of the dispersion approach [26] is that the subleading power "hadronic" photon correction is taken into account effectively by modifying the spectral function in the real-photon channel at the price of introducing two nonperturbative parameters (i.e., the vector meson mass m ρ and the effective threshold parameter s 0 ). This effective method allows continuous improvement of the theoretical accuracy for predicting the pion-photon form factor by including the next-to-next-to-leading order (NNLO) QCD correction to the twist-2 contribution and the finite-width effect of the unstable vector mesons in the hadronic dispersion relation [27][28][29][30]. Further applications of this technique were pursued in radiative leptonic B-meson decay [31,32] and electro-production of the pseudoscalar eta mesons [33] and of tensor mesons [34] in an attempt to "overcome" the difficulty of rapidity divergences emerged in the direct QCD calculations of the subleading power contributions. It is then in demand to provide an independent QCD approach to compute the above-mentioned power corrections for the sake of boosting our confidence on the reliability of both theoretical tools. Another objective of this paper is to construct the light-cone sum rules (LCSRs) for the hadronic photon effect

Factorization of the leading power contribution
The purpose of this section is to compute the leading power contribution to the pion-photon form factor at one loop π(p)|j em µ |γ(p ) = g 2 em µναβ q α p β ν (p )F γ * γ→π 0 (Q 2 ) , with both the NDR and HV schemes for the γ 5 matrix in D dimensions, where q = p − p , p refers to the four-momentum of the pion, the on-shell photon carries the four-momentum p and j em µ = q g em Q qq γ µ q , 0123 = −1 .
We further introduce a light-cone vectorn µ parallel to the photon momentum p , define another light-cone vector n µ along the direction of the momentum p in the massless pion limit, and employ the following power counting scheme at large momentum transfer n · p ∼ n · p ∼ O( Q 2 ) , n · p ∼ O(Λ 2 / Q 2 ) .
2.1 QCD factorization of F γ * γ→π 0 (Q 2 ) at tree level QCD factorization for the leading-twist contribution to the γ * γ → π 0 form factor at tree level can be established by inspecting the four-point QCD matrix element at leading order (LO) in α s , where x indicates the momentum fraction carried by the collinear quark of the pion andx ≡ 1 − x. Computing the two diagrams displayed in figure 1 yields Here, O A, µν (0) and O B, µν (0) denote the tree-level partonic matrix elements of the collinear operators O A, µν and O B, µν in soft-collinear effective theory (SCET) and the convolution integration is represented by an asterisk. The manifest definition of the SCET operator O j, µν in the momentum space are given by with the collinear Wilson line and To facilitate the determination of the hard function entering the leading power factorization formula of F γ * γ→π 0 (Q 2 ), we employ the SCET operator basis where O E, µν is an evanescent operator vanishing in four dimensions and It is evident that the effective operator O 1, µν cannot couple with a collinear pion state due to the parity conservation. Taking advantages of the operator identities we observe that the two tree-level diagrams in figure 1 give rise to the identical contribution to the pion-photon transition form factor and such observation can be further generalized to all orders in QCD applying the charge-conjugation transformation. We now employ the operator matching equation with the evanescent operator Figure 2: Diagrammatical representation of the one-loop contribution to the partonic amplitude γ γ * → qq induced by two electromagnetic currents. The corresponding symmetric diagrams obtained by exchanging the two photon states are not shown. and expand all quantities to the tree level, yielding Utilizing the definition of the leading twist pion DA on the light cone it is straightforward to derive the tree-level factorization formula of the π 0 γ * γ form factor

QCD factorization of
We proceed to compute the NLO QCD correction to the four-point partonic amplitude Π (1) µ at leading power in 1/Q 2 for the determination of the hard function T 2 at O(α s ). It needs to be stressed that the QCD matrix element Π µ defined by two electromagnetic currents is independent of the prescription of γ 5 in the D-dimensional space and the renormalization scheme dependence of the perturbative matching coefficient T 2 comes solely from the radiative correction to the twist-2 pion DA φ π (x, µ), whose definition depends on the precise treatment of the Dirac matrix γ 5 in dimensional regularization, due to the infrared subtraction. Evaluating the hard contribution to the one-loop diagrams displayed in figure 2 with the method of regions [55] immediately leads to where the ellipses represent terms proportional to O 1, µν (x, x ) (0) and O E, µν (x, x ) (0) . Adding up different pieces together we can readily obtain the QCD matrix element Π µ at O(α s ) where the γ 5 -prescription independent amplitude A 2,hard reads Expanding the matching equation (13) to the one-loop order yields Now we are in a position to derive the master formula for the one-loop perturbative matching coefficient T (1) i by implementing both the ultraviolet (UV) renormalization and the infrared (IR) subtraction. Following the strategy presented in [18] the UV renormalized matrix element of the SCET operator O i, µν at O(α s ) is given by where the bare matrix element M (1) ij,bare depends on the IR regularization scheme R. Applying the dimensional regularization for both the UV and IR divergences, the bare matrix element M (1) ij,bare vanishes due to scaleless integrals entering the relevant one-loop computation. Inserting (24) into (23) and comparing the coefficient of O 2, µν (0) give rise to It is evident that the SCET operators O 1, µν and O 2, µν cannot mix into each other under QCD renormalization due to the parity conservation, hence Z 12 = 0. In addition, the IR subtraction term T (0) 2 * Z (1) 22 will remove the collinear contribution to the QCD amplitude Π µ at one loop so that the matching coefficient T only encodes the information of strong interaction dynamics at the hard scale. Technically, the collinear subtraction has been automatically implemented in the above computation of the QCD matrix element Π µ , since only the hard contribution computed with the expansion by regions enters the expression of A (1) 2,hard displayed in (22). We are now ready to discuss the renormalization constant Z E2 of the evanescent operator O E, µν for the derivation of the final result of the matching coefficient T (1) 2 . Applying the renormalization prescription that the IR finite matrix element of the evanescent operator O E, µν vanishes with dimensional regularization applied only to the UV divergences and with the IR singularities regularized by any parameter other than the dimensions of spacetime [19,20] and making use of the identity (24) yield The one-loop matching coefficient of the physical operator O 2, µν can be readily obtained by substituting (26) into (25) The one-loop contribution to the matrix element of the evanescent operator O E, µν depends on the renormalization prescription of γ 5 in the D-dimensional space. We will employ both the NDR and HV schemes of γ 5 below for the illustration of the prescription independence of the factorization formula of F γ * γ→π 0 (Q 2 ) at O(α s ) and at leading power in 1/Q 2 . This amounts to showing that the renormalization scheme dependence of the short-distance coefficient function cancels against that of the twist-2 pion DA precisely. Evaluating the SCET diagrams displayed in figure 3 with the Wilson-line Feynman rules, we find that only a single diagram 3(a) can generate a nonvanishing contribution to M (1)off E2 using the NDR scheme of γ 5 , which turns out to be proportional to the spin-dependent term of the Brodsky-Lepage evolution kernel [1,2]. The manifest expression of the infrared subtraction term T (0) is then given by However, computing the one-loop matrix element of the evanescent operator O E, µν with the HV scheme of γ 5 leads to Inserting (22), (28) and (29) into the master formula (27), we obtain where the renormalization scheme dependent parameter δ is given by Our result of T in the NDR scheme is identical to that presented in [6] using the trace formalism.
We now aim at demonstrating the renormalization prescription independence of the oneloop factorization formula for the pion-photon form factor where the superscript "∆" indicates the γ 5 -scheme in dimensional regularization. Taking advantage of the relation of the twist-2 pion DA between the NDR and HV schemes where the finite renormalization kernel Z −1 HV is given by [56] It is then straightforward to show that which cancels against the renormalization scheme dependence of the NLO hard function T (1), ∆ 2 (x, µ) as displayed in (30). We emphasize again that the γ 5 -prescription independence of the leading power factorization formula for F γ * γ→π 0 (Q 2 ) stems from the fact that the QCD matrix element q(x p)q(x p)|j em µ |γ(p ) itself is free of the γ 5 ambiguity in dimensional regularization and both the NDR and HV prescriptions can be employed to construct QCD factorization theorems for hard processes provided that the corresponding matching coefficients are computed in an appropriate way without overlooking the potential evanescent operators.
The renormalization scale independence of the factorization formula (32) can be readily verified by making use of the evolution equation of the pion DA φ π (x, µ) where the evolution kernel V (x, y) can be expanded perturbatively in QCD with the "plus" function defined as and the LO Brodsky-Lepage kernel given by [1,2] V 0 (x, y) It is appropriate to point out that the one-loop evolution kernel V 0 (x, y) is independent of the γ 5 prescription in the complex D-dimensional space, however, the two-loop evolution kernel V 1 (x, y) does depend on the renormalization scheme. Applying the renormalization-group (RG) equation (36) then leads to We further turn to sum the parametrically large logarithms of Q 2 /µ 2 in the short-distance function at next-to-leading-logarithmic (NLL) accuracy employing the standard RG approach in the momentum space. Technically, the desired NLL resummation can be readily achieved by setting the factorization scale of order µ ∼ Q 2 and evolving the leading twist pion DA up to that scale at two loops. The NLO evolution kernel V 1 (x, y) was first obtained with the diagrammatic approach in the light-cone gauge [57,58], then in the Feynman gauge [59,60], and was further reconstructed [61] based upon the knowledge of the conformal anomalies and the available forward DGLAP splitting functions at O(α 2 s ). The two-loop evolution potential V 1 (x, y) can be organized as where N f is the number of the active quark flavours. The explicit expressions of the kernels V N , V G and V F are given by [62] where we have introduced the functions F (x, y) and H(x, y) as follows In order to perform the NLL QCD resummation, it turns out to be convenient to adopt the Gegenbauer expansion of the pion DA where the LO coefficient a 0 (µ) = 1 is renormalization invariant due to the normalization condition. The exact solution to the two-loop RG equation (36) can be constructed from the forward anomalous dimensions and the special conformal anomaly matrix in the Gegenbauer moment space [63,64], and we obtain (see also [22]) where both n and k are non-negative even integers and the explicit expressions of E NLO V,n and d k V,n are collected in Appendix A. Inserting (47) into (32) and employing the technique developed in [65], we obtain where the hard coefficient C n (Q 2 , µ) in the NDR scheme of γ 5 is given by with the harmonic number defined as H n = Σ k=n k=1 1/k.

The subleading-power correction from photon DAs
In this section we aim at evaluating the power suppressed contribution to the pion-photon form factor due to the hadronic photon effect at O(α s ) with the LCSR approach. To this end, we construct the following vacuum-to-photon correlation function defined with an electromagnetic current (2) carrying a four-momentum q µ and a pion interpolating current j π whose explicit structure is as follows Here we have introduced the convention ⊥ µναβ ≡ g ρ ⊥ µ ρναβ . Following the standard strategy, the primary task for the sum-rule construction consists in the demonstration of QCD factorization for the considered correlation function (51) at space-like interpolating momentum p = p +q. In contrast to the factorization proof of the leading power contribution presented in Section 2, the QCD matrix element (51) itself depends on the γ 5 prescription in D-dimensional space manifestly. We will employ both the NDR and HV schemes of the Dirac matrix γ 5 to establish the QCD factorization formula of the transition amplitude (51) at O(α s ) and then derive the NLL resummation improved LCSR for the hadronic photon correction to the π 0 γ * γ form factor. Furthermore, the power counting rule for the external momenta will be adopted to determine the perturbative matching coefficient entering the factorization formula of G µ (p , q) to the one-loop order.

The hadronic photon effect at tree level
QCD factorization for the correlation function (51) at tree level can be established by investigating the following four-point QCD amplitude at LO in α s . Evaluating the two diagrams in figure 4 leads to where r = −p 2 /Q 2 , η u = 1 and η d = −1. The partonic matrix element of the (anti)-collinear SCET operator O A,µ at tree level is given by The explicit definition of the (anti)-collinear operator O j,µ in the momentum space is where we have suppressed the flavour indices of O j,µ for brevity, Γ A, µ = γ 5 n γ µ,⊥ and the corresponding Wilson line is defined as To achieve the hard-collinear factorization for the correlation function (51), we introduce the SCET operator basis where O E,µ is evidently an evanescent operator. Applying the operator matching equation including the evanescent operator and expanding all quantities to the tree level, we can readily find that Making use of the leading twist DA of the photon defined in [35] 0|χ(0) Wc(0, y) σ αβ χ(y)|γ(p ) the tree-level factorization formula of the form factor G(p 2 , Q 2 ) can be written as where the magnetic susceptibility of the quark condensate χ(µ) encodes the dynamical information of the QCD vacuum [66].
Applying the standard definition for the pion decay constant we can write down the hadronic dispersion relation of G(p 2 , Q 2 ) The final expression of the LCSR for the hadronic photon correction to the pion-photon form factor can then be derived by implementing the continuum subtraction and the Borel transformation with the aid of the parton-hadron duality with u 0 = Q 2 /(s 0 + Q 2 ). Employing the power counting scheme for the sum rule parameters we can readily obtain the scaling behaviour of the hadronic photon effect at large Q 2

The hadronic photon effect at one loop
To construct the NLL LCSR for the hadronic photon effect, we first need to establish the one-loop factorization formula for the correlation function (51) at the leading power in 1/Q 2 . Following the strategy for demonstrating QCD factorization of the leading power contribution presented in Section 2, the perturbative matching coefficient entering the factorization formula of the form factor G(p 2 , Q 2 ) can be determined by evaluating the one-loop diagrams for the QCD amplitude Π µ (54) in figure 5. We will compute the hard contributions from these diagrams one-by-one with both the NDR and HV schemes applying the strategy of regions. The one-loop QCD correction to the electromagnetic vertex diagram displayed in figure  5(a) is obviously free of the γ 5 ambiguity and a straightforward calculation yields where the term proportional to O E,µ (x, x ) is not shown explicitly. Due to the appearance of γ 5 in the pion interpolating current, the one-loop QCD correction to the pion vertex diagram depends on the γ 5 prescription employed in the reduction of the Dirac algebra. Computing the hard contribution from the one-loop diagram 5(b) in both the NDR and HV schemes gives where Π (0a) µ represents the tree-level contribution to the diagram 4(a) and can be obtained from (55) We finally turn to compute the hard contribution from the one-loop box diagram shown in figure 5(d), which depends on the actual prescription of γ 5 adopted in the calculation of the corresponding QCD amplitude. Evaluating the contribution from the box diagram with both the NDR and HV schemes, we find that the corresponding hard coefficients only contribute at O( ), vanishing in four-dimensional space. Explicitly, Collecting different pieces together, it is straightforward to derive the NLO QCD correction to the four-point amplitude Π µ where the renormalization prescription dependent hard amplitude A 1,hard is given by Applying the strategy to implement the IR subtraction for the four-point QCD amplitude Π µ discussed in Section 2, the master formula for the one-loop hard coefficient of the physical SCET operator O 1,µ can be written as where the bare matrix element M represents the QCD mixing of the evanescent operator O E,µ into O 1,µ at one loop. It is evident that the infrared subtraction term T (0) suffers from the γ 5 ambiguity in dimensional regularization. The corresponding SCET diagrams at one loop are in analogy to that displayed in figure 3, but with the vertex "⊗" indicating an insertion of O E,µ . Computing these effective diagrams with dimensional regularization applied to the UV divergences and with the IR singularities regularized by the fictitious gluon mass, we find that M (1)off E1 vanishes at one loop with the NDR scheme of γ 5 and it receives a nonvanishing contribution of O( ) with the HV scheme of γ 5 from the effective diagram with a collinear gluon exchange between two external quarks. We are then led to conclude that Inserting (78) into (77) immediately yields for both the NDR and HV schemes of the γ 5 matrix, with A 1,hard presented in (75) and (76). We mention in passing that the γ 5 scheme dependence of the short-distance function T (1) 1 will not be cancelled by the one-loop QCD correction to the twist-2 photon DA defined by the light-cone matrix element of the tensor current, which is clearly free of the γ 5 ambiguity in dimensional regularization, and the γ 5 ambiguity of A 1,hard can be traced back to the renormalization prescription dependence of the QCD amplitude (54) itself.
To preserve the one-loop character of the axial anomaly, an additional finite counterterm must be introduced [11] when performing the UV renormalization of the pseudoscalar current in the HV scheme. Making use of (61), (75), (76) and (79), it is then straightforward to verify that at one loop, which provides a nontrivial check to justify the obtained one-loop hard amplitude T 1 . The NLO factorization formula for the vacuum-to-photon correlation function can be further derived as follows With the NLO hard coefficient function T (1) 1 at hand, we can also obtain the one-loop shortdistance function entering the factorization formula of the H → J/ψ γ form factor at leading power in 1/m 2 H by taking the r → ∞ limit of T 1 and by performing the analytical continuation in the variable p 2 , which reproduces the expression displayed in (3.17) of [67] (see also [68,69]) computed from an alternative approach precisely.
We are now in a position to demonstrate the factorization-scale independence of (82) by employing the RG equation of the leading twist photon DA with the perturbative expansion of the evolution kernel The one-loop renormalization kernel V 0 (x, y) is given by [70] V 0 (x, y) Taking into account the factorization scale dependence of T 1 (x, µ), we can further deduce The residual µ-dependence of G(p 2 , Q 2 ) evidently originates from the UV renormalization of the QCD pseudoscalar current defining the correlation function (51). Taking advantage of the evolution equation of the QCD renormalization constant for the pseudoscalar current [71,72] and distinguishing the renormalization scale of the QCD current from the factorization scale due to the IR subtraction (see [73] for more details), we can then find that the expression for the form factor G(p 2 , Q 2 ) (82) is indeed factorization-scale invariant at O(α s ). We proceed to perform the NLL resummation for the parametrically large logarithms in the short-distance function T 1 , which can be achieved alternatively by fixing the factorization scale as µ ∼ Q 2 and by evolving the twist-2 photon DA from the hadronic scale to that scale. To this end, we need the two-loop coefficient of the evolution kernel V (x, y) [61,74,75] where the explicit expressions of the kernel functions are [74] V with Applying the Gegenbauer expansion of the twist-2 photon DA [35] φ and implementing the conformal consistency relation discussed in [63], the two-loop evolution of the Gegenbauer moment b n (µ 0 ) can be constructed as follows with even k, n ≥ 0. The detailed expressions the evolution functions E NLO T,n and d k T,n can be found in Appendix A. Combining everything together we arrive at the NLL resummation improved factorization formula where the perturbative matching coefficient C n (Q 2 , µ) is defined by 1 (x, µ) We will not present the analytical result of C n (Q 2 , µ) by evaluating the appeared convolution integral explicitly, since the continuum subtraction needs to be performed for the dispersion representation of (97) in order to construct the desired LCSRs for the hadronic photon correction to the pion-photon form factor. Employing the spectral representations of the convolution integrals displayed in Appendix B, it is straightforward to derive the dispersion form of the NLL factorization formula where we have exploited the symmetric property of the photon DA φ γ (x, µ) = φ γ (x, µ) due to the charge-parity conservation. The resulting QCD spectral densities ρ (i) (s, Q 2 ) (i = 0, 1) can be written as where P indicates the principle-value prescription. The NLL LCSRs for the subleading power contribution to the π 0 γ * γ form factor can be further derived as Collecting different contributions together, we now present the final expression for the pion-photon form factor including the twist-4 correction computed in [22,26] where the manifest expressions of F LP γ * γ→π 0 and F NLP γ * γ→π 0 are displayed in (49) and (101), respectively. The obtained factorization formula of F tw−4 γ * γ→π 0 from both the two-particle and the three-particle pion DAs at tree level reads where the definition of the twist-4 pion DA F π can be found in (38) of [22] and keeping only the leading conformal spin (i.e., "S"-wave) contribution we obtain The nonperturbative parameter δ 2 π is defined by the local QCD matrix element 0|g sqGµν γ ν q|π(p) = i f π δ 2 π (µ) p µ , with the renormalization-scale evolution at one loop Several comments on the general structure of the π 0 γ * γ form factor (102) are in order.
• It is apparent that the twist-four correction to the π 0 γ * γ form factor is suppressed by a factor of δ 2 π /Q 2 compared with that of the leading twist contribution. Such subleading power contribution turns out to be numerically significant at Q 2 ≤ 5 GeV 2 due to the large prefactor "80/3" entering the asymptotic expression of F π (x, µ), however, it is still far from sufficient to generate the scaling violation at Q 2 ∼ 40 GeV 2 indicated by the BaBar measurement [21]. Furthermore, it is of high interest to compute the NLO correction to the twist-four contribution in order to develop a better understanding of factorization properties of the high twist effects, where the infrared subtraction for constructing the factorization formula is complicated by the mixing of different twist-four pion DAs under the QCD renormalization.
• The twist-six correction to the pion-photon form factor computed from the dispersion approach [22] is partially absorbed into the hadronic photon effect F NLP γ * γ→π 0 (Q 2 ) displayed in (101). The precise correspondence of distinct contributions in two frameworks cannot be established without identifying the operator definitions of the "soft" corrections in [22], which originate from the nonperturbative modification of the QCD spectral density appeared in the dispersion form of the π 0 γ * γ * form factor.
• The subleading power corrections from the yet higher twist pion/photon DAs, which are not taken into account in this work, are conjectured to be suppressed by only one power of Λ 2 /Q 2 due to the absent correspondence between the twist counting and the largemomentum expansion [22]. A manifest calculation of the two-particle and three-particle corrections to the pion-photon form factor from the twist-three and twist-four photon DAs based upon the LCSR approach is in demand to verify this interesting hypothesis.
• We do not include the NNLO QCD correction to the leading power contribution in the large β 0 approximation [8] on account of the absence of a complete NNLO contribution, which also necessitates the three-loop evolution equation of the twist-two pion DA [76] to obtain the factorization formula at the next-to-next-to-leading-logarithmic accuracy. A recent discussion of the NNLO radiative corrections in the framework of the dispersion approach can be found in [30].

Numerical analysis
We are now ready to explore the phenomenological consequences of the hadronic photon correction to the pion-photon form factor applying the master formula (102). In doing so, we will first need to specify the non-perturbative models for the twist-2 pion DA, the magnetic susceptibility χ(µ), the Gegenbauer moments of the photon DA, and to determine the "internal" sum rule parameters entering the expression (101).

Theory input parameters
The fundamental ingredients entering the NLL factorization formula of the leading power contribution are the Gegenbauer moments of the twist-2 pion DA. Tremendous efforts have been devoted to the determinations of the lowest moment a 2 (µ) from the direct calculations with the QCD sum rules pioneered by Chernyak and Zhitnitsky (CZ) [77] and with the lattice simulations, and from the indirect calculations by matching the LCSR predictions with the experimental data. To quantify the systematic uncertainty from the Gegenbauer moments, we will consider the following four models for the leading twist pion DA The obtained Gegenbauer coefficients in the Bakulev-Mikhailov-Stefanis (BMS) model [30,78] are computed from the QCD sum rules with non-local condensates absorbing the high-order terms in the operator-product-expansion (OPE) partially (see, however, [79]). The first and second nontrivial Gegenbauer moments of the KMOW model [80] are determined by comparing the LCSR predictions for the pion electromagnetic form factor, including the NLO correction to the twist-2 effect and the subleading terms up to twist-6, with the intermediate-Q 2 data from the JLab experiment. The holographic model of the twist-2 pion DA [81] φ Hol π (x, µ 0 ) = is motivated by the correspondence between the string theory in the five-dimensional antide Sitter space and conformal field theories in the physical space-time (see also [82] for a similar end-point behaviour of the pion DA) and implementing the Gegenbauer expansion of φ Hol π (x, µ 0 ) leads to the expression of a n displayed in (107). For the phenomenological analysis of the π 0 γ * γ form factor, we will truncate the expansion of the "holographic" model at n = 12, which was demonstrated to be a good approximation in [22]. It needs to point our that the values of the second Gegenbauer coefficient in the first three models of (107) are in line with the recent lattice determinations [83] within the theory uncertainties and the CZ model is introduced for the illustration purpose to understand the model dependence of the predictions for the pion-photon form factor.
The normalization parameter for the twist-four pion DAs will be taken as δ 2 π (1 GeV) = (0.2 ± 0.04) GeV 2 computed from the QCD sum rules [84] (see also [85]). We further adopt the value of the quark condensate density qq (1 GeV) = − 256 +14 −16 MeV 3 determined in [80]. A key nonperturbative quantity appearing in the twist-2 photon DA is the magnetic susceptibility of the quark condensate χ(µ) describing a response of the QCD vacuum in the presence of an external photon field. Different QCD-based approaches have been proposed to evaluate χ(µ) (see, e.g., [35,86,87]) with the aid of the resonance information from the experimental data and the interval χ(1 GeV) = (3.15 ± 0.3) GeV −2 [35] will be employed in the numerical calculations. In contrast, our understanding of the higher Gegenbauer moments of the leading twist photon DA is rather limited, even for the leading non-asymptotic correction due to b 2 (µ 0 ). The available information of the second Gegenbauer coefficient mainly comes from the QCD sum rules constructed from the correlation function with a light-ray tensor operator and a local vector current, which unfortunately give rise to the theory predictions sensitive to the choice of the input parameters. The crude estimate b 2 (1 GeV) = 0.07 ± 0.07 from [85] will be used in our numerical analysis and an independent determination from the lattice QCD calculation will be very welcome in the future.
A natural choice of the factorization scale in the leading power factorization formula (49) is µ 2 = x Q 2 with 1/4 ≤ x ≤ 3/4 corresponding to the characteristic virtuality of the intermediate quark displayed in figure 1(a), and it will be frozen at µ = 1 GeV for x Q 2 < 1 GeV 2 at low Q 2 in order not to run into the nonperturbative QCD regime (see [8] for the discussion about the BLM proposal). Along similar lines, the factorization scale entering the NLL LCSRs for the hadronic photon effect (101) will be taken as µ 2 = x M 2 + x Q 2 as widely employed in the sum rule calculations [22]. Finally, the determination of the Borel mass M 2 and the threshold parameter s 0 can be achieved by applying the standard strategies described in [88,89], and we can readily obtain in agreement with the intervals adopted in [90].
4.2 Predictions for the π 0 γ * γ form factor Now we will turn to investigate the phenomenological significance of distinct terms contributing to the pion-photon form factor. Taking the BMS model for the twist-two pion DA as an example, it is evident from figure 6 that the twist-four correction and the hadronic photon contribution generate the destructive and constructive interference with the leading power effect (a similar observation for the high twist corrections already made in [22]) and there appears to be a strong cancellation between these two mechanisms in the whole Q 2 ≤ 40 GeV 2 region. However, both subleading power effects become rapidly suppressed with the growing of the momentum transfer squared in contrast to the numerically sizeable soft power correction estimated from the dispersion approach [22]. Such discrepancy may be ascribed to the very definition of the "soft" effect in the formalism of [26], roughly corresponding to the ρ-resonance contribution to the π 0 γ * γ form factor with the parton-hadron duality approximation, which has no transparent counterpart in the framework of perturbative QCD factorization. In addition, the NLL radiative corrections are observed to give rise to approximately O(15 %) (almost Q 2 -independent) shift to the LL predictions for both the leading power contribution and the hadronic photon effect.
To understand the phenomenological impact of the QCD resummation for the large logarithms appearing in the factorization formula for the leading power contribution and in the LCSRs for the hadronic photon correction, we further present in figure 7 our predictions for the π 0 γ * γ form factor, at LL, NLO and NLL accuracy, with the BMS model. The NLO QCD corrections are found to induce O (25 %) reduction of the tree-level results at 10 GeV 2 ≤ Q 2 ≤ 40 GeV 2 , however, the NLL resummation effect will enhance the NLO predictions by an amount of approximately O (10 %), in accordance with the pattern for the perturbative QCD corrections observed in [32,88]. Inspecting the model dependence of pionphoton form factor on the leading twist pion DA displayed in figure 7 implies that the theory predictions with both the holographic and KMOW models can reasonably balance the BaBar and Belle data at high Q 2 without resorting to the "exotic" end-point behaviour as advocated in [51,52]. In fact, we have checked that the predicted π 0 γ * γ form factor with the flat pion DA will overshoot both the BaBar and Belle data, in most Q 2 region of interest, at least in our framework. Given the fact that the end-point behaviour of the twist-two pion DA in the  [21] (orange circles) and Belle [47] (brown spades) are also displayed here. holographic model differs from the standard postulation, motivated by the conformal expansion analysis, as employed for the KMOW model, we conclude that the local information of the pion DA cannot be extracted from the experimental measurements of the pion-photon form factor even in the leading power approximation. It needs further to point out that our predictions with the holographic and KMOW models do not match the experimental data at 2 GeV 2 ≤ Q 2 ≤ 8 GeV 2 well, where the power suppressed contributions from the yet highertwist pion and photon DAs will become more pronounced and actually the large-momentum expansion applied for the construction of the factorization formula also becomes questionable. By contrast, the theory predictions from the dispersion approach [22,23] can result in a satisfactory description of the BaBar and Belle data in the whole Q 2 region by introducing the nonperturbative modification of the QCD spectral density function. Moreover, it becomes apparent that the computed pion-photon form factor with the BMS model and the asymptotic pion DA are less favorable by the experimental measurements at high Q 2 , albeit with the reasonable agreement achieved at low Q 2 . Also, confronting the theory predictions from the CZ model with the BaBar and Belle data indicates a large value of the second Gegenbauber moment a 2 (µ 0 ) is not favored, in agreement with the recent lattice QCD calculations [83,92].
We present our final predictions for the π 0 γ * γ form factor from the expression (102) with three different models of the twist-two pion DA in figure 8, including the theory uncertainties due to the variations of the input parameters discussed before. We already assigned 20% uncertainty for the first six nontrivial Gegenbauer coefficients of the holographic model in the numerical estimation for the illustration purpose. It turns out that the dominant theory uncertainties originate from the shape parameters of the leading twist pion and photon DAs instead of the variations of the factorization scales. Precision determinations of the higher Gegenbauer coefficients for both two DAs along the lines of [83,92] will be essential to pin down the presently sizeable theory uncertainty in order to meet the challenge of the (potentially) more accurate experimental measurements at the BEPCII collider [93] and the SuperKEKB accelerator.

Conclusion
Applying the standard OPE technique with the evanescent operator(s) we revisited the demonstration of QCD factorization for the pion-photon transition form factor at leading power in 1/Q 2 with both the NDR and HV schemes for γ 5 in the D-dimensional space. It has been shown explicitly at one loop that the renormalization scheme dependence of the short-distance matching coefficient and the twist-two pion DA are cancelled out precisely rendering the γ 5prescription independence of the factorization formula for the leading power contribution to F γ * γ→π 0 (Q 2 ). This can be readily understood from the fact that the QCD matrix element defined by two electromagnetic currents are free of the γ 5 ambiguity and the renormalization scheme dependence of the hard function arises from the infrared subtraction term completely. In the same vein, we established QCD factorization of the desired correlation function at one loop for the construction of the LCSRs for the hadronic photon contribution to the pionphoton form factor. By contrast, the corresponding QCD matrix element defined with an interpolating current for the pion and an electromagnetic current suffers from the γ 5 ambi-guity and the leading twist photon DA is independent of the γ 5 prescription in dimensional regularization. The finite renormalization term introduced in the HV scheme to restore the appropriate Ward-Takahashi identities was found to provide the very transformation function to construct the hard matching coefficient in the NDR scheme. The NLL resummation of the parametrically large logarithms was also implemented by solving the relevant two-loop evolution equations in momentum space.
Taking into account the leading power contribution and the hadronic photon effect at NLL and the twist-four pion DA correction at tree level, we further explored the phenomenological consequence of the perturbative QCD corrections and the subleading power contributions. Interestingly, the observed strong cancellation between the two power suppressed mechanisms leads to the insignificant correction to the leading power contribution (almost) in the whole Q 2 region accessible at the current experiments. In addition, we paid a particular attention to the model dependence of the theory predictions for F γ * γ→π 0 (Q 2 ) on the twist-two pion DA. Both the holographic and KMOW models turned out to balance the BaBar and Belle data reasonably well at high Q 2 , despite the visible discrepancy at low Q 2 which could be compensated by the unaccounted subleading power corrections of both perturbative and nonperturbative origins. It was also demonstrated that the end-point behaviour of the pion DA cannot be extracted by matching the theory predictions for the pion-photon form factor with the experimental measurements.
Aiming at a better confrontation with the BaBar and Belle data, further improvements of our calculations can be made by first carrying out the perturbative correction to the twist-4 contribution from both the two-particle and three-particle pion DAs, which is also of conceptual interest in the framework of perturbative QCD factorization, and then by evaluating the high twist contributions from the photon DAs with the LCSR approach. Phenomenological applications of the techniques discussed in this work can be also pursued in the context of the γ * γ → (η ( ) , η c ) transition form factors [33,94] for understanding the quark-gluon structure of eta mesons and heavy quarkonium states, the radiative leptonic B-meson decays for the determination of the inverse moment λ B , the radiative penguin decays of B-mesons for the precision test of the quark-flavour structure of the Standard Model, and the radiative heavyhadron decays for constraining the magnetic susceptibility of the quark condensate [86]. To conclude, the anatomy of the subleading power contributions for the exclusive hadronic reactions is of high interest for understanding the general structures of the large momentum/mass expansion in QCD and for hunting new physics in the quark-flavour sector as indicated by the various flavour "anomalies" observed at the ongoing experiments.
The manifest expressions of the RG functions E NLO T,n and d k T,n can be obtained from that of E NLO V,n and d k V,n given above with the replacement rule γ

B Spectral representations
We present the dispersion representations of convolution integrals entering the NLL QCD factorization formula (82) in order to construct the sum rules for the hadronic photon correction to the pion-photon form factor. We have verified the spectral representations in what follows numerically by checking the corresponding dispersion integrals.