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$

The cross section values of $J/\psi\ \eta_c$ and $J/\psi\ J/\psi$ production in $e^+e^-$ annihilation are estimated within one-loop approximation near $Z$ pole, as well as at higher energies. Both intermediate bosons, $\gamma$ and $Z$, are taken into account. It is shown that at $Z$ mass the NLO contribution increases the cross-section values by 3.5 times.


INTRODUCTION
An associative production of J/ψ η c in e + e − annihilation was studied experimentally in details by Belle and BaBar Collaborations at the energy around 10.6 GeV [1,2].It turned out that the cross-section value predicted within the LO order approximation was by a factor of more than 5 smaller than the experimental measurement.This significant difference has initiated an intensive study of various corrections to the LO mechanism.
Two sources of corrections were mainly concerned.The first one is the internal motion of the quarks inside quarkonium (see the Ref. [3], which initiated the discussion on this contribution and then been followed by Refs.[4][5][6][7][8][9][10][11][12][13][14][15]).The second source of enhancement is due to QCD loop corrections and the one-loop corrections have already been determined [16,17].Nowadays, the corrections are known up to two-loops accuracy at the energy of B-factories [18].There are results, where both types of corrections were considered and applied [19,20].Evaluating the theoretical results in this field, one tends to believe that indeed both mechanisms are required to correctly describe the data.
Future charmonia studies are in plans of two big projects, International Linear Collider (ILC) and Future Circular Collider (FCC).Both projects will make available e + e − collisions at Z mass and above, specifically the energy range announced for FCC extends from √ s = 90 GeV to 400 GeV, while the collision energy of ILC has to be tuned to 250 GeV.Furthermore, the studies of Z-boson decays into two charmonia are certainly of an interest for the running LHC experiments.The J/ψ η c and J/ψ J/ψ pair production near Z-boson pole were calculated in LO approximation [21,22].Moreover, Z-boson decay mode to two charmonia was studied in the framework of a light cone formalism, what allowed to include the internal motion of quarks inside charmonium [23].Complementing the aforementioned results, another type of the oneloop correction is introduced in this report, specifically the J/ψ η c and J/ψ J/ψ pair production in the e + e − annihilation is computed up to the one-loop accuracy including both Z boson and photon, The theoretical study on J/ψ J/ψ pair production at e + e − collisions via a double photon exchange has been performed [24,25].Unfortunately, the attempts to make the calculations of the process at B-factories had no progress [1,2,26].In the present report, the e + e − annihilation via a single boson only is concerned.
The report follows up the previous research on B c -pair production in e + e − annihilation [27].
Meanwhile, the B c -pair production via γγ-fusion has also been covered in Ref. [28].

THE METHOD
The discussed production of the charmonium pair via a single boson exchange is affected by several selection rules.First of all, neither photon nor Z may decay to two identical η c mesons, since the η c pair must be in a P-wave state with a symmetric wave function, what is impossible.
Second, J/ψ J/ψ pair can not be produced by a single photon exchange due to a charge parity conservation.Similarly to the photon case, the vector part of the Z-boson vertex does not contribute to the J/ψ J/ψ production.Third, due to the charge parity conservation the axial part of the Z-boson vertex does not contribute to the J/ψ η c amplitude.The listed above selection rules are explicitly reproduced in the calculations presented below, hence providing the additional verification of the procedure.
The production of double heavy bound states is effectively described by the NRQCD factorization [29].The factorization formalism is introduced to factor out the perturbative degrees of freedom, therefore to separate the production mechanism into hard (short distance) and soft (long distance) subprocesses.Given the fact that m c >> m c v, where v is the velocity of c-quark in charmonium, the short distance interaction corresponds to the perturbative part of cc-pair production, whereas the long distance interaction describes the bound state formation and dynamics.
In our computations of the J/ψ η c production matrix elements, we start from the matrix element of four heavy quark production e + e − → c(p c )c(p c)c(q c )c(q c) with heavy quarks and antiquarks defined on their mass shells: As we assign v = 0 before the projection onto the bound states, the momentum P of the vector charmonium and the momentum Q of the pseudoscalar charmonium are related with the heavy quark momenta as follows below, To construct the bound states, we replace the spinor products v(p c)ū(p c ) and v(q c)ū(q c ) by the appropriate covariant projectors for color-singlet, spin-singlet and spin-triplet states as per where m = 2m c , ε is the polarization of the J/ψ meson, satisfying the following constraints: In exactly the same way, we express the matrix element of the two vector charmonia, denoting their momenta by P 1 and P 2 , and their polarizations by ε 1 and ε 2 .
The operators (2) close the fermion lines into traces.The contributions of LO diagrams to the amplitude always contain only one trace, while the NLO contributions contain one or two traces as illustrated at Figure 1.
The factorized matrix elements have the forms specified below in Eqs. 3 and 4, A e + e − → J/ψ(P 1 ) J/ψ(P where M µ (P, Q) ε µ and M µν (P 1 , P 2 ) ε 1 µ ε 2 ν are the hard production matrix elements of the two quark-antiquark pairs, projected on the quark-antiquark states with zero relative velocities and the appropriate quantum numbers using projectors (2); and R J/ψ (0), R ηc (0) are the radial wave function values at origin.The added caveat is that a real gluon radiation does not contribute to the NLO corrections in the studied processes, since within the applied approximation both heavy quark pairs are produced in the color singlet states.Thereby the NLO QCD corrections include only the contribution of the interference between the LO amplitudes and the one loop amplitudes.The total squared amplitude is of the order of O(α 2 α 3 S ).It contains the following seven terms, specifically The so-called "on shell" scheme is used for renormalization of masses and spinors and M S scheme is adopted for coupling constant renormalization as per where C = 4πµ 2 m 2 e −γ E and γ E is the Euler constant.The counter-terms are obtained from the leading order diagrams as follows below, The NLO amplitude ÃNLO has been calculated using the physical spinors and masses, as well as the physical value of coupling constant.The isolated singularities are further canceled with the singular parts of A CT so that A N LO = ÃNLO +A CT remains finite for the renormalized amplitude.

CALCULATION DETAILS
For technical reasons, we calculate separately the amplitudes of the e + e − fusion into the virtual Z boson and photon, and consequently the amplitudes of the Z-boson and photon transitions into charmonia.
There are 4 LO diagrams and 86 one loop diagrams for the Z * decay and the same number of diagrams for the γ * decay according to Figure 2. The diagrams and the corresponding analytic expressions are generated with the FeynArts-package [30] in Wolfram Mathematica.
The following toolchain is used for the computations: FeynArts → FeynCalc [31] (FeynCalcFormLink [32], TIDL) → Apart only.The $Apart function does the extra simplification of the amplitudes by partial fractioning of IR-divergent integrals.Finally, the FIRE package provides the complete IBP reduction of the amplitudes to master integrals, using the strategy mostly based on the Laporta algorithm [36].
The master integrals are then evaluated by substitution of their analytical expressions using the X-package.The computations are performed analytically, and the numerical values of the parameters are substituted only at the last step.
The conventional dimensional regularization (CDR) scheme with D-dimensional momenta (loop and external) and Dirac matrices are adopted for the calculations.The so-called naive interpretation of γ 5 is applied: γ 5 matrices anticommute with all other matrices and therefore, they are cancelled out in traces with an even number of γ 5 .In traces with an odd number of γ 5 matrices the remaining γ 5 are moved to the right and replaced per Eq. 10, where ε αβσρ is either 4-or D-dimensional Levi-Civita tensor.It is checked, that results of calculations do not depend on the dimension of Levi-Civita tensor in (10).The choice of its dimension has slightly affected the traces evaluation process, but it has no effect on the renormalized amplitudes.
It should be noticed that the diagrams with triangle quark loops, e.g. the diagram 4 in Figure 2, do not contribute to the J/ψ η c production due to the C-parity conservation, and we can directly verify this condition in our calculations.
The diagrams with two heavy quark traces, e.g. the diagram 3 in Figure 2, add about 3% to the total cross section of J/ψ η c production and do not contribute to the J/ψ J/ψ production cross section at all.
As explained earlier, the η c η c pair can not be produced in the photon or Z decays, indeed this selection rule is directly reproduced in our calculations at both LO and NLO levels.
After the FIRE reduction only one-, two-and three-point integrals of types A 0 , B 0 , and C 0 are left in the amplitudes.Some integrals of types A 0 and B 0 do contribute to the amplitude with the singular coefficient 1 D−4 .Therefore, in general case, one should carefully evaluate these specific integrals: terms proportional to O(ε) in the master integral expansion may contribute to the finite part of the amplitude.However in the considered processes the ∼ 1 D−4 terms cancel each other contrary to the case of B c pair production [27].The infinite amplitude parts coming from the divergent master integrals A 0 and B 0 contain only O(1/ε) poles.Working with automatic tools we do not distinguish ε IR and ε U V :  In the presented calculations we use the strong coupling constant evaluated with two loops accuracy as per 1179.N f = 5 is chosen for the interaction energies above 30 GeV, and N f = 4 is chosen for the interaction energies below 30 GeV.
The same value for the renormalization scale and for the coupling scale is used, Calculating the loop amplitudes we assume that u-, d-and s-quarks are massless.The fine structure constant is fixed in the Thomson limit α = 1/137.The numerical values of other parameters are oulined in Table I.
Since the analytical expressions for NLO cross sections are too cumbersome, we avoid showing ones in the text while presenting only the numerical values of the cross sections at several energies.The results are encapsulated in the Table II below.
On the contrary, the analytical expressions for leading order cross sections are quite simple and we present them below, highlighting the contributions of the γ * decay, the Z * decay and the interference between γ * and Z * as per Eqs.11-12, where Table II: The cross section values within the NLO approximation for different collision energies and renormalization scales.Our results for J/ψ η c production at low energies do reproduce the results of earlier studies [16,17,19,20], what is truly essential.
It can be seen in Figures 3 and 4, that the cross sections of both studied processes have maximum at √ s ∼ 7 ÷ 8 GeV.
As J/ψ J/ψ production proceeds only through the virtual Z, it is expected that near the threshold such a process will be strongly suppressed in comparison with the J/ψ η c production process.Our calculations are in agreement with this expectation: the discussed suppression is of the order 10 −6 at the energies below 10 GeV and it lessens with the increasing energy.Since Z-boson exchange dominates the area around Z pole, the cross sections of J/ψ J/ψ and J/ψ η c production at the corresponding energy area could be similar, in contrast to the e + e − → J/ψ J/ψ (blue curves) as a function of interaction energy: the NLO approach (solid curves) versus the LO approach (dashed curves).
The production cross section values are shown for the interaction energies above 30 GeV. case of the production near the threshold.This consideration is confirmed by our results: as seen in Figure 5 and Table II the cross section values are quite close to each other near Z pole.Both LO and NLO computations of the ratio r = σ J/ψ ηc /σ J/ψ J/ψ at the energy √ s = M Z are in a good agreement with the earlier results [23], where the width ratio R = Γ (Z → J/ψ η c ) /Γ (Z → J/ψ J/ψ) was calculated within the LO accuracy: Also, it is interesting to note, that as it is seen from Figures 4 and 5, the cross sections for the J/ψ J/ψ production at the energies of B-factories and at the energies near Z pole are comparable in a magnitude.
The NLO calculations of the exclusive production of quarkonium states are known to encounter the problems related to the double logarithmic terms at high energies.In this study we confirm the result of the previous studies [37,38] done for the process e + e − γ − → J/ψ η c .We  have obtained the double logarithmic terms in the expansion at √ s >> m as per Also we demonstrate for the first time the same behaviour both for the processes e + e − Z − → J/ψη c and e + e − Z − → J/ψ J/ψ.The presented case is different from the B c -pair production one.As it is shown in [27] the relative contribution of NLO QCD corrections to the e + e − → B does not increase with the increasing energy.This fact requires an additional study, which we plan to perform in the future works.
Asymptotically the cross sections fall off with the increase of the energy: the LO contributions decrease as O (1/s 4 ) while the NLO decrease as O(ln 2 s/s 4 ). Figure 6 demonstrates that the ratio σ N LO /σ LO increases with energy.It seems that the specific behaviour of the NLO corrections can not be compensated by the scale choice at very high energies, as explained in Ref. [37,38].As a result, at the energies about 2M Z the NLO contribution is responsible for a fivefold increase in the cross section.
Obviously, the Z-boson exchange dominates in the J/ψ η c production around the Z mass, as it is seen in Figures 11 and 12.At the Z mass, the ratio σ(Z * + γ * )/σ(γ * ) amounts ≈ 60.
As one moves away from the Z pole, the contribution of the Z-boson exchange diminishes in such a way that the ratio σ It is worth to mention that the P -symmetry is not violated in the considered processes, because the V − A interference, which could cause such a violation, does not contribute either to the J/ψ η c production, or to the J/ψ J/ψ production.
Discussing the process e + e − Z − → J/ψ J/ψ, it is natural to suggest a significance for the search of Z → J/ψ J/ψ decays in LHC detectors.Currently, the studies of Z decays to double quarkonia states are motivated by the CMS study [39], where the search for Higgs and Z decays to J/ψ and Υ pairs was performed for the first time.In this way our work complements the   predictions of [23] and predicts that the width Γ(Z → J/ψ J/ψ) and Γ(Z → J/ψ η c ) are approximately 3.5 times larger at the NLO approximation.

CONCLUSIONS
The cross sections of J/ψ η c pair and J/ψ J/ψ pair production in the e + e − single boson annihilation are calculated within the QCD one loop approximation with the γ exchange, the Z-boson exchange and the γ − Z interference considered.It is found that the one loop QCD corrections are responsible for a significant, up to a fivefold, increase of the cross section values at all investigated energies.It is obtained, that σ N LO /σ LO ≈ 3.5 at Z pole for both investigated processes.Obviously, the same enhancement by a factor of 3.5 applies to the widths of decays Z → J/ψ J/ψ and Z → J/ψ η c .
The results obtained in the paper might be useful for future studies of charmonia physics at ILC and FCC colliders.Furthermore, the results are directly related to the searches of rare

Figure 1 :
Figure 1: The examples of diagrams for e + e − → J/ψ η c process with one and two traces: the LO diagram (a); the NLO diagrams (b) and (c).

Figures 3 -
Figures 3-6 demonstrate the energy dependence of the production cross sections calculated at the scale µ = √ s.As clearly seen in the presented Figures, the NLO contributions do significantly increase the cross section values.To estimate the theoretical uncertainties due to the scale choice we vary the µ value from √ s/2 to 2 √ s and present the results in Figures 7-10.

Figure 3 :
Figure 3: The production cross sections for the process e + e − → J/ψ η c calculated within the LO approximation (dashed curve) and within the NLO approximation (solid curve) as a function of interaction energy.The production cross section values are shown for the interaction energies below 30 GeV.

Figure 4 :
Figure4: The same as in Figure3, but for the process e + e − → J/ψ J/ψ.

Figure 5 :
Figure 5: The production cross sections for the processes e + e − → J/ψ η c (red curves) and

Figure 6 :
Figure 6: The ratio σ N LO /σ LO for the process e + e − → J/ψ η c (red curve) and e + e − → J/ψJ/ψ (blue curve) as a function of interaction energy.

s / 2 < μ < 2 Figure 7 :
Figure 7: The dependence of NLO cross section on the interaction energy for the process e + e − → J/ψ η c at different scale values: √ s < µ < 2 √ s.The production cross section values are shown for the interaction energies below 30 GeV.

Figure 8 :
Figure8: The same as in Figure7, but for the process e + e − → J/ψ J/ψ.

Figure 9 :Figure 10 :
Figure 9: The dependence of NLO cross section on the interaction energy for the processes e + e − → J/ψ η c (purple) and e + e − → J/ψ J/ψ (pink) at different scale values: √ s < µ < 2 √ s.The production cross section values are shown for the interaction energies above 30 GeV.

Figure 12 :
Figure 12: The cross sections ratios for the process e + e − → J/ψ η c at NLO as a function of interaction energy.

Z
-boson decays into double quarkonia states in LHC detectors.Authors would like to thank I. Gorelov, A. Onishchenko and A. Luchinsky for help and constructive discussions.The work was supported by the RFBR (grant No. 20-02-00154 A).I. Belov acknowledges the support from the "Basis" Foundation (grant No. 20-

Table I :
The parameter values.