Next-to-next-to-leading-order QCD corrections to pion electromagnetic form factors

We investigate the next-to-next-to-leading order (NNLO) QCD radiative corrections to the pion electromagnetic form factor with large momentum transfer. We explicitly verify the validity of the collinear factorization to two-loop order for this observable, and obtain the respective IR-finite two-loop hard-scattering kernel in the closed form. The NNLO QCD correction turns to be positive and significant. Incorporating this new ingredient of correction, we then make a comprehensive comparison between the finest theoretical predictions and numerous pion form factor measurements in both space-like and time-like regions. Our phenomenological analysis provides strong constraint on the second Gegenbauer moment of the pion light-cone distribution amplitude (LCDA) obtained from recent lattice QCD studies.

Introduction.Originally proposed by Yukawa as the strong nuclear force carrier in 1935 [1], subsequently discovered in the cosmic rays in 1947 [2], the π mesons have always occupied the central stage throughout the historic advancement of the strong interaction.As the lightest particles in the hadronic world (hence the highlyrelativistic bound systems composed of light quark and gluons), π mesons entail extremely rich QCD dynamics, exemplified by the color confinement and chiral symmetry breaking.Notwithstanding extensive explorations during the past decades, there still remain some great myths about the internal structure of the π mesons.
A classic example of probing the internal structure of the charged π is the pion electromagnetic (EM) form factor: with Q 2 ≡ −(P − P ′ ) 2 .The electromagnetic current is defined by J µ em = f e f f γ µ f , with e u = 2/3 and e d = −1/3 indicating the electric charges of the u and d quarks.
During the past half century, the pion EM form factor has been intensively studied experimentally .From the theoretical perspective, the pion EM form factor at small Q 2 can be investigated in chiral perturbation theory [30] and lattice QCD [31][32][33][34][35]. On the other hand, at large momentum transfer, the F π (Q 2 ) is expected to be adequately described by perturbative QCD.Within the collinear factorization framework tailored for hard exclusive reactions [36][37][38][39][40][41][42] (for a review, see [43]), at the lowest order in 1/Q, the pion EM form factor can be expressed in the following form: where T (x, y) signifies the perturbatively calculable hardscattering kernel, and Φ π (x, µ F ) represents the nonperturbative yet universal leading-twist pion light-cone distribution amplitude (LCDA), i.e., the probability amplitude of finding the valence u and d quark inside π + carrying the fractional momenta x and x ≡ 1 − x, respectively.The leading-twist pion LCDA assumes the following operator definition: with W signifies the light-like gauge link to ensure the gauge invariance.Conducting the UV renormalization for (3), one is led to the celebrated Efremov-Radyushkin-Brodsky-Lepage (ERBL) evolution equation [38,40]: with V (x, y) referring to the perturbatively calculable ERBL kernel.Eq. ( 2) is expected to hold to all orders in perturbative expansion.The hard-scattering kernel can thus be FIG. 1: Sample parton-level Feynman diagrams for the reaction γ * π(P ) → π(P ′ ) at various perturbative orders.
expanded in the power series: , and N c = 3 is the number of colors.The leading order (LO) result was known shortly after the advent of QCD [37,38,40,42,[44][45][46]. The next-toleading order (NLO) correction was originally computed by three groups in early 80s [47][48][49].Scrutinizing the previous calculations, in 1987 Braaten and Tse traced the origin of the discrepancies among the earlier work and presented the correct expression of the NLO hardscattering kernel [50].In 1998, Melić, Nizić, and Passek conducted a comprehensive phenomenological study by incorporating the NLO correction as well as the evolution effect of pion LCDA [51].The central goal of this work is to compute the next-to-next-to-leading order (NNLO) perturbative correction to pion EM form factor, and critically examine its phenomenological impact 1 .Setup of perturbative matching.The strategy of deducing the short-distance coefficients is through the standard matching procedure.Since the hard-scattering kernel is insensitive to the long-distance physics, it is legitimate to replace the physical π + by a free quark-antiquark pair u d, and compute both sides of (2) in perturbation theory, order by order in α s .To make things simpler, we neglect the transverse motion and assign the momenta of the u and d in the incoming "pion" to be uP and ūP , and assign the momenta of the u and d in the outgoing "pion" to be vP ′ and vP ′ , with u, v ranging from 0 to 1.
In the left-hand side of (2), we extract the scalar form factor F (u, v) through the partonic reaction γ * + u(uP ) d(ūP ) → u(vP ′ ) d(vP ′ ).Some typical Feynman diagrams through two-loop order are depicted in Fig. 1.It is subject to a perturbative expansion: 1 It is worth mentioning that, historically there has been debate on the applicability of collinear factorization at experimentally accessible moderate energy range, say, Q < 6 GeV [52][53][54][55].The central issue is whether the nonfactorizable soft mechanism can dominate the perturbative QCD contribution or not.
Once beyond the tree level, the UV and IR divergences inevitably arise and we use the dimensional regularization (DR) to regularize both types of divergences.Nevertheless, the bare "pion" LCDA remains intact since the scaleless integrals vanish in DR.The renormalized "pion" LCDA is related to the bare one via which is solely comprised of various IR poles.Z(x, y) in ( 9) signifies the renormalization function in the MS scheme, which can be cast into the following Laurent-expanded form in ǫ: Note that the prefactor of single pole in (10) is related to the ERBL kernel V (x, y) in (4) via V (x, y) = −α s ∂Z 1 /∂α s [56].Note that the two-loop [49,[57][58][59][60] and three-loop corrections [61] to V (x, y) have been available.
The two-loop renormalized "pion" LCDA Φ (2) also contains double IR pole.The Z 2 can be obtained through the recursive relation [62] where With the aid of ( 9) and ( 10), we then determine the O(α s ) and O(α 2 s ) corrections to the renormalized "pion" LCDA in (7).
At one-loop order, the matching equation for a fictitious pion state becomes where ⊗ x signifies the convolution over x.Note that the renormalized scalar form factor F (1) (u, v) still contains single collinear pole.However, the renormalized Φ (1) (x|u) and Φ (1) (y|v) also contains the same IR poles.
Upon solving this matching equation, one ends with both UV and IR-finite T (1) (x, y).Our expressions agree with the known NLO result [51].
To the desired two-loop order, the following matching equation descends from (2): +Φ (1) (x|u) ⊗ x More Severe IR divergences are expected to arise in both F (2) (u, v) and Φ (2) (x|u).Clearly one also needs compute Description of the calculation.We use HepLib [63] and FeynArts [64] to generate Feynman diagrams and the corresponding amplitudes for the reaction γ * + u(uP ) d(ūP ) → u(vP ′ ) d(vP ′ ).We employ the covariant projector technique to enforce each u d pair to bear zero helicity.For our purpose it suffices to adopting the naive anticommutation relation to handle γ 5 in DR.There are about 1600 two-loop diagrams, one of which is depicted in Fig. 1c).We employ the package Apart [65] to conduct partial fraction, and FIRE [66] for integration-bypart reduction.We end up with 116 independent master integrals (MIs).The MIs are calculated by utilizing the differential equations method [67][68][69].Note that these MIs are considerably more involved than than what are encountered in the two-loop corrections for the π − γ transition form factor [70,71].We have attempted two independent ways to construct and solve the differentialequation systems, one of which is based on the method developed in [72][73][74][75].The analytic results are expressed in terms of the Goncharov Polylogarithms (GPLs) [76].Two independent calculations yield the identical answer.We also numerically check our results against the package AMFLOW [77] and found perfect agreement.Technical details will be included in the future long write-up.Upon renormalizing the QCD coupling in MS scheme, we end up with a rather lengthy expression for F (2) (u, v).
Being UV finite, it still contains severe IR divergences which start at order-1/ǫ 2 IR .Inspecting the matching equation ( 13), piecing all the known ingredients together, we are able to solve for the intended two-loop hardscattering kernel.Hearteningly, T (2) (x, y) is indeed IR finite.Therefore, our explicit calculation verifies that the collinear factorization does hold at two-loop level for the pion EM form factor.The analytical expression of T (2) (x, y) is too lengthy to be reproduced here.For the sake of clarity, in the suppletory material we provide the asymptotic expressions of T (1,2) (x, y) near the endpoint regions.
Master formula for pion EM form factor at NNLO.Given a certain parametrized form of pion LCDA, the two-fold integration in (2) turns out to be difficult to conduct numerically, mainly due to numerical instability caused by the spurious singularity as x → y/x → ȳ in T (2) (x, y).Our recipe to circumvent this technical challenge is to present the two-loop results in an analytical manner, which enables us to achieve exquisite precision.
The leading-twist pion LCDA is conveniently expanded in the Gegenbauer polynomial basis: where the pion decay constant f π = 0.131 GeV, and ′ signifies the sum over even integers.Note all the nonperturbative dynamics is encoded in the Gegenbauer moments a n (µ F ).
Substituting ( 14) into (2), we reexpress the pion EM form factor as with mn defined by (16) For simplicity, we will set µ R = µ F = µ and n L = 3 from now on.The two-fold integrations in ( 16) can be readily worked out at tree and one-loop levels.For instance, we have with Remarkably, the two-loop coefficients T mn can also be computed analytically, thanks to the fact that T (2) can be expressed in terms of the GPLs.Although the integrand in ( 16) contains about O(10 5 ) individual terms, the final result after two-fold integration becomes exceedingly compact, which can be expressed in terms of the rational numbers and Riemann zeta function.For instance, the expression of T With the input from Table I, Eq. ( 15) constitutes our master formula for yielding phenomenological predictions through the two-loop accuracy.Compared with the original factorization formula (2), we have simplified an integration task into an algebraic one.
It is straightforward to adapt our master formula from the space-like region to the time-like one, provided that one makes the replacement L µ → L µ + iπ in Table I, with Q 2 now indicating the squared invariant mass of the π + π − pair.
Input parameters.As the key nonperturative input, our knowledge on the pion LCDA is still not confirmative enough.In the early days, it is popular to assume asymptotic form, CZ parametrization [43] and the BSM parameterizations [78,79].In recent years there have emerged extensive investigations on the profile of the pion DA from different methodologies, including QCD light-cone sum rule [80] with nonlocal condensate [78,81] or fitted from dispersion relation [82] or Platykurtic [83], Dyson-Schwinger equation [84,85], basis light-front quantization [86], light-front quark model [87], holographic QCD [88], and very recently, from the lattice simulation [89,90].The predicted values of various Gegenbauer moments are scattered in a wide range.
Since lattice QCD provides the first-principle predictions, in this Letter we will take the most recent lattice results as inputs.In 2019 RQCD Collaboration has presented a precise prediction for the second Gegenbauer moment of pion LCDA in scheme, with a 2 (2 GeV) = 0.116 +0.019 −0.020 [89].
An important progress in lattice QCD is expedited by the advent of the Large-momentum effective theory (LaMET) a decade ago [91,92], which allows one to access the light-cone distributions in Euclidean lattice directly in the x space.Very recently, the LPC Collaboration has presented the whole profile of the pion LCDA [90], from which various Gegenbauer moments can be inferred: a 2 (2 GeV) = 0.258 ± 0.087, a 4 (2 GeV) = 0.122 ± 0.056, a 6 (2 GeV) = 0.068 ± 0.038.It is curious that the value of a 2 reported by LPC Collaboration is about twice greater than that reported by the RQCD Collaboration.This discrepancy might be attributed to the fact that the LaMET approach receives large power correction in the endpoint region.On the other hand, it is very challenging for the local operator matrix element approach [89] to compute the higher Gegenbauber moments, thus difficult to reconstruct the whole profile of the LCDA.Phenomenological exploration.We use the three-loop evolution equation [61,93] to evolve each a n evaluated at 2 GeV by lattice simulation to any intended scale µ.We only retain those Gegenbauer moments with n up to 6.We also use the package FAPT [94] to evaluate the running QCD coupling constant to three-loop accuracy.
For the sake of comparison, we take the pion EM form factor data in the spacetime region from NA7 collaboration [11], Cornell data compiled by Bebek et al. [5], and the reanalyzed JLab data [16], and take the timelike pion EM form factor data entirely from the BaBar experiment [27].We discard many irrelevant small-Q 2 data.In Fig. 2 and Fig. 3, we confront our predictions at various perturbative accuracy with the available data, including both space-like and time-like regime.One clearly visualizes that NNLO correction is positive and substantial.In Fig. 2, we set the various Gegenbauer moments of pion LCDA to the central values given by LPC Collaboration [90].It appears the NNLO predictions are significantly overshooting the experimental data at large Q 2 (> 5 GeV 2 ), especially for the time-like regime with high statistics data.This symptom is mainly due to the large value of a 2 .
In Fig. 3 we present our predictions with a 2 taken from RQCD while setting the values of a 4 and a 6 to zero.We find satisfactory agreement between our NNLO predictions and the data, both in space-like and time-like regimes.This may indicates that the value of a 2 given by RQCD might be more trustworthy.It is of utmost importance for RQCD and LPC collaborations to settle the discrepancy in the value of a 2 in the future.
The prospective Electron-Ion Collider (EIC) program plans to measure the space-like pion EM form factor with Q 2 as large as 30 GeV 2 [95], where perturbative QCD should be definitely reliable.We are eagerly awaiting to confronting our NNLO predictions with the future EIC data.Summary.In this work we report the first calculation of the NNLO QCD corrections to the pion electromagnetic form factor.We have explicitly verified the validity of the collinear factorization to two-loop order for this observable, and obtain the UV, IR-finite two-loop hard-scattering kernel in closed form.The NNLO QCD correction turns to be positive and important.We then confront our finest theoretical predictions with various space-like and time-like pion form factor data.Our phenomenological study reveals that adopting the second Gegenbauer moment computed by RQCD can yield a decent agreement with large-Q 2 data (above the resonance region in the time-like case).Nevertheless, to make a definite conclusion, it seems imperative to resolve the discrepancy between LPC and RQCD Collaboration on the value of a 2 in the future study.Furthermore, we look forward to the future high-statistics larger-Q 2 pion EM form factor data for critically testing our NNLO predictions.It is also very interesting to confront our NNLO predictions with the available high-quality kaon EM form factor data.
The limiting behavior of T (1) (x, y, µ) in two other corners, x → 1, y → 1 and x → 1, y → 0, can be obtained from the above formulas by making the substitutions x → x, y → ȳ and e u ↔ −e d .
For the two-loop hard-scattering kernel, we have lim x→0 y→0 Similar to the one-loop case, the limiting behavior of T (2) (x, y, µ) in two other corners, x → 1, y → 1 and x → 1, y → 0, can also be deduced by making the substitutions x → x, y → ȳ and e u ↔ −e d .
tions from the last two regions cancel with each other, and the first two regions yield identical results once the replacement e d ↔ −e u is made.For the sake of clarity, we also tabulate the numerical results of T (k) mn asy in Table II.In comparison with Table I in the main text, we observe that the agreement between T (1,2) mn asy and the exact results becomes increasingly satisfactory as m, n increase.

4 FIG. 2 :
FIG. 2: Theoretical predictions vs. data for Q 2 Fπ(Q 2 ) in the space-like (left panel) and time-like (right panel) regions.We take the central values of a2, a4 and a6 determined by LPC.The red, green and blue curves correspond to the LO, NLO and NNLO results, and the respective bands are obtained by sliding µ from Q/2 to Q. Experimental data points are taken from NA7 [11], Bebek et al. [5], Huber et al. [16] and BaBar [27].

FIG. 3 :
FIG.3: Same as Fig.2, except the predictions are made by taking the central value of a2 determined by RQCD, with a4 and a6 set to zero.

TABLE I :
The numerical values for T