Next-to-Next-to-Leading-Order QCD Prediction for the Photon-Pion Form Factor

We accomplish the complete two-loop computation of the leading-twist contribution to the photon-pion transition form factor $\gamma \, \gamma^{\ast} \to \pi^0$ by applying the hard-collinear factorization theorem together with modern multi-loop techniques. The resulting predictions for the form factor indicate that the two-loop perturbative correction is numerically important. We also demonstrate that our results will play a key role in disentangling various models of the twist-two pion distribution amplitude thanks to the envisaged precision at Belle II.


INTRODUCTION
It is widely accepted that the exclusive two-photon production of a light pseudoscalar meson π 0 is of utmost importance for probing the partonic landscape of composite hadrons with unprecedented precision and for deepening our understanding of the factorization properties of hard exclusive QCD reactions in general. Historically, the photon-pion transition form factor (TFF) had been explored intensively (see, for instance [1][2][3]) even before the advent of QCD due to its apparent connection to the axial anomaly [4][5][6] when taking the vanishing virtuality limits for the two photons. Likewise, for the energetic photo-production at the light-like distance, the dynamical behavior of the TFF γ * γ ( * ) → π 0 can be predicted by applying the operator-product-expansion (OPE) technique for the time-ordered product of two electromagnetic currents [7,8]. The appearing correlation function with the distinct kinematic set-up concurrently provides the theory description of a variety of two-photon processes including deeply inelastic lepton-hadron scattering (DIS) [9][10][11] and deeply virtual Compton scattering (DVCS) [12][13][14], whose importance in establishing QCD as the theory of strong interactions and in accessing the transverse-distance structure of the target hadron cannot be exaggerated (see, for instance [15][16][17][18]). In particular, the double-virtual photon-to-pion form factor has been demonstrated to be an indispensable ingredient for determining the hadronic light-by-light scattering (HLbL) contribution to the anomalous magnetic moment of the muon in the dispersive framework [19]. The comprehensive data-driven strategy of determining the pion-pole contribution to the HLBL scattering further illustrates that the yielding hadronic uncertainties of (g − 2) µ can be reduced substantially with the improved determination of the single-virtual photon-pion form factor [20,21]. Moreover, the collinear QCD dynamics dictating the TFF at large momentum transfer is encapsulated in the universal light-cone distribution amplitude (DA), which is undoubtedly of decisive importance for the modelindependent description of semileptonic and nonleptonic B-meson decays, such as B → π −ν , B → π + − and B → π π, on the basis of QCD factorization demanding the detailed information on the very pion DA as the fundamental non-perturbative input [22][23][24][25]. Additionally, the lattice QCD technique has been applied to address the photon-pion form factor at low photon virtualities with very encouraging predictions [26,27].
Experimentally, the pion TFF γ γ * → π 0 with one on-shell and one off-shell photon turns out to be more accessible from the "single-tagged" measurements of the differential e + e − → e + e − π 0 cross section [28][29][30]. The unexpected scaling violation of the TFF reported by the BaBar measurements [29] triggered a storm of enthusiasm in the theory community with diverse speculations from both perturbative and non-perturbative QCD perspectives. The subsequent Belle data [30] covering the same kinematical region Q 2 ∈ [4, 40] GeV 2 (where Q 2 stands for the virtuality of the off-shell photon), however, did not reveal the pronounced increase of the pion TFF in the high-Q 2 region thus creating a moderate tension at the level of 2σ in comparison with the BaBar results. The substantial improvement on the integrated luminosity and the trigger efficiency of the Belle II experiment at SuperKEKB [31] will evidently enable a clarification of the observed discrepancy between the BaBar and the Belle data and urgently demand the high-precision theory computation of the γ γ * → π 0 TFF in QCD.
At leading power (LP) in Λ 2 QCD /Q 2 the photon-pion TFF can be expressed in terms of the perturbatively calculable coefficient function (CF) and the twist-two pion DA [8] (see also [32][33][34]). The short-distance CF at the next-to-leading order (NLO) in α s had been determined more than three decades ago [35][36][37]. The next-to-nextto-leading-order (NNLO) computation of the perturbative matching coefficient was further carried out in the large-β 0 approximation [38] and in the so-called conformal scheme [39]. Consequently, the factorization-scale invariance of the resulting expression for the pion TFF was not achieved in the modified minimal subtraction (MS) scheme at two loops. Accomplishing the full NNLO QCD prediction to the TFF γ γ * → π 0 is therefore of the top priority, on the one hand, for developing the perturbative factorization formalism of hard exclusive reactions; and on the other hand, for elevating our capability to confront the yielding theory predictions with the forthcoming precision of Belle II measurements. In this Letter, we will take advantage of the QCD factorization program to obtain analytically the two-loop hard function by evaluating an appropriate bare QCD matrix element at O(α 2 s ) by means of dedicated multi-loop computational strategies and by implementing the ultraviolet (UV) renormalization and the infrared (IR) subtractions, with an exploratory phenomenological analysis regarding the numerical significance of the predicted TFF.

THE PHOTON-PION FORM FACTOR
We first set up the theory framework for establishing the hard-collinear factorization formula of the TFF F γπ (Q 2 ). The latter is defined in terms of the matrix element of the electromagnetic current between an on-shell photon with momentum p and a pion with momentum p π(p)|j em µ |γ(p ) = g 2 em µναβ q α p β ν (p )F γπ (Q 2 ) , (2) where q = p − p refers to the transfer momentum, ν (p ) is the polarization vector of the on-shell photon, and e q denotes the electric charge of the quark field in units of g em (the positron charge). Moreover, we have employed the convention 0123 = −1 and the notation Q 2 = −q 2 . Introducing further two light-cone vectors n µ andn µ with n 2 =n 2 = 0 and n ·n = 2 allows for the decomposition p µ = (n · p)/2 n µ and p µ = (n · p)/2n µ . Applying the hard-collinear factorization theorem results in the LP contribution to the γ γ * → π 0 form factor where the pion decay constant f π = (130.50±0.02±0.03± 0.13) MeV is determined from the experimental measurement [40] and µ F represents the factorization scale for which we tacitly adopt µ F = µ UV ≡ µ in (3). Here the π 0 wavefunction |π 0 = 1 √ 2 (|uū −|dd ) with exact isospin symmetry is employed. The CF T 2 can be expanded perturbatively in the form (similarly for any other partonic quantity) The twist-two pion DA φ π (x, µ F ) is defined by the renormalized matrix element on the light-cone where [zn, 0] is the Wilson line to ensure gauge invariance. We also note that φ π (x, µ F ) = φ π (x, µ F ) due to charge symmetry, wherex ≡ 1 − x. The corresponding renormalization group (RG) equation at one loop [8,41] implies the series expansion of φ π (x, µ F ) in terms of Gegenbauer polynomials, with the vanishing odd-moments a 1,3,... (µ F ). It is convenient to evaluate the CF by inspecting the following correlation function which can be parameterized by the two form factors for the bilinear quark currents with the spin structures [42] Γ µν thanks to the QED Ward identity and charge symmetry.
The perpendicular component of any four-vector p µ is defined by p µ = (np)n µ /2 + (np)n µ /2 + p µ ⊥ . We will dedicate the next section to the two-loop calculation of the bare quantity Π (2) µν and subsequently derive the master formula for the CF T 2 in (3).

TWO-LOOP CALCULATION
The techniques entering the calculation of the bare two-loop amplitude have become standard in contemporary multi-loop computations. We first generate the Feynman diagrams with both Feynarts [43] and by means of an in-house routine. Sample diagrams are shown in Fig. 1. After taking into account that certain diagrams have color factor zero, vanish due to the Furry theorem and/or represent the flavor-singlet contributions, a set of 42 diagrams (plus the ones with the two photons exchanged) has to be computed. The entire calculation is carried out in dimensional regularization with d = 4 − 2 . The Dirac and tensor reductions are performed with inhouse routines, yielding two-loop scalar integrals, which get further processed with FIRE [44], an implementation of Laporta's algorithm [45,46] relying on integrationby-parts identities [47,48]. In addition, we exploit the fact that for the quark momenta p 1 ≡ xp ∝xp ≡ p 2 , which yields relations between integrals based on momentum conservation and enables us to arrive at the minimal set of master integrals. At this stage, additional Dirac structures besides Γ µν A,B disappear from the sum of all diagrams, making the QED Ward identity and charge symmetry manifest and providing a valuable check of our calculation.
In total, we get 12 independent master integrals, the analytic -expansion of which can be written in terms of harmonic polylogarithms (HPLs) [49][50][51]. We relegate the explicit analytic results of all master integrals to a forthcoming longer write-up.

ULTRAVIOLET RENORMALIZATION AND INFRARED SUBTRACTIONS
We are now prepared to extract the two-loop contribution to T 2 by first introducing the operator basis where Γ µν with σ µν ⊥ = 1 2 [γ µ ⊥ , γ ν ⊥ ] and then by exploiting the matching equation It is evident from the definitions g µν ⊥ = g µν − n µnν /2 − n νnµ /2 and µν ⊥ ≡ µναβn α n β /2 that O µν E is an evanescent operator vanishing at d = 4, and the CP-even operator O µν 1 cannot couple to the pseudoscalar π 0 state. We therefore encounter a unique physical operator O µν 2 [42]. The correlator Π µν assumes the following form to all orders in α s in terms of the tree-level matrix elements of the effective operators O µν where Z α = 1 − a s β 0 / + O(a 2 s ). Note that, hereafter, we disregard O µν 1 completely due to parity. The expansion of T 2 was already given in Eq. (4), and due to scaleless integrals in dimensional regularization the matrix elements of the light-cone operators are expanded as Sums over repeated indices are understood to run over {2, E}, O µν a (0) is the tree-level matrix element of O µν a , and the renormalization constants Z ( ) 22 can be extracted from the Efremov-Radyushkin-Brodsky-Lepage (ERBL) kernel with the one-and two-loop result given in [8,41] and [52][53][54][55][56], respectively. Inserting Eqs. (4), (12), and (13) into Eq. (11) and comparing coefficients of O µν 2,E (0) at O(a 2 s ) leads to the following "master formula" for the CF at NNLO All quantities multiplied by a divergent term in (14) must be expanded beyond O( 0 ) to correctly capture all finite terms. It remains important to emphasize the role of the evanescent operator in the master formula (14). Despite the vanishing of evanescent operators in four dimensions, their non-vanishing bare matrix element and evanescentto-physical mixing (captured by Z E2 [57]) are crucial for correctly implementing the IR subtraction and establishing the hard-collinear factorization in our d-dimensional framework. These steps are necessary in view of the fact that taking the limit d → 4 does not commute with the hard-collinear factorization for the correlation function Π µν (see also [58,59]).
Substituting the newly derived two-loop correction to the hard function, whose explicit expression is given in the supplemental material, into the factorization formula (3) and adopting further the asymptotic pion DA enable us to write down Generalizing (15) to account for the non-asymptotic correction will give rise to lengthy expressions due to the intricacy of the NNLO CF. Consequently, we will only evaluate such higher moment effect numerically in the phenomenological exploration. In particular, it has been verified explicitly that the photon-pion form factor computed from the factorization formula (3) is independent of the factorization scale at O(α 2 s ), with the aid of the twoloop evolution equation of the twist-two pion DA [52][53][54][55][56]. Subsequently, the complete next-to-next-to-leadinglogarithmic (NNLL) resummation of log(Q 2 /Λ 2 QCD ) appearing in (3) is achieved by virtue of the three-loop evolution equation of φ π (x, µ F ) in the naive dimensional regularization (NDR) scheme [60], which is implemented throughout our computation.

NUMERICAL ANALYSIS
We now turn to explore the numerical implications of the achieved NNLL computation of F γπ (Q 2 ) with three phenomenological models of φ π (x, µ 0 ) at µ 0 = 1 GeV, Model I is inspired from the AdS/QCD correspondence [61] with modifications from the most recent lattice result a 2 (2 GeV) = 0.116 +0.019 −0.020 [62], which corresponds to a 2 (µ 0 ) = 0.176 +0.027 −0.029 . The construction of Model II is achieved by matching the light-cone QCD sum rules of the pion electromagnetic form factor in the space-like region onto the modified dispersion integral of the modulus of the time-like form factor [63]. The intervals of a 2 and a 4 in Model III [64][65][66] are determined from the method of QCD sum rules [67]. To facilitate the implementation of the three-loop evolution of the leading-twist pion DA [60,68], we will apply the Gegenbauer expansion for the "holographic"-type Model I by discarding a 2n (µ 0 ) with n ≥ 7 (see also [69]).  (16), where the red solid curve is obtained by adding further the subleading power contributions evaluated in [42,70].
Adopting Model I as our default choice, we present in Fig. 2 the resulting impacts of QCD radiative corrections to the TFF at LP as well as the numerical features of the power suppressed terms from the twist-four collinear dynamics of the pion system [70] and from the long-distance photon correction [42]. Inspecting the numerical patterns of the various dynamical mechanisms displayed in Fig. 2 implies that the NNLL twist-two correction can reduce the corresponding NLL prediction by an amount of (4 ∼ 7) % at Q 2 ∈ [3, 40] GeV 2 (in agreement with the previous estimate [65]), while the genuine one-loop perturbative correction is responsible for a constant ∼ 10% reduction of the LL contribution in the same kinematic domain. By contrast, the above-mentioned subleading contributions are observed to bring about the negligible corrections [42]. The full two-loop computation of the LP effect in the TFF is therefore indeed of extraordinary phenomenological significance. Additionally, the three-loop QCD correction estimated with the "naïve non-Abelianization" prescription [71] is expected to be approximately O(10 − 20)% of the one-loop radiative effect without adding further the systematic uncertainty.
FIG. 3. Theory predictions of the TFF γγ * → π 0 for the three models presented in (16). The color bands are due to variation of the factorization/renormalization scale µ. For a comparison, we also display the experimental measurements from the CLEO [28] (purple squares), BaBar [29] (orange circles) and Belle [30] (brown spades) Collaborations.
We proceed to present in Fig. 3 the theory predictions of the TFF for the three phenomenological models of the leading-twist pion DA. As the dominating nonperturbative uncertainties of our numerical analysis from φ π (x, µ 0 ) can be inferred from the discrepancies of the obtained theoretical results for the distinct sample models, the quoted uncertainties in (16) are not included in the three shaded (red, green, gray) bands in Fig. 3 to highlight the perturbative uncertainty from the variation of µ 2 = x Q 2 with 1/4 ≤ x ≤ 3/4 [69]. The distinctive snapshot of the well-separated uncertainty bands for the three models is particularly encouraging to differentiate numerous speculations on the intricate behaviors of the twist-two pion DA. Furthermore, the theory predictions from the incomplete two-loop computation of F γπ (Q 2 ) in the large β 0 approximation differ from our full O(α 2 s ) results by almost a factor of two at high-Q 2 . Additionally, the NNLL improved hard-collinear factorization tends to validate the scaling behaviour of the photon-pion form factor at Q 2 ≥ 20 GeV 2 with an impressive precision of O(1 %), which will be confronted with the future Belle II data. As expected, such scaling behaviour is broken down at low momentum transfer, numerically Q 2 ≤ 5 GeV 2 , due to the pronounced subleading power contributions. A more detailed investigation of the asymptotic behaviour of F γπ (Q 2 ) in the hard-collinear factorization framework can be pushed forward by introducing the scaling-rate quantity Ω(Q 2 ) as proposed in [66,72].
We finally address the actual impact of the NNLL correction to the photon-pion form factor on the model- FIG. 4. Illustration of the a4(µ0) dependence of the photonpion form factor at Q 2 = 27.30 GeV 2 with the second Gegenbauer moment determined from [62] and with the vanishing higher Gegenbauer moments (i.e., a n≥6 (µ0) = 0). independent determination of the twist-two pion DA in anticipation of the upcoming precision Belle II measurements. To this end, we investigate the intrinsic sensitivity of the γ γ * → π 0 form factor with regard to the poorly constrained Gegenbauer moment a 4 (µ 0 ), by employing the improved lattice result for a 2 (µ 0 ) [62] and by further leaving out the higher moments a n≥6 (µ 0 ) for illustration purpose. The numerical implication displayed in Fig. 4 reveals that the two-loop computation of the photon-pion TFF is indeed highly beneficial for pinning down the theory uncertainty of the extracted value of a 4 (µ 0 ) (up to a factor of two improvement approximately when compared with the corresponding one-loop result).

CONCLUSIONS
In summary, we have endeavored to carry out, for the first time, the complete two-loop computation of the leading-twist contribution to the photon-pion form factor in the hard-collinear factorization framework. The phenomenological significance of such higher-order perturbative effect was further addressed by employing the three non-perturbative models of the twist-two pion DA. The genuine O(α 2 s ) correction to the hard matching coefficient was demonstrated to be numerically important at Q 2 ∈ [3, 40] GeV 2 , and in particular, the previously evaluated radiative correction at O(α 2 s β 0 ) only accounts for a moderate portion (approximately O(50 %)) of the full NNLL QCD effect. Moreover, our improved theory predictions of the form factor γ γ * → π 0 appear to be promising for better constraining the shape parameters of the pion DA. The established two-loop QCD factorization formula of the photon-pion form factor is then naturally expected to be of notable importance for exploiting the many facets of the strong interaction dynamics embedded in hard exclusive reactions.

Supplemental Materials to Next-to-Next-to-Leading-Order QCD Prediction for the Photon-Pion Form Factor
Analytic expressions in the MS-scheme The master formula (14) was derived under the assumption that dimensional regularization is used on both sides of Eq. (11) for both UV and IR divergences. However, to determine the UV-renormalization constants Z where Z E2 is IR finite as expected albeit both M E2 and M 22 are IR divergent. We now present the explicit expressions of the two-loop hard coefficient function T 2 (x) = β 0 C F K (2) β (x)/x + K (2) β (x)/x + C 2 F K F (x)/x + K F (x)/x + C F /N c K N (x)/x + K H 0 (x)