Study on $\eta_{c2}(\eta_{b2})$ electromagnetic decay into double photons

Within the framework of nonrelativistic QCD (NRQCD) factorization formalism, we compute the helicity amplitude as well as the decay width of $\eta_{Q2}$ ($Q=c,b$) electromagnetic decay into two photons up to next-to-next-to-leading order (NNLO) in $\alpha_s$ expansion. For the first time, we verify the validity of NRQCD factorization for the D-wave quarkonium decay at NNLO. We find that the $\mathcal{O}(\alpha_s)$ and $\mathcal{O}(\alpha_s^2)$ corrections to the helicity amplitude are negative and moderate, nevertheless both corrections combine to suppress the leading-order prediction for the decay width significantly. By approximating the total decay width of $\eta_{Q2}$ as the sum of those for the hadronic decay and the electric $E1$ transition, we obtain the branching ratios ${\rm Br}(\eta_{c2}\to 2\gamma)\approx 5\times10^{-6}$ and ${\rm Br}(\eta_{b2}\to 2\gamma)\approx 4\times10^{-7}$. To explore the potential measurement on $\eta_{Q2}$, we further evaluate the production cross section of $\eta_{Q2}$ at LHCb at the lowest order in $\alpha_s$ expansion. With the kinematic constraint on the longitudinal rapidity $4.5>y>2$ and transverse momentum $P_T>(2-4)m_Q$ for $\eta_{Q2}$, we find the cross section can reach $2-50$ nb for $\eta_{c2}$, and $1-22$ pb for $\eta_{b2}$. Considering the integrated luminosity $\mathcal{L}=10\, {\rm fb}^{-1}$ at $\sqrt{s}=7$ TeV and $\sqrt{s}=13$ TeV, we estimate that there are several hundreds events of $pp\to \eta_{c2}\to 2\gamma$. Since the background is relatively clean, it is promising to reconstruct $\eta_{c2}$ through its electromagnetic decay. On the contrary, due to small branching ratio and production cross section, it is quite challenging to detect $\eta_{b2}\to 2\gamma$ at LHCb.


I. INTRODUCTION
Heavy quarkonium, as a multiscale system, is an ideal laboratory for testing the interplay between perturbative and nonperturbative QCD.Its mass spectrum has been predicted by various potential models, and most of the low-lying quarkonium states have been probed by the experiment.However there still exit some undiscovered states.Among the missing states in charmonium family, η c2 ( 1 D 2 ) is the only spin-singlet low-lying D-wave state.A full understanding on η c2 in both theory and experiment can help to illuminate the interquark force and reveal the nature of the strong interaction.
The mass of η c2 is predicted to range from 3.80 to 3.88 GeV [1][2][3][4][5][6][7], which well lies between the D D and the D * D thresholds.Quite different from ψ ′′ , the decay of η c2 into D D is forbidden, which is accounted for by the parity conservation.Thus η c2 is a narrow resonance, and its main decay modes are considered to be hadronic decay and electric E1 transition, which have been well investigated in the references.
The C = +1 charmonia bear a considerable branching ratio of electromagnetic decay into double photons, e.g., Br(η c → 2γ) ≈ 1.57 × 10 −4 , Br(χ c0 → 2γ) ≈ 2.04 × 10 −4 , and Br(χ c2 → 2γ) ≈ 2.85 × 10 −4 [24], which have been well measured by the experiment.Although the electric E1 transition η c2 → h c γ comprises one of the main decay channel, it is anticipated that the branching ratio for η c2 electromagnetic decay is of the same magnitude as those in η c and χ c .Compared with the hadronic decay, the electromagnetic decay has the advantage of bearing a clean background.Thus we propose to detect η c2 through η c2 → 2γ.In this work, we will evaluate the partial width of η Q2 → 2γ (the heavy quark flavor Q = c, b) within framework of the well established NRQCD factorization formalism.
NRQCD factorization formalism is widely employed to tackle heavy quarkonium decay and production.Within this framework, the production cross section or decay width can be systematically disentangled the short-distance and long-distance effect, formalized by a double expansion in powers of heavy quark velocity v Q and strong coupling constant α s .The perturbative contributions with the scale larger than the heavy quark m Q is encoded into the short-distance coefficients (SDCs), while the nonperturbative effects are depicted by the NRQCD long-distance matrix elements (LDMEs).
In addition, to provide aid for experimental search for η Q2 through η Q2 → 2γ, we will explore the η Q2 production at colliders, and estimate the corresponding number of events.
The η c2 associated production with a photon at B factory has been obtained in Ref. [42].The η c2 + J/ψ production at B factory can be found in Ref. [43].Unfortunately, the cross sections of both channels are too small to probe η c2 → 2γ.Considering the significant luminosity as well as the sizeable cross sections for quarkonia production at LHC, exemplified by σ(pp → η c +X) ∼ 0.5 µb with the transverse momentum cut p T > 4m c [44], we anticipate the cross section for η c2 production is also considerable.Actually, the cross section of η c2 at LHC can be simply estimated by σ(η c2 ) ∼ c × 0.5 ∼ 5 nb.Thus, we expect there are amount number of events for pp → η c2 → 2γ, which renders η c2 → 2γ a promising channel to probe η c2 .
The remainder of this paper is organized as follows.In Sec.II, we make Lorentz decomposition for the amplitude of η Q2 → 2γ, and present the decay width in term of the helicity amplitudes.In Sec.III, we outline the NRQCD factorization formalism for the helicity amplitude, and briefly describe the theoretical framework to deduce the NRQCD SDC.In Sec.IV, we first introduce the technicalities encountered in performing loop calculation, and then present our final results for the SDC.Sec.V is devoted to the phenomenological analysis and discussion.A theoretical prediction on the production cross section of η Q2 at LHCb is also contained in this section.We present our summary in Sec.VI.

II. THEORETICAL FORMULA FOR THE DECAY WIDTH
The partial width of η Q2 → 2γ can be expressed in term of helicity amplitude where J = 2 denotes the spin of η Q2 ,1 2! accounts for the indistinguishability of the two identical photons, 1 8π corresponds to the phase space factor, and A λ 1 ,λ 2 signify the helicity amplitudes of η Q2 → γ(λ 1 )γ(λ 2 ) with λ 1,2 = ±1 being the helicity of the photons.By invoking the parity invariance, we only enumerate the independent helicity amplitudes 1 .
To extract the helicity amplitudes, we first decompose the amplitude A of η Q2 → 2γ by Lorentz invariance.By Bose symmetry, transversality for the polarization of photon, together with parity conservation, the amplitude can be generically expressed as where p denotes half of the momentum of η Q2 , k 1 and k 2 signify the momenta of the two outgoing photons, c i (i = 1, 2, 3) are Lorentz invariant and refer to form factors of the corresponding Lorentz structure, ǫ 1 and ǫ 2 represent the polarization vectors of the photons, and ǫ H denotes the polarization tensor of η Q2 .
It is straightforward to deduce the helicity amplitudes A λ 1 ,λ 2 from Eq. (2).To carry out the calculation, it is convenient to construct the explicit expressions of the polarization tensor ǫ H and polarization vectors ǫ 1 and ǫ 2 [46].We define the polarization vectors The five polarization tensors of η Q2 are readily expressed as If assuming the photon with momentum k 1 is outgoing in the positive z direction, its polarization vector ǫ 1 equals to ǫ + for helicity +1, and ǫ − for helicity −1, while the helicity polarization vector of the backward photon is just reversed.Substituting the explicit expressions Eq. ( 3) and Eq. ( 4) into Eq.( 2), we obtain the helicity amplitudes in term of the three form factors where A 1,−1 explicitly vanishes.

III. NRQCD FACTORIZATION FORMALISM FOR THE HELICITY AMPLI-TUDE
Owing to the strong interaction inner the hadron η Q2 , the helicity amplitude A 1,1 is nonperturbative.Fortunately, we can employ the NRQCD factorization formalism to factorize the helicity amplitude into [18] A where C 1,1 represents the perturbative SDC which depicts a heavy quark pair annihilation into double photons, µ F signifies the factorization scale, and with ǫ H being the polarization tensor of η Q2 .LDME 0|χ † K1 D 2 ψ(µ F )|η c2 deciphering the nonperturbative effect in the hadron is process-independent, and can be related to the second derivative of the radial wave function at the origin through To get the helicity amplitude, we must determine the SDC.Since the SDC is irrelevant to the nonperturbative hadronization effect, it can be computed through the standard matching technique.Concretely, we can replace the physical meson η Q2 with a heavy quark pair Q Q, carrying the same quantum number as 1 D 2 .The factorization formalism is also valid to the free quark state Q Q( 1 D 2 ), therefore after the replacement, Eq. ( 6) becomes to where we have suppressed the factorization scale and the high-order relativistic corrections.
The SDC in Eq. ( 9) is exactly the same as in Eq. ( 6).Since the amplitude A 1,1 and matrix element are now perturbative, both sides of Eq. ( 9) are calculable.In principle, one can solve the SDC C 1,1 in any prescribed α s order.
In the following, we briefly describe the procedure to evaluate the perturbative helicity amplitude A 1,1 ( 1 D 2 ).We assign the momenta of the Q and Q quarks to be where p and q represent half of the total momentum and relative momentum of Q Q pair, respectively.The on-shell condition enforces that In our calculation, we first evaluate the amplitude A of Q Q → 2γ, then employ the covariant spin-projector to extract the spin-singlet pattern of Q Q.To be consistent with the decay width formula (1), we utilize the nonrelativistically normalized spin-singlet/colorsinglet projector [47], which reads The L = 2 orbital partial wave can be projected out by differentiating the color-singlet/spinsinglet quark amplitude with respect to the relative momentum q, followed by setting q to zero Now, we have collected all the necessary ingredients to calculate the amplitude of Q Q( 1 D 2 ) → 2γ.Subsequently, we can pick up the Lorentz invariant form factor c i through Eq. ( 2), and obtain the perturbative helicity amplitude with the aid of Eq. ( 5).Meanwhile, the perturbative NRQCD matrix element 0|χ † K1 D 2 ψ|Q Q( 1 D 2 ) can also be carried out at desired α s order.At lowest order in α s , we have Finally, it is straightforward to determine the SDC C 1,1 at prescribed α s order through Eq. ( 9).For the future convenience, we re-express the partial width of η Q2 → 2γ in term of SDC IV. SDC UP TO NNLO In this section, we first describe the computational technicalities utilized to evaluate the perturbative amplitude in detail, then present our main results for the SDC C 1,1 .
LO NLO NNLO ("light by light") We employ FeynArts [48] to generate the Feynman diagrams and the corresponding amplitude for Q Q → 2γ.The representative Feynman diagrams are illustrated in Fig. 1.We employ the spin-singlet/color-singlet projector (12) and apply the recipe as specified in (13) to project out the intended amplitude for Q Q( 1 D 2 ) → 2γ up to two-loop level.Subsequently, the packages FeynCalc [49,50] and FormLink [51,52] are employed to perform the Dirac trace and Lorentz contraction.
For the NLO and NNLO corrections, we carry out the derivative of the amplitude with regard to the relative momentum q prior to perform loop integration, which amounts to directly extract the contribution from the hard region.We employ the package Apart [53] and FIRE [54] to conduct partial fraction and the corresponding integration-by-part (IBP) reduction.Finally, we have 3 one-loop master integrals (MIs) and 80 two loop MIs.There exist some complex-valued two-loop integrals originating from the light-by-light (lbl) diagrams, which are relatively hard to carry out numerically.It deserves mentioning that although the MIs encountered here are almost the same sets as in the processes η c → 2γ and χ c → 2γ, the computation complexity is much more involved.Since we have taken the second derivative of the amplitude with respect to the relative momentum q, some coefficients of the MIs are relatively larger as well as more divergent in 1/ǫ expansion compared with those in η c and χ c decay.Thus, to reach the desired precision, we must perform numerical integration over the MIs to higher accuracy.
For the real-valued MIs, we directly use CubPack/HCubature [55,56] to carry out the integration.In contrast to the application of SD to the Euclidean region, the singularities encountered in the physical region lie inside, rather than sit on, the integration boundary, which render the integration hard to be numerically evaluated.To overcome this difficulty, we conduct integration contour deformation via the variable transformation prior to decomposing the sectors [57], and determine the integration contour through optimizing a set of contour parameters [38].For more technical details, we refer the readers to Refs.[38,41].
The lbl Feynman diagrams are both gauge invariant and free of any UV and IR divergences in sum.On the contrary, there exist UV divergence at one loop, and both UV and IR divergences at two loop.The UV divergence originates from the integration over the loop momentum, which can be eliminated through the standard renormalization procedure.We implement the on-shell renormalization for the heavy quark wave function and mass up to O(α 2 s ) [58][59][60], and MS renormalization for the strong coupling constant.Thus, the ultimate amplitude is completely UV finite.There remains a piece of unremoved IR divergence, nevertheless this IR pole can be factored into the NRQCD LDME, so that the NRQCD SDC becomes IR finite.As a consequence, both of the LDME and the corresponding twoloop SDC C 1,1 bear ln µ F dependence, nevertheless their product must be independent of factorization scale.
After some hard work, we finally obtain the SDC C 1,1 expanded in power of the strong coupling constant α s where α denotes the electromagnetic coupling constant, e Q signifies the electric charge of the heavy quark, β 0 = 11 3 C A − 2 3 (n L + n H ) corresponds to the one loop coefficient of the QCD β function, where n H = 1, and n L signifies the number of the active quark flavor (n L = 3 for η c2 , and n L = 4 for η b2 ).The exact occurrence of the ln µ 2 R is demanded by the renormalization group invariance.
∆ (1) and ∆ (2) correspond to the O(α s ) and O(α 2 s ) corrections to the SDC.The expression of ∆ (1) is analytically obtained ∆ (2) can be expressed as The coefficient of the factorization scale dependence ln µ F term corresponds to the anomalous dimension of the NRQCD operator in Eq. ( 6), which is consistent with the result in Ref. [61] 2 , and thereby the NRQCD factorization is verified in η Q2 electromagnetic decay.The terms of ∆ (2) reg and ∆ (2) lbl in (18) represent the contributions from the regular and 'light-by-light' Feynman diagrams, which are illustrated in Fig. 1.
Furthermore, we can organize the ∆ reg and ∆ (2) lbl according to the color structure, where 2 In Ref. [61], the authors computed the anomalous dimensions of spin-single and spin-triplet currents for heavy quark pair with arbitrary orbital angular momentum.The anomalous dimension of 1 D 2 can be obtained by utilizing Eq. ( 40) in Ref. [61] with color-singlet Wilson coefficients of the potentials given in Refs.[62][63][64][65].

and ∆
(2 By setting the renormalization scale µ R = m Q and factorization scale µ F = 1 GeV, we get the radiative corrections to C 1,1 at various perturbative orders, where r = 2.12 + 0.0005 i for η c2 and r = 1.11 + 0.005 i for η b2 .We find that both the O(α s ) and O(α 2 s ) corrections to the helicity amplitude are negative as well as moderate.It seems that the perturbative expansion for η Q2 → 2γ exhibits a decent convergence, however as will be found, the radiative corrections accurate up to O(α 2 s ) change the LO decay width considerably.
For completeness, it is necessary to deduce the explicit expression of decay width.Applying the formula Eq. ( 15) and the expression of helicity amplitude in Eq. ( 16), we readily obtain the decay width of where the symbol Re signifies the real part of the argument.

A. predictions for the decay width
To make concrete prediction, we first choose the input parameters.We take the heavy quark mass to be m c = 1.68 GeV and m b = 4.78 GeV, which correspond to the twoloop charm quark and bottom quark pole masses converted from the corresponding MS masses [66].We evaluate the electromagnetic coupling constant as α(2m c ) ≈ 1 132 and α(2m b ) ≈ 1 131 by the formulas in Ref. [67], and evaluate α s at each energy scale by RunDec [66].
The NRQCD LDME is related to the second derivative of the 1D radial wave function at the origin in Eq. ( 8), which are well determined by the nonrelativistic potential model.In this work, we take the R ′′ D (0) from the Cornell potential model [68,69], and determine the LDMEs as With these input parameters, we present our predictions for the decay widths of η c2 /η b2 → 2γ at various levels of accuracy in α s in Tab.I.The uncertainties affiliated with the decay TABLE I: NRQCD predictions for the decay width of η Q2 → 2γ at various levels of accuracy in α s .We take the two-loop quark pole masses m c = 1.68 GeV and m b = 4.78 GeV.The NRQCD LDMEs are evaluated by the Cornell potential model.The errors are estimated by sliding the renormalization scale µ R from m Q to 2m Q with center value µ R = √ 2m Q .By taking the total decay width as Γ total (η c ) ≈ 445.1 keV and Γ total (η b ) ≈ 29.1 keV, we also present the branching ratio Br(η Q2 → 2γ).width are estimated by varying µ R from m Q to 2m Q with the central values evaluated at √ 2m Q . 3From the table, we have several observations.First, the NNLO decay width is much smaller than the LO one for the channel η c2 → 2γ, which is accounted for by the sizeable and negative O(α s ) and O(α 2 s ) radiative corrections.Second, the decay width is insensitive to the factorization scale µ F .Third, the decay width of η b2 → 2γ is considerably smaller than the case of η c2 , which is mainly caused by the heavier quark mass and smaller electric charge for bottom quark.
We can also predict the branching ratio of η Q2 → 2γ.According to current theoretical computation, η c2 decay predominately through the electric E1 transition and hadronic decay.If assuming η Q2 decay is saturated by these two decay patterns, we can approximate the total decay width through [11] Γ total (η Q2 ) ≈ Γ(η Q2 → LH) + Γ(η Q2 → γh Q ), (25) where LH denotes the abbreviation for light hadrons.The hadronic decay width of η Q2 up to NLO has been known for a long time [11], and the prediction for electromagnetic E1 transition of 1 D 2 → 1 P 1 from Cornell potential model can be found in Refs.[2,8].Thus, we readily obtain the total decay width Γ total (η c2 ) = 142.1 + 303.0 = 445.1 keV and Γ total (η b2 ) = 3.8 + 25.3 = 29.1 keV, where we have reevaluated the hadronic decay width with the parameters selected in this work.Consequently, the branching ratio of η Q2 → 2γ is illustrated in Tab.I.It is remarkable that Br(η b2 → 2γ) is much smaller than Br(η c → 2γ), which renders the search for η b2 through its electromagnetic decay quite challenging.
B. η Q2 production at colliders In the following, we will evaluate the production cross section of η Q2 at B factory and LHC.The η c2 associated production with a photon at B factory has been studied in Ref. [42], where √ s = 10.58GeV is the center-of-mass (CM) energy at B factory.With the aforementioned input parameters, we immediately arrive at σ(e + e − → η c2 + γ) = 1.49fb.Due to the small cross section and branching ratio, it seems impossible to detect η c2 through its electromagnetic decay at B factory.Now we turn to the hadron collider LHC, where amount number of quarkonia can be produced due to the considerable production cross section, e.g., the cross section of η c can reach around 0.5µb through single parton scattering (SPS) [44].For pp → η Q2 + X, the differential cross section through SPS can be factorized as where we have neglected the color-octet contribution, f i/p (x) represents the parton distribution function of a proton, and dσ denotes the partonic cross section.Since the gluon distribution is overwhelming in the proton at small momentum fraction x, we expect that the gluon scattering will denominate the cross section.Thus, it is reasonable for us to consider gluon-gluon partonic scattering to estimate the cross section of η Q2 .In addition, we will carry out the partonic cross section dσ at lowest order in α s .For concreteness, we consider the production of η Q2 at LHCb detector, where a kinematic constraint on the longitudinal rapidity of η Q2 4.5 > y > 2 is implemented.We further take a transverse momentum P T cut for the η Q2 to guarantee the validity of the factorization formula (27).In our computation, we employ CTEQ14 PDF sets [74] for the proton PDF.
In Tab II, we present the theoretical predictions for the cross section of η Q2 at two benchmark CM energy √ s = 7 TeV and 13 TeV with various transverse momentum cutoffs for η Q2 .From the table, we find the cross section of η b2 is smaller than that of η c2 by roughly three order of magnitude.Combing the luminosity at LHCb, we can estimate the number of events for η Q2 production.With the integrated luminosity L = 10 fb −1 at each CM energy, there are 10 7 − 10 8 η c2 and 10 4 − 10 5 η b2 event produced at LHCb.Therefore, LHCb will be an ideal platform to probe η Q2 .Furthermore, multiplying the branching ratio of η Q2 → 2γ, we predict that there are several hundreds of double-photon events through pp → η c2 → 2γ, which is a promising channel to probe this undiscovered charmonium.In contrast, the η b2 is proved to be hard to detect through its electromagnetic decay at LHCb.

VI. SUMMARY
Applying the NRQCD factorization formalism, we evaluate the η Q2 electromagnetic decay into double photons up to O(α 2 s ) radiative corrections.For the first time, we scrutinize the validity of the NRQCD factorization for D-wave quakonium decay at NNLO.Both the O(α s ) and O(α 2 s ) corrections to the decay width of η Q2 → 2γ are negative.Although the radiative corrections to the helicity amplitude are moderate, the corrections change the LO decay width significantly, especially for η c2 .By assuming η Q2 decay is saturated by the electric E1 transition and the hadronic decay, we obtain the branching ratios Br(η c2 → 2γ) ≈ 5 × 10 −6 and Br(η b2 → 2γ) ≈ 4 × 10 −7 .We have also studied the η Q2 production at LHCb.By imposing kinematic restriction on the longitudinal rapidity and transverse momentum of η Q2 , we predict the cross sections to be 2−50 nb for η c2 and 1−22 pb for η b2 for various transverse momentum cutoffs.Thus, it is promising to observe η c2 through its electromagnetic decay at LHCb, while quite challenging to detect η b2 at current integrated luminosity.

TABLE II :
The cross sections of η c2 and η b2 at the lowest order in α s expansion at LHCb.By taking the CM colliding energy √ s = 7 TeV and 13 TeV, we evaluate the integrated cross sections with η Q2 longitudinal rapidity constraint 4.5 > y > 2 and transverse momentum P T cut.