Leading-twist parton distribution amplitudes of S-wave heavy-quarkonia

The leading-twist parton distribution amplitudes (PDAs) of ground-state $^1S_0$ and $^3S_1$ $c\bar c$- and $b\bar b$-quarkonia are calculated using a symmetry-preserving continuum treatment of the meson bound-state problem which unifies the properties of these heavy-quark systems with those of light-quark bound-states, including QCD's Goldstone modes. Analysing the evolution of $^1S_0$ and $^3S_1$ PDAs with current-quark mass, $\hat m_q$, increasing away from the chiral limit, it is found that in all cases there is a value of $\hat m_q$ for which the PDA matches the asymptotic form appropriate to QCD's conformal limit and hence is insensitive to changes in renormalisation scale, $\zeta$. This mass lies just above that associated with the $s$-quark. At current-quark masses associated with heavy-quarkonia, on the other hand, the PDAs are piecewise convex-concave-convex. They are much narrower than the asymptotic distribution on a large domain of $\zeta$; but nonetheless deviate noticeably from $\varphi_{Q\bar Q}(x) = \delta(x-1/2)$, which is the result in the static-quark limit. There are also material differences between $^1S_0$ and $^3S_1$ PDAs, and between the PDAs for different vector-meson polarisations, which vanish slowly with increasing $\zeta$. An analysis of moments of the root-mean-square relative-velocity, $\langle v^{2m}\rangle$, in $^1S_0$ and $^3S_1$ systems reveals that $\langle v^4\rangle$-contributions may be needed in order to obtain a reliable estimate of matrix elements using such an expansion, especially for processes involving heavy pseudoscalar quarkonia.


Introduction.
In studying hard exclusive processes within the Standard Model there are many instances in which one may appeal to factorisation theorems so that, at leading-order in a systematic expansion, the amplitude involved can be written as a convolution of a hard-scattering kernel, calculable in perturbation theory, and the so-called leading-twist PDA of the hadron involved, ϕ(x), where x is the light-front fraction of the hadron's total momentum carried by the struck parton.
For mesons, the PDA is a light-front projection of the system's Bethe-Salpeter wave-function onto the light-front. It is therefore process independent and hence plays a crucial role in explaining and understanding a wide range of a given meson's properties and interactions. The PDA is also essentially nonperturbative, i.e. it cannot be calculated using perturbation theory. The last two decades have witnessed significant progress toward the computation of realistic meson Bethe-Salpeter amplitudes [5][6][7][8]; and, critically, the last two years have seen the development of novel techniques which enable the reliable calculation of meson PDAs from such Bethe-Salpeter amplitudes [9,10]. Predictions are now available for the (leading) twist-Email addresses: corresponding authors -yxliu@pku.edu.cn (Yu-Xin Liu ), cdroberts@anl.gov (Craig D. Roberts) two PDAs of the π-, K-, ρ-and φ-mesons [11][12][13], and for twist-three pion and kaon PDAs [14,15].
Given that the last fifteen years have seen a dramatic expansion of interest in heavy-quark systems, owing to advances in both theoretical methods, and experimental activity and discoveries [16,17], it is an opportune moment to use the continuum approach indicated above in order to compute the twist-two PDAs of S -wave heavy-quarkonia. The theoretical interest is plain: one thereby arrives at a unified description and explanation of the leading-twist PDAs for almost all empirically accessible pseudoscalar and vector mesons. This means, e.g. that one simultaneously obtains an understanding of the structure of QCD's Goldstone modes and the η b -meson, and can track the structural rearrangements which take place as a growth in current-quark mass drives an evolution between them. There is also a phenomenological imperative: heavy-quarkonia PDAs appear in the analysis of numerous hard exclusive processes, e.g. quarkonia production at high-energies [18,19]; J/Ψ + η c pair production in e + e − annihilation [19,20]; B c → η c transitions [21]; decays of heavy S -wave quarkonia into lighter vector mesons [22]; deeply virtual quarkonia production, which can be used to probe the gluon distribution in the proton [23]; and Higgs boson decays into quarkonia [24].
Notably, it is often supposed that such and kindred processes may be treated accurately using non-relativistic QCD (NRQCD) [25], wherewith the leading-twist quarkonia PDAs are approximated as ϕ QQ (x, ζ) = δ(x − 1/2), where ζ is the scale of the momentum transfer involved, and effects of nonzero Q-Q relative velocity, v 2 > 0, are treated perturbatively. However, whereas this might be true for processes that only involve bquarks, corrections as large as a factor of two or more have been found when c-quarks are involved [26][27][28][29]. A calculation of ϕ QQ (x, ζ) in a framework that is capable of unifying this PDA with those of light-mesons can therefore also serve as valuable check on the fidelity of the NRQCD approximation. The Dyson-Schwinger equations (DSEs) [5][6][7][8] have this feature and we use that framework herein.

Parton distribution amplitudes.
In calculating the leadingtwist PDAs of heavy quarkonia, we follow Refs. [9,13] and consider projections of their Bethe-Salpeter wave functions onto the light-front: where: P, V denote, respectively, pseudoscalar ( 1 S 0 ) and vector ( 3 S 1 ) quarkonia; the trace is over color and spinor indices; Λ dq is a Poincaré-invariant regularisation of the four-dimensional integral, with Λ the ultraviolet regularisation mass-scale; Z 2,T (ζ, Λ) are, respectively, the renormalisation constants for the quark wave-function and the tensor vertex, which we compute in the chiral limit at the renormalisation scale ζ; n is a lightlike fourvector; P is the meson's four-momentum, with P 2 = −m 2 P,V , and n · P = −m P,V , with m P,V being the meson's mass.
In writing Eqs. (1) we have used the fact that there are only two independent vector-meson PDAs at leading-twist [30]: ϕ V (x), ϕ ⊥ V (x) describe, respectively, the light-front fraction of the meson's total momentum carried by the quark in a longitudinally or transversely polarised bound-state. We have also adopted the convention 1 0 dx ϕ(x) = 1, so that f P , f V , f ⊥ V are decay constants. The first two are measurable; but, whilst the tensor couplings f T V are gauge-and Poincaré-invariant, they depend on the renormalisation scale: f T V (ζ) → 0 as ζ → ∞. (Further details are available in Appendix A of Ref. [13].) The Bethe-Salpeter wave functions in Eqs. (1) can be written with Γ the relevant meson's Bethe-Salpeter amplitude and S the dressed propagator for the quark in that bound-state. We have defined q + = q + ηP, q − = q − (1 − η)P, η ∈ [0, 1]. Owing to Poincaré invariance, no observable can legitimately depend on η, i.e. the definition of the relative momentum. On the other hand, the choice η = 1/2 is computationally convenient when solving the Bethe-Salpeter equation that describes a bound-system constituted from equal-mass valence partons. It is appropriate to remark here that in order to produce gauge invariant results, Eqs. (1) should contain contributions that derive from a Wilson line, W[−xn/2, xn/2], drawn between the bound-state's valence constituents. The Wilson line vanishes in light-cone gauge and hence does not contribute when this choice is employed. On the other hand, light-cone gauge is seldom practicable in either model calculations or quantitative nonperturbative analyses in continuum QCD. Herein, as is typical in nonperturbative DSE studies, we employ Landau gauge because, inter alia [31][32][33]: it is a fixed point of the renormalisation group; that gauge for which sensitivity to model-dependent differences between Ansätze for the fermion-gauge-boson vertex are least noticeable; and a covariant gauge, which is readily implemented in numerical simulations of lattice-regularised QCD. It is therefore significant that W[−xz/2, xz/2] is not quantitatively important in the calculation of leading-twist PDAs [34].
With realistic meson Bethe-Salpeter amplitudes in hand, it is straightforward to follow Refs. [9,13] and obtain PDAs for heavy quarkonium bound-states from Eqs. (1). The first step is to compute these moments: In our Poincaré-covariant framework, arbitrarily many moments can be calculated, in principle and practice.
Having computed a sufficient number of the moments for a light-quark system, one could then reconstruct the associated PDA using the "Gegenbauer-α" procedure introduced in Refs. [9,35], which is ideal for representing the broad, concave amplitudes that are characteristic of such systems. For heavy quarkonia, on the other hand, one expects the PDAs to be piecewise convex-concave-convex on x ∈ [0, 1], as is typical of finite-width representations of δ(x − 1/2); and hence we proceed by assuming that can serve as an efficacious replacement for the expansion of heavy-quarkonia PDAs in terms of order-α Gegenbauer polynomials. This function should be viewed as an informed assessment of the likely prior-distribution, in the sense of a Bayesian analysis of the reconstruction problem.
At this point, with 2m max moments computed for a given heavy-quarkonium state using the appropriate formula in Eq. (3), one determines a in Eq. (4) by minimising Heavy quarkonia Bethe-Salpeter amplitudes. To continue with our calculation of S -wave heavy-quarkonia valence-quark PDAs, the dressed-quark propagators and Bethe-Salpeter amplitudes associated with these bound states are needed. We compute these quantities using gap and Bethe-Salpeter equation solutions obtained using the rainbow-ladder (RL) truncation of QCD's DSEs [5][6][7][8] and the interaction introduced in Ref. [36]. The RL truncation is the leading order in a systematic, symmetry-preserving procedure that enables a tractable formulation of the continuum bound-state problem [37,38]. It is widely used in hadron physics and known to be accurate for light-quark ground-state vector-and isospin-nonzeropseudoscalar-mesons [5][6][7][8], and properties of the nucleon and ∆-baryon [39][40][41][42], because corrections in these channels largely cancel owing to parameter-free preservation of relevant Ward-Green-Takahashi identities (WGTIs).
The RL truncation has also been explored in connection with heavy-light mesons and heavy-quarkonia [43][44][45][46][47][48]. Those studies reveal that beyond-RL corrections to the dressed-quarkgluon vertex and hence the Bethe-Salpeter kernel are critical in heavy-light systems; and an interaction strength for the RL kernel fitted to pion properties alone is not optimal in the treatment of heavy quarkonia. Both observations are readily understood; but we focus on the latter because it is relevant herein.
The interaction in Ref. [36] is deliberately consistent with that determined in studies of QCD's gauge sector, which indicate that the gluon propagator is a bounded, regular function of spacelike momenta, q 2 , that achieves its maximum value on this domain at q 2 = 0 [49][50][51][52], and the dressed-quark-gluon vertex does not possess any structure which can qualitatively alter these features [53,54]. It also preserves the one-loop renormalisation group behaviour of QCD so that, e.g. the quark massfunction is independent of the renormalisation point, and the infrared behaviour is determined by a single parameter, conventionally expressed as (ς G ) 3 := Dω. Computations [55][56][57] show that observable properties of light-quark ground-state vectorand isospin-nonzero pseudoscalar-mesons are practically insensitive to variations of ω ∈ [0.4, 0.6] GeV, so long as (The midpoint ω = 0.5 GeV is usually employed in calculations.) This feature also extends to numerous properties of the nucleon and ∆-baryon [58]. The value of ς G is chosen so as to obtain the measured value of the pion's leptonic decay constant, f π ; and in RL truncation this requires ς RL G = 0.87 GeV. Following Ref. [59], however, it has become possible to employ far more sophisticated kernels for the gap and Bethe-Salpeter equations, which overcome the weaknesses of RL truncation in all channels studied thus far. This new technique, too, is symmetry preserving; but it has an additional strength, i.e. the capacity to express dynamical chiral symmetry breaking (DCSB) nonperturbatively in the integral equations connected with bound-states. That is a crucial advance because DCSB is an important emergent phenomena within the Standard Model: it is the origin of more than 98% of the visible mass in the Universe [60]. Owing to this feature, the new scheme is described as the "DCSB-improved" or "DB" truncation. It preserves successes of the RL truncation; but has also enabled elucidation of many novel nonperturbative features of QCD [9,[61][62][63].
In a realistic DB truncation, ς DB G = 0.55 GeV; a value which coincides with that predicted by solutions of gauge-sector gap equations in QCD [64]. Since all dressing of the quark-gluon vertex vanishes in the heavy-quark limit, so that RL truncation must become valid, then the aforementioned agreement entails both that ς DB G should be the infrared mass-scale appropriate for the RL analysis of truly heavy-heavy systems and provide more realistic results in such treatments of empirically accessible heavy-quarkonia. Herein we use both ς RL G and ς DB G ; and, as will become clear from those comparisons with experiment which are possible, the expectations described here are confirmed.
With the kernels of the gap and Bethe-Salpeter equations specified, one can employ standard algorithms and obtain numerical solutions for the dressed-quark propagators and Bethe-Salpeter amplitudes we require. Those solutions yield the predictions for quarkonia static properties listed in Table 1 The leptonic decay constants in Table 1 computed using ς DB G differ with experiment by a 10% root-mean-square relativeerror (rms-re) and by the same amount when compared with the lQCD values, which themselves differ from experiment by 11%. Measured this way, they are similar to the DSE results in Refs. [45,46]. On the other hand, the rms-re between the ς RL G results and experiment is 43%. The CQM values listed in the last row emphasise the difficulty, understood by practitioners [73], that such approaches encounter when attempting to simultaneously describe light-and heavy-quarkonia in the absence of a veracious expression of relevant WGTIs [74,75]. (N.B. Our numerically determined Bethe-Salpeter wave-functions for heavy 1 S 0 quarkonia satisfy Eq. (18) in Ref. [74], a corollary of the axial-vector WGTI, with an accuracy of 98.9 ± 0.4%.) It is interesting to note that in the limit of infinitely-heavy quarks, when all other mass-scales can be neglected, one has [76][77][78]: m P = 2M Q , where one may choose M Q = M Q (0), i.e. the value of the appropriate dressed-quark mass-function at the origin, since this function no longer runs; and 2m Q (ζ) f V = m V f ⊥ V (ζ). Deviations from these predictions are one crude measure of the impact of essentially dynamical effects in heavyquarkonia. Working from Table 1, ς RL G results differ from these expectations with a 7% rms-re, whilst the ς DB G rms-re is 4%.     (5) can now be used to compute the twist-two PDAs. For systems composed of light-quarks, i.e. those with current masses 0.1 GeV, which is roughly the s-quark mass, the (n · q + ) m factor in Eqs. (3) produces a highly-oscillatory integrand and thus reliable values for the moments cannot be obtained using a direct approach to computing the integrals. In these cases, the procedure described in Ref. [9], based on generalised spectral representations of the light-quark propagators and bound-state amplitudes, is necessary and efficacious. With increasing current-quark mass, however, owing to a damping influence from the large quark mass, this problem is shifted to progressively higher moments, which are also of diminishing magnitude and hence have little real impact, so that a "bruteforce" approach is feasible for heavier quarkonia. That is how we proceed with Eqs. (3)-(5) herein, viz. direct integration using interpolations of numerical solutions for the propagators and Bethe-Salpeter amplitudes. In order to eliminate dependence on the upper-bound of the momentum integration, which is a remnant of the oscillation problem just described, we introduced a factor 1/(1 + k 2 r 2 ) m for each moment ξ 2m ; computed the moment as a function of r; and extrapolated to r = 0. This procedure produced the values in Table 2.
It is evident from Table 2 that in heavy-quarkonium systems one obtains a reasonable nonzero signal for moments m ≤ m max = 4. Using these tabulated moments, Eq. (5) yields twist-two PDAs of the form in Eq. (4) with the widthparameters "a" in Table 3. In performing the least-squares fit we found ǫ RL I = 20 ± 14% and ǫ DB I = 16 ± 9%. These numbers can serve as an estimate of the errors in our moments and the as-   Table 3: Computed values of the width parameter "a" in Eq. (4), which characterises the twist-two PDA of each heavy-quarkonium system.  Table 3 are depicted in Fig. 1. Consistent with the pattern that has already been established, we consider the results obtained using ς DB G , depicted in the upper panel, to be the more realistic; and, notably, peakheights and widths in this case show a natural ordering: where "< N " means "'narrower than".
The PDAs obtained with ς RL G , drawn in the lower panel of Fig. 1, are anomalous in a number of ways, e.g. the -PDAs are unnaturally sharply peaked and hence ϕ J/Ψ < N ϕ Υ ⊥ . Such behaviour can be traced to an over-concentration of interactionstrength in the far-infrared when one requires a good description of light-meson observables using RL truncation [9][10][11][12]15]. This is corrected when DB kernels are employed [64]; and that improvement is mimicked in RL-studies of heavy-quarkonia which broaden the interaction by increasing ω in Eq. (6) [48].
Consideration of Fig. 1 reveals a curious interplay between PDA evolution with current-quark mass and renormalisation scale. Plainly, as Λ QCD /ζ → 0, where Λ QCD is QCD's fixed renormalisation-induced scale, all twist-two meson PDAs must approach ϕ asy (x). On the other hand, for any fixed ζ, ϕ QQ (x; ζ) → δ(x − 1/2) as Λ QCD /m Q (ζ) → 0. Given that each light-quarkonium PDA is a concave function of unit area which is broader than ϕ asy (x) and ERBL evolution [1][2][3] is smooth, uniform and area-preserving, then with increasing current-quark mass there must be a critical value, m c q (ζ), cor-4 responding to a particular value of the renormalisation-groupinvariant current-quark massm c q , at which ϕ qq (x; ζ) = ϕ asy (x). So long as the current-quark mass remains fixed bym c q , then this light-quarkonim PDA cannot change under ERBL evolution, viz. it is a fixed point. This occurs at the following masses: Whilst the actual values may change modestly in response to improvements of our framework, one may reliably conclude that the critical value typically lies just above the s-quark mass. The moments in Table 2 can be used to compute the rms relative velocity of valence-constituents in the quarkonium system under consideration. At leading-order of an expansion in this velocity [18,22]: v 2n = (2n + 1) ξ 2n . Using ς DB G , we find: Evidently, the generalised Gremm-Kapustin relation [79]: v 2n ≈ v 2 n , valid in nonrelativistic potential models, fails at every order for each system owing to the persistent x = 0, 1 endpoint tails of our computed PDAs. Despite this, Eq. (9) suggests that computation of physical observables using an expansion in rms relative velocity should converge reasonably quickly, but with v 4 -contributions that are sizeable for 1 S 0 systems. Figure 2: Comparison between our predictions for heavy-quarkonia PDAs and those reconstructed from moments computed using sum rules (SR) [18,19,21] or a light-front CQM (QM) [20]. Our predictions are the curves labelled η c , η b , (J/Ψ) ⊥ , (J/Ψ) .
Pointwise forms for the twist-two PDAs of cc-quarkonia have been estimated using sum rules [18,19,21] and a lightfront CQM [20]. Within their errors, those analyses report ϕ η c ≈ ϕ J/Ψ ⊥ ≈ ϕ J/Ψ . Our framework does not suffer from this drawback. There is little information on the PDAs of bbquarkonia; but a result for ϕ η b can be inferred from Ref. [21].
We have taken the PDA moments reported in Refs. [18][19][20][21] and reconstructed pointwise forms for the associated heavyquarkonia PDAs using the method described in connection with Eqs. (4), (5) herein. Figure 2 displays the resulting comparison. The form of ϕ η b determined from Ref. [21] is almost identical to our prediction. However, whilst the forms of ϕ η c ≈ ϕ J/Ψ ⊥ ≈ ϕ J/Ψ reconstructed from Refs. [20,21] are quite similar to each other, they are markedly different from our predictions, especially insofar as we find material differences between these three systems. (N.B. Although Refs. [18,19] used a function for ϕ I (ξ, a) in Eqs. (4) that is somewhat different from ours, in comparing the results we found only immaterial pointwise differences.)

Conclusion.
We computed the leading-twist PDAs of 1 S 0 and 3 S 1 cc-and bb-quarkonia using a framework that has already provided explicit forms for PDAs of light-quark mesons, and hence arrived at a unified picture of systems ranging from QCD's Goldstone modes, whose properties are greatly influenced by dynamical chiral symmetry breaking, to those in which explicit chiral symmetry breaking is the dominant effect.
In this connection we examined the evolution of meson PDAs with current-quark mass,m q , and found that the broad, concave twist-two PDAs of light-quark systems change smoothly with increasingm q . Consequently, there is always a value of m q =:m c q for which a given quarkonium PDA matches the asymptotic form, ϕ asy (x) = 6x(1 − x), and hence no longer evolves with renormalisation scale, ζ. This value ofm q lies just above that associated with the s-quark. Form q >m c q the PDAs are piecewise convex-concave-convex; but naturally evolve to the fixed point ϕ asy (x) as ζ is increased.
Heavy-quarkonia involvem q ≫m c q , in which case the PDAs are piecewise convex-concave-convex and much narrower than ϕ asy on a large domain of ζ. Nevertheless, for realistic values ofm q the PDAs deviate noticeably from ϕ QQ (x) = δ(x − 1/2), which is the limiting form in the case of static quarks. There are also material differences between 1 S 0 and 3 S 1 PDAs and between the PDAs for different vector-meson polarisations, although these vanish slowly with increasing ζ. Considering these features, we computed moments of the rms relativevelocity, v 2m , in 1 S 0 and 3 S 1 systems and found that in order to obtain accurate estimates for observables using an expansion in these terms it is likely that v 4 -corrections will require calculation, especially for processes involving 1 S 0 systems.