Improved analysis of double $J/\psi$ production in $Z$-boson decay

In this paper, we present an improved calculation for the decay rate of the rare $Z$-boson decay into $J/\psi + J/\psi$. This decay is dominated by the photon fragmentation mechanism, i.e., the transition $Z\to J/\psi + \gamma^{*}$ followed by the fragmentation $\gamma^{*}\to J/\psi$. In our calculation, the amplitude of $\gamma^{*}\to J/\psi$ is extracted from the measured value of $\Gamma(J/\psi \to e^+ e^-)$, and the amplitude of $Z\to J/\psi + \gamma^{*}$ is calculate through the light-cone approach. The higher-order QCD and relativistic corrections in the amplitude of $\gamma^{*}\to J/\psi$ and the large logarithms of $m_{_Z}^2/m_c^2$ that appear in the amplitude of $Z\to J/\psi + \gamma^{*}$ are resummed in our calculation. Besides, the non-fragmentation amplitude is calculated based on the NRQCD factorization, and the next-to-leading order QCD and relativistic corrections are included. The obtained branching fraction for this $Z$ decay channel is $8.66 ^{+1.48} _{-0.69}\times 10^{-11}$.


I. INTRODUCTION
The rare Z-boson decay into double J/ψ provides a good platform for studying the heavy quarkonium production and for searching of new physics beyond the standard model (SM).In 2019, the CMS collaboration presented the first search for the decay Z → J/ψ+J/ψ based on the data from pp collisions with an integrated luminosity of 37.5 fb −1 at √ s = 13 TeV at the large hadron collider (LHC).The observed upper limit on the branching fraction for the decay is Br(Z → J/ψ+J/ψ) < 2.2×10 −6 [1].In 2022, the CMS collaboration updated the upper limit on this branching fraction by using a data set of pp collisions corresponding to an integrated luminosity of 138 fb −1 .Their updated upper limit is Br(Z → J/ψ +J/ψ) < 1.4×10 −6 [2].Although this decay process has not been observed up to now, the prospect of observing this decay at future colliders is optimistic [3].As is well-known, there are several lepton colliders, e.g., ILC [4], CEPC [5], FCC-ee [6], and super Z factory [7], are under consideration.These lepton colliders are planed to operate at the Z pole for a period of time, and a large number of Z bosons are expected to be produced, e.g., at the FCC-ee, about 5 × 10 12 Z bosons will be produced [6].Therefore, these lepton colliders provide new opportunities for studying the rare decay Z → J/ψ + J/ψ.
On the theoretical side, the decay Z → J/ψ + J/ψ was first studied in Ref. [8] in 1990.In Ref. [9], this decay was reanalysed at leading order (LO) in α s and v (v stands for the typical velocity of the c quark or the c quark in the rest frame of J/ψ) based on the nonrelativistic QCD (NRQCD) factorization.In these references, only the contribution of the LO QCD Feynman diagrams were considered.In Ref. [10], Gao and Gong pointed out that besides the LO QCD diagrams, the other transition via Z → J/ψ + γ * followed by the fragmentation γ * → J/ψ can also bring significant contribution.Their calculation shows that the branching fraction of Z → J/ψ + J/ψ can be enhanced to 1.1 × 10 −10 from 1.5 × 10 −13 after including the single photon fragmentation contribution.Recently, the next-to-leading order (NLO) QCD correction to the decay rate of Z → J/ψ + J/ψ has been calculated in Refs.[11,12].The results there show that the NLO QCD correction can diminish the decay rate significantly [11].
In this paper, we devote ourselves to presenting an improved calculation for the decay width of Z → J/ψ+J/ψ.As the photon fragmentation contribution dominates this decay process, we first focus on the amplitude for the photon fragmentation diagrams.We express the amplitude of the photon fragmentation diagrams into two parts: one for Z → J/ψ + γ * , and the other for γ * → J/ψ.For the amplitude of γ * → J/ψ, we directly extract it from the experimentally measured leptonic decay width of the J/ψ, rather than calculating it order by order in α s and v. Using this approach, the QCD and relativistic corrections to the amplitude of γ * → J/ψ are resummed up to all orders.For the amplitude of Z → J/ψ + γ * , we adopt the light-cone approach suggested by Refs.[13,14], in which the amplitude is expressed as the convolution of the hard-scattering kernel with the J/ψ light-cone distribution amplitude (LCDA).The hard-scattering kernel can be calculated with the perturbative QCD (pQCD), and the LCDA for J/ψ can be calculated through the NRQCD factorization [15].Under the light-cone approach, logarithms of m 2 Z /m 2 c appearing in the amplitude can be resummed through solving the evolution equation of the LCDA.In addition to the contribution of the photon fragmentation diagrams, we also calculate the contribution from the non-fragmentation diagrams.We shall adopt the fixed-order approach under the NRQCD factorization, which is adopted in previous calculations on the decay rate of Z → J/ψ + J/ψ, to calculate the nonfragmentation amplitude.
The rest of the paper is organized as follows.In Sec.II, we give the calculation formalism for the photon fragmentation amplitude, the non-fragmentation amplitude, and the decay width of Z → J/ψ + J/ψ.In Sec.III, we present the numerical results and discussions.Section IV is reserved as a summary.The dominant mechanism of the decay Z → J/ψ+J/ψ is via the photon fragmentation, i.e., the transition Z → J/ψ + γ * followed by γ * → J/ψ.To understand this mechanism, we present two typical LO Feynman diagrams for the photon fragmentation mechanism in Fig. 1 1 .The amplitude for the photon fragmentation mechanism can be written as where the factor −i/m 2 J/ψ arises from the photon propagator that will fragment into a J/ψ.iM γ * →J/ψ,µ denotes the amplitude for the transition of a virtual photon into the J/ψ, and it can be expressed as [16,17] where J µ (x) is the electromagnetic current, 1 It is noted that we shall adopt an improved approach where some important higher-order corrections are resummed, rather than the conventional fixed-order approach.Hence, Fig. 1 is only used to illustrate the photon fragmentation mechanism, we do not use it in our calculation.
Here, the sum extends over all quark flavors.The J/ψphoton effective coupling constant g J/ψγ can be determined using the decay width of the J/ψ to an electronpositron pair [16], i.e., Γ(J/ψ → e + e − ) = 4πα 2 (m J/ψ )g 2 J/ψγ 3m 3 The method of calculating the amplitude iM γ * →J/ψ,µ based on Eqs.( 2), ( 3) and ( 4) has an important advantage compared to the fixed-order calculation based on the NRQCD factorization.Under the NRQCD factorization, an amplitude is expressed as a double expansion in α s and v.The calculations of higher-order QCD and relativistic corrections are very difficult.On the other hand, the calculation of the amplitude iM γ * →J/ψ,µ based on Eqs.( 2), ( 3) and ( 4) already contains the QCD and relativistic corrections to all orders.
Next, let us demonstrate the calculation of the amplitude M µ Z→J/ψ+γ * .We shall adopt the approximation M µ Z→J/ψ+γ * ≈ M µ Z→J/ψ+γ , which will only lead to an error of O(m 2 J/ψ /m 2 Z ).The decay Z → J/ψ + γ has been studied extensively based on the NRQCD and light-cone approaches [18][19][20][21][22][23].In this work, we adopt the light-cone approach to calculate the amplitude M µ Z→J/ψ+γ .Under the light-cone approach, the amplitude of Z → J/ψ + γ is expressed by [21] iM Z→J/ψ+γ,µ = iA ǫ ξµνρ ǫ ξ Z ǫ * ν J/ψ p ρ γ , where Here, e c = +2/3 is the electric charge of the charm quark, f J/ψ is the decay constant of the longitudinally polarized J/ψ, T H (x, µ) is the hard-scattering kernel, φ J/ψ is the longitudinally polarized LCDA of the J/ψ, ǫ Z is the polarization vector of the Z boson, ǫ * J/ψ is the polarization vector of J/ψ, and p γ is the momentum of the photon, µ is the factorization scale, and x is the momentum fraction of the c quark in the J/ψ.The couplings g Z and g c A are defined as , where G F is the Fermi constant and t c 3L = 1/2 is the weak isospin for the left-handed c quark.
The hard-scattering kernel T H (x, µ) can be calculated through pQCD, and it has been calculated up to order α s in Ref. [20].The expression of the hard-scattering kernel, up to order α s , is where Here, C F = (N 2 c − 1)/(2N c ) with N c = 3 is the quadratic Casimir operator.In order to avoid large logarithms appearing in the hard-scattering kernel, we set the factorization scale as µ = m Z in Eq. (6).
The longitudinally polarized LCDA φ J/ψ is defined as [24,25]: where p is the quarkonium momentum, the LCDA is normalized as The coordinate z lies along the plus light-cone direction, and the gauge link [z, 0] that makes the definition of the LCDA gauge invariant is where P denotes the path ordering.Under the NRQCD factorization, the LCDA φ J/ψ is expressed as a double expansion in α s and v [20,26], i.e., Here, in order to avoid large logarithms appearing in the LCDA, we set the initial factorization scale of the LCDA as where the covariant derivative operator ψ † ↔ Dχ ≡ ψ † Dχ− (Dψ) † χ, ψ and χ are Pauli spinor fields that describe c quark annihilation and c quark creation, respectively.
To solve the ERBL equation, i.e., Eq.( 18), it is convenient to expand the LCDA in terms of the Gegenbauer polynomials C (3/2) n (2x − 1), which are the eigenfunctions of the LO evolution kernel.The Gegenbauer expansion of the LCDA is (20) where a n (µ) is the n th Gegenbauer moment of the LCDA φ J/ψ (x, µ), and where the normalization factor Solving the ERBL equation in the Gegenbauer moment space, the obtained the evolved Gegenbauer moments have the following form: The expressions for U nk (µ, µ 0 ) up to next-leadinglogarithmic (NLL) accuracy can be found in Ref. [33].
Similarly, the hard-scattering kernel can also be expanded in terms of the Gegenbauer polynomials as The n th Gegenbauer moment Making use of the Gegenbauer expansions of the LCDA and the hard-scattering kernel and the orthogonality property of the Gegenbauer polynomials, we can express the convolution of the hard-scattering kernel and the LCDA in Eq.( 6) as As stated below Eq.( 8), we have set µ = m Z so that the hard-scattering kernel T H (x, µ) does not contain large logarithms.The large logarithms of m 2 Z /m 2 c will appear in the LCDA φ J/ψ (x, µ = m Z ).In our calculation, these large logarithms are resummed up to NLL accuracy through the evolution of the LCDA from the initial scale µ 0 = m c to the scale µ = m Z .
It is found that the sum over n in Eq.( 26) is divergent [21].In order to address the nonconvergent problem in the Gegenbauer moment summation, the authors of Ref. [34] introduced the Abel-Padé method, which combines the Abel summation and the Padé approximant.
Under this method, the summation in Eq.( 26) is first written as an Abel summation as Then, we keep 20 nonzero terms in the second line of Eq.( 27), and generate a 10 × 10 Padé approximant.It is worth noting that the limit of z → 1 should be taken after the generation of the Padé approximant.
In addition to calculating the convolution in Eq.( 6), we also need to determine the value of the decay constant f J/ψ .Setting z to 0 in Eq.( 9) and using the normalization With this relation, we can relate the decay constant f J/ψ with the leptonic decay width of J/ψ [21], i.e., The decay constant f J/ψ can be determined through this relation straightforwardly.

FIG. 2. Six typical tree-level and one-loop non-fragmentation
Feynman diagrams that contribute to the decay Z → J/ψ + J/ψ.
In this subsection, we present the key formulas for calculating the non-fragmentation amplitude.Compared to the fragmentation amplitude, the non-fragmentation amplitude is significantly suppressed [10,11].We will calculate it directly based on the NRQCD factorization.
Under the NRQCD factorization [15], the nonfragmentation amplitude up to order-v 2 for the decay Z → J/ψ + J/ψ can be written as where c 0 and c 2 are short-distance coefficients (SDCs) that can be expanded in powers of α s .In this work, we calculate the SDC c 0 up to NLO in α s , and calculate the SDC c 2 at LO in α s .J/ψ|ψ † σ • ǫχ|0 is the LDME for J/ψ.v 2 J/ψ is a ratio of LDMEs, which has been defined in Eq. (12).The factor 2m J/ψ arises from that we use the relativistic normalization for the J/ψ in M nfr Z→J/ψ+J/ψ on the left-hand side of Eq.( 30), while we use conventional nonrelativistic normalization for the NRQCD matrix element on the right-hand side of Eq. (30).
The SDCs can be determined through the matching method.According to this method, we first replace the J/ψ state in Eq.( 30) by an on-shell (cc) pair with a small relative velocity v, i.e., where the amplitude M nfr can be calculated using perturbative theory in QCD, and the LDME (cc) 1 ]|ψ † σ • ǫχ|0 for an on-shell (cc) pair can be calculated using perturbative theory in NRQCD.Then the SDCs can be determined by comparing both sides of Eq.( 31) order by order.The SDCs are insensitive to the long-distance physical effects, thus the SDCs extracted from the on-shell (cc) pair production can be applied to the J/ψ production.
To calculate the amplitude M nfr up to order α s and order v 2 , we assign the momenta of the c and c quarks as p i = 1 2 P i + q i and pi = 1 2 P i − q i , where i = 1, 2. Then we have the relations In the rest frame of a (cc) pair, the total and relative momenta of the (cc) pair are P i = (2E q , 0) and q i = (0, q i ) where E q = q 2 i + m 2 c .In deriving these relations, we have implicitly adopted the relation q 2 1 = q 2 2 = m 2 c v 2 .It is convenient to adopt the covariant projection method to enforce the produced (cc) pairs to be the corresponding spin and color state.The spin projector for the spin-triplet state is [35 where ǫ * is the polarization vector of the spin-triplet state.The color projector for color-singlet state is where 1 stands for the unit matrix of the SU (3) c group.
To obtain the amplitude in terms of powers of v, we first expand the amplitude in powers of q i to the required order, where is the non-fragmentation amplitude, which has been abbreviated as M on right hand side of Eq.( 35).To project out the S-wave contribution, we make the replacement where The odd-power terms of q i vanish for the S-wave contribution.Then we obtain the amplitude up to order v 2 , In the calculation, we expand the amplitude in powers of the relative momenta q i before carrying out the loop integration.This amounts to calculating the contributions arising from the hard region in the language of the strategy of region [36].Consequently, the Coulomb divergences arising from the potential region do not appear in our calculation.
There are ultraviolet (UV) and infrared (IR) divergences in the loop-diagram calculation.We utilize dimensional regularization to regularize both the UV and IR divergences, with spacetime dimension D = 4−2ǫ.The UV divergences are removed through renormalization.We define the renormalized charm-quark field, charm-quark mass, and gluon field in the on-shell (OS) scheme, while define the renormalized strong coupling constant in the modified minimal subtraction (MS) scheme.The corresponding renormalization constants are Z i = 1 + δZ i , and where γ E is the Euler constant, µ r is the renormalization scale, β 0 = 11/3 − 2n f /3, β ′ 0 = 11/3 − 2n lf /3, n f is the number of active quark flavors and n lf = 3 is the number of light quark flavors.
In dimensional regularization, the γ 5 matrix should be noted.In this work, we employ the Larin scheme [37], where the product of a γ µ matrix and a γ 5 matrix is expreesed as: In this γ 5 scheme, a finite renormalization is required to restore the axial current Ward identity for the axial vector vertex.The corresponding finite renormalization constant Z 5 up to order α s for the axial vector vertex is given by: In order to obtain the SDCs using Eq.( 31), in addition to the amplitude M nfr , we also need to calculate the matrix element (cc) where λ denotes the helicity of the (cc) pair.
In the calculation of the non-framentation amplitude, the packages FeynArts [38] is adopted to generate Feynman diagrams, the package Feyncalc [39,40] is adopted to carry out the traces over the Dirac and color matrices, the package $Apart [41] is adopted to do partial fraction, and the package FIRE [42] is adopted to do integrationby-parts (IBP) reduction for the loop integrals.After the IBP reduction, the one-loop integral are reduced to mater integrals, and these master integrals are calculated numerically using the package LoopTools [43].

C. Decay width
The total amplitude for the decay Z → J/ψ + J/ψ is the sum of the fragmentation and non-fragmentation amplitudes, i.e., M Z→J/ψ+J/ψ = M fr Z→J/ψ+J/ψ + M nfr Z→J/ψ+J/ψ .(43) Then the decay width can be calculated through where the factor 1/2! comes from that there are two identical particles in the final state, 1/3 comes from the average over the polarizations of the Z boson, denotes that the polarizations of the initial and final states are summed over.

III. NUMERICAL RESULTS AND DISCUSSIONS
In the numerical calculation, the necessary input parameters are taken as follows: where m c stands for the pole mass of the charm quark, and its value is taken from Ref. [44].The values of m Z , m J/ψ , and G F are taken from the Particle Data Group (PDG) [45].The value of the running electromagnetic coupling α(m J/ψ ) is taken from Ref. [46].For the running strong coupling constant, we use the two-loop formula where β 1 = 102 − 38n f /3 is the two-loop coefficient of the QCD β function.According to α s (m Z ) = 0.118 [45], we obtain Λ With the determined values of Λ QCD , the strong coupling constant at an arbitrary scale (µ r ≫ Λ QCD ) can be directly calculated through Eq.( 46).
The values of g J/ψγ and f J/ψ can be extracted from the lepton decay width of J/ψ through Eq.( 4) and Eq.( 29), respectively.
Using the measured value Γ(J/ψ → e + e − ) = 5.55keV [47], we obtain The coupling constant g J/ψγ has a relative minus sign compared to f J/ψ , which can be realized by comparing Eq.( 2) and Eq.( 28).
We take the values for the LDME O 1 J/ψ = and the ratio of LDMEs q 2 J/ψ from Ref. [44], i.e., The value of v 2 J/ψ can be derived through v 2 J/ψ = q 2 J/ψ /m 2 c directly.

A. Basic results
this subsection, we present the numerical results for the decay width of the decay Z J/ψ + J/ψ.In orto have a glance on the magnitude of the photon fragmentation, non-fragmentation, and interference contributions to the decay width of Z → J/ψ + J/ψ, we first present the numerical results with the input parameters taken as their central values.A detailed analysis of the uncertainty of the decay width will be presented in the next subsection.Numerical results for the decay width are shown in Table I, where the photon fragmentation, nonfragmentation, interference, and total contributions are given explicitly.In the calculation, we have set the renormalization scale in the non-fragmentation amplitude as µ r = m Z .
From Table I, we can find that the largest contribution comes from the photon fragmentation, which accounts for 74.1% of the total contribution.The interference of the photon fragmentation amplitude and the non-fragmentation amplitude accounts for 24.8% of the total contribution, while the non-fragmentation contribution only accounts for 1.1% of the total contribution.Although the photon contribution is suppressed by the (α/α s )2 compared to the non-fragmentation contribution, the photon fragmentation contribution is two orders of magnitude lager than the non-fragmentation contribution.This is because the invariant mass of the photon propagator in the fragmentation diagrams is m J/ψ , while the invariant mass of the gluon propagator in the LO non-fragmentation diagrams is m Z /2, i.e., the fragmentation contribution is enhanced by a factor of (m Z /2m J/ψ ) 4 from the gauge boson propagator compared to the non-fragmentation contribution.

B. Uncertainty analysis
In this subsection, we present an estimate of the theoretical uncertainties for the decay width of Z → J/ψ + J/ψ.The main uncertainty sources include the factorization/renormalization scales, the charm quark mass, the leptonic decay width Γ(J/ψ → e + e − ), the LDME O 1 J/ψ , and the ratio q 2 J/ψ .In the following estimation, when we discuss the uncertainty from one parameter, other parameters will be kept to be their central values.
There are several factorization/renormalization scales involved in the calculation: the initial factorization scale µ 0 and the final factorization scale µ in the calculation of the photon fragmentation amplitude, as well as the renormalization scale µ r in the calculation of the nonfragmentation amplitude.We estimate the uncertainties by varying these scales in the ranges µ 0 ∈ [1GeV, 2m c ], µ ∈ [m Z /2, 2m Z ], and µ r ∈ [m Z /2, 2m Z ].For the uncertainty caused by the charm quark mass, we estimate it by taking m c = 1.4 ± 0.2 GeV.It is noted that the values for the LDME O 1 J/ψ and the ratio q 2 J/ψ obtained in Ref. [44] depend on the value of the charm quark pole mass.Therefore, when considering the uncertainty caused by the charm quark mass, we not only consider the dependence of the SDCs on the charm quark mass, but also the dependence of the LDME O 1 J/ψ and the ratio q 2 J/ψ on the charm quark mass.The dependence of the LDME and the ratio q 2 J/ψ on the charm quark mass has been given in the Ref. [44].In our calculation, the parameters g J/ψγ and f J/ψ are extracted from the measured value of Γ(J/ψ → e + e − ), and the LDME O 1 J/ψ and the ratio q 2 J/ψ obtained in Ref. [44] also depend on the measured value of Γ(J/ψ → e + e − ).Thus, Γ(J/ψ → e + e − ) is an important uncertainty source.For the uncertainty caused by the measured value of Γ(J/ψ → e + e − ), we take the uncertainties of the leptonic width as Γ(J/ψ → e + e − ) = 5.55 ± 0.14 ± 0.02 keV [47].Since the uncertainties caused by the charm quark mass and the measured value of Γ(J/ψ → e + e − ) have already been considered, we will omit the uncertainties caused by charm quark mass and the measured value of Γ(J/ψ → e + e − ) when considering the uncertainties caused by the LDME O 1 J/ψ and the ratio q 2 J/ψ .The uncertainties for the LDME O 1 J/ψ and the ratio q 2 J/ψ with omitting the uncertainties caused by the charm quark mass and the measured value of Γ(J/ψ → e + e − ) [44] are q 2 J/ψ = 0.441 +0.139 −0.138 GeV 2 . (52) Then we obtain the uncertainties arising from these uncertainty sources Γ Z→J/ψ+J/ψ = 2.16 +0.06+0.32+0.12+0.10+0.07−0.07−0.05−0.10−0.08−0.06× 10 −10 GeV, (53) where the first uncertainty is caused by the factorization/renormalization scales 2 , the second uncertainty is caused by the charm quark mass, the third uncertainty is caused by the measured value of Γ(J/ψ → e + e − ), the fourth uncertainty is cased by the LDME O 1 J/ψ and the last uncertainty is caused by the ratio q 2 J/ψ .Adding these uncertainties in quadrature, we obtain the total theoretical uncertainty for the decay width Γ Z→J/ψ+J/ψ = 2.16 +0.37 −0.17 × 10 −10 GeV.
C. Comparison with other theoretical calculation and current experiment limit In this subsection, we present a comparison of our calculation with previous theoretical calculation and the experimental upper limit on the decay Z → J/ψ + J/ψ.TABLE II.The branching fraction of the decay Z → J/ψ + J/ψ.Our result and the result from Ref. [11] are theoretical predictions, while the limit in the last column is the current experimental upper limit given by Ref. [2]  The theoretical predictions and the current experimental upper limit on the branching fraction of the decay Z → J/ψ + J/ψ are shown in Table II.From the table, we can see that our branching fraction for Z → J/ψ+J/ψ is smaller than that in Ref. [11] by about 22%, which is about −1.0σ in the uncertainties of Ref. [11].
In Ref. [11], the photon fragmentation, nonfragmentation, and the interference contributions were calculated up to NLO in α s based on the NRQCD factorization.It was found in Ref. [11] that the NLO QCD corrections are significant, e.g., the NLO QCD correction diminishes the photon fragmentation contribution by 20 − 40% with different choices of the renormalization scale.In this work, the photon fragmentation amplitude is calculated with an improved approach.More explicitly, the photon fragmentation amplitude is written as the product of the amplitude of Z → J/ψ + γ * and the amplitude of γ * → J/ψ.The amplitude of γ * → J/ψ is extracted from the measured value of Γ(J/ψ → e + e − ), and the amplitude of Z → J/ψ + γ * is calculated through the light-cone approach where the large logarithms of m 2 Z /m 2 c are resummed through the evolution of the J/ψ LCDA.Using this improved approach, the higher-order QCD and relativistic corrections in the amplitude of γ * → J/ψ and the large logarithms of m 2 Z /m 2 c in the amplitude of Z → J/ψ + γ * are resummed in our calculation.Moreover, in our calculation of the non-fragmentation amplitude, besides the NLO QCD corrections, the NLO relativistic corrections are also included.
From Table II, we can see that our branching fraction for Z → J/ψ + J/ψ is compatible with the current experimental upper limit established by the CMS Collaboration.

IV. SUMMARY
In this paper, we have presented a calculation for the decay rate of Z → J/ψ + J/ψ.The decay process Z → J/ψ + J/ψ is dominated by the single photon fragmentation mechanism, i.e., the transition Z → J/ψ + γ * followed by the fragmentation γ * → J/ψ.The amplitude of γ * → J/ψ has already been calculated up to O(α 3 s ) under the NRQCD factorization in Refs.[48,49].It seems that the perturbative expansion of the amplitude of γ * → J/ψ does not converge in those calculations.Moreover, the higher-order corrections of the amplitude of Z → J/ψ + γ * involve the large logarithms of m 2 Z /m 2 c .Considering these problems, we have presented an improved calculation for the decay rate of Z → J/ψ + J/ψ compared to the fixed-order calculation.In our calculation, the amplitude for the photon fragmentation is expressed as the product of the amplitude of Z → J/ψ + γ * and the amplitude of γ * → J/ψ.Then, the amplitude of γ * → J/ψ is extracted from the measured value of Γ(J/ψ → e + e − ), and the amplitude of Z → J/ψ + γ * is calculated by the light-cone approach where the large logarithms of m 2 Z /m 2 c are resummed through solving the evolution equation of the J/ψ LCDA.The advantage of our calculation is that the higher-order QCD and relativistic corrections in the amplitude of γ * → J/ψ and the large logarithms of m 2 Z /m 2 c in the amplitude of Z → J/ψ + γ * are resummed.For the non-fragmentation amplitude, we adopt the fixed-order approach under the NRQCD factorization.In our calculation of the nonfragmentation amplitude, the NLO QCD and NLO relativistic corrections are included simultaneously.
The obtained branching fraction for Z → J/ψ + J/ψ is 8.66 +1. 48 −0.69 × 10 −11 , which is compatible with the current experimental upper limit Br(Z → J/ψ + J/ψ) = 1.4 × 10 −6 given by the CMS collaboration [2].In recent years, several high-luminosity e + e − colliders, such as ILC, CEPC, FCC-ee, and Super Z factory are proposed.A large number of Z boson events are expected to be generated at these future e + e − colliders.For instance, about 5 × 10 12 Z bosons can be generated at the FCC-ee [6].According to the branching fraction obtained in this work, about 400 Z → J/ψ + J/ψ events are expected to be generated at the FCC-ee.One may expect that this decay can be studied at these future e + e − colliders.

TABLE I .
The decay width (in unit: 10 −12 GeV) of the decay Z → J/ψ + J/ψ.The photon fragmentation, nonfragmentation, interference, and total contributions are given explicitly.