$Z$ boson radiative decays to a $P$-wave quarkonium at NNLO and LL accuracy

In this work, we study the radiative decay of $Z$ boson to a $P$-wave quarkonium $H$ in association with a photon, where $H$ can be $\chi_{QJ}$, $h_Q$ with $Q=c,b$ and $J=0,1,2$. The helicity amplitudes and the unpolarized decay widths are evaluated up to QCD next-to-next-to-leading order (NNLO) within the framework of nonrelativistic QCD (NRQCD). For the first time, we check the NRQCD factorization for $h_Q$ exclusive production at two-loop order. The leading logarithms (LL) of $m_Z^2/m_Q^2$ in the leading-twist short-distance coefficients, which may potentially ruin the perturbative convergence, are resummed to all orders of $\alpha_s$ by employing the light-cone factorization. We find the radiative corrections are considerable for $\chi_{Q2}$ and $h_Q$ productions, while are moderate or even minor for other channels. We also notice that the LL resummation can change the leading-order predictions for decay widths by more than 25\% for $\chi_{c0,2}$ and $h_c$ productions, and by around 50\% for $\chi_{c1}$ production. However, effects of the LL resummation on the next-to-leading-order and NNLO predictions are notably mitigated. Some phenomenological explorations are also performed.


Q
in the leading-twist short-distance coefficients, which may potentially ruin the perturbative convergence, are resummed to all orders of α s by employing the light-cone factorization.We find the radiative corrections are considerable for χ Q2 and h Q productions, while are moderate or even minor for other channels.We also notice that the LL resummation can change the leading-order predictions for decay widths by more than 25% for χ c0,2 and h c productions, and by around 50% for χ c1 production.However, effects of the LL resummation on the next-to-leading-order and NNLO predictions are notably mitigated.Some phenomenological explorations are also performed.

I. INTRODUCTION
The radiative decay of Z boson to a quarkonium serves to an ideal platform to study the interplay of the perturbative and nonperturbative nature of QCD.To data, experimentalists have made many endeavors to search for such processes [1][2][3], yet failed to find any signals.In recent years, several high-luminosity lepton colliders are proposed, such as ILC [4], FCC-ee [5] and CEPC [6], which are planed to run at Z mass pole for a period of time.Undoubtedly, tremendous Z bosons will be accumulated.Thus it will provide more opportunities to probe these rare decay processes.
The exclusive processes Z → quarkonium + γ have been extensively studied on the theoretical side.The computation on these processes can date back to the earlier 1980s by the authors in Ref. [7].In Ref. [8], these processes have also been studied at lowest order in α s and v 2 in both the nonrelativisitic 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], the analytic expressions of the amplitudes for Z → quarkonium + γ were obtained in the leading-power LC approximation at NLO in α s .In Ref. [13], calculations of the rates for Z → V + γ, where V signifies a vector quarkonium J/ψ or Υ, were presented.The calculations were accurate up to the leading-power LC approximation at NLO in α s and v. Shortly afterwards, the decay rates for Z → V + γ were restudied in Ref. [14], where the resummation of the leading logarithms (LL) of m 2 Z /m 2 Q , with m Z and m Q being the masses of the Z boson and heavy quark Q respectively, were carried out.In Ref. [15]. the authors have 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.Very recently, the decay rates for Z → Υ(nS)+γ have been calculated up to NLO in α s based on the NRQCD, which are proposed to determine the Zb b coupling [16].As relevant studies, the cross sections of e + e − → charmonium + γ at Z factories have been computed at LO and NLO in α s in Ref. [17] and in Ref. [18] respectively, and the cross sections of e + e − → bottomonium + γ at Z factories have been studied in Ref. [19].
In this work, we study the processes of Z boson radiative decays to a P -wave quarkonium, i.e., Z → H + γ, where H can be χ QJ , h Q with Q = c, b and J = 0, 1, 2. Based on the NRQCD factorization and helicity formulas, we compute the various helicity amplitudes at next-to-next-to-leading order (NNLO) in α s and leading order (LO) in v. Since the two typical energy scales m Q and m Z involved in these processes are widely separated, the NRQCD short-distance coefficients (SDCs) receive contribution from large logarithms of m 2 Z /m 2 Q , which may potentially ruin the perturbative convergence in α s .Fortunately, it was pointed out [20] that the NRQCD SDCs can be refactorized in the framework of LC formalism, in which the large logarithms can be resummed by employing the celebrated Efremov-Radyushkin-Brodsky-Lepage (ERBL) equation [10,21].We will carry out the LL resummation for the leading-twist helicity amplitudes.Thus our computation for Z → H +γ will be at NNLO in α s at fixed-order accuracy, meanwhile at all orders in α s at LL accuracy.
The rest of the paper is organized as follows.In Sec.II, we employ the helicity amplitude formalism to analyze the Z → H + γ processes, and build the polarized and unpolarized decay rates out of various helicity amplitudes.In Sec.III, we factorize the helicity amplitudes by employing the NRQCD factorization formalism, and parameterize the corresponding SDCs through NNLO in α s .The key technical ingredients of extracting the SDCs affiliated with each helicity amplitude through α 2 s are sketched, and values of the SDCs at various perturbative levels are presented.The corresponding details about constructions of all the helicity projectors are presented in Appendix A. We devote Sec.IV to the LC factorization for the leading-twist helicity SDCs.The resummation of the LL is formulated and explicitly carried out.In Sec.V, a detailed phenomenological analysis is performed.Finally, we summarize in Sec.VI.

II. THE GENERAL FORMULA
It is convenient to employ the helicity amplitude formalism to analyze the hard exclusive production process.The differential decay width of Z boson with polarization (along the zaxis) S z into a quarkonium H and a photon, the helicities of which are λ 1 and λ 2 , respectively, can be expressed as [22,23] where P denotes the spatial components of the H momentum, A H λ 1 ,λ 2 represents the amplitude corresponding to the helicity configuration (λ 1 , λ 2 ), and d 1 Sz,λ 1 −λ 2 (θ) is the Wigner function.Here, θ is the angle between the direction of P and the z-axis.Note that the constraint, λ 1 − λ 2 ≤ 1, is guaranteed by the angular momentum conservation.The magnitude of the spatial momentum |P| is readily determined via where m H denotes the mass of the quarkonium H, and the Källen function is defined via λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz.
Integrating over the polar angle θ and averaging over the polarization of Z, we finally obtain the integrated decay width of Z → H + γ for the helicity configuration (λ 1 ,λ 2 ) as (3) Thanks to the parity invariance [22], we have the following relations, Thus the number of independent helicity amplitudes for χ Q0 , χ Q1 , χ Q2 and h Q production can be reduced to one, two, three, and two, respectively.In the limit of m Z ≫ m Q , the helicity amplitude A H λ 1 ,λ 2 satisfies the asymptotic behavior where r is defined via r = m Q /m Z .One power of r in Eq. ( 5) 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 [24,25].
In terms of the independent helicity amplitudes, the unpolarized decay widths can be explicitly written as The main task of this work remains to compute the helicity amplitudes.Z boson interacts with quark-anti-quark pair through the tree-level weak interaction as where g is the weak coupling in SU(2) L × U(1) Y electro-weak 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 signfies the Weinberg angle.
Z boson can decay to χ QJ + γ only through the vectorial interaction, and decay to h Q + γ only through the axial-vectorial interaction.Therefore, for simplicity, it is convenient to explicitly extract the electro-weak coupling from the helicity amplitudes as

III. FRAMEWORK OF NRQCD COMPUTATION A. NRQCD factorization
According to the NRQCD formalism [9], the helicity amplitude A H λ 1 ,λ 2 can be factorized into where C H λ 1 ,λ 2 signify the dimensionless SDCs, N c = 3 is the number of the color, and the NRQCD long-distance matrix elements (LDMEs) are defined via where ψ † and χ denote the Pauli spinor fields creating a heavy quark and antiquark in NRQCD, respectively, and with ǫ χ Q1 , ǫ h Q , and ǫ χ Q2 representing the polarization vector/tensor of χ Q1 , h Q , and χ Q2 respectively.Exploiting the heavy quark spin symmetry, we can make the following approximations In Eq. ( 9), the factor √ 2m H appears on the right side because we adopt relativistic normalization for the quarkonium H, but we use conventional nonrelativistic normalization for the LDMEs.In this work, we will not compute the relativistic corrections, therefore it is reasonable to take the approximation m H ≈ 2m Q .
Through Eq. ( 5) and Eq. ( 9), we can readily deduce the helicity selection rule for the SDCs Q .The SDCs are insensitive to the nonperturbative hadronization effects, therefore they can be determined with the aid of the standard perturbative matching technique.That is, by replacing the physical H meson with a fictitious onium composed of a free Q Q pair, carrying the quantum number 3 P J for χ QJ and 1 P 1 for h Q , we compute both sides of Eq. ( 9) .After that, we are able to solve for the desired SDCs.For more details, we refer the readers to Ref. [26].
It is convenient to parametrize the SDCs in powers of α s where µ R and µ Λ signify the renormalization scale and factorization scale, respectively, where n f is the number of active quark flavors.The explicit ln µ 2 R term is deduced from the renormalizationgroup invariance.γ H represent the anomalous dimensions associated with the NRQCD bilinear currents carrying the quantum numbers 3 P J or 1 P 1 , the expressions of which read [27] The occurrence of ln µ Λ is required by the NRQCD factorization.According to the factorization, the µ Λ dependence in the SDCs should be thoroughly canceled by that in the LDMEs.As illustrated in Fig. 1, we classify the Feynman diagrams into 'regular' part and 'nonregular' part at O(α 2 s ).Correspondingly, C H,( 2) nonreg,λ 1 ,λ 2 in Eq. ( 14) represent contributions from the 'regular' Feynman diagrams and the 'nonregular' Feynman diagrams respectively.
The quark-level Feynman diagrams and Feynman amplitudes are generated using FeynArts [28].Employing the color and spin projectors followed by enforcing spin-orbit coupling, we obtain the hadron-level amplitudes order by order in α s with the aid of the packages FeynCalc [29] and FormLink [30].To evaluate the helicity amplitudes, we employ the technique of helicity projection.The concrete expressions of all the helicity projectors are presented in Appendix.A.
It is well known that, in dimensional regularization, the anticommutation relation {γ µ , γ 5 } and the cyclicity of Dirac trace can not be satisfied simultaneously.In practical computation, the naive-γ 5 scheme [31], 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.Due to the lack of the cyclicity of the trace, one must fix a reading point for a fermion loop with odd number of γ 5 .In our work, we will select the vertex of Z boson as the reading point.
Then, it is straightforward to obtain the LO helicity SDCs: Once beyond the LO, we adopt the standard shortcut to directly extract the SDCs, i.e., compute the hard region in the context of strategy of region [32].Utilizing the packages Apart [33] and FIRE [34], we can further reduce the loop integrals into linear combinations of master integrals (MIs).Finally, we end up with 6 one-loop MIs, which are computed using Package-X [35], and roughly 320 two-loop MIs, the evalutaion of which is a challenging work.Fortunately, a powerful new algorithm, dubbed Auxiliary Mass Flow (AMF), has recently been pioneered by Liu and Ma [36][37][38][39].Its main idea is to set up differential equations with respect to an auxiliary mass variable, with the vacuum bubble diagrams as the boundary conditions.Remarkably, these differential equations can be solved iteratively with very high numerical precision.In this work, we utilize the newly released package AMFlow [40] to compute all the two-loop MIs.After implementing the on-shell renormalization scheme for the heavy quark mass and field strength [41], and the MS renormalization scheme for the QCD coupling, the UV poles are exactly canceled, while an uncancelled single IR pole still remains.This symptom is a common feature specific to the NRQCD factorization, which have been encountered many times in NNLO perturbative calculations involving quarkonium.This IR pole can be factored into the NRQCD LDME, so that the NRQCD SDCs become IR finite.We have numerically verified that the coefficient of the remaining IR pole equal to one-quarter of the anomalous dimension in Eq. ( 15) with high precision, as required by the NRQCD factorization.
The analytic expressions of λ 1 ,λ 2 can be readily obtained.Instead of presenting the cumbersome expressions, here we merely present their asymptotic expansions in r → 0: where the real parts of the results in Eqs.(17a-f) are consistent with these in Ref. [42] 1 , and the results in Eqs.(17a,c,f,h) are consistent with these in Ref. [12].It is rather challenging for us to analytically compute λ 1 ,λ 2 , so we turn to numerically evaluate their values.To perform the numerical computation, we take m Z = 91.1876GeV from the latest particle data group (PDG) [43], and the pole masses of the charm quark and bottom quark to be m c = 1.69 GeV and m b = 4.80 GeV, which are converted from the MS masses mc ( mc ) = 1.28 GeV and mb ( mb ) = 4.18 GeV [43] at two-loop level by use of the package RunDec [44].
We tabulate the results of the SDCs C H,(1) nonreg,λ 1 ,λ 2 in Tab.I for charmonium production, and in Tab.II for bottomonium production 2 .For the sake of reference, we explicitly remain the n l , n c and n b dependence in the SDCs, where n l denotes the number of the light quarks, and n c = 1 and n b = 1 signify the numbers of the charm quark and bottom quark respectively.

IV. LC FACTORIZATION FOR THE LEADING-TWIST SDCS A. The LC factorization
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 1 In Ref. [42], only the cross sections of e + e − → χ cJ + γ are presented, where the imaginary parts of the amplitudes are ignored. 2Since weak interaction in Standard Model (SM) is a chiral gauge theory, the gauge anomaly should be avoided for physical processes.To satisfy the condition of anomaly free, when evaluating C hQ,(2) nonreg,λ1,λ2 which corresponds to contribution from Fig. 1d, we have included all six flavor quark loops.TABLE I: NRQCD predictions to the various helicity SDCs defined in Eq.( 14) for charmonium production.For simplicity, we define the symbols f 1 ≡ where g u V and g d V correspond to the values of g V for up-type quark and down-type quark respectively.the spirit of Ref. [20], the LC factorization formula for the SDCs is written as where C H,LL 0 0,1 represents the asymptotic expansion of C H,(0) 0,1 in r → 0, 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.At LO in α s , we have We can reproduce the asymptotic expansions of the SDCs C H,(1) 0,1 in Eq. ( 17) exactly with the corresponding NLO corrections to T H and φH which have been calculated in [12].More importantly, we can employ the LC factorization to resum the LL terms (α s ln r 2 ) n in C H 0,1 .

B. Resummation of the LL with the ERBL equation
The leading twist LCDAs φH obey the celebrated ERBL equation [21,45] with the Brodsky-Lepage (BL) kernel where the subscript "+" implies the familiar "plus" prescription.Solving the ERBL equation, we can obtain the SDCs with all the LL terms (α s ln r 2 ) n resummed.Formally, we have with where with "⋆" stands for the convolution It is well known that the BL kernel has the eigen-system that where and the eigen-functions with C (3/2) n (2x − 1) being order 3/2 Gegenbauer polynomials.We can decompose the LCDA φ(0) where Hence, we can get exp (κC Employing and the decompositions We get for H = χ Q1 , and for With the explicit expansion of the formal solution in Eq. ( 24) in κ and the formulas we can obtain the expansion of K LL H in α s , for H = χ Q1 , and for H = χ Q0 , χ Q2 and h Q .
C. SDCs by combining the fixed-order results and the LL resummation In the following, we will improve the leading-twist SDCs by combining the NRQCD fixedorder results given in Sec.III and the LL resummation sketched in Sec.IV A and Sec.IV B.
To avoid double counting, it is necessary to subtract the terms of O(α n s ln n r) from the fixed-order SDCs.Formally, we have where the superscripts 'LO', 'NLO' and 'NNLO' denote the fixed-order SDCs accurate up to O(α 0 s ), O(α 1 s ) and O(α 2 s ) respectively, and the superscripts 'LL 0 ', 'LL 1 ' and 'LL 2 ' signify the C H,LL 0,1 truncated at O(α 0 s ), O(α 1 s ) and O(α 2 s ) respectively, which have been explicitly computed in Sec.IV B.
In Tab.III, we present the theoretical predictions on the squared leading-twist SDCs |C H 0,1 | 2 at various levels of accuracy.In the table, we also enumerate the values of K LL H , which are computed by applying the Eq. ( 35) and Eq. ( 36) with α s (m Z ) = 0.1181, taken from the PDG, and α s (m c ) = 0.3240 and α s (m b ) = 0.2151, evaluated through the renormalization group running at two loop with the aid of the package RunDec.It is worth noting that, in order to accelerate the convergence, we have used the so-called Abel-Padé method [46] to sum the series in Eq. ( 35) and Eq.(36).From Tab.III, we have several observations.Firstly, we find the fixed-order predictions for charmonium production are close to these for bottomonium production, however the LL resummation can give arise to some difference.Secondly, we notice that the effect of the LL resummation can considerably change the LO results, especially for charmonium production, for instance, it can change the LO results by more than 25% for χ c0,2 and h c production, and by around 50% for χ c1 production.The magnitude of the LL resummation for bottomium production is roughly half of that for charmonium case.Finally, it is worth noting that the effects of the LL resummation on the NLO and NNLO predictions are continuously becoming milder.It can be explained by the fact that some of the LL resummation have been included in the radiative corrections.Concretely, the O(α s ln r) contribution has been included in the NLO prediction, while both the O(α s ln r) and O(α 2 s ln 2 r) contributions have been included in NNLO prediction.As a consequence, the remaining contribution from the LL resummation can correspondingly reduce its effect on the NLO and NNLO predictions.

V. PHENOMENOLOGY
To make concrete phenomenological predictions, we need to fix the various input parameters.We take m Z = 91.1876GeV from PDG, and the heavy quark pole masses m c = 1.69 GeV and m b = 4.80 GeV, as mentioned which are converted from their MS masses.We take the running QED coupling constant evaluate at the mass of m Z as α(m Z ) = 1/128.943[47].We take the NRQCD factorization scale µ Λ = 1 GeV.The NRQCD LDMEs for χ QJ and h Q are approximated by the first derivative of the Schrödinger radial wave function at origin through where the 1P radial wave functions at origin, evaluated in Buchmüller-Tye (BT) potential model, are taken from Ref. [48].In addition, we take s 2 W = 0.231, and the value of the total decay width of Z boson Γ Z = 2.4952 GeV from the PDG [43].
With all ingredients in hand, we can compute the decay widths of Z → H + γ.The unpolarized decay widths at various levels of accuracy for charmonium production and bottomonium production are separately tabulated in Tab.IV and Tab.V.In the tables, we have included the uncertainties from ambiguity of renormalization scale and QCD higher-order corrections.We should emphasize that the values of the Schrödinger wave functions may largely affect the theoretical predictions on the decay rates, i.e., |R ′ 1P,cc (0)| 2 can range from 0.07 to 0.13 GeV 5 , |R ′ 1P,b b(0)| 2 can range from 0.93 to 2.07 GeV 5 in Refs.[48,49], which may change the central values of the decay rates of quarkonium production by roughly a factor of 2.
It is interesting to note that both the O(α s ) and O(α 2 s ) corrections to Z → χ Q2 /h Q + γ are sizable and negative.The LL resummation turns out to further decrease the decay widths.Incorporating all the perturbative corrections and LL resummation reduces the LO prediction by roughly half of magnitude.In contrast, both the O(α s ) and O(α 2 s ) corrections are moderate or even minor for other channels.The situation is quite similar to the case in Ref. [26], where the radiative corrections are significant for e + e − → χ c2 + γ at B factory, however inconsiderable for e + e − → χ c0,1 + γ.
It is enlightening to compare the strengths of the decay widths for different quarkonium production.For charmoinum production, we find that h c + γ production has the biggest branching fraction, followed by χ c1 +γ production.Although the branching fraction of χ c2 +γ is two times larger than that of χ c0 +γ at LO, their branching fractions are nearly the same at 'NNLO+LL' accuracy.For bottomonium production, we notice that the branching fractions of χ b1 + γ, h b + γ, χ c2 + γ, χ b0 + γ make the most, second, third, and last strengths.
It is worth noting that, for the same quantum number of quarkonium, the branching fraction of charmonium production is larger than that of bottomium production.Finally, we estimate number of the quarkonium production at the proposed super Z factories, such as the Z-factory mode in CEPC, where Z boson yield will reach 7 × 10 11 [6].Thus it is expected that there will be several hundreds and thousands of charmonia and bottomonia production through Z → H + γ.The signal for Z → χ QJ + γ production can be measured by probing χ QJ with a recoiling hard photon, where χ QJ can be reconstructed from their transition to γ + J/ψ (γ + Υ) with J/ψ(Υ) → ℓℓ.Due to the low multiplicative branching ratio, and extra event-selection rules to suppress the backgrounds, it will be rather difficult to measure Z → χ QJ + γ in experiment.Alternatively, the χ QJ may be reconstructed through its hadronic decays, however, the experimental measurements on χ QJ + γ are still challenging.The condition for Z → h Q + γ is even worse.

VI. SUMMARY
In summary, we study the exclusive decay processes of Z → χ QJ /h Q + γ in the NRQCD framework.The amplitudes of all the helicity configurations and the unpolarized decay widths are evaluated up to O(α 2 s ).It is the first time that the NRQCD factorization for h Q exclusive production at two loop is verified explicitly.The LL of m 2 Z /m 2 Q in the leadingtwist SDCs are resummed to all orders of α s by employing the LC factorization.We find the radiative corrections are considerable for χ Q2 and h Q productions, while are moderate or even minor for other channels.We also notice that the LL resummation can change the LO results by more than 25% for χ c0,2 and h c production, and by around 50% for χ c1 production.However, effects of the LL resummation on the NLO and NNLO predictions are notably mitigated.We expect that several hundreds and thousands of charmonia and bottomonia will be produced through Z → H + γ at the proposed super-Z factories.

TABLE II :
NRQCD predictions to the various helicity SDCs for bottomonium production.For simplicity, we define the symbols f1 ≡

TABLE IV :
Unpolarized (total) decay widths for Z → charmonium+γ at various levels of accuracy.The predictions are obtained by setting µ R = m Z / √ 2. The first error at accuracy of 'NNLO' and 'NNLO+LL' is estimated by varying µ R from m Z /2 to m Z , while the second error is from the QCD higher-order corrections, which are estimated through α 3 s ∼ 0.002.

TABLE V :
Unpolarized (total) decay widths for Z → bottomonium+γ at various levels of accuracy.The source of the theoretical uncertainties is the same as that in Tab IV.