$Z$ boson radiative decays to a $S$-wave quarkonium at NNLO and NLL accuracy

Within the framework of nonrelativistic QCD (NRQCD) factorization formalism, we compute QCD next-to-next-to-leading order (NNLO) corrections to the helicity amplitudes as well as the decay width of $Z\to H+\gamma$, where $H$ can be $\eta_Q (Q=c,b), J/\psi$, or $\Upsilon$. In addition, we resum the next-to-leading logarithms (NLL) of ${m_Z^2}/{m_Q^2}$ to all orders of $\alpha_s$ for the leading-twist helicity amplitude by employing the light-cone factorization approach. It is worth mentioning that we obtain the analytic expressions of the truncated NLL at $\alpha_s^2$. We find that the $\mathcal{O}(\alpha_s)$ corrections are around 10\% for $\eta_c$ and $\Upsilon$ productions, however insignificant for $J/\psi$ and $\eta_b$ productions. The $\mathcal{O}(\alpha_s^2)$ corrections are moderate for charmonium production, while very small for bottomonium production. Moreover, it is found that the NLL resummation can considerably alter the NRQCD prediction, especially for $J/\psi$ production. Combining the NRQCD and light-cone computation, we make phenomenological predictions on the decay widths and branching fractions. In addition, we investigate the dependence of the theoretical results on the heavy quark mass, and find the branching fraction of $Z\to H+\gamma$ monotonically decreases as $m_Q$ increases.


I. INTRODUCTION
It is an ideal platform to study the interplay between perturbative and nonperturbative QCD through the radiative decay of the Z boson to a quarkonium.To date, experimentalists have made many attempts to search for such processes [1][2][3], yet have failed to find any signals.In recent years, several high-luminosity lepton colliders, such as ILC [4], FCC-ee [5] and CEPC [6], have been proposed to run at the Z pole mass for a period of time.It will undoubtedly provide an opportunity to accumulate a large number of Z bosons, thus increasing the chances to probe rare decay processes.
The exclusive processes Z → quarkonium + γ have been extensively studied on the theoretical side, with the earliest computations dating back to the 1980s [7].In Ref. [8], Luchinsky studied these processes at lowest order in α s and v 2 in both the nonrelativistic QCD (NRQCD) [9] and light-cone (LC) factorization formalisms [10,11], where v represents the typical velocity of the heavy quark in the quarkonium rest frame.In Ref. [12], Wang et al. obtained the analytic expressions of the amplitudes for Z → quarkonium + γ in the leading-power LC approximation at next-to-leading order (NLO) in α s .Furthermore, Huang et al. presented calculations of the rates for Z → V + γ accurate up to the leading-power LC approximation at NLO both in α s and v [13], where V denotes a vector quarkonium.Shortly afterwards, the resummation of the leading logarithms (LL) of m 2 Z /m 2 Q for the decay width of Z → V + γ was carried out [14].Bodwin et al. further considered the resummation of logarithms of m 2 Z /m 2 Q for the O(α s ) corrections as well as the O(v 2 ) corrections [15].A combination of the next-to-leading logarithms (NLL) resummation and NLO fixed-order results was carried out for η Q + γ production in Ref. [16].The decay widths for Z → Υ(nS) + γ have been calculated up to NLO in α s based on the NRQCD, which are proposed to determine the Zb b coupling [17].Very recently, the decay widths of Z boson radiative decays to a P -wave quarkonium have been computed accurately up to next-to-next-toleading order (NNLO) in α s based on the NRQCD and LL resummation based on the LC approach [18].Moreover, the cross sections of e + e − → quarkonium + γ at Z factories have been computed in Refs.[19][20][21].It is worth mentioning that some efforts toward the NNLO perturbative corrections to e + e − → quarkonium + γ at B factories have been made in Refs.[22][23][24] in recent years.
In this work, we investigate the radiative decay of the Z boson to a quarkonium H (H can be η Q , J/ψ or Υ) by including both the NNLO perturbative corrections and the NLL resummation.We first compute the helicity amplitudes at NNLO in α s and leading order (LO) in v within the framework of NRQCD.To reduce the ambiguity in choosing the energy scale and uncertainty from the higher-order corrections arising from the large logarithms of m 2 Z /m 2 Q , we employ the LC formalism [25] to refactorize the NRQCD shortdistance coefficients (SDCs) and utilize the celebrated Efremov-Radyushkin-Brodsky-Lepage (ERBL) equation [10,26] to resum the large logarithms of m 2 Z /m 2 Q .Concretely, we will perform the NLL resummation for the leading-twist helicity amplitudes, i.e., resuming both the ) to all orders of α s .The paper is organized as follows.In Sec.II, we present the theoretical framework to compute the decay widths of Z → H + γ.In Sec.III, we employ the NRQCD formalism to factorize the helicity amplitudes, introduce the procedure and techniques to compute the SDCs, and present the results of the helicity SDCs at various perturbative levels.Section IV is devoted to the LC factorization for the leading-twist helicity SDCs.In addition, resummation of the NLL is formulated and explicitly carried out.The analytic expressions of the truncated NLL at α 2 s are also obtained.A detailed phenomenological analysis is performed in Sec.V. Finally, we summarize in Sec.VI.In Appendix A, we construct the helicity projectors.In Appendix B, the explicit expressions for the Brodsky-Lepage (BL) kernels are given.In Appendix C, we present some useful convolution formulas.

II. THEORETICAL FRAMEWORK FOR DECAY WIDTH
Applying the helicity amplitude formalism to analyze the hard exclusive production process proves to be convenient.The unpolarized decay widths of Z → H + γ can be expressed in terms of helicity amplitudes where |P| denotes the magnitude of the H spatial momentum: where m H refers to the mass of the quarkonium H, and the Källen function is defined via represents the helicity amplitude of Z → H(λ 1 ) + γ(λ 2 ) with λ 1 and λ 2 being the helicities of the H and outgoing photon respectively.To deduce (1), we have applied the parity invariance [27] to relate different helicity amplitude (3) Obviously, there are one independent helicity amplitude for η Q production, and two for J/ψ or Υ.
In the limit of m Q ≪ m Z , the helicity amplitude A H λ 1 ,λ 2 satisfies the asymptotic behavior where r = m Q /m Z .In (4), one power of r originates from the large momentum transfer which are required for the heavy-quark pair to form the heavy quarkonium with small relative momentum, and the other powers arise from the helicity selection rule in perturbative QCD [28,29].
To obtain the decay width, it is crucial to work out each helicity amplitude, which is the chief task of this work.Z boson interacts with quark-antiquark pair through the tree-level weak interaction as where g is the weak coupling in SU(2) L × U(1) Y electroweak gauge theory, g V = 1 − 8s 2 W /3 and g A = 1 for the up-type quark, and g V = −1 + 4s 2 W /3 and g A = −1 for the down-type quark.Here we have defined s W ≡ sin θ W , and c W ≡ cos θ W , where θ W signifies the Weinberg angle.
The Z boson can decay to η Q + γ through the vectorial interaction, while decay to J/ψ(Υ) + γ through the axial-vectorial interaction.For simplicity, it is convenient to explicitly extract the electroweak coupling from the helicity amplitudes:

III. NRQCD COMPUTATION A. The NRQCD factorization
According to the NRQCD factorization formalism [9], the helicity amplitude A H λ 1 ,λ 2 can be factorized into The nonrelativistically normalized long-distance matrix elements (LDMEs) are for η Q , and for J/ψ or Υ, where ψ † and χ denote the Pauli spinor fields creating a heavy quark and antiquark in NRQCD respectively, and ǫ H represents the polarization vector of J/ψ or Υ.The dimensionless SDC C H λ 1 ,λ 2 , signifying the perturbative contribution, can be evaluated either through the standard matching procedure or the method of region [30].In this work, we use the latter to compute these SDCs.The asymptotic behavior of the helicity SDCs can be straightforwardly deduced from (4) and ( 7) where µ R and µ Λ indicate the renormalization scale and factorization scale respectively, β 0 = (11/3)C A −(4/3)T F n f is the one-loop coefficient of the QCD β function, where n f is the number of active quark flavors.The explicit ln µ 2 R term is deduced from the renormalizationgroup invariance.γ H represents the anomalous dimension associated with the NRQCD bilinear currents carrying the quantum number 1 S 0 or 3 S 1 [31]: The occurrence of ln µ 2 Λ is demanded by the NRQCD factorization.Note that the γ H ln µ 2 Λ terms in (11) exactly cancel the µ Λ dependence of the NRQCD matrix element, so that the helicity amplitudes/decay widths are independent of µ Λ .For convenience, we have classified the two-loop Feynman diagrams into two groups, the 'regular' and the 'nonregular'.Some representative Feynman diagrams are illustrated in Fig. 1.Correspondingly, C H,( 2) 11) represent the contributions from the 'regular' part and 'nonregular' part respectively.
We briefly outline the calculation.We begin with the quark-level process for Z → Q Q + γ.The package FeynArts [32] is employed to generate the Feynman diagrams and the corresponding Feynman amplitudes through order O(α 2 s ).We utilize the spin projectors to enforce Q Q in 1 S 0 for η Q and 3 S 1 for J/ψ or Υ, and employ the packages FeynCalc [33] and FormLink [34] to obtain the hadron-level amplitudes order by order in α s .The helicity amplitudes are evaluated with the aid of the helicity projectors, which are constructed in Appendix.A.
It is well known that, in dimensional regularization, the anticommutation relation {γ µ , γ 5 } and the cyclicity of Dirac trace cannot be satisfied simultaneously.In practical computation, the naive-γ 5 scheme [35], which keeps the anticommutation relation {γ µ , γ 5 }, is frequently applied.In this scheme, spurious anomaly, which spoils chiral symmetry and hence gauge invariance, can be avoided.Because of the lack of the cyclicity of the trace, one must fix a reading point for a fermion loop with an odd number of γ 5 .In this work, we will select the vertex of the Z boson as the reading point for the Feynman diagram in Fig. 1(d), and the final state quarkonium as the reading point for the other Feynman diagrams.For more details, we refer the readers to Ref. [18].
At lowest order in v, we neglect the relative momentum in Q Q pair prior to carrying out the loop integration, which amounts to directly extracting the SDCs from the hard loop region [30].With the aid of the packages Apart [36] and FIRE [37], we can reduce the loop integrals into linear combinations of master integrals (MIs).Finally, we end up with 6 one-loop MIs, and around 320 two-loop MIs, which are evaluated by the powerful package AMFlow [38][39][40][41][42].
After implementing the on-shell renormalization scheme for the heavy quark mass and field strength [43,44], and the MS renormalization scheme for the QCD coupling, the UV pole is exactly eliminated, while an uncanceled single IR pole still remains, which can be factored into the NRQCD LDME, so that the SDC becomes IR finite.
It is straightforward to obtain the LO helicity SDCs: The analytic expressions of C H,(1) λ 1 ,λ 2 can also be readily obtained.Instead of presenting the cumbersome expressions, we list their asymptotic expansions in the limit of r → 0: where the real part of (14a) is consistent with that in Refs.[12,45], and the expression of (14c) is consistent with that in Ref. [12].It becomes much more challenging to deduce the analytical expressions for all the encountered two-loop MIs.In this work, we are content with high-precision numerical results.To perform the numerical computation, we take m Z = 91.1876GeV from the particle data group (PDG) [46], and the charm quark and bottom quark pole masses to be m c = 1.69 GeV and m b = 4.80 GeV, which are converted from the MS masses m c (m c ) = 1.28 GeV and m b (m b ) = 4.18 GeV [46] at two-loop level by use of the package RunDec [47].The numerical values of the various helicity SDCs are tabulated in Table I.For reference, we explicitly keep the n L , n c , and n b dependence for the SDCs, where n L denotes the number of light quark flavors, and n c = 1 and n b = 1 signify the numbers of the charm and bottom quarks respectively.The dependence of the theoretical results on the heavy quark mass will be investigated in Sec.V.
where g u V and g d V correspond to the values of g V for up-type quark and down-type quark respectively.Besides the NRQCD factorization formalism, we can also employ the LC factorization framework to calculate the decay amplitude for Z → H +γ at the leading twist.By following the spirit of Ref. [25], the LC factorization formula for the SDC is written as where the hard-kernel T H and the leading-twist LC distribution amplitude (LCDA) φH are perturbatively calculable around the scale m Z and m Q , respectively.Up to O(α s ), we have the expansions The explicit hard-kernels and LCDAs up to O(α s ) are given in [12] as1 T (0) T φ(0) φ( 1) for H = J/ψ and Υ.Here x = 1 − x and ∆ = 0 for the NDR scheme and ∆ = 1 for the HV scheme.Note that the terms proportional to δ(x − 1/2) in φ(1) H actually contribute to the one-loop corrections to decay constants f P,V .

B. The NLL resummation with the ERBL equation
The leading twist LCDAs of quarkonia obey the celebrated ERBL equation [26,48] with the BL kernel expanded in α s H (x, y) + . . . .
The BL kernel at the lowest order of α s given in [26,48] for a pseudoscalar meson is the same as that for a longitudinally polarized vector meson while their BL kernels at O(α 2 s ) given in [49][50][51][52][53] differ from each other by [54] V where ∆ = 0 for the NDR scheme and ∆ = 1 for the HV scheme.The explicit expressions of these BL kernels are summarized in Appendix B.
The formal solution can be where the evolution kernel is with P standing for the ordering on α s which can be expanded in form of the Dyson series and the " * " denotes the appropriate convolution over the light-fractions.Therefore, we have the renormalization group improved SDCs for Z → η Q + γ and Z → J/ψ(Υ) + γ at the leading power of expansion in NRQCD factorization as 1.The truncated NLL resumed SDCs The explicit truncation of the perturbative expansion of Here we have used the RG evolution of the strong coupling α s (µ) at the NLL level and its perturbative expansion With the useful convolutions listed in Appendix C, we have the truncation of the NLL resumed SDCs up to O(α 2 s ) explicitly for H = η Q , and for H = J/Ψ(Υ).Note that the γ 5 -scheme dependences of the hard-kernel T P and LCDA φP cancel with each other eventually in the above expansion of K P,NLL as it should.

The NLL resummation to all orders with the Gegenbauer polynomial expansion
The resummation of the LL to all orders of α s are commonly done with the assists of the Gegenbauer polynomials by noticing such polynomials are the eigenfuctions of the BL kernel at the lowest order of α s .The NLL resummation to all orders of α s can be also done in a similar way by considering the nondiagonal part in V (1) with C (3/2) n (x) being the order 3/2 Gegenbauer polynomials, and the Gegenbauer moments φH,n (µ) = 4(2n + 3) (n + 1)(n + 2) φH,n (µ) can have the similar perturbative expansion in α s around the scale m Q as φH (x, m Q , µ).
Solving the ERBL equation in the Gegenbauer moments space, we have formally in which the matrix-elements of the NLL evolution kernel in the Gegenbauer moments space U H n,k are given in Ref. [55].
Similarly, we have the Gegenbauer expansion for the hard-kernels with T H,n (µ) can have the similar perturbative expansion in α s around the scale m Z as T H (x, m Z , µ).Hence, we can get Specifically, we have where It is worth mentioning that K P is indeed independent of the γ 5 -scheme used in calculations of the hard-kernel T P and LCDA φP as well as V P , as the authors of Ref. [16] indicated.

C. The SDCs by combining the NRQCD prediction and the LC resummation
We proceed to compute the leading-twist SDCs by combining the NRQCD prediction and the LC resummation.To avoid double counting, one should subtract α n s ln n r terms from the NRQCD prediction for the LL resummation, and subtract both the α n s ln n r and α n+1 s ln n r terms for the NLL resummation.It is convenient to introduce the following symbols where K H,NLL can be found in ( 30) and ( 31), and K H,LL signify the sum of α n s ln n r terms in K H,NLL .Thus, we formally have  II.Note that the strong coupling constant α s (m Q ) is evaluated through the running formula (28).To accelerate the convergence, we employ the Abel-Padé approach [56] to sum the series in (39).From Table II, we find that the LL resummation can significantly improve the LO NRQCD prediction, however only slightly alter the higher order predictions.It can be explained by that some dominant α n s ln n r contributions have already been included in the higher order NRQCD SDCs, i.e., the O(α s ln r) contribution has been included in the O(α s ) NRQCD SDC, and the O(α 2 s ln 2 r) contribution have been included in the O(α 2 s ) NRQCD SDC.To be contrary, the effect of the NLL resummation does not become weaker as the perturbative order increases, e.g., the difference between 'NNLO+NLL' and 'NNLO' is as large as that between 'NLO+NLL' and 'NLO' .Furthermore, unlike the case for LL resummation, we find the difference between the NLL resummation K H and its truncated expansion K H,NLL actually does not decrease from one-loop to two-loop accuracy.The crucial reason is that we evaluate the value of α s (m Q ) in (38) by the running formula (28), while obtain its truncated expression K H,NLL by employing (29).The values of α s (m Q ) are quite different when using (28) and (29), for example α s (m c ) = 0.297 by (28) and α s (m c ) = 0.228 by (29). 2 By taking m c = 1.69 GeV, m Z = 91.1876GeV, and α s (m Z ) = 0.1181, we obtain the values of K J/ψ(0,0) , K J/ψ(1,0) and K J/ψ(0,1) to be 1.217, −10.555 − 4.227i and −9.673, respectively.As a comparison, if expanding K J/ψ in powers of α s , and truncating K J/ψ(0,0) to O(α 2 s ), K J/ψ(1,0) and K J/ψ(0,1) to O(α s ), we obtain the values of the three quantities to be 1.203, −10.658 − 4.909i, and −8.657 respectively.
Since contribution from the NLL resummation is considerable even at two-loop order, particularly for J/ψ production, the NLL resummation is important to improve the theoretical predictions.

V. PHENOMENOLOGY
In phenomenological analysis, we take s 2 W = 0.231, m Z = 91.1876GeV, the total decay width of the Z boson Γ Z = 2.4952 GeV [46], and fix the running QED coupling α(m Z ) = 1/128.943[57].The default value of µ R is chosen µ R = m Z / √ 2 and we have varied µ R from m Z /2 to m Z to estimate the theoretical uncertainties in computing the NNLO perturbative corrections.In addition, we approximate the NRQCD LDMEs at µ Λ = 1 GeV by the Schrödinger radial wave function at the origin where the radial wave functions at the origin are evaluated from Buchmüller-Tye (BT) potential model [58].By taking the heavy quark pole masses m c = 1.69 GeV and m b = 4.80 GeV, we tabulate the unpolarized decay widths and branching fractions for Z → H + γ at various levels of accuracy in Table III.Since the axial-vectorial interaction of Zcc is roughly three times larger than the vectorial interaction, the branching fraction for Z → J/ψ + γ is much larger than that for Z → η c + γ.In addition, the O(α s ) corrections are negligible for J/ψ and η b production, while can reach 10% for the other two channels.The O(α 2 s ) corrections are moderate for the charmonium production, however are small for the bottomonium production.Moreover, we find the NLL resummation can considerably alter the NRQCD predictions, particularly for Z → J/ψ + γ.It is worth noting that the uncertainty from the renormalization scale µ R is inconsiderable.To be honest, we must emphasize that different choice of the values of the LDMEs [58][59][60][61][62][63] may largely affect the theoretical predictions.
It is intriguing to compare our theoretical prediction on the branching fraction with that from the LC models in literature.In Table III, we list the LC predictions from Ref. [8] and Ref. [14] in the fifth column and sixth column respectively.In Ref. [8], the authors take the same LCDA for η c and J/ψ, which was obtained at µ = m c in Ref. [64] inspired by the QCD sum rule.In Ref. [14], the authors assumed the LCDA for J/ψ and Υ at µ = 1 GeV to be of Gaussian form.Different LCDA corresponds to the different internal quark motion assumption.In spite of having very different nonperturbative inputs, we find our 'NNLO+NLL' prediction for η c production is consistent with the result from Ref. [8], and our prediction for Υ production roughly agrees with the result from Ref. [14].On the other hand, our prediction for J/ψ production is a bit smaller than the values from both references.The distinct LCDA input and the considerable NLL resummation effect mainly account for the difference.There is no doubt that the ambiguity in choosing LCDA can cause large uncertainties.
We proceed to investigate the dependence of the theoretical results on the heavy quark mass.In Fig. 2, we plot the branching fraction for Z → H + γ as a function of m Q at various levels of accuracy.We find that the branching fraction monotonically decreases as m Q increases.[8] and Ref. [14] are also listed in the fifth column and sixth column respectively.The two uncertainties in the fifth column arise from uncertainties of the leptonic decay constant of charmonium and the LCDA parameters respectively.The three uncertainties in the sixth column originate from the factorization scale dependence, the quarkonium decay constants, and the LCDA parameters, respectively.

Channel
Order Γ total (eV) Br(×10 −9 ) Br(×10 Finally, we estimate signal events at the future super Z factories, such as the Z-factory mode in CEPC, where the Z boson yield can reach 7 × 10 11 [6].It is expected that there will be around 5 × 10 4 J/ψ or Υ and 10 4 η Q events produced through Z → H + γ.The J/ψ or Υ can be reconstructed through their leptonic decays, thus providing several thousands of γℓℓ events.It is promising to search these two channels.However, due to the lack of a clean decay mode for η Q , the experimental measurement of Z → η Q + γ would be quite challenging.

VI. SUMMARY
In summary, we have computed the O(α s ) and O(α 2 s ) corrections to the helicity amplitudes and decay widths for Z → H + γ by applying the NRQCD factorization approach.It is found that the O(α s ) corrections are moderate for η c and Υ productions, however tiny for J/ψ and η b productions.The O(α 2 s ) corrections are considerable for charmonium pro- duction, while small for bottomonium production.In addition, we find that the branching fraction for J/ψ production is much larger than that for η c production.The NLL of m 2 Z /m 2 Q in the leading-twist SDCs are resummed to all orders of α s by employing the celebrated ERBL equation.We find that the NLL resummation can considerably alter the NRQCD prediction, particularly for J/ψ production, so that the NLL resummation is important to improve the theoretical prediction.We also find the branch fraction for Z → H + γ monotonically decreases as m Q increases.In addition, we compare our predictions on the branching fractions with these from the LC models in literature.We find our prediction for η c production is consistent with Ref. [8], our prediction for Υ production is slightly smaller than Ref. [14], and our prediction for J/ψ is a bit smaller than Ref. [14].
It is expected that there will be around 5 × 10 4 J/ψ or Υ events produced through Z boson radiative decay at the future super Z factories.Therefore it seems that the observation prospects of Z → H + γ are promising in the future.at the level of NLLs, where ∆ = 0 for the NDR scheme, and ∆ = 1 for the HV scheme.One can see in the truncated expansion of the resumed SDC in Eq. ( 30), the ∆-dependence vanishes at each order of α s as it should.The corresponding convolutions for H = J/Ψ(Υ) can be obtained similarly.

- 9 )FIG. 2 :
FIG. 2: Branching fraction for Z → H + γ as a function of m Q .The green band denotes the uncertainty from µ R .

TABLE I :
NRQCD predictions to the various helicity SDCs.For simplicity, we define the symbols

TABLE III :
Unpolarized decay widths and branching fractions for Z → H + γ at various levels of accuracy.For comparison, the LC predictions from Ref.