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

It is generally accepted that the pion electromagnetic form factor (EMFF) is of paramount importance for probing the infrared structure of QCD scattering amplitudes of hard exclusive reactions at leading power and beyond and for unraveling the precise mechanisms that dictate the intricate nature of composite hadron systems. Advancing our understanding towards this gold-plated form factor is especially suitable for delivering deep and far-reaching insights into the delicate interplay between the emergent and Higgs mass generation mechanisms for such lightest pseudoscalar mesons as the Nambu-Goldstone modes of QCD. Historically, the systematic analysis of the pion form factor played an indispensable role in establishing the perturbative factorization formalism for the entire domain of exclusive hadronic processes with large momentum transfers Lepage and Brodsky (1979, 1980); Efremov and Radyushkin (1980a, b); Duncan and Mueller (1980a, b); Brodsky and Lepage (1989). As a consequence, the dimensional-counting rules for the power-law behaviour of numerous hard exclusive reactions on the basis of the parton model Brodsky and Farrar (1973); Matveev et al. (1973) can be placed on a solid footing within this field-theoretical framework (see Gross et al. (2023) for an overview). Additionally, the developed perturbative factorization technique has been extended to the first-principles calculations of the deeply virtual Compton scattering Ji (1997); Radyushkin (1997); Diehl et al. (1999a, b), the hard diffraction of light mesons at high energies Collins et al. (1997); Brodsky et al. (1994); Pire et al. (2021), and the exclusive heavy-to-light $B$-meson decays at large hadronic recoil Beneke and Feldmann (2001, 2004); Bauer et al. (2001, 2003). Moreover, the time-like scalar-current pion form factors are intimately connected with the chirally-enhanced weak annihilation contributions to the flagship charmless two-body hadronic $B\to\pi\pi$ decays Beneke et al. (2001); Duraisamy and Kagan (2010); Lü et al. (2023), which provide an essential testing ground for CP violation in the Standard Model (SM) with controllable theoretical uncertainties.

Experimentally, the pion EMFF at low momentum transfers ($Q^{2}\in[0.015,\,0.253]\,{\rm GeV^{2}}$) has been measured from the elastic scattering of high-energy pions off the atomic electrons in a liquid-hydrogen target at Fermilab Dally et al. (1981, 1982) and CERN Amendolia et al. (1984, 1986). The challenging measurement of the pion form factor at intermediate momentum transfers, up to $Q^{2}=10\,{\rm GeV^{2}}$, has been further carried out at Cornell Bebek et al. (1976a, b, 1978), DESY Brauel et al. (1979); Ackermann et al. (1978) and JLab Volmer et al. (2001); Tadevosyan et al. (2007); Horn et al. (2006); Blok et al. (2008); Huber et al. (2008), by exploiting the so-called Sullivan mechanism and by extrapolating the determined longitudinal cross section at negative values of the Mandelstam variable $t$ to the pion pole at $t=m_{\pi}^{2}$ Sullivan (1972). Extracting the charged pion form factor over a wide range of the moderate momentum transfers ($Q^{2}\in[2.0,\,6.0]\,{\rm GeV^{2}}$) can be also anticipated from the E12-06-101 experiment with the upgraded JLab accelerator Arrington et al. (2022). Accessing the space-like pion form factor at yet higher momentum transfers ($10\,{\rm GeV^{2}}\leq Q^{2}\leq 40\,{\rm GeV^{2}}$) with high precision will be one of the key targets for the forthcoming EIC experiment at BNL Abdul Khalek et al. (2022).

According to the hard-collinear factorization theorem Lepage and Brodsky (1979, 1980); Efremov and Radyushkin (1980a, b); Duncan and Mueller (1980a, b), the leading-power contribution to the pion EMFF in the large-momentum expansion can be expressed in terms of the short-distance matching coefficient and the twist-two light-cone distribution amplitude (LCDA) for the (anti)-collinear pion state. The next-to-leading order (NLO) QCD correction to the perturbatively calculable hard function had been explicitly determined with the diagrammatic factorization approach more than forty years ago Field et al. (1981); Dittes and Radyushkin (1981); Sarmadi (1982); Khalmuradov and Radyushkin (1985); Braaten and Tse (1987) (see Melic et al. (1999) for further verification). Including the complete NLO radiative correction in the factorization analysis of the pion EMFF turns out to be extraordinarily advantageous to pin down the intrinsic theory uncertainty from varying the renormalization and factorization scales in comparison with the counterpart tree-level prediction Melic et al. (1999). Accomplishing the next-to-next-to-leading-order (NNLO) computation of the short-distance coefficient function rigorously is therefore in high demand, on the one hand, for enriching and developing the perturbative factorization formalism at an unprecedented level, and on the other hand, for improving further the resulting theory prediction of the charged pion form factor in order to match the ever-growing precision of the dedicated experimental measurements at JLab and EIC. It is our primary objective to fill such an important and longstanding gap in this Letter, by evaluating an appropriate $5$-point partonic matrix element at ${\cal O}(\alpha_{s}^{3})$ with the contemporary advanced multi-loop computational strategies and by carrying out the ultraviolet (UV) renormalization and infrared (IR) subtractions in a meticulous and factorization-compatible manner CCF . Phenomenological implications of the thus determined two-loop QCD correction to the space-like pion form factor at large momentum transfers will be then explored comprehensively, by adopting four sample models for the leading-twist pion distribution amplitude at a low reference scale.

We first lay out the theoretical framework for constructing the hard-collinear factorization formula of the pion EMFF at leading power in an expansion in powers of $\Lambda_{\rm QCD}^{2}/Q^{2}$, where $\Lambda_{\rm QCD}$ stands for the strong interaction scale. Applying the general decomposition for the hadronic matrix element of the quark electromagnetic current enables us to write down

where $p$ and $p^{\prime}$ correspond to the four-momenta carried by the initial and final pion states, respectively. Employing the vector-current conservation condition immediately leads to $\tilde{F}_{\pi}(Q^{2})=0$, thus leaving us with a single form factor $F_{\pi}(Q^{2})$. Implementing further the charge conjugation transformation for the matrix element on the left-hand side of (1) indicates an isospin relation between the two charged pion form factors $F_{\pi^{-}}(Q^{2})=-F_{\pi}(Q^{2})$ for an arbitrary value of $Q^{2}$. In addition, the electric charge conservation determines the normalization condition for the pion form factor in the forward-scattering limit $F_{\pi}(0)=1$ (see, for instance, Khodjamirian (2020) for an elementary discussion). Here the customary notation $Q^{2}=-(p-p^{\prime})^{2}$ has been employed and the quark electromagnetic current is given by

with $e_{u}=2/3$ and $e_{d}=-1/3$ for the up and down quarks. Introducing two light-like reference vectors $n_{\mu}$ and $\bar{n}_{\mu}$ with the constraints $n^{2}=\bar{n}^{2}=0$ and $n\cdot\bar{n}=2$ then allows for the decomposition $p_{\mu}=({n\cdot p}/2)\,\bar{n}_{\mu}$ and $p^{\prime}_{\mu}=({\bar{n}\cdot p^{\prime}}/2)\,n_{\mu}$ at the leading-power accuracy.

Taking advantage of the modern effective theory field technique Bauer et al. (2002); Rothstein (2004), the hard-collinear factorization formula for the pion form factor at large momentum transfers can be cast in the desired form (see Li and Sterman (1992); Li et al. (2011, 2013, 2014); Li and Wang (2015) for an alternative formalism)

which is valid to all orders in perturbation theory and at leading power in the $\Lambda_{\rm QCD}^{2}/Q^{2}$ expansion. We take the charged pion decay constant from the three-flavour FLAG average $f_{\pi}=(130.2\pm 1.2)\,{\rm MeV}$ Aoki et al. (2022) with an increased uncertainty from adding an approximate $0.7\%$ charm sea-quark contribution. Apparently, the hard-scattering kernel $T_{1}$ depends on both the renormalization scale $\nu$ and the factorization scale $\mu$, the latter of which corresponds to the resolution with which the microscopic structure of the $\pi$-meson is being probed. This short-distance coefficient function can be expanded perturbatively in terms of the strong coupling constant (similarly for any other QCD quantity)

The leading-twist pion distribution amplitude $\phi_{\pi}$ in the factorized expression (3) can be defined by the renormalized QCD matrix element on the light-cone Braun and Filyanov (1989, 1990)

where $[\tau\,\bar{n},0]$ is the finite-length collinear Wilson line ensuring gauge invariance. The one-loop renormalization-group (RG) equation for the light-ray operator on the left-hand side of (5) implies the conformal partial wave expansion of $\phi_{\pi}(x,\mu)$ in terms of the Gegenbauer polynomials with multiplicatively renormalizable coefficients (see Braun et al. (2003) for an excellent overview)

The normalization condition $\int_{0}^{1}dx\,\phi_{\pi}(x,\mu)=1$ has been adopted throughout this work and the odd moments $a_{1,3,...}(\mu)$ vanish due to the ${\rm G}$-parity symmetry.

The short-distance matching coefficient $T_{1}$ can be routinely determined by investigating the $5$-point QCD matrix element below

where the external momenta can be restricted to their leading components $p_{1}=x\,p$, $p_{2}=\bar{x}\,p$, $p_{1}^{\prime}=y\,p^{\prime}$ and $p_{2}^{\prime}=\bar{y}\,p^{\prime}$ with the “bar notation” $\bar{x}\equiv 1-x$, $\bar{y}\equiv 1-y$. We will perform the computation of the bare amplitude in dimensional regularization with $D=4-2\,\epsilon$, where UV and IR divergences manifest themselves as poles up to the second order in $\epsilon$. The former divergences are evidently cancelled by the UV renormalization for the strong coupling constant $\alpha_{s}$ in the $\overline{\rm MS}$ scheme, while the latter disappear after executing the nowadays standard IR subtraction procedure.

We will dedicate this section to a brief technical description of the two-loop QCD calculation for the considered partonic quantity $\Pi_{\mu}$. We start with generating the NNLO Feynman diagrams with FeynArts Hahn (2001) and, independently, by means of an in-house routine. Taking into account the observation that a subset of the two-loop diagrams yield the vanishing contribution, due to the Furry theorem and/or the zero colour (electric) charge factors, we eventually encounter $1066$ non-vanishing diagrams of our interest. The two sample Feynman diagrams are depicted in Figure 1 explicitly.

All tensor integrals can be expressed in terms of scalar integrals and algebraic tensorial structures using the Passarino-Veltman decomposition Passarino and Veltman (1979). We then perform the Dirac and tensor reductions with in-house Mathematica routines, by applying further the QCD equations of motion and on-shell conditions. The large amount of two-loop scalar integrals are further reduced to an irreducible set of master integrals with the packages Apart Feng (2012) and FIRE Smirnov (2008). In total, we obtain $57$ two-loop master integrals that appear in our computation, $11/16/30$ of which appear to depend on $1/2/3$ physical scale(s). Unsurprisingly, we encounter the entire list of the $12$ one- and two-scale master integrals entering the NNLO computation of the photon-pion form factor Gao et al. (2022a, b). The evaluation of the remaining master integrals can be accomplished in an analytic fashion by employing the method of differential equations Kotikov (1991); Gehrmann and Remiddi (2000). In order to identify a suitable linear transformation such that the master integrals in the new basis satisfy a set of differential equations in an $\epsilon$-factorized form (i.e., canonical form) Henn (2013), we take advantage of Lee’s algorithm Lee (2015) as implemented in the program Libra Lee (2021). The boundary conditions to the differential equations are determined by capitalizing on the package AMFlow Liu and Ma (2023) based upon Liu et al. (2018); Liu and Ma (2022). The yielding numerical results of the master integrals at three distinct kinematic points $(x,y)\in\{(1/2,1/5),\,(1/3,1/6),(1/4,1/8)\}$, with approximately $100$ digits, allow us to construct the desired boundary constants of the differential equations in terms of transcendental numbers using the PSLQ algorithm Ferguson et al. (1999, 1998). It is thus straightforward to express the master integrals through Goncharov polylogarithms (GPLs) Goncharov (1998, 2001) up to the required order in $\epsilon$. We will present the analytic expressions for all master integrals in the forthcoming write-up.

We proceed to deduce the master formula for the hard-scattering kernel $T_{1}$ by performing the UV renormalization and IR subtractions in the presence of evanescent operators, whose significance in the perturbative treatment of quantum field theory has been explored at length ever since the pioneering era of dimensional regularization Breitenlohner and Maison (1977); Bonneau (1980a, b); Collins (2023). To achieve this goal, we exploit the matching equation for the interested QCD amplitude $\Pi_{\mu}$ onto the effective matrix elements

with the following choice of the collinear operator basis

One can readily verify that ${\cal O}_{1}$ is the only physical operator in our problem and the other two operators are evanescent, i.e., vanish algebraically in four dimensions. To reduce our notation to the essentials, we strip off the collinear Wilson lines and the position arguments of quark fields from ${\cal O}_{k}$, and represent them merely by their flavour and Dirac structure. Following the conventions of Beneke and Jager (2006), we label the collinear and anti-collinear effective fields moving into the directions of $\bar{n}$ and $n$ by $\xi$ and $\chi$, respectively, which evidently satisfy $\not{\bar{n}}\,\xi=0$ and $\bar{\chi}\,\not{n}=0$.

We are now in a position to organize the perturbative expansion of the renormalized QCD correlation function $\Pi_{\mu}$ in terms of the partonic tree-level matrix elements of the effective operators ${\cal O}_{k}$

where the renormalization constant of the strong coupling in the standard $\overline{\rm MS}$ scheme takes the form of $Z_{\alpha}=1-{\alpha_{s}\over 4\pi}\,{1\over\epsilon}\,\beta_{0}\,+\left({\alpha_{s}\over 4\pi}\right)^{2}\left({1\over\epsilon^{2}}\,\beta_{0}^{2}-{1\over 2\epsilon}\,\beta_{1}\right)+{\cal O}(\alpha_{s}^{3})$ (see Herzog et al. (2017); Luthe et al. (2017); Chetyrkin et al. (2017) for the five-loop expression). The scalar quantities $A_{k}^{(\ell)}$ in (10) represent the bare $\ell$-loop on-shell QCD amplitudes. In the same vein, the UV-renormalized matrix elements of the collinear operators ${\cal O}_{k}$ can be expanded according to

further taking into account that scaleless integrals vanish in dimensional regularization. Since the collinear fields for distinct directions already decouple at the hard scale, the renormalization constant $Z_{11}^{(\ell)}$ can be therefore determined from the celebrated Efremov-Radyushkin-Brodske-Lepage (ERBL) kernel with the one- and two-loop results obtained in Lepage and Brodsky (1980); Efremov and Radyushkin (1980b) and Sarmadi (1984); Dittes and Radyushkin (1984); Katz (1985); Mikhailov and Radyushkin (1985); Belitsky et al. (1999). Substituting Eqs. (4), (10) and (11) into the matching equation (8) leads us to derive the master formulae for the short-distance coefficient function

by comparing the coefficient of the tree-level matrix element of the physical operator $\langle{\cal O}_{1}\rangle^{(0)}$. Following the prescription proposed in Dugan and Grinstein (1991); Herrlich and Nierste (1995); Buras and Weisz (1990); Jamin and Pich (1994); Buras (2023), the renormalization constants for the evanescent operators are adjusted to ensure that the IR-finite matrix elements $\langle{\cal O}_{k}\rangle$ ($k=2,3$) vanish at an arbitrary scale. Adopting the preferred collinear operator basis (9) yields the peculiar evanescent-to-physical operator mixing under the RG evolution, which turns out to be adequately captured by the finite renormalization constant $Z_{21}^{(2)}$ at two-loop order. This justifies the essential role of introducing evanescent operators in constructing the QCD factorization formulae with dimensional regularization beyond the leading logarithmic approximation Beneke and Jager (2006, 2007); Becher and Hill (2004); Hill et al. (2004); Beneke and Yang (2006); Buchalla et al. (1996); Buras et al. (2000); Wang and Shen (2015, 2017); Gao et al. (2020, 2022a); Huang et al. (2024); Li et al. (2020); Cui et al. (2023).

Inserting the newly obtained two-loop expression of the hard-scattering kernel $T_{1}$ into the collinear factorization formula (3) and then performing the two-fold integration over $x\in(0,1)$ and $y\in(0,1)$ analytically with the Mathematica package PolyLogTools Duhr and Dulat (2019) in the asymptotic approximation (namely, $\phi_{\pi}^{\rm Asy}(x,\mu)=6\,x\,(1-x)$) results in

where $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ and $C_{A}=N_{c}$ denote the Casimir operators of the fundamental and adjoint representations of the ${\rm SU}(N_{c})$ gauge group with the standard normalization $T_{F}=1/2$. In addition, $n_{\ell}$ stands for the number of active quark flavours, and $\zeta_{n}$ represents the Riemann zeta function with $\zeta_{2}=\pi^{2}/6$, $\zeta_{3}\cong 1.202056903$, $\zeta_{4}=\pi^{4}/90$, and $\zeta_{5}\cong 1.036927755$ Vermaseren (1999); Smirnov (2012). Including the higher conformal spin contributions from the twist-two pion LCDA $\phi_{\pi}$ leads to the lengthy results for the non-asymptotic corrections to the charged pion form factor, whose explicit expressions with the $m=12$ truncation of the Gegenbauer expansion are collected in the Supplemental Material for completeness. In particular, we have verified that the thus achieved NNLO QCD computation of $F_{\pi}(Q^{2})$ with the perturbative factorization formula (3) is truly independent of the renormalization/factorization scale at the ${\cal O}(\alpha_{s}^{3})$ accuracy, by applying the two-loop evolution equation of $\phi_{\pi}(x,\mu)$ Sarmadi (1984); Dittes and Radyushkin (1984); Katz (1985); Mikhailov and Radyushkin (1985); Belitsky et al. (1999). Taking advantage of the momentum-space RG formalism enables us to accomplish an all-order summation of the enhanced logarithms of $Q^{2}/\Lambda_{\rm QCD}^{2}$ entering the factorized expression of $F_{\pi}(Q^{2})$ in the next-to-next-to-leading-logarithmic (NNLL) approximation, which necessitates an implementation of the three-loop evolution of the leading-twist pion distribution amplitude Braun et al. (2017); Strohmaier (2018).

We are now prepared to explore the phenomenological implication of the newly determined two-loop QCD correction to the pion EMFF, with an emphasis on the detailed comparison with both the available experimental measurements and the state-of-the-art lattice QCD results. To achieve this goal, we proceed by first discussing our choice for the phenomenological models of the twist-two pion LCDA appearing in the hard-collinear factorization formula (3). The first model $\phi_{\pi}^{\rm Model\,\,I}(x,\mu_{0})=\left[\Gamma(2+2\alpha_{\pi})/\Gamma^{2}(1+\alpha_{\pi})\right]\,(x\bar{x})^{\alpha_{\pi}}$ is motivated from the anti-de Sitter-QCD correspondence Brodsky and de Teramond (2008) with the non-perturbative parameter $\alpha_{\pi}(\mu_{0})=0.585^{+0.061}_{-0.055}$ Khodjamirian et al. (2021) determined by matching to the updated lattice result of the second Gegenbauer moment $a_{2}(\mu_{0})=0.116^{+0.019}_{-0.020}$ at the reference scale $\mu_{0}=2.0\,{\rm GeV}$ Bali et al. (2019). Our second model $\{a_{2},\,a_{4},a_{6},a_{8}\}\,(\mu_{0})=\{0.181(32),\,0.107(36),\,0.073(50),\,0.022(55)\}$ Cheng et al. (2020) is obtained from the comparison of the light-cone sum rule for the pion EMFF Braun et al. (2000) with the experimental measurements (see Bijnens and Khodjamirian (2002); Khodjamirian et al. (2011) for the earlier construction along this line). By contrast, the numerical intervals for the two lowest conformal coefficients in model III $\{a_{2},\,a_{4}\}\,(\mu_{0})=\{0.149^{+0.052}_{-0.043},\,-0.096^{+0.063}_{-0.058}\}$ Bakulev et al. (2001); Mikhailov et al. (2016); Stefanis (2020) are extracted from the method of QCD sum rules with non-local condensates Mikhailov and Radyushkin (1992). Furthermore, we employ an alternative model (hereafter labelled as model IV) with the first three non-vanishing Gegenbauer coefficients $\{a_{2},\,a_{4},\,a_{6}\}\,(\mu_{0})=\{0.196(32),\,0.085(26),\,0.056(15)\}$ Cloet et al. (2024) determined from the lattice prediction for the dynamic shape of the pion distribution amplitude in the range of $x\in[0.25,0.75]$ using large momentum effective theory Ji (2013, 2014); Ji et al. (2021) and from a simple power-law parametrization of the corresponding end-point behaviour. Additionally, we will vary the renormalization scale for the strong coupling $\alpha_{s}$ in the interval $\nu^{2}\in[Q^{2}/2,\,2\,Q^{2}]$ with the central value $Q^{2}$. Following Agaev et al. (2011), the factorization scale $\mu$ characterizing the virtualities of quark and gluon propagators in the hard-scattering partonic process will be taken as $\mu^{2}=\langle x\rangle\,Q^{2}$ with $1/4\leq\langle x\rangle\leq 3/4$ (see Melic et al. (1999, 2002); Stefanis (1999); Mojaza et al. (2013); He et al. (2006); Lu et al. (2009); Wang and Shen (2016, 2017); Di Giustino et al. (2024) for an elaborate discussion on the renormalization/factorization scale setting).

In order to facilitate an in-depth exploration of the dynamical pattern dictating the pion form factor, we display explicitly in Figure 2 the obtained theory predictions for the leading-power hard-gluon-exchange contributions at leading-logarithmic (LL), NLL, and NNLL accuracy accomplished in this Letter as well as the subleading power corrections arising from I) an atypical “end-point” configuration of twist-two in the one-loop approximation and II) the higher-twist pion distribution amplitudes up to the twist-six accuracy with the dispersion technique Braun et al. (2000), by adopting model I of the pion distribution amplitude $\phi_{\pi}$ as our default choice. Inspecting the distinctive feature for a variety of higher-order corrections in Figure 2 implies that the newly determined two-loop QCD correction to the hard-scattering contribution based upon the hard-collinear factorization theorem can substantially enhance the corresponding NLL prediction of the pion form factor at intermediate and large momentum transfers: numerically at the level of $(30-50)\%$. In analogy to the QCD anatomy of the radiative leptonic $\bar{B}\to\gamma\ell\bar{\nu}_{\ell}$ decay with an energetic photon Braun and Khodjamirian (2013); Wang (2016); Wang and Shen (2018); Beneke et al. (2018); Khodjamirian et al. (2024), the soft non-factorizable correction to the charged pion form factor can constantly shift the NNLL resummation improved leading-power contribution by approximately an amount of ${\cal O}(25\%)$ in the kinematic domain $5.0\,{\rm GeV^{2}}\leq Q^{2}\leq 20.0\,{\rm GeV^{2}}$. We are then led to conclude that an ironclad and fully analytical extraction of the two-loop short-distance coefficient function $T_{1}^{(2)}$ in the factorized expression (3) is evidently vital for obtaining the robust and accurate theory prediction of the flagship hadron form factor $F_{\pi}$ and for advancing further the QCD factorization programme targeting at the high-precision computation of hard exclusive reactions. It remains interesting to observe that the very inclusion of the NNLO QCD correction to the pion EMFF turns out to be especially beneficial for better accommodating the benchmark lattice simulation results at intermediate momentum transfers Ding et al. (2024).

We now turn to present in Figure 3 our final theory predictions for the pion form factor $F_{\pi}$ with the four different phenomenological models of the twist-two pion LCDA discussed above, confronting with the two experimental measurements Bebek 78 Bebek et al. (1978) and JLab 08 Huber et al. (2008) at four distinct kinematic points of $Q^{2}=\{2.45,\,3.33,\,6.30,\,9.77\}\,{\rm GeV^{2}}$. It can then be observed that the yielding numerical prediction with the anti-de Sitter-QCD inspired model fulfilling an additional constraint from the lattice determination of $a_{2}(\mu_{0})$ Bali et al. (2019) (i.e., our model I) provides us with the most optimized description of the available experimental data points and the model-independent lattice QCD results Ding et al. (2024) simultaneously. The extraordinary snapshot of the well-separated uncertainty bands for the first three sample models due to the variations of the renormalization/factorization scales in the preferred intervals is particularly encouraging to acquire new insights on the intricate behaviour of the leading-twist pion distribution amplitude, in combination with the envisaged precision measurements at JLab and EIC. In contrast with the NNLL QCD computation for the photon-pion transition form factor Gao et al. (2022a); Braun et al. (2021), our two-loop prediction of the pion EMFF based upon the hard-collinear factorization prescription approaches the desired scaling behaviour in the formal $Q^{2}\rightarrow\infty$ limit rather slowly. This intriguing pattern can be attributed to the fact that the determined hierarchy between the subleading conformal spin effect and the counterpart asymptotic contribution grows steadily with the increasing loop order $\ell$ at realistic momentum transfers accessible in the current and forthcoming experimental facilities, thus postponing the onset of the asymptotic regime to an enormously higher value of $Q^{2}$ (see Braun et al. (2002, 2006); Anikin et al. (2013); Huang et al. (2024) for further discussions in the context of the nucleon electromagnetic form factors).

We finally address the genuine benefit of the full two-loop QCD calculation of the charged pion form factor on the model-independent extraction of the shape parameters dictating the twist-two pion distribution amplitude on the light-cone. To this end, we proceed to explore the intrinsic sensitivity of the exclusive $\pi^{+}\,\gamma^{\ast}\to\pi^{+}$ form factor on the currently poorly constrained conformal coefficient $a_{4}(\mu_{0})$ by combining the achieved NNLL prediction of the leading-twist contribution with the previously determined subleading power corrections in the $\Lambda_{\rm QCD}^{2}/Q^{2}$ expansion. It is evident from Figure 4 that taking into account the newly obtained two-loop QCD correction to the pion EMFF form factor will indeed be advantageous to enhance the sensitivity of extracting the key shape parameter $a_{4}(\mu_{0})$ in anticipation of the prospective precision EIC measurements Abdul Khalek et al. (2022), complementary to an alternative strategy on the basis of the factorization analysis of the photon-pion form factor Gao et al. (2022a).

In conclusion, we have endeavored to accomplish for the first time the rigorous two-loop QCD computation of the pion form factor in an analytical fashion, by applying the modern effective field theory formalism that enables us to implement the UV renormalization and IR subtractions of evanescent operators systematically. Crucially, we have demonstrated that the thus determined NNLO QCD correction to the short-distance coefficient function can bring about an enormous impact on the pion form factor over a wide range of momentum transfers, by employing the four phenomenologically acceptable models of the leading-twist pion distribution amplitude. In particular, the very inclusion of the two-loop radiative correction allowed for an improved extraction of the essential shape parameters dictating the intricate profile of the pion distribution amplitude, when confronted with the encouraging measurements from the upcoming EIC experiment. Extending and developing further our factorization prescription to the charged kaon form factor and to the more challenging $B_{c}\to\eta_{c}\,\ell\,\bar{\nu}_{\ell}$ transition form factors Böer (2018); Bell and Feldmann (2007); Bell (2006); Böer et al. (2024) will be highly beneficial for deepening our understanding towards the diverse facets of the strong interaction dynamics encoded in the hard-scattering processes.

It is our pleasure to thank Yong-Kang Huang for illuminating discussions. The research of Y.J. was supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Sino-German Collaborative Research Center TRR110 “Symmetries and the Emergence of Structure in QCD” (DFG Project-ID 196253076, NSFC Grant No. 12070131001, - TRR 110). J.W. is supported by the National Natural Science Foundation of China with Grant No. 12005117, No. 12321005, and No. 12375076, and the Taishan Scholar Foundation of Shandong province with Grant No. TSQN201909011. The work of Y.F.W. is supported in part by the National Natural Science Foundation of China with Grant No. 12405117. Y.M.W. acknowledges support from the National Natural Science Foundation of China with Grant No. 12075125 and No. 2475097.