One-loop corrections to the processes $e^+e^- \to \gamma, Z \to J/\psi~\eta_c$ and $e^+e^- \to Z \to J/\psi~J/\psi$

e+e− → Z0 → J/ψ J/ψ A. V. Berezhnoy,1, ∗ I. N. Belov,1, 2, † S. V. Poslavsky,3, ‡ and A. K. Likhoded3, § 1SINP MSU, Moscow, Russia 2Physics department of MSU, Moscow, Russia 3NRC “Kurchatov Institute”–IHEP, Protvino, Russia The cross section values of J/ψ ηc and J/ψ J/ψ production in e+e− annihilation are estimated within one-loop approximation near Z0 pole, as well as at higher energies. Both intermediate bosons, γ and Z0, are taken into account. It is shown that at Z0 mass the NLO contribution increases the cross-section values by 3.5 times.

Study of charmonia production in the future is encouraged by two big projects: ILC and FCC. Both of them propose e + e − collisions at energies of order of Z-boson's mass: energy range announced for FCC is √ s = 90 ÷ 400 GeV and √ s = 250 GeV is proposed for ILC.
Also it is worth to note, the Z 0 -boson decays to two charmonia may be of some interest for experiments at the LHC. The J/ψ η c and J/ψ J/ψ pair production near Z 0 -boson pole were researched within LO in the papers [21,22]. Recently in the work [23] Z 0 -boson decays to two charmonia were studied in the framework of lightcone formalism, which allowed to take into account internal motion of quark inside charmonium. Complementing these works, we take Discussing J/ψJ/ψ production in e + e − annihilation we should mention the theoretical study of this final state production via double photon exchange [24,25] and so far unsuccessful attempts to find such a process at B-factories [1,2,26]. However, in this paper we restrict ourselves by studying e + e − annihilation into one boson only.
In this paper we consider double charmonia production in close relation to the B c pair production in e + e − -annihilation, which was studied in our previous work [27] (see also [28], where the B c pair production was studied in γγ-fusion).

THE METHOD
The discussed production of the charmonium pair through single boson exchange is affected by several selection rules. First of all neither photon nor Z-boson can decay to η c η c -pair becаuse the pair of η c mesons in final state must simultaneously have the total angular momentum 1 and the symmetric wave function, which is impossible. Also J/ψ J/ψ pair can not be produced by single photon exchange due to charge parity conservation. Equally to photon the vector part of Z-boson vertex can not contribute to J/ψ J/ψ production. Also due to the charge parity conservation the axial part of Z-boson vertex does not contribute to J/ψ η c amplitude.
These selection rules were explicitly reproduced by our calculations and therefore provide the additional checks of our work.
Production of double heavy bound states is effectively described by NRQCD factorization.
The factorization formalism is designed to factor out the perturbative degrees of freedom and therefore separate the production mechanism into hard and soft subprocesses. The physics reasoning lies in separation of short distance and long distance interactions. If we stand by scales hierarchy m c >> m c v, where v velocity of c-quark in charmonium then short distance interaction corresponds to perturbative production of cc-pairs and long distance interaction can describe the bound state dynamics separately.
In this study we follow so-called δ-approximation neglecting the internal motion of quark in quarkonium and put equal the velocities of c-andc-quarks before the projection onto the bound state Ψ cc . Four final c-quarks are fixed two by two and form either two vectors (P 1 and P 2 ) or vector and pseudoscalar (P and Q correspondingly): The bound states momenta P, Q are kept on mass shells with equal masses: The factorized matrix elements have the following form: where A J/ψ ηc = A µ (P, Q) ε µ (P ) , A J/ψ J/ψ = A µν (P 1 , P 2 ) ε µ (P 1 ) ε ν (P 2 ) and R J/ψ, ηc (0) are radial wave functions in origin.
It should be clarified that we describe bound states in colour singlets only. We apply the projection technique to form the bound states. At the lowest order by v spin sums v λ 1 (pc)u λ 2 (p c ) ( 1 2 , λ 1 ), ( 1 2 , λ 2 )|S, s z over λ 1 , λ 2 can be expressed in terms of projector operators: where m = 2m c mass of charmonium. These operators are placed in the vertices of bound states and close the fermion lines into traces. Figure 1 shows the diagrams for e + e − → J/ψ η c process with the corresponding projectors. At next-to-leading order each of projectors can merge quark-antiquark pair either of one fermion line or from two different fermion lines.  An important feature of the studied challenge is absence of corrections for real gluon radiation (since we hold two colour-singlet final states). Thereby complete analysis of QCD corrections involves interference between LO and NLO amplitudes as well as interference between intermediate γ and Z. The full square of amplitude is of order of O(α 2 α 3 S ) and comprises the following seven terms: The renormalization procedure is organized by building counter-terms from the leading order diagrams. The renormalization constants are listed in (7)- (9). "On shell" scheme is fixed for mass and spinors renormalization and M S scheme is fixed for coupling constant. At next-toleading order we start the automatic calculation ofÃ N LO with physical spinors, masses and coupling constant. Final c-quarks are kept on mass shells: p 2 c = m 2 c . The isolated singularities are further cancelled with singular parts of A CT so that A N LO =Ã N LO + A CT remains a finite expression for renormalized amplitude:

CALCULATION DETAILS
Working with Feynman diagrams one can factorize the annihilation process and firstly consider the Z-boson's decay into two cc bound states (see Figure 2) as well as the same decay of virtual photon and then square the matrix element along with electric current. There are 4 diagrams at leading-order and 86 diagrams with one loop for Z * decay and the same number of diagrams for γ * decay. The diagrams and the corresponding analytic expressions are generated with FeynArts-package [29] in Mathematica. The computation strategy is based on the following toolchain: FeynArts → FeynCalc [30] (FeynCalcFormLink [31], TIDL) → Apart where ε αβσρ is either 4-or D-dimensional. It is checked that calculation result does not depend on ε αβσρ dimension in (10). The choice of ε αβσρ dimension slightly affects the traces evaluation process but it has no effect on the renormalized amplitudes.
The set of evaluated diagrams contains ones with triangle loops (see diagram 4 in the Figure 2). These diagrams do not contribute to any of examined processes because of P -parity violation and assume the verification test for calculations. The contribution from diagrams with two distinct traces projection (see diagram 3 in the Figure 2) is relevant only for J/ψ η c production. However it is pretty small (about 3% of total amplitude). We treat bubble loops (see diagram 8 in the Figure 2) massive for t-, b-and c-quarks while masses of light quarks u, d, s are neglected. One of the additional checks is that we explicitly obtain zero for prohibited process of η c η c production at both LO and NLO levels.
At next-to-leading order loop integration is carried out after projection onto the bound states. The δ-approximation used implies that momenta of relative motion of quark-antiquark pairs do not appear in loop integrals. This approach significantly reduces the computation cost.
It is worth mentioning that usage of FIRE reduction raises terms ∼ 1 D−4 in amplitudes. Therefore one should carefully execute the processing of master integrals the terms O(ε) in masters asymptotics might contribute to the finite part of amplitude. However in this study the amplitude terms ∼ 1 D−4 cancel each other (unlike the study [27]). After FIRE we treat only one-, two-and three-point integrals A 0 , B 0 , C 0 . The infinite part coming from divergent masters A 0 , B 0 carries poles O(1/ε) only. Working with automatic tools we do not distinguish ε IR and ε U V rather handle them together: ε IR = ε U V = ε.  In the presented calculations the strong coupling constant is taken with two loops accuracy: where L = ln (Q 2 /Λ 2 ); β 0 = 11 − 2 3 N f and β 1 = 102 − 38 3 N f with N f = 6; the reference value is α S (M Z ) = 0.1185. We have chosen the same scale for renormalization and coupling scales: The fine structure constant is fixed in Thomson limit α = 1/137. Numerical values for all the rest parameters are shown in Table I.

RESULTS
We present the analytic expressions for leading order cross sections in Equations (12) to (15) and highlight the contributions from γ annihilation, Z annihilation and interference between γ and Z. At next-to-leading order expressions are quite large so it is more reasonable to introduce the reference values (see Table II). where: a γZ = tan 2 θ w (3 csc 4 θ w − 20 csc 2 θ w + 32) 16  It can be seen that NLO contribution significantly enhances the leading order values. The maximums for cross sections correspond to √ s ≈ 7 GeV for J/ψ η c and to √ s ≈ 9 GeV for J/ψ J/ψ. One can make sure that our results for J/ψ η c production at low energies reproduce ones of the earlier works [16,17,19,20]. Particularly we obtain σ J/ψ ηc ≈ 15.5 fb at B-factories  good agreement with results of the paper [23], where R was calculated within LO accuracy: It is interesting to note, that as it is seen from the comparison of the NLO study of exclusive production of quarkonium states is known to face problems in double logarithmic terms that specify the corrections magnitude at high energies. In this study we confirm the result of the previous researches [36,37] for the process e + e − γ − → J/ψ η c and obtain the double logarithmic terms in the expansion at √ s >> m: Also we demonstrate for the first time the same behaviour both for the processes e + e − Z 0 − → J/ψ η c and e + e − Z 0 − → J/ψ J/ψ.
In this context pair charmonia production differs from pair B c production. As it is shown in paper [27] NLO QCD corrections to e + e − → B are stable with energy increasing. This fact requires separate consideration, which will be carried out in the next study.
We provide the Figures 7 and 8 to demonstrate the importance of account for Z 0 exchange in J/ψ η c production. At Z 0 mass σ(Z * + γ * )/σ(γ * ) ≈ 60. At energies away from Z 0 mass deviation of this ratio from one is related to the sign of γ −Z interference member (see Figure 8).
It is worth to mention that P -symmetry is not violated in the studied challenge. Although P -asymmetry could be induced by V − A interference in Z-quarks coupling it has no place in the particular case due to selection rules. Charge parity conservation lefts only one part of V − A structure to act (either vector or axial vector). Both of them are not able to violate P -symmetry alone.  for e + e − → J/ψ η c at next-to-leading-order. Besides the interest in studying double charmonia production in e + e − collisions we can not fail to refer to searches for Z → J/ψ J/ψ decays at LHC. Now the study of Z decays to double quarkonia states is motivated by the first CMS search [38] where Higgs and Z boson decays to J/ψ and Υ pairs are examined in the four-muon final state. In this way our work complements predictions of [23] showing that Γ(Z → J/ψ J/ψ) and Γ(Z → J/ψ η c ) should be approximately 3.5 times larger at next-to-leading order of NRQCD.

CONCLUSIONS
We calculate cross sections for associative J/ψ η c and J/ψ J/ψ production in e + e − annihilation at next-to-leading order of NRQCD. We consider single boson annihilation including γ exchange, Z-boson exchange and γ − Z interference. It is shown that at energies near Z-bosons mass as well as at higher energies annihilation with photon only becomes insufficient. QCD corrections O(α S ) significantly enhance the leading-order cross sections. At energies below 2M Z cross sections are enhanced up to 5 times. In Z-pole we obtain σ N LO ≈ 3.5 σ LO . The same enhancement in 3.5 times applies to widths of decays Z → J/ψ J/ψ and Z → J/ψ η c .
The results performed in the paper might be relevant for future study of charmonia physics at ILC or FCC colliders. In addition they are directly related to searches for rare decays of Z to double quarkonia states at LHC.