# Small Batch Sizes Improve Training of Low-Resource Neural MT

Àlex R. Atrio1,2 Andrei Popescu-Belis1,2
1HEIG-VD / HES-SO 2EPFL Yverdon-les-Bains Lausanne Switzerland Switzerland
{alejandro.ramirezatrio<EMAIL_ADDRESS>
###### Abstract

We study the role of an essential hyper-parameter that governs the training of
Transformers for neural machine translation in a low-resource setting: the
batch size. Using theoretical insights and experimental evidence, we argue
against the widespread belief that batch size should be set as large as
allowed by the memory of the GPUs. We show that in a low-resource setting, a
smaller batch size leads to higher scores in a shorter training time, and
argue that this is due to better regularization of the gradients during
training.

##  1 Introduction

Training Transformers for low-resource neural machine translation (NMT), i.e.
when only small parallel corpora are available, raises the challenge of
finding optimal hyper-parameters. While several fixed configurations of the
Transformer (Vaswani et al., 2017) have been empirically validated by the
community, such as ‘Base’ or ‘Big’, the settings of many other hyper-
parameters rely on tips from practitioners. However, these values are not
always suitable to low-resource settings, and systematic studies in these
settings are rare (Araabi and Monz, 2020; Van Biljon et al., 2020).

In this paper, we show that the best values of a hyper-parameter that is
essential for training, namely _batch size_ , differ in low-resource settings
from those commonly accepted when larger data sets are available. We analyze
the role of small batch sizes, inspired by studies in computer vision (Keskar
et al., 2016), and then pinpoint empirically the optimal trade-off between a
high batch size (for efficiency) and a small one (for regularization).
Although large batch sizes were found to lead to higher-quality models in
experiments with high-resource NMT (Popel and Bojar, 2018; Xu et al., 2020),
we show here that smaller batch sizes can outperform the latter, likely due to
a regularizing effect in the gradient update. Moreover, we show that this
finding is invariant to changes in tokenization methods.

The paper is organized as follows. In Section 2, we discuss batch size from a
machine learning perspective, showing why smaller values of batch size may act
as regularizers. Then, in Section 3, we review studies of hyper-parameters in
NMT. In Section 4, we present the parameters of our Transformer and the data
from the WMT 2020 Low-resource task (Fraser, 2020) and other sources that we
use in our experiments. In Section 5, we provide empirical evidence that
smaller batch sizes are preferable in low-resource settings.

##  2 ML Perspective on Batch Size

Machine learning theory argues that performing back-propagation with large
batch sizes leads to better optimization, because the estimates of the
gradients are more accurate. Conversely, using small batches during training
leads to noisier gradient estimations, i.e. with a larger variance in
comparison to the gradient computed over the entire training set. Still, one
advantage of small batch sizes is that they are more likely to make parameters
converge towards flatter minima of the loss (Goodfellow et al., 2016, Chapter
8.1.3), as explained below. Such flatter minima have better generalization
capacities, i.e. they maintain performance when presented with a new test set.

Keskar et al. (2016) define a _flat minimizer_ – as opposed to a _sharp_ one –
as a point in the parameter space that is a local minimum of the loss
function, and where this function varies slowly in a relatively large
neighborhood. Keskar et al. (2016) point to the following _generalization gap_
: training with large batch sizes tends to converge towards sharp minimizers,
which offer poorer generalization capacities. Conversely, _small batch sizes
allow convergence towards flat minimizers_ , which are likely to generalize
better. Thus, smaller batch sizes have _exploration abilities_ : the search is
more likely to exit the basins of sharp minimizers, and to tend towards flat
minimizers, from where noise will not cause them to exit.

Since a sharp minimizer requires high precision to be described, unlike a flat
one, the more noise there is in the gradient, the more unlikely it is that the
parameters will converge towards a sharp minimizer. This is precisely the
contribution of a smaller batch size: introduce noise in the gradient
estimation. According to this theoretical view, above a certain threshold of
the batch size, the generalization capacities of a model deteriorate. The
threshold depends on several hyper-parameters, including the batch size. Its
role has not been fully settled yet, with observations and conclusions varying
widely across studies (Dinh et al., 2017; Hoffer et al., 2017; Goyal et al.,
2017; Li et al., 2017; Kawaguchi et al., 2017). Moreover, these studies are on
image data sets, with fully connected or with convolutional NNs, which differ
substantially from NMT settings.

##  3 The Role of Batch Size in Neural MT

Several recent studies in NMT have considered batch size among other hyper-
parameters, but they have either been in high-resource settings (Popel and
Bojar, 2018; Xu et al., 2020) or have given only marginal attention to batch
size (Sennrich and Zhang, 2019; Araabi and Monz, 2020).

Popel and Bojar (2018) reported that BLEU scores increased with batch size
(including when using more GPUs) in a Transformer-based NMT system, although
with diminishing returns, recommending in particular that “batch size should
be set as high as possible”. Their experiments were performed using mainly two
datasets, with respectively 58M and 15M sentence pairs. It thus remains an
open question whether their findings regarding batch size also apply when much
less training data is available.

Sennrich and Zhang (2019) experimented with a recurrent network in a low-
resource setting and found that smaller batch sizes were beneficial, along
with other forms of regularization. They experimented with two batch sizes of
4,000 and 1,000 tokens, and observed improvements with the latter of 0.30 and
0.04 BLEU points on data sets with 5k and 160k sentence pairs, respectively.
It is difficult to predict from these results what the optimal batch size is
for Transformer-based NMT.

Araabi and Monz (2020) studied the role of 15 hyper-parameters of the
Transformer, with several sizes of low-resource datasets. For the largest
training sizes tested (80k and 165k sentence pairs), larger batch sizes
improved performance, with respectively 8,192 and 12,288 versus 4,096 for the
other sizes. For lower training sizes, smaller batch sizes did not improve
performance, which the authors explain by Transformer’s need for larger
batches. In our view, an alternative explanation is the order of optimization
of the hyper-parameters (a grid search in which they optimize one hyper-
parameter at a time): batch size is #12 out of 15, so by the time several
sizes are compared, regularization has already been introduced in the model by
dropouts on words, activation, and layers. Late optimization of batch size, of
warmup steps (#14) or of learning rate (#15) cannot properly determine their
regularizing effects.

Xu et al. (2020) proposed to compute gradients while accumulating minibatches,
and observed that increasing batch size stabilizes gradient direction up to a
certain point, after which it starts to fluctuate. They used this criterion to
dynamically adjust batch sizes while training. In their experiments with large
training sets (4.5M and 36M sentence pairs), their average batch size was
around 26k on two GPUs, and never lower than 7k. Their observations on the
gradient direction as more minibatches are accumulated are consistent with the
findings of Popel and Bojar (2018) who see diminishing returns when increasing
batch size.

##  4 Datasets and Systems

We train NMT systems with two low-resource parallel corpora, listed in the
first two lines of Table 1: the Upper Sorbian (HSB) to German (DE) training
data of the WMT 2020 Low-Resource Translation Task (Fraser, 2020) and a low-
size excerpt of the German to English News Commentary v13 (Bojar et al.,
2018), from which we randomly sampled 60k parallel lines. For the HSB-DE
models, we also use the development and test sets provided by the WMT 2020 and
2021 Low-Resource Translation Task (Libovický and Fraser, 2021), each
consisting of 2k sentences, and for DE-EN we sample a development set and a
test set from the original corpus, with 2k sentences each as well. We apply a
common filtering process for all data used: we delete from all our data the
sentences that are not between 2 and 300 words long, with resulting numbers of
lines shown in Table 1.

Dataset | Lang. | Orig. | Filt. |  $\Delta$%
---|---|---|---|---
WMT20 Low-res. | HSB-DE | 60k | 59.8k | 0.29
News Comm. v13 | DE-EN | 60k | 59.9k | 0.20
Sorbian Institute | HSB | 339k | 339k | 0.00
Witaj | HSB | 222k | 220k | 0.84
Web | HSB | 134k | 121k | 9.98
Europarl v8 | DE | 2.2M | 2.2M | 0.79
News Comm. v15 | DE | 422k | 411k | 2.58
JW300 | DE | 2.3M | 2.2M | 4.44
Europarl v3 | DE | 790k | 785k | 0.69
Europarl v3 | EN | 790k | 782k | 1.07
Table 1: Numbers of lines in the original and filtered corpora used in our
experiments. HSB stands for Upper Sorbian and $\Delta$% for the proportion of
lines filtered out. The only parallel corpora used for training NMT are the
first two ones; the other corpora are only used to train the SentencePiece
model.

We build subword vocabularies using the Unigram LM model (Sennrich et al.,
2016; Kudo, 2018) as implemented in
SentencePiece111https://github.com/google/sentencepiece, with the monolingual
corpora from Table 1. We train a shared model for HSB-DE with a vocabulary of
32k pieces, character coverage of 0.98, nbest = 1 and alpha = 0. The HSB data
adds up to 740k sentences, and we sample 680k sentences from three DE corpora,
and add them to the 60k sentences from the DE side of the parallel HSB-DE
corpus. To train the SentencePiece model for the DE-EN, for comparison
purposes, we treat German as a low-resource language, and sample 680k lines of
English and German from Europarl v3 (Tiedemann, 2012), which we combine
respectively with the 60k lines extracted from the DE-EN parallel corpus.

We use the Transformer-Base (Vaswani et al., 2017) in the implementation
provided by OpenNMT (Klein et al., 2017, 2020), with the parameters given in
Appendix A. Unless otherwise specified, we follow OpenNMT-py’s recommended
values for the hyper-parameters.222https://opennmt.net/OpenNMT-
py/examples/Translation.html

When using several GPUs with gradient accumulation, each GPU processes several
batches, which are then accumulated across all GPUs and used to update the
model at each step. Therefore, the effective batch size is $B\times G\times
A$, where $B$ is the individual batch size, $G$ is the number of GPUs and $A$
the number of accumulated batches, and differs from the batch_size hyper-
parameter $B$. We train all models on two GeForce RTX 1080Ti GPUs with 11 GB
of memory each and accumulate gradients over two minibatches ($A=2$),
following OpenNMT-py’s recommendation. Therefore, the batch_size parameter is
not our effective batch size, which is four times larger. Throughout this
work, we will refer to batch size $B$ as the batch_size parameter, and report
true _epochs_ , which we define as computed with the effective batch size as
$S\times B_{\mathit{eff}}/N$, for $S$ training steps, $B_{\mathit{eff}}$
effective batch size, and $N$ number of source tokens in the training set.

Following OpenNMT-py’s recommendations, we set the Adam hyper-parameters at
$\beta_{1}=0.9,\ \beta_{2}=0.998,\ \epsilon=10^{-8}$ and apply at each step a
scaling factor of two to Noam’s learning rate schedule, setting warmup steps
to 8k. Translations are generated with a beam width of seven, with an ensemble
of the last four saved checkpoints. We report BLEU scores (Papineni et al.,
2002) obtained with SacreBLEU (Post, 2018) on detokenized text.

##  5 Experimental Results

Figure 1: BLEU scores on the test set for HSB-DE models trained with different
batch sizes.

To study the impact of batch sizes in a low-resource setting, we train various
HSB-DE and DE-EN models for 700 epochs with the following batch sizes: 100,
250, 500, 1,000, 2,500, 5,000, 7,500, 10,000, and 10,240 (this is the largest
one that fits in our GPU memory).

###  5.1 NMT Performance

Batch | HSB-DE | DE-EN
---|---|---
Size | dev | test | dev | test
|  Xent |  BLEU |  chrF |  TER |  Xent |  BLEU |  chrF |  TER |  Xent |  BLEU |  chrF |  TER |  Xent |  BLEU |  chrF |  TER
500 |  0.03 |  48.12 |  71.13 |  37.35 |  0.03 |  41.53 |  67.34 |  43.84 |  0.11 |  37.35 |  58.04 |  54.54 |  0.11 |  37.72 |  58.35 |  54.60
1,000 |  0.02 |  49.23 |  72.07 |  36.35 |  0.02 |  42.26 |  67.93 |  43.16 |  0.05 |  38.03 |  59.39 |  52.91 |  0.05 |  38.67 |  59.68 |  52.71
2,500 |  0.03 |  48.28 |  71.63 |  37.02 |  0.03 |  41.18 |  67.36 |  44.02 |  0.04 |  33.83 |  56.70 |  56.27 |  0.04 |  35.51 |  57.76 |  55.47
5,000 |  0.03 |  46.99 |  70.74 |  38.05 |  0.03 |  40.28 |  66.62 |  45.24 |  0.05 |  32.47 |  55.20 |  57.88 |  0.05 |  33.97 |  56.16 |  57.08
7,500 |  0.03 |  46.05 |  70.29 |  38.87 |  0.03 |  39.10 |  65.94 |  46.18 |  0.05 |  32.67 |  55.99 |  57.67 |  0.05 |  33.80 |  56.72 |  57.21
10,000 |  0.04 |  44.61 |  69.19 |  40.00 |  0.04 |  38.41 |  65.67 |  46.45 |  0.05 |  31.84 |  55.20 |  58.35 |  0.05 |  33.50 |  56.14 |  57.63
10,240 |  0.04 |  45.59 |  70.12 |  39.26 |  0.04 |  38.19 |  65.39 |  46.79 |  0.06 |  31.49 |  55.00 |  58.65 |  0.06 |  33.03 |  55.78 |  58.07
Table 2: Loss and scores for models trained for 700 epochs with various batch
sizes for HSB-DE and DE-EN directions. All differences in BLEU on the dev and
test sets are statistically significant at the 95% level, except for the pairs
in similar colors.

NMT performance on the HSB-DE test set throughout the training is shown in
Figure 1, with BLEU scores depending on the number of epochs. The evolution
depending on real training time (wall time) is similar in terms of rankings.
Thus, the following analysis holds whether we train the models for the same
amount of epochs or of hours.

The final scores on the development and test sets are given in Table 2, sorted
by batch sizes. We provide first the actual loss of the model (‘Xent’ for
cross-entropy), and then three typical NMT scores: BLEU (Papineni et al.,
2002), chrF (Popović, 2015) and Translation Error Rate (Snover et al., 2006).
The 100 and 250 batch size models did not reach BLEU scores significantly
above zero, and are not included among the results in the table.

We test the statistical significance of the differences between each score and
the others, with 95% confidence, using the paired bootstrap resampling tool
from SacreBLEU (Post, 2018).333github.com/mjpost/sacrebleu with the signature
nrefs:1$|$bs:1000$|$seed:12345$|$case:mixed$|$eff:no
$|$tok:13a$|$smooth:exp$|$version:2.0.0. All differences between higher and
lower BLEU scores are statistically significant, except the pairs highlighted
in similar colors in Table 2.444The difference in BLEU between the following
pairs is not significant. For HSB-DE, 2,500 vs. 500, and 10,240 vs. 10,000, on
the test set; and 2,500 vs. 500, and 7,500 vs. 10,240 on the dev set. For DE-
EN these are 7,500 vs. 5,000, and 10,000 vs. 7,500, on the test set; and 1,000
vs. 500, 5,000 vs. 7,500, and 10,000 vs 10,240, on the dev set. The best NMT
scores, which are always obtained with a batch size of 1,000, are
significantly higher than all the other ones, including those obtained with
the largest possible batch sizes for our GPU (10,000 or 10,240). We thus
select two values for further experiments: a batch size of 1,000 as our
highest-scoring model, and one of 10,000 as the maximum allowed by our GPU
memory. A simple ratio of 10 holds between the two values.

These empirical results are contrary to those from Popel and Bojar (2018), who
observe that increasing the batch size for Transformer-Base produces higher
scores, although with diminishing returns after a certain threshold. We
hypothesize that the main explanation is the difference between the amounts of
training data: in our low-resource setting, we use 60k sentences, while Popel
and Bojar (2018) use 57M sentences. Our findings are consistent with those of
Keskar et al. (2016), who also observe that the optimal batch size is at the
lower end of the range, on a computer vision task with convolutional and
fully-connected NNs.

###  5.2 Asymptotic Performance

An alternative explanation for the previous results is that the learning rate
is too small for the larger batch sizes, which require more time to converge.
To test whether the differences observed above between small and large batch
sizes depend on the actual training time, we continue training the 1,000 and
10,000 batch size models for HSB-DE and DE-EN for twice as many epochs as
above (1400). The BLEU scores and their increases with respect to training for
700 epochs are given in Table 3. The performance gap (from +3.85 to +3.25
BLEU) between small and large batch sizes is not overturned by training the
models for much longer.

The scores from our best system (1,000 batch size, 42.81 BLEU on the test set)
are similar to scores obtained by _baselines_ of the five highest-scoring
teams at the WMT20 Low-resource shared task on HSB-DE (Fraser, 2020). While
the scores of Scherrer et al. (2020) and Li et al. (2020) are not comparable
due to a different architecture or the use of unsupervised pre-training, the
baseline scores of Knowles et al. (2020), Libovický et al. (2020) and
Kvapilíková et al. (2020) are respectively 44.1, 43.4, and 38.7. The first one
is higher than our best BLEU by 1.29, likely due to the use of 43M lines of CS
and DE data for the subword vocabulary, vs. 700k in our case.

Batch | HSB-DE |  | DE-EN
---|---|---|---
size |  dev |  test |  |  dev |  test
1,000 |  49.52 (+0.29) |  42.81 (+0.55) |  |  38.67 (+0.64) |  39.24 (+0.57)
10,000 |  46.44 (+1.83) |  39.56 (+1.15) |  |  33.19 (+1.35) |  34.42 (+0.92)
Table 3: BLEU scores for models trained for 1,400 epochs. The scores for 1,000
are significantly higher (at 95%) than those for 10,000. In parenthesis, the
absolute difference with BLEU scores after 700 epochs.

###  5.3 Invariance with respect to Vocabulary

We additionally perform two comparisons that show that the above results hold
regardless of the tokenizer and the vocabulary size. First, we test whether
the score difference is preserved with an unshared SentencePiece vocabulary,
i.e. when not sharing the source (HSB and DE) and the target (DE and EN)
vocabularies.

Second, we train two NMT models for HSB-DE using a Byte Pair Encoding (BPE)
vocabulary (Sennrich et al., 2016), which we generate using the learn_bpe.py
tool from OpenNMT-py, with 32k merge operations and the remaining parameters
at default values. Table 4 shows BLEU scores on the development sets for batch
sizes of 1,000 and 10,000. The previously observed differences in score
between the batch sizes still hold, and we see that a shared SentencePiece
vocabulary leads to a better NMT system than an unshared or a BPE one.

Batch size | SP unshared |  BPE
---|---|---
|  HSB-DE |  DE-EN |  HSB-DE
1,000 |  46.80 (-2.43) |  35.90 (-2.13) |  46.21 (-3.02)
10,000 |  41.99 (-2.62) |  30.09 (-1.75) |  43.35 (-1.26)
Table 4: BLEU scores on the dev set for HSB-DE and DE-EN models trained with
SentencePiece (SP) vocabularies not shared between source and target (left)
and BPE subwords (right). The scores for 1,000 are significantly higher (at
95%) than those for 10,000. In parenthesis, the difference with BLEU scores
obtained with the SP shared vocabulary.

##  6 Conclusion and Future Work

In this work, we have shown that insights from computer vision on the
regularizing effect of small batch sizes are also applicable to NMT. Our
results, focused on a low-resource setting, challenge those of previous NMT
studies with large amounts of training data, and the general belief that batch
sizes should be as large as they fit in the GPU memory. We have shown that
training with small batch sizes leads to models that generalize better, and
found the optimal batch size below which performance degrades.

Future work should explore how the learning rate must be adjusted depending on
the batch size, and whether a dynamically scheduled combination of batch size
and learning rate can provide an even better regularizer. For instance, it
should be tested if dynamic batch sizes as proposed by Xu et al. (2020) can
also improve performance in a low-resource setting, with batch size thresholds
changed to measure an optimal level of noise.

## Acknowledgments

We thank for their financial support the Swiss National Science Foundation
(grant n. 175693 for the DOMAT project: On-demand Knowledge for Document-level
Machine Translation) and Armasuisse (FamilyMT project). We especially thank
Dr. Ljiljana Dolamic (Armasuisse) for her support in the FamilyMT project. We
are also grateful to the anonymous reviewers and to Giorgos Vernikos for their
helpful suggestions.

## References

  * Araabi and Monz (2020) Ali Araabi and Christof Monz. 2020.  Optimizing Transformer for low-resource neural machine translation.  In _Proceedings of the 28th International Conference on Computational Linguistics_ , pages 3429–3435, Barcelona, Spain (Online). International Committee on Computational Linguistics.
  * Bojar et al. (2018) Ondřej Bojar, Christian Federmann, Mark Fishel, Yvette Graham, Barry Haddow, Philipp Koehn, and Christof Monz. 2018.  Findings of the 2018 conference on machine translation (WMT18).  In _Proceedings of the Third Conference on Machine Translation: Shared Task Papers_ , pages 272–303, Belgium, Brussels. Association for Computational Linguistics.
  * Dinh et al. (2017) Laurent Dinh, Razvan Pascanu, Samy Bengio, and Yoshua Bengio. 2017.  Sharp minima can generalize for deep nets.  In _International Conference on Machine Learning_ , pages 1019–1028. PMLR.
  * Fraser (2020) Alexander Fraser. 2020.  Findings of the WMT 2020 shared tasks in unsupervised MT and very low resource supervised MT.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 765–771, Online. Association for Computational Linguistics.
  * Goodfellow et al. (2016) Ian Goodfellow, Yoshua Bengio, and Aaron Courville. 2016.  _Deep Learning_.  MIT Press.
  * Goyal et al. (2017) Priya Goyal, Piotr Dollár, Ross Girshick, Pieter Noordhuis, Lukasz Wesolowski, Aapo Kyrola, Andrew Tulloch, Yangqing Jia, and Kaiming He. 2017.  Accurate, large minibatch SGD: Training Imagenet in 1 hour.  _arXiv preprint arXiv:1706.02677_.
  * Hoffer et al. (2017) Elad Hoffer, Itay Hubara, and Daniel Soudry. 2017.  Train longer, generalize better: Closing the generalization gap in large batch training of neural networks.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17, page 1729–1739, Red Hook, NY, USA. Curran Associates Inc.
  * Kawaguchi et al. (2017) Kenji Kawaguchi, Leslie Pack Kaelbling, and Yoshua Bengio. 2017.  Generalization in deep learning.  _arXiv preprint arXiv:1710.05468_.
  * Keskar et al. (2016) Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy, and Ping Tak Peter Tang. 2016.  On large-batch training for deep learning: Generalization gap and sharp minima.  _arXiv preprint arXiv:1609.04836_.
  * Klein et al. (2020) Guillaume Klein, François Hernandez, Vincent Nguyen, and Jean Senellart. 2020\.  The OpenNMT neural machine translation toolkit: 2020 edition.  In _Proceedings of the 14th Conference of the Association for Machine Translation in the Americas (Volume 1: Research Track)_ , pages 102–109, Virtual. Association for Machine Translation in the Americas.
  * Klein et al. (2017) Guillaume Klein, Yoon Kim, Yuntian Deng, Jean Senellart, and Alexander Rush. 2017\.  OpenNMT: Open-source toolkit for NMT.  In _Proceedings of ACL 2017, System Demonstrations_ , pages 67–72. Association for Computational Linguistics.
  * Knowles et al. (2020) Rebecca Knowles, Samuel Larkin, Darlene Stewart, and Patrick Littell. 2020.  NRC systems for low resource German-Upper Sorbian machine translation 2020: Transfer learning with lexical modifications.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 1112–1122, Online. Association for Computational Linguistics.
  * Kudo (2018) Taku Kudo. 2018.  Subword regularization: Improving neural network translation models with multiple subword candidates.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 66–75.
  * Kvapilíková et al. (2020) Ivana Kvapilíková, Tom Kocmi, and Ondřej Bojar. 2020.  CUNI systems for the unsupervised and very low resource translation task in WMT20.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 1123–1128, Online. Association for Computational Linguistics.
  * Li et al. (2017) Hao Li, Zheng Xu, Gavin Taylor, Christoph Studer, and Tom Goldstein. 2017.  Visualizing the loss landscape of neural nets.  _arXiv preprint arXiv:1712.09913_.
  * Li et al. (2020) Zuchao Li, Hai Zhao, Rui Wang, Kehai Chen, Masao Utiyama, and Eiichiro Sumita. 2020\.  SJTU-NICT’s supervised and unsupervised neural machine translation systems for the WMT20 news translation task.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 218–229, Online. Association for Computational Linguistics.
  * Libovický et al. (2020) Jindřich Libovický, Viktor Hangya, Helmut Schmid, and Alexander Fraser. 2020.  The LMU Munich system for the WMT20 very low resource supervised MT task.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 1104–1111, Online. Association for Computational Linguistics.
  * Libovický and Fraser (2021) Jindřich Libovický and Alexander Fraser. 2021.  Findings of the WMT 2021 shared tasks in unsupervised MT and very low resource supervised MT.  In _Proceedings of the Sixth Conference on Machine Translation_ , Online. Association for Computational Linguistics.
  * Papineni et al. (2002) Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu. 2002.  Bleu: a method for automatic evaluation of machine translation.  In _Proceedings of the 40th Annual Meeting of the Association for Computational Linguistics_ , pages 311–318, Philadelphia, Pennsylvania, USA. Association for Computational Linguistics.
  * Popel and Bojar (2018) Martin Popel and Ondřej Bojar. 2018.  Training tips for the Transformer model.  _The Prague Bulletin of Mathematical Linguistics_ , 110(1):43–70.
  * Popović (2015) Maja Popović. 2015.  chrF: character n-gram f-score for automatic MT evaluation.  In _Proceedings of the Tenth Workshop on Statistical Machine Translation_ , pages 392–395, Lisbon, Portugal. Association for Computational Linguistics.
  * Post (2018) Matt Post. 2018.  A call for clarity in reporting BLEU scores.  In _Proceedings of the Third Conference on Machine Translation: Research Papers_ , pages 186–191, Belgium, Brussels. Association for Computational Linguistics.
  * Scherrer et al. (2020) Yves Scherrer, Stig-Arne Grönroos, and Sami Virpioja. 2020.  The University of Helsinki and Aalto university submissions to the WMT 2020 news and low-resource translation tasks.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 1129–1138, Online. Association for Computational Linguistics.
  * Sennrich et al. (2016) Rico Sennrich, Barry Haddow, and Alexandra Birch. 2016.  NMT of rare words with subword units.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1715–1725.
  * Sennrich and Zhang (2019) Rico Sennrich and Biao Zhang. 2019.  Revisiting low-resource neural machine translation: A case study.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 211–221, Florence, Italy. Association for Computational Linguistics.
  * Snover et al. (2006) Matthew Snover, Bonnie Dorr, Richard Schwartz, Linnea Micciulla, and John Makhoul. 2006.  A study of translation edit rate with targeted human annotations.  In _Proceedings of the Conference of the Association for Machine Translation in the Americas (AMTA 2006)_.
  * Tiedemann (2012) Jörg Tiedemann. 2012.  Parallel data, tools and interfaces in OPUS.  In _Proceedings of the Eight International Conference on Language Resources and Evaluation (LREC’12)_ , Istanbul, Turkey. European Language Resources Association (ELRA).
  * Van Biljon et al. (2020) Elan Van Biljon, Arnu Pretorius, and Julia Kreutzer. 2020.  On optimal Transformer depth for low-resource language translation.  _arXiv preprint arXiv:2004.04418_.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in Neural Information Processing Systems 30_ , pages 5998–6008.
  * Xu et al. (2020) Hongfei Xu, Josef van Genabith, Deyi Xiong, and Qiuhui Liu. 2020.  Dynamically adjusting Transformer batch size by monitoring gradient direction change.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 3519–3524, Online. Association for Computational Linguistics.

##  Appendix A Appendix

The hyper-parameters used to train our models are the following ones:

src_words_min_frequency: 2
tgt_words_min_frequency: 2
valid_batch_size: 200
max_generator_batches: 2
optim: adam
learning_rate: 2.0
adam_beta2: 0.998
decay_method: noam
accum_count: 2
warmup_steps: 8000
label_smoothing: 0.1
max_grad_norm: 0
param_init: 0
param_init_glorot: true
normalization: tokens
encoder_type: transformer
decoder_type: transformer
position_encoding: true
layers: 6
heads: 8
rnn_size: 512
word_vec_size: 512
transformer_ff: 2048
dropout: 0.1
batch_size: 1000
batch_type: tokens


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

Guang-Yu Wanga<EMAIL_ADDRESS>Xing-Gang Wua<EMAIL_ADDRESS>Xu-Chang Zhenga
<EMAIL_ADDRESS>Jiang Yana<EMAIL_ADDRESS>Jia-Wei Zhangb
<EMAIL_ADDRESS>a Department of Physics, Chongqing Key Laboratory for
Strongly Coupled Physics, Chongqing University, Chongqing 401331, People’s
Republic of China
b Department of Physics, Chongqing University of Science and Technology,
Chongqing, 401331, People’s Republic of China

###### Abstract

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/\psi$ 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\to J/\psi+J/\psi$ based on the data from $pp$
collisions with an integrated luminosity of $37.5\,{\rm fb}^{-1}$ at
$\sqrt{s}=13\,{\rm TeV}$ at the large hadron collider (LHC). The observed
upper limit on the branching fraction for the decay is ${\rm Br}(Z\to
J/\psi+J/\psi)<2.2\times 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\,{\rm fb}^{-1}$.
Their updated upper limit is ${\rm Br}(Z\to J/\psi+J/\psi)<1.4\times 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\times
10^{12}$ $Z$ bosons will be produced [6]. Therefore, these lepton colliders
provide new opportunities for studying the rare decay $Z\to J/\psi+J/\psi$.

On the theoretical side, the decay $Z\to J/\psi+J/\psi$ was first studied in
Ref.[8] in 1990. In Ref.[9], this decay was reanalysed at leading order (LO)
in $\alpha_{s}$ and $v$ ($v$ stands for the typical velocity of the $c$ quark
or the $\bar{c}$ quark in the rest frame of $J/\psi$) 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\to J/\psi+\gamma^{*}$ followed by the fragmentation $\gamma^{*}\to
J/\psi$ can also bring significant contribution. Their calculation shows that
the branching fraction of $Z\to J/\psi+J/\psi$ can be enhanced to $1.1\times
10^{-10}$ from $1.5\times 10^{-13}$ after including the single photon
fragmentation contribution. Recently, the next-to-leading order (NLO) QCD
correction to the decay rate of $Z\to J/\psi+J/\psi$ 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\to J/\psi+J/\psi$. 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\to J/\psi+\gamma^{*}$, and
the other for $\gamma^{*}\to J/\psi$. For the amplitude of $\gamma^{*}\to
J/\psi$, we directly extract it from the experimentally measured leptonic
decay width of the $J/\psi$, rather than calculating it order by order in
$\alpha_{s}$ and $v$. Using this approach, the QCD and relativistic
corrections to the amplitude of $\gamma^{*}\to J/\psi$ are resummed up to all
orders. For the amplitude of $Z\to J/\psi+\gamma^{*}$, 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/\psi$ light-cone
distribution amplitude (LCDA). The hard-scattering kernel can be calculated
with the perturbative QCD (pQCD), and the LCDA for $J/\psi$ can be calculated
through the NRQCD factorization [15]. Under the light-cone approach,
logarithms of $m_{{}_{Z}}^{2}/m_{c}^{2}$ 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\to J/\psi+J/\psi$, to calculate the non-
fragmentation 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\to J/\psi+J/\psi$. In
Sec.III, we present the numerical results and discussions. Section IV is
reserved as a summary.

##  II Calculation formalism

###  II.1 Fragmentation amplitude

Figure 1: Two typical LO Feynman diagrams for the photon fragmentation
mechanism of the decay $Z\to J/\psi+J/\psi$.

The dominant mechanism of the decay $Z\to J/\psi+J/\psi$ is via the photon
fragmentation, i.e., the transition $Z\to J/\psi+\gamma^{*}$ followed by
$\gamma^{*}\to J/\psi$. To understand this mechanism, we present two typical
LO Feynman diagrams for the photon fragmentation mechanism in Fig.1 111It 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.. The amplitude for the photon
fragmentation mechanism can be written as

$\displaystyle i\mathcal{M}^{\rm fr}_{Z\to J/\psi+J/\psi}$
$\displaystyle=i\mathcal{M}^{\mu}_{Z\to
J/\psi(P_{1})+\gamma^{*}(P_{2})}\frac{-i}{m_{J/\psi}^{2}}i{\cal
M}_{\gamma^{*}\to J/\psi,\mu}$ $\displaystyle~{}~{}~{}+(1\leftrightarrow 2),$
(1)

where the factor $-i/m_{J/\psi}^{2}$ arises from the photon propagator that
will fragment into a $J/\psi$. $i{\cal M}_{\gamma^{*}\to J/\psi,\mu}$ denotes
the amplitude for the transition of a virtual photon into the $J/\psi$, and it
can be expressed as [16, 17]

$\displaystyle i{\cal M}^{\mu}_{\gamma^{*}\to J/\psi}$
$\displaystyle=-ie\langle
J/\psi|J^{\mu}(x=0)|0\rangle=-ieg_{J/\psi\gamma}\epsilon_{J/\psi}^{*\mu},$ (2)

where $J^{\mu}(x)$ is the electromagnetic current,

$\displaystyle J^{\mu}(x)=\sum_{q}e_{q}\bar{q}(x)\gamma^{\mu}q(x).$ (3)

Here, the sum extends over all quark flavors. The $J/\psi$-photon effective
coupling constant $g_{J/\psi\gamma}$ can be determined using the decay width
of the $J/\psi$ to an electron-positron pair [16], i.e.,

$\displaystyle\Gamma(J/\psi\to
e^{+}e^{-})=\frac{4\pi\alpha^{2}(m_{J/\psi})g_{J/\psi\gamma}^{2}}{3m_{J/\psi}^{3}}.$
(4)

The method of calculating the amplitude $i{\cal M}_{\gamma^{*}\to J/\psi,\mu}$
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 $\alpha_{s}$
and $v$. The calculations of higher-order QCD and relativistic corrections are
very difficult. On the other hand, the calculation of the amplitude $i{\cal
M}_{\gamma^{*}\to J/\psi,\mu}$ 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
$\mathcal{M}^{\mu}_{Z\to J/\psi+\gamma^{*}}$. We shall adopt the approximation
$\mathcal{M}^{\mu}_{Z\to J/\psi+\gamma^{*}}\approx\mathcal{M}^{\mu}_{Z\to
J/\psi+\gamma}$, which will only lead to an error of ${\cal
O}(m_{J/\psi}^{2}/m_{{}_{Z}}^{2})$. The decay $Z\to J/\psi+\gamma$ 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 $\mathcal{M}^{\mu}_{Z\to J/\psi+\gamma}$.

Under the light-cone approach, the amplitude of $Z\to J/\psi+\gamma$ is
expressed by [21]

$\displaystyle i\mathcal{M}_{Z\to
J/\psi+\gamma,\mu}=i\mathcal{A}\,\epsilon_{\xi\mu\nu\rho}\epsilon_{Z}^{\xi}\epsilon_{J/\psi}^{*\nu}p_{\gamma}^{\rho},$
(5)

where

$\displaystyle
i\mathcal{A}=-\frac{ee_{c}g_{{}_{Z}}g_{A}^{c}m_{J/\psi}}{m_{{}_{Z}}^{2}}f_{J/\psi}^{\parallel}\int_{0}^{1}dx\,T_{H}(x,\mu)\phi_{J/\psi}^{\parallel}(x,\mu).$
(6)

Here, $e_{c}=+2/3$ is the electric charge of the charm quark,
$f_{J/\psi}^{\parallel}$ is the decay constant of the longitudinally polarized
$J/\psi$, $T_{H}(x,\mu)$ is the hard-scattering kernel,
$\phi_{J/\psi}^{\parallel}$ is the longitudinally polarized LCDA of the
$J/\psi$, $\epsilon_{Z}$ is the polarization vector of the $Z$ boson,
$\epsilon^{*}_{J/\psi}$ is the polarization vector of $J/\psi$, and
$p_{\gamma}$ is the momentum of the photon, $\mu$ is the factorization scale,
and $x$ is the momentum fraction of the $c$ quark in the $J/\psi$. The
couplings $g_{{}_{Z}}$ and $g_{A}^{c}$ are defined as
$g_{{}_{Z}}=2(\sqrt{2}G_{F})^{1/2}m_{{}_{Z}}$,
$g_{A}^{c}=\frac{1}{2}t_{3L}^{c}$, where $G_{F}$ is the Fermi constant and
$t_{3L}^{c}=1/2$ is the weak isospin for the left-handed $c$ quark.

The hard-scattering kernel $T_{H}(x,\mu)$ can be calculated through pQCD, and
it has been calculated up to order $\alpha_{s}$ in Ref.[20]. The expression of
the hard-scattering kernel, up to order $\alpha_{s}$, is

$T_{H}(x,\mu)=T_{H}^{(0)}(x,\mu)+\frac{\alpha_{s}(\mu)}{4\pi}T_{H}^{(1)}(x,\mu),$
(7)

where

$\displaystyle T_{H}^{(0)}(x,\mu)=$ $\displaystyle\,\frac{1}{x(1-x)},$
$\displaystyle T_{H}^{(1)}(x,\mu)=$
$\displaystyle\,C_{F}\frac{1}{x(1-x)}\bigg{\\{}\big{[}3+2x\,{\rm ln}(1-x)$
$\displaystyle\,+2(1-x)\,{\rm ln}x\big{]}\left({\rm
ln}\frac{m_{{}_{Z}}^{2}}{\mu^{2}}-i\pi\right)$ $\displaystyle\,+x\,{\rm
ln}^{2}(1-x)+(1-x)\,{\rm ln}^{2}x$ $\displaystyle\,-(1-x)\,{\rm
ln}(1-x)-x\,{\rm ln}x-9\bigg{\\}}.$ (8)

Here, $C_{F}=(N_{c}^{2}-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 $\mu=m_{{}_{Z}}$ in Eq.(6).

The longitudinally polarized LCDA $\phi_{J/\psi}^{\parallel}$ is defined as
[24, 25]:

$\displaystyle\langle J/\psi(p)|\bar{Q}(z)\gamma^{\mu}[z,0]Q(0)|0\rangle$
$\displaystyle=$ $\displaystyle-p^{\mu}\frac{\epsilon_{J/\psi}^{*}\cdot
z}{p\cdot z}f_{J/\psi}^{\parallel}m_{J/\psi}\int_{0}^{1}dx\,e^{ip\cdot
zx}\phi_{J/\psi}^{\parallel}(x,\mu),$ (9)

where $p$ is the quarkonium momentum, the LCDA is normalized as
$\int_{0}^{1}dx\,\phi_{J/\psi}^{\parallel}(x,\mu)=1$. 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

$\displaystyle[z,0]=P\,{\rm
exp}\left[ig_{s}\int_{0}^{z}dxA_{a}^{+}(x)T^{a}\right],$ (10)

where $P$ denotes the path ordering.

Under the NRQCD factorization, the LCDA $\phi_{J/\psi}^{\parallel}$ is
expressed as a double expansion in $\alpha_{s}$ and $v$ [20, 26], i.e.,

$\displaystyle\phi_{J/\psi}^{\parallel}(x,\mu_{0})$
$\displaystyle=\phi_{J/\psi}^{\parallel(0)}(x,\mu_{0})+\langle
v^{2}\rangle_{J/\psi}\phi_{J/\psi}^{\parallel(v^{2})}(x,\mu_{0})$
$\displaystyle+\frac{\alpha_{s}(\mu_{0})}{4\pi}\phi_{J/\psi}^{\parallel(1)}(x,\mu_{0})+{\cal
O}(\alpha_{s}^{2},\alpha_{s}v^{2},v^{4}).$

Here, in order to avoid large logarithms appearing in the LCDA, we set the
initial factorization scale of the LCDA as $\mu_{0}=m_{c}$. The quantity
$\langle v^{2}\rangle_{J/\psi}$ is a ratio of NRQCD LDMEs, whose definition is

$\langle
v^{2}\rangle_{J/\psi}=\frac{\langle\bm{q}^{2}\rangle_{J/\psi}}{m_{c}^{2}}=\frac{\langle
J/\psi(\lambda)|\psi^{\dagger}(-\frac{i}{2}\tensor{\bm{D}})^{2}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle}{m_{c}^{2}\,\langle
J/\psi(\lambda)|\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle},$
(12)

where the covariant derivative operator $\psi^{\dagger}\tensor{\bf
D}\chi\equiv\psi^{\dagger}{\bf D}\chi-({\bf D}\psi)^{\dagger}\chi$, $\psi$ and
$\chi$ are Pauli spinor fields that describe $c$ quark annihilation and
$\bar{c}$ quark creation, respectively.

The expression for the LO LCDA is [20]

$\displaystyle\phi_{J/\psi}^{\parallel(0)}(x,\mu_{0})=$
$\displaystyle\,\delta(x-\tfrac{1}{2}).$ (13)

The order-$v^{2}$ LCDA $\phi_{J/\psi}^{\parallel(v^{2})}(x,\mu_{0})$ was
calculated in Ref.[27], and the order-$\alpha_{s}$ LCDA
$\phi_{J/\psi}^{\parallel(1)}(x,\mu_{0})$ was calculated in Ref.[20]. Their
expressions are

$\displaystyle\phi_{J/\psi}^{\parallel(v^{2})}(x,\mu_{0})=$
$\displaystyle\,\frac{\delta^{(2)}(x-\frac{1}{2})}{24},$ (14)
$\displaystyle\phi_{J/\psi}^{\parallel(1)}(x,\mu_{0})=$
$\displaystyle\,C_{F}\theta(1-2x)\Bigg{\\{}\Bigg{[}\left(4x+\frac{8x}{1-2x}\right)$
$\displaystyle\,\cdot\left(\log\frac{\mu_{0}^{2}}{m_{c}^{2}(1-2x)^{2}}-1\right)\Bigg{]}_{+}-\left[8x\right]_{+}$
$\displaystyle\,+\left[\frac{16x(1-x)}{(1-2x)^{2}}\right]_{++}\Bigg{\\}}+(x\leftrightarrow
1-x),$ (15)

where $\delta^{(2)}$ stands for the second derivative of the delta function.
The $+$ and $++$ functions are defined as

$\displaystyle\int_{0}^{1}{\rm d}x\,[f(x)]_{+}g(x)=$
$\displaystyle\,\int_{0}^{1}{\rm d}x\,f(x)\left[g(x)-g(1/2)\right],$ (16)
$\displaystyle\int_{0}^{1}{\rm d}x\,[f(x)]_{++}g(x)=$
$\displaystyle\,\int_{0}^{1}{\rm d}x\,f(x)[g(x)-g(1/2)$
$\displaystyle\,-g^{\prime}(1/2)(x-1/2)].$ (17)

To calculate the amplitude based on Eq.(6), we need the LCDA at the scale
$\mu=m_{{}_{Z}}$. We calculate the LCDA
$\phi_{J/\psi}^{\parallel}(x,\mu=m_{{}_{Z}})$ through solving the evolution
equation of the LCDA, where the initial LCDA
$\phi_{J/\psi}^{\parallel}(x,\mu_{0}=m_{c})$ is used as the boundary
condition. The scale evolution of the LCDA $\phi_{J/\psi}^{\parallel}(x,\mu)$
obeys the Efremov-Radyushkin-Brodsky-Lepage (ERBL) equation [13, 28, 29],
i.e.,

$\displaystyle\frac{\partial}{\partial\,{\rm
ln}\,\mu^{2}}\phi_{J/\psi}^{\parallel}(x,\mu)=\int_{0}^{1}dy\,V_{\parallel}[x,y;\alpha_{s}(\mu)]\,\phi_{J/\psi}^{\parallel}(y,\mu).$
(18)

where the evolution kernel can be expanded in $\alpha_{s}$ as

$\displaystyle V_{\parallel}[x,y;\alpha_{s}(\mu)]$
$\displaystyle=\frac{\alpha_{s}(\mu)}{4\pi}V^{(0)}(x,y)+\left(\frac{\alpha_{s}(\mu)}{4\pi}\right)^{2}V^{(1)}(x,y)$
$\displaystyle~{}~{}~{}+{\cal O}(\alpha_{s}^{2}),$ (19)

At present, the order-$\alpha_{s}$ evolution kernel $V^{(0)}(x,y)$ is given in
Refs.[13, 28], and the order-$\alpha_{s}^{2}$ evolution kernel $V^{(1)}(x,y)$
is given in Refs.[30, 31, 32].

To solve the ERBL equation, i.e., Eq.(18), it is convenient to expand the LCDA
in terms of the Gegenbauer polynomials $C_{n}^{(3/2)}(2x-1)$, which are the
eigenfunctions of the LO evolution kernel. The Gegenbauer expansion of the
LCDA is

$\displaystyle\phi_{J/\psi}^{\parallel}(x,\mu)=\sum_{n=0}^{\infty}a_{n}^{\parallel}(\mu)\,x(1-x)\,C_{n}^{(3/2)}(2x-1),$
(20)

where $a_{n}^{\parallel}(\mu)$ is the $n_{\rm th}$ Gegenbauer moment of the
LCDA $\phi_{J/\psi}^{\parallel}(x,\mu)$, and

$\displaystyle
a_{n}^{\parallel}(\mu)=N_{n}\int_{0}^{1}dx\,C_{n}^{(3/2)}(2x-1)\phi_{J/\psi}^{\parallel}(x,\mu),$
(21)

where the normalization factor

$\displaystyle N_{n}=\frac{4(2n+3)}{(n+1)(n+2)}.$ (22)

Solving the ERBL equation in the Gegenbauer moment space, the obtained the
evolved Gegenbauer moments have the following form:

$a_{n}^{\parallel}(\mu)=\sum_{k=0}^{n}U_{nk}(\mu,\mu_{0})a_{k}^{\parallel}(\mu_{0}).$
(23)

The expressions for $U_{nk}(\mu,\mu_{0})$ up to next-leading-logarithmic (NLL)
accuracy can be found in Ref.[33].

Similarly, the hard-scattering kernel can also be expanded in terms of the
Gegenbauer polynomials as

$\displaystyle
T_{H}(x,\mu)=\sum_{n=0}^{\infty}N_{n}b_{n}(\mu)C_{n}^{(3/2)}(2x-1).$ (24)

The $n_{\rm th}$ Gegenbauer moment

$\displaystyle
b_{n}(\mu)=\int_{0}^{1}dx\,x(1-x)C_{n}^{(3/2)}(2x-1)T_{H}(x,\mu).$ (25)

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

$\displaystyle\int_{0}^{1}dx\,T_{H}(x,\mu)\phi_{J/\psi}^{\parallel}(x,\mu)$
$\displaystyle=\sum_{n=0}^{\infty}b_{n}(\mu)a_{n}(\mu)$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{k=0}^{n}b_{n}(\mu)U_{nk}(\mu,\mu_{0})a_{k}(\mu_{0}),$
(26)

As stated below Eq.(II.1), we have set $\mu=m_{{}_{Z}}$ so that the hard-
scattering kernel $T_{H}(x,\mu)$ does not contain large logarithms. The large
logarithms of $m_{{}_{Z}}^{2}/m_{c}^{2}$ will appear in the LCDA
$\phi_{J/\psi}^{\parallel}(x,\mu=m_{{}_{Z}})$. In our calculation, these large
logarithms are resummed up to NLL accuracy through the evolution of the LCDA
from the initial scale $\mu_{0}=m_{c}$ to the scale $\mu=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

$\displaystyle\sum_{n=0}^{\infty}\sum_{k=0}^{n}b_{n}(\mu)U_{nk}(\mu,\mu_{0})a_{k}(\mu_{0})$
$\displaystyle=\lim_{z\to
1}\sum_{n=0}^{\infty}\sum_{k=0}^{n}b_{n}(\mu)z^{n}U_{nk}(\mu,\mu_{0})a_{k}(\mu_{0}).$
(27)

Then, we keep 20 nonzero terms in the second line of Eq.(27), and generate a
$10\times 10$ Padé approximant. It is worth noting that the limit of $z\to 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/\psi}^{\parallel}$. Setting
$z$ to 0 in Eq.(9) and using the normalization
$\int_{0}^{1}dx\,\phi_{J/\psi}^{\parallel}(x,\mu)=1$, we obtain

$\langle
J/\psi|\bar{c}(0)\gamma^{\mu}c(0)|0\rangle=-f_{J/\psi}^{\parallel}m_{J/\psi}\epsilon_{J/\psi}^{*\mu}.$
(28)

With this relation, we can relate the decay constant $f_{J/\psi}^{\parallel}$
with the leptonic decay width of $J/\psi$ [21], i.e.,

$\Gamma(J/\psi\to
e^{+}e^{-})=\frac{4\pi}{3m_{J/\psi}}\alpha^{2}(m_{J/\psi})e_{c}^{2}f_{J/\psi}^{\parallel
2}.$ (29)

The decay constant $f_{J/\psi}^{\parallel}$ can be determined through this
relation straightforwardly.

###  II.2 Non-fragmentation amplitude

Figure 2: Six typical tree-level and one-loop non-fragmentation Feynman
diagrams that contribute to the decay $Z\to J/\psi+J/\psi$.

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 non-fragmentation amplitude up to
order-$v^{2}$ for the decay $Z\to J/\psi+J/\psi$ can be written as

$\displaystyle{\cal M}^{\rm nfr}_{Z\to J/\psi+J/\psi}$
$\displaystyle=2m_{J/\psi}\left(c_{0}+c_{2}\langle
v^{2}\rangle_{J/\psi}\right)\langle
J/\psi|\psi^{\dagger}{\bm{\sigma}}\cdot{\bm{\epsilon}}\chi|0\rangle^{2},$ (30)

where $c_{0}$ and $c_{2}$ are short-distance coefficients (SDCs) that can be
expanded in powers of $\alpha_{s}$. In this work, we calculate the SDC $c_{0}$
up to NLO in $\alpha_{s}$, and calculate the SDC $c_{2}$ at LO in
$\alpha_{s}$. $\langle
J/\psi|\psi^{\dagger}{\bm{\sigma}}\cdot{\bm{\epsilon}}\chi|0\rangle$ is the
LDME for $J/\psi$. $\langle v^{2}\rangle_{J/\psi}$ is a ratio of LDMEs, which
has been defined in Eq.(12). The factor $2m_{J/\psi}$ arises from that we use
the relativistic normalization for the $J/\psi$ in ${\cal M}^{\rm nfr}_{Z\to
J/\psi+J/\psi}$ 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/\psi$ state in Eq.(30) by an on-shell
$(c\bar{c})$ pair with a small relative velocity $v$, i.e.,

$\displaystyle{\cal M}^{\rm
nfr}_{Z\to(c\bar{c})[^{3}S_{1}^{[1]}]+(c\bar{c})[^{3}S_{1}^{[1]}]}$
$\displaystyle=\left(c_{0}+c_{2}\,v^{2}\right)\langle(c\bar{c})[^{3}S_{1}^{[1]}]|\psi^{\dagger}{\bm{\sigma}}\cdot{\bm{\epsilon}}\chi|0\rangle^{2},$
(31)

where the amplitude ${\cal M}^{\rm
nfr}_{Z\to(c\bar{c})[^{3}S_{1}^{[1]}]+(c\bar{c})[^{3}S_{1}^{[1]}]}$ can be
calculated using perturbative theory in QCD, and the LDME
$\langle(c\bar{c})[^{3}S_{1}^{[1]}]|\psi^{\dagger}{\bm{\sigma}}\cdot{\bm{\epsilon}}\chi|0\rangle$
for an on-shell $(c\bar{c})$ 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 $(c\bar{c})$ pair
production can be applied to the $J/\psi$ production.

To calculate the amplitude ${\cal M}^{\rm
nfr}_{Z\to(c\bar{c})[^{3}S_{1}^{[1]}]+(c\bar{c})[^{3}S_{1}^{[1]}]}$ up to
order $\alpha_{s}$ and order $v^{2}$, we assign the momenta of the $c$ and
$\bar{c}$ quarks as $p_{i}=\frac{1}{2}P_{i}+q_{i}$ and
$\bar{p}_{i}=\frac{1}{2}P_{i}-q_{i}$, where $i=1,2$. Then we have the
relations

$\displaystyle
p_{i}^{2}=\bar{p}_{i}^{2}=m_{c}^{2},\,P_{i}^{2}=4E_{q}^{2},\,P_{i}\cdot
q_{i}=0,$ $\displaystyle q_{i}^{2}=m_{c}^{2}-E_{q}^{2}=-m_{c}^{2}v^{2}.$ (32)

In the rest frame of a $(c\bar{c})$ pair, the total and relative momenta of
the $(c\bar{c})$ pair are $P_{i}=(2E_{q},\bm{0})$ and $q_{i}=(0,\bm{q}_{i})$
where $E_{q}=\sqrt{\bm{q}_{i}^{2}+m_{c}^{2}}$. In deriving these relations, we
have implicitly adopted the relation
$\bm{q}_{1}^{2}=\bm{q}_{2}^{2}=m_{c}^{2}v^{2}$.

It is convenient to adopt the covariant projection method to enforce the
produced $(c\bar{c})$ pairs to be the corresponding spin and color state. The
spin projector for the spin-triplet state is [35]

$\displaystyle\Pi=$
$\displaystyle\frac{1}{4\sqrt{2}E_{q}(E_{q}+m_{c})}(\not{\bar{p}}_{i}-m_{c})\not{\epsilon}^{*}(\not{P}_{i}+2E_{q})$
(33) $\displaystyle\cdot(\not{p}_{i}+m_{c}),$

where $\epsilon^{*}$ is the polarization vector of the spin-triplet state. The
color projector for color-singlet state is

$\displaystyle\Lambda_{1}=\frac{\bf 1}{\sqrt{3}},$ (34)

where ${\bf 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,

$\displaystyle{\cal M}^{\rm
nfr}_{Z\to(c\bar{c})[^{3}S_{1}^{[1]}]+(c\bar{c})[^{3}S_{1}^{[1]}]}$
$\displaystyle={\cal
M}\Big{|}_{q_{1}=q_{2}=0}+\sum_{i=1}^{2}q_{i}^{\alpha}\frac{\partial{\cal
M}}{\partial q_{i}^{\alpha}}\Big{|}_{q_{1}=q_{2}=0}$
$\displaystyle~{}~{}~{}+\frac{1}{2!}\sum_{i,j=1}^{2}q_{i}^{\alpha}q_{j}^{\beta}\frac{\partial^{2}{\cal
M}}{\partial q_{i}^{\alpha}\partial
q_{j}^{\beta}}\Big{|}_{q_{1}=q_{2}=0}+\cdots,$ (35)

where ${\cal M}^{\rm
nfr}_{Z\to(c\bar{c})[^{3}S_{1}^{[1]}]+(c\bar{c})[^{3}S_{1}^{[1]}]}$ is the
non-fragmentation amplitude, which has been abbreviated as ${\cal M}$ on right
hand side of Eq.(35). To project out the $S$-wave contribution, we make the
replacement

$\displaystyle
q_{i}^{\mu}q_{i}^{\nu}\to\frac{{\bm{q}}_{i}^{2}}{(D-1)}\Pi^{\mu\nu}.$ (36)

where

$\displaystyle\Pi^{\mu\nu}=-g^{\mu\nu}+\frac{P_{i}^{\mu}P_{i}^{\nu}}{P_{i}^{2}},$
(37)

The odd-power terms of $q_{i}$ vanish for the S-wave contribution. Then we
obtain the amplitude up to order $v^{2}$,

$\displaystyle{\cal M}^{\rm
nfr}_{Z\to(c\bar{c})[^{3}S_{1}^{[1]}]+(c\bar{c})[^{3}S_{1}^{[1]}]}$
$\displaystyle={\cal
M}\Big{|}_{q_{1}=q_{2}=0}+\frac{1}{2!}\sum_{i=1}^{2}\frac{{\bm{q}}_{i}^{2}\Pi_{\alpha\beta}}{(D-1)}\frac{\partial^{2}{\cal
M}}{\partial q_{i}^{\alpha}\partial q_{i}^{\beta}}\Big{|}_{q_{1}=q_{2}=0}.$
(38)

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\epsilon$. 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 ($\overline{\rm MS}$) scheme. The corresponding
renormalization constants are $Z_{i}=1+\delta Z_{i}$, and

$\displaystyle\delta Z_{2}^{\rm OS}=$ $\displaystyle-
C_{F}\dfrac{\alpha_{s}}{4\pi}\left[\dfrac{1}{\epsilon_{UV}}+\frac{2}{\epsilon_{IR}}-3\gamma_{E}+3\ln\frac{4\pi\mu_{r}^{2}}{m_{c}^{2}}+4\right],$
$\displaystyle\delta Z_{m}^{\rm OS}=$
$\displaystyle-3C_{F}\dfrac{\alpha_{s}}{4\pi}\left[\dfrac{1}{\epsilon_{UV}}-\gamma_{E}+\ln\dfrac{4\pi\mu_{r}^{2}}{m_{c}^{2}}+\frac{4}{3}\right],$
$\displaystyle\delta Z_{3}^{\rm OS}=$
$\displaystyle\frac{\alpha_{s}}{4\pi}\bigg{[}(\beta^{\prime}_{0}-2C_{A})\left(\dfrac{1}{\epsilon_{UV}}-\frac{1}{\epsilon_{IR}}\right)$
$\displaystyle-\frac{4}{3}T_{F}\left(\dfrac{1}{\epsilon_{UV}}-\gamma_{E}+\ln\dfrac{4\pi\mu_{r}^{2}}{m_{c}^{2}}\right)\bigg{]},$
$\displaystyle\delta Z_{g}^{\overline{\mathrm{MS}}}=$
$\displaystyle-\dfrac{\beta_{0}}{2}\dfrac{\alpha_{s}}{4\pi}\left[\dfrac{1}{\epsilon_{UV}}-\gamma_{E}+\ln(4\pi)\right],$
(39)

where $\gamma_{E}$ is the Euler constant, $\mu_{r}$ is the renormalization
scale, $\beta_{0}=11/3-2n_{f}/3$, $\beta^{\prime}_{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 $\gamma_{5}$ matrix should be noted. In
this work, we employ the Larin scheme [37], where the product of a
$\gamma_{\mu}$ matrix and a $\gamma_{5}$ matrix is expreesed as:

$\displaystyle\gamma_{\mu}\gamma_{5}=\frac{i}{6}\epsilon_{\mu\rho\sigma\tau}\gamma^{\rho}\gamma^{\sigma}\gamma^{\tau}.$
(40)

In this $\gamma_{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 $\alpha_{s}$ for the axial
vector vertex is given by:

$\displaystyle Z_{5}=1-\frac{\alpha_{s}}{\pi}C_{F}.$ (41)

In order to obtain the SDCs using Eq.(31), in addition to the amplitude ${\cal
M}^{\rm nfr}_{Z\to(c\bar{c})[^{3}S_{1}^{[1]}]+(c\bar{c})[^{3}S_{1}^{[1]}]}$,
we also need to calculate the matrix element
$\langle(c\bar{c})[^{3}S_{1}^{[1]}]|\psi^{\dagger}{\bm{\sigma}}\cdot{\bm{\epsilon}}\chi|0\rangle$.
The matrix element up to order $\alpha_{s}$ is

$\displaystyle\sum_{\lambda}|\langle(c\bar{c})[^{3}S_{1}]|\psi^{\dagger}{\bm{\sigma}}\cdot{\bm{\epsilon}}(\lambda)\chi|0\rangle|^{2}=2N_{c}(D-1)(2E_{q})^{2},$
(42)

where $\lambda$ denotes the helicity of the $(c\bar{c})$ 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 integration-by-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].

###  II.3 Decay width

The total amplitude for the decay $Z\to J/\psi+J/\psi$ is the sum of the
fragmentation and non-fragmentation amplitudes, i.e.,

$\displaystyle\mathcal{M}_{Z\to J/\psi+J/\psi}=$
$\displaystyle\,\mathcal{M}^{\rm fr}_{Z\to J/\psi+J/\psi}+\mathcal{M}^{\rm
nfr}_{Z\to J/\psi+J/\psi}.$ (43)

Then the decay width can be calculated through

$\displaystyle\Gamma_{Z\to J/\psi+J/\psi}$
$\displaystyle=\frac{1}{2!}\frac{1}{3}\frac{\sqrt{m_{{}_{Z}}^{2}/4-m_{J/\psi}^{2}}}{8\pi
m_{{}_{Z}}^{2}}\sum\left|\mathcal{M}_{Z\to J/\psi+J/\psi}\right|^{2},$ (44)

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, $\sum$ 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:

$\displaystyle m_{{}_{Z}}=91.1876\,{\rm GeV},\;\;\;m_{c}=1.4\pm 0.2\,{\rm
GeV},$ $\displaystyle m_{J/\psi}=3.0969\,{\rm
GeV},\;\;\alpha(m_{J/\psi})=1/132.6,$ $\displaystyle G_{F}=1.16638\times
10^{-5}{\rm GeV}^{-2},$ (45)

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/\psi}$, and $G_{F}$ are
taken from the Particle Data Group (PDG) [45]. The value of the running
electromagnetic coupling $\alpha(m_{J/\psi})$ is taken from Ref.[46]. For the
running strong coupling constant, we use the two-loop formula

$\displaystyle\alpha_{s}(\mu_{r})=\frac{4\pi}{\beta_{0}\,{\rm
ln}(\mu_{r}^{2}/\Lambda_{\rm QCD}^{2})}\left[1-\frac{\beta_{1}{\rm ln}\,{\rm
ln}(\mu_{r}^{2}/\Lambda_{\rm QCD}^{2})}{\beta_{0}^{2}\,{\rm
ln}(\mu_{r}^{2}/\Lambda_{\rm QCD}^{2})}\right],$ (46)

where $\beta_{1}=102-38n_{f}/3$ is the two-loop coefficient of the QCD $\beta$
function. According to $\alpha_{s}(m_{{}_{Z}})=0.118$ [45], we obtain
$\Lambda^{n_{f}=5}_{\rm QCD}=0.226\,{\rm GeV}$ and $\Lambda^{n_{f}=4}_{\rm
QCD}=0.328\,{\rm GeV}$. With the determined values of $\Lambda_{\rm QCD}$, the
strong coupling constant at an arbitrary scale ($\mu_{r}\gg\Lambda_{\rm QCD}$)
can be directly calculated through Eq.(46).

The values of $g_{J/\psi\gamma}$ and $f_{J/\psi}^{\parallel}$ can be extracted
from the lepton decay width of $J/\psi$ through Eq.(4) and Eq.(29),
respectively. Using the measured value $\Gamma(J/\psi\to e^{+}e^{-})=5.55{\rm
keV}$ [47], we obtain

$\displaystyle g_{J/\psi\gamma}=-0.832\,{\rm GeV}^{2},$ (47) $\displaystyle
f_{J/\psi}^{\parallel}=403\,{\rm MeV}.$ (48)

The coupling constant $g_{J/\psi\gamma}$ has a relative minus sign compared to
$f_{J/\psi}^{\parallel}$, which can be realized by comparing Eq.(2) and
Eq.(28).

We take the values for the LDME
$\langle\mathcal{O}_{1}\rangle_{J/\psi}=\frac{1}{3}\sum_{\lambda}|\langle
J/\psi((\lambda)|\psi^{\dagger}{\bm{\sigma}}\cdot{\bm{\epsilon}}(\lambda)\chi|0\rangle|^{2}$
and the ratio of LDMEs $\langle{\bm{q}}^{2}\rangle_{J/\psi}$ from Ref.[44],
i.e.,

$\displaystyle\langle\mathcal{O}_{1}\rangle_{J/\psi}=0.440\,{\rm GeV}^{3},$
(49) $\displaystyle\langle{\bm{q}}^{2}\rangle_{J/\psi}=0.441\,{\rm GeV}^{2}.$
(50)

The value of $\langle v^{2}\rangle_{J/\psi}$ can be derived through $\langle
v^{2}\rangle_{J/\psi}=\langle\bm{q}^{2}\rangle_{J/\psi}/m_{c}^{2}$ directly.

###  III.1 Basic results

In this subsection, we present the numerical results for the decay width of
the decay $Z\to J/\psi+J/\psi$. In order to have a glance on the magnitude of
the photon fragmentation, non-fragmentation, and interference contributions to
the decay width of $Z\to J/\psi+J/\psi$, 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.

Table 1: The decay width (in unit: $10^{-12}\,{\rm GeV}$) of the decay $Z\to J/\psi+J/\psi$. The photon fragmentation, non-fragmentation, interference, and total contributions are given explicitly. | Fragmentation | Non-fragmentation | Interference | Total
---|---|---|---|---
$\Gamma$ | $154.0$ | $2.4$ | $60.2$ | $216.6$

Numerical results for the decay width are shown in Table 1, where the photon
fragmentation, non-fragmentation, interference, and total contributions are
given explicitly. In the calculation, we have set the renormalization scale in
the non-fragmentation amplitude as $\mu_{r}=m_{{}_{Z}}$.

From Table 1, 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 fragmentation contribution is suppressed by
the factor $(\alpha/\alpha_{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/\psi}$,
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/\psi})^{4}$ from the gauge boson propagator
compared to the non-fragmentation contribution.

###  III.2 Uncertainty analysis

In this subsection, we present an estimate of the theoretical uncertainties
for the decay width of $Z\to J/\psi+J/\psi$. The main uncertainty sources
include the factorization/renormalization scales, the charm quark mass, the
leptonic decay width $\Gamma(J/\psi\to e^{+}e^{-})$, the LDME
$\langle\mathcal{O}_{1}\rangle_{J/\psi}$, and the ratio
$\langle{\bm{q}}^{2}\rangle_{J/\psi}$. 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 $\mu_{0}$ and the final
factorization scale $\mu$ in the calculation of the photon fragmentation
amplitude, as well as the renormalization scale $\mu_{r}$ in the calculation
of the nonfragmentation amplitude. We estimate the uncertainties by varying
these scales in the ranges $\mu_{0}\in[1{\rm GeV},2m_{c}]$,
$\mu\in[m_{{}_{Z}}/2,2m_{{}_{Z}}]$, and
$\mu_{r}\in[m_{{}_{Z}}/2,2m_{{}_{Z}}]$. For the uncertainty caused by the
charm quark mass, we estimate it by taking $m_{c}=1.4\pm 0.2\,{\rm GeV}$. It
is noted that the values for the LDME $\langle\mathcal{O}_{1}\rangle_{J/\psi}$
and the ratio $\langle{\bm{q}}^{2}\rangle_{J/\psi}$ 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 $\langle\mathcal{O}_{1}\rangle_{J/\psi}$ and the ratio
$\langle{\bm{q}}^{2}\rangle_{J/\psi}$ on the charm quark mass. The dependence
of the LDME and the ratio $\langle{\bm{q}}^{2}\rangle_{J/\psi}$ on the charm
quark mass has been given in the Ref.[44]. In our calculation, the parameters
$g_{J/\psi\gamma}$ and $f_{J/\psi}^{\parallel}$ are extracted from the
measured value of $\Gamma(J/\psi\to e^{+}e^{-})$, and the LDME
$\langle\mathcal{O}_{1}\rangle_{J/\psi}$ and the ratio
$\langle{\bm{q}}^{2}\rangle_{J/\psi}$ obtained in Ref.[44] also depend on the
measured value of $\Gamma(J/\psi\to e^{+}e^{-})$. Thus, $\Gamma(J/\psi\to
e^{+}e^{-})$ is an important uncertainty source. For the uncertainty caused by
the measured value of $\Gamma(J/\psi\to e^{+}e^{-})$, we take the
uncertainties of the leptonic width as $\Gamma(J/\psi\to e^{+}e^{-})=5.55\pm
0.14\pm 0.02\,{\rm keV}$ [47]. Since the uncertainties caused by the charm
quark mass and the measured value of $\Gamma(J/\psi\to e^{+}e^{-})$ have
already been considered, we will omit the uncertainties caused by charm quark
mass and the measured value of $\Gamma(J/\psi\to e^{+}e^{-})$ when considering
the uncertainties caused by the LDME $\langle\mathcal{O}_{1}\rangle_{J/\psi}$
and the ratio $\langle{\bm{q}}^{2}\rangle_{J/\psi}$. The uncertainties for the
LDME $\langle\mathcal{O}_{1}\rangle_{J/\psi}$ and the ratio
$\langle{\bm{q}}^{2}\rangle_{J/\psi}$ with omitting the uncertainties caused
by the charm quark mass and the measured value of $\Gamma(J/\psi\to
e^{+}e^{-})$ [44] are

$\displaystyle\langle\mathcal{O}_{1}\rangle_{J/\psi}=0.440^{+0.065}_{-0.054}\,{\rm
GeV}^{3},$ (51)
$\displaystyle\langle{\bm{q}}^{2}\rangle_{J/\psi}=0.441^{+0.139}_{-0.138}\,{\rm
GeV}^{2}.$ (52)

Then we obtain the uncertainties arising from these uncertainty sources

$\displaystyle\Gamma_{Z\to J/\psi+J/\psi}$
$\displaystyle=2.16^{+0.06+0.32+0.12+0.10+0.07}_{-0.07-0.05-0.10-0.08-0.06}\times
10^{-10}\ {\rm GeV},$ (53)

where the first uncertainty is caused by the factorization/renormalization
scales 222The uncertainties arising from the scales $\mu_{0}$, $\mu$, and
$\mu_{r}$ have been added in quadrature., the second uncertainty is caused by
the charm quark mass, the third uncertainty is caused by the measured value of
$\Gamma(J/\psi\to e^{+}e^{-})$, the fourth uncertainty is cased by the LDME
$\langle\mathcal{O}_{1}\rangle_{J/\psi}$ and the last uncertainty is caused by
the ratio $\langle{\bm{q}}^{2}\rangle_{J/\psi}$. Adding these uncertainties in
quadrature, we obtain the total theoretical uncertainty for the decay width

$\displaystyle\Gamma_{Z\to J/\psi+J/\psi}=2.16^{+0.37}_{-0.17}\times 10^{-10}\
{\rm GeV}.$ (54)

Using the decay width obtained in the present paper and the total $Z$-boson
decay width $\Gamma_{{}_{Z}}=2.4955\,{\rm GeV}$ from PDG [45], we obtain the
branching fraction for the decay $Z\to J/\psi+J/\psi$

$\displaystyle{\rm Br}(Z\to J/\psi+J/\psi)=8.66^{+1.48}_{-0.69}\times
10^{-11}.$ (55)

###  III.3 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\to
J/\psi+J/\psi$.

Table 2: The branching fraction of the decay $Z\to J/\psi+J/\psi$. 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] on this branching fraction. Br(this work) |  ${\rm Br}$(Ref.[11]) |  ${\rm Br}$(Ref.[2])
---|---|---
$8.66^{+1.48}_{-0.69}\times 10^{-11}$ | $1.110^{+0.334+0.054}_{-0.241-0.001}\times 10^{-10}$ | $<1.4\times 10^{-6}$

The theoretical predictions and the current experimental upper limit on the
branching fraction of the decay $Z\to J/\psi+J/\psi$ are shown in Table 2.
From the table, we can see that our branching fraction for $Z\to
J/\psi+J/\psi$ is smaller than that in Ref.[11] by about $22\%$, which is
about $-1.0\sigma$ in the uncertainties of Ref.[11].

In Ref.[11], the photon fragmentation, non-fragmentation, and the interference
contributions were calculated up to NLO in $\alpha_{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\to J/\psi+\gamma^{*}$ and the
amplitude of $\gamma^{*}\to J/\psi$. 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 calculated through the light-cone
approach where the large logarithms of $m_{{}_{Z}}^{2}/m_{c}^{2}$ are resummed
through the evolution of the $J/\psi$ LCDA. Using this improved 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}$
in the amplitude of $Z\to J/\psi+\gamma^{*}$ 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 2, we can see that our branching fraction for $Z\to J/\psi+J/\psi$
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\to
J/\psi+J/\psi$. The decay process $Z\to J/\psi+J/\psi$ is dominated by the
single photon fragmentation mechanism, i.e., the transition $Z\to
J/\psi+\gamma^{*}$ followed by the fragmentation $\gamma^{*}\to J/\psi$. The
amplitude of $\gamma^{*}\to J/\psi$ has already been calculated up to ${\cal
O}(\alpha_{s}^{3})$ under the NRQCD factorization in Refs.[48, 49]. It seems
that the perturbative expansion of the amplitude of $\gamma^{*}\to J/\psi$
does not converge in those calculations. Moreover, the higher-order
corrections of the amplitude of $Z\to J/\psi+\gamma^{*}$ involve the large
logarithms of $m_{{}_{Z}}^{2}/m_{c}^{2}$. Considering these problems, we have
presented an improved calculation for the decay rate of $Z\to J/\psi+J/\psi$
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\to
J/\psi+\gamma^{*}$ and the amplitude of $\gamma^{*}\to J/\psi$. Then, 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 calculated by the light-cone approach where the large logarithms of
$m_{{}_{Z}}^{2}/m_{c}^{2}$ are resummed through solving the evolution equation
of the $J/\psi$ LCDA. The advantage of our calculation is that 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}$ in the
amplitude of $Z\to J/\psi+\gamma^{*}$ are resummed. For the non-fragmentation
amplitude, we adopt the fixed-order approach under the NRQCD factorization. In
our calculation of the non-fragmentation amplitude, the NLO QCD and NLO
relativistic corrections are included simultaneously.

The obtained branching fraction for $Z\to J/\psi+J/\psi$ is
$8.66^{+1.48}_{-0.69}\times 10^{-11}$, which is compatible with the current
experimental upper limit ${\rm Br}(Z\to J/\psi+J/\psi)=1.4\times 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\times 10^{12}$ $Z$
bosons can be generated at the FCC-ee [6]. According to the branching fraction
obtained in this work, about $400$ $Z\to J/\psi+J/\psi$ 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.

Acknowledgments: We thank Hai-Bing Fu for helpful discussions on the evolution
of the light-cone distribution amplitude of $J/\psi$. This work was supported
in part by the Natural Science Foundation of China under Grants No. 12005028,
No. 12175025, No. 12275036 and No. 12347101, by the Fundamental Research Funds
for the Central Universities under Grant No. 2020CQJQY-Z003, by the Chongqing
Natural Science Foundation under Grants No. CSTB2022NSCQ-MSX0415 and No.
cstc2021jcyj-msxmX0681, and by the Chongqing Graduate Research and Innovation
Foundation under Grant No. ydstd1912.

## References

  * [1] A. M. Sirunyan et al. [CMS], Search for Higgs and $Z$ boson decays to J/$\psi$ or Y pairs in the four-muon final state in proton-proton collisions at s=13TeV, Phys. Lett. B 797, 134811 (2019).
  * [2] A. Tumasyan et al. [CMS], Search for Higgs boson decays into Z and J/$\psi$ and for Higgs and $Z$ boson decays into J/$\psi$ or Y pairs in pp collisions at s=13 TeV, Phys. Lett. B 842, 137534 (2023).
  * [3] D. d’Enterria and V. Le, Rare and exclusive few-body decays of the Higgs, Z, W bosons, and the top quark, arXiv:2312.11211 [hep-ph].
  * [4] H. Baer et al. [ILC], The International Linear Collider Technical Design Report - Volume 2: Physics, arXiv:1306.6352 [hep-ph].
  * [5] J. B. Guimarães da Costa et al. [CEPC Study Group], CEPC Conceptual Design Report: Volume 2 - Physics & Detector, arXiv:1811.10545 [hep-ex].
  * [6] A. Abada et al. [FCC], FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2, Eur. Phys. J. ST 228, 261-623 (2019).
  * [7] J. P. Ma and Z. X. Zhang (The super Z-factory group), Preface, Sci. China: Phys., Mech. Astron. 53, 1947 (2010).
  * [8] L. Bergstrom and R. W. Robinett, ON THE RARE DECAYS $Z\to VV$ AND $Z\to VP$, Phys. Rev. D 41, 3513 (1990).
  * [9] A. K. Likhoded and A. V. Luchinsky, Double Charmonia Production in Exclusive $Z$ Boson Decays, Mod. Phys. Lett. A 33, 1850078 (2018).
  * [10] D. N. Gao and X. Gong, Note on rare Z-boson decays to double heavy quarkonia, Chin. Phys. C 47, 043106 (2023).
  * [11] C. Li, Z. Sun and G. Y. Zhang, Analysis of double-J/$\psi$ production in Z decay at next-to-leading-order QCD accuracy, JHEP 10, 120 (2023).
  * [12] X. Luo, H. B. Fu, H. J. Tian and C. Li, Next-to-leading-order QCD correction to the exclusive double charmonium production via $Z$ decays, arXiv:2209.08802 [hep-ph].
  * [13] G. P. Lepage and S. J. Brodsky, Exclusive Processes in Perturbative Quantum Chromodynamics, Phys. Rev. D 22, 2157 (1980).
  * [14] V. L. Chernyak and A. R. Zhitnitsky, Asymptotic Behavior of Exclusive Processes in QCD, Phys. Rept. 112, 173 (1984).
  * [15] G. T. Bodwin, E. Braaten and G. P. Lepage, Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium, Phys. Rev. D 51, 1125-1171 (1995).
  * [16] G. T. Bodwin, E. Braaten, J. Lee and C. Yu, Exclusive two-vector-meson production from e+ e- annihilation, Phys. Rev. D 74, 074014 (2006).
  * [17] G. T. Bodwin, F. Petriello, S. Stoynev and M. Velasco, Higgs boson decays to quarkonia and the $H\bar{c}c$ coupling, Phys. Rev. D 88, 053003 (2013).
  * [18] B. Guberina, J. H. Kuhn, R. D. Peccei and R. Ruckl, Rare Decays of the Z0, Nucl. Phys. B 174, 317 (1980).
  * [19] A. V. Luchinsky, Leading order NRQCD and Light-Cone Analysis of Exclusive Charmonia Production in Radiative $Z$-boson Decays, arXiv:1706.04091.
  * [20] X. P. Wang and D. Yang, The leading twist light-cone distribution amplitudes for the S-wave and P-wave quarkonia and their applications in single quarkonium exclusive productions, JHEP 06, 121 (2014).
  * [21] G. T. Bodwin, H. S. Chung, J. H. Ee and J. Lee, $Z$-boson decays to a vector quarkonium plus a photon, Phys. Rev. D 97, 016009 (2018).
  * [22] W. L. Sang, D. S. Yang and Y. D. Zhang, $Z$-boson radiative decays to an S-wave quarkonium at NNLO and NLL accuracy, Phys. Rev. D 108, 014021 (2023).
  * [23] G. Y. Wang, X. C. Zheng, X. G. Wu and G. Z. Xu, $Z$-boson decays into $S$-wave quarkonium plus a photon up to ${\cal O}(\alpha_{s}v^{2})$ corrections, Phys. Rev. D 109, 074004 (2024).
  * [24] P. Ball and V. M. Braun, The Rho meson light cone distribution amplitudes of leading twist revisited, Phys. Rev. D 54, 2182-2193 (1996).
  * [25] P. Ball, V. M. Braun, Y. Koike and K. Tanaka, Higher twist distribution amplitudes of vector mesons in QCD: Formalism and twist - three distributions, Nucl. Phys. B 529, 323-382 (1998).
  * [26] J. P. Ma and Z. G. Si, NRQCD Factorization for Twist-2 Light-Cone Wave-Functions of Charmonia, Phys. Lett. B 647, 419-426 (2007).
  * [27] W. Wang, J. Xu, D. Yang and S. Zhao, Relativistic corrections to light-cone distribution amplitudes of S-wave Bc mesons and heavy quarkonia, JHEP 12, 012 (2017).
  * [28] A. V. Efremov and A. V. Radyushkin, Factorization and Asymptotical Behavior of Pion Form-Factor in QCD, Phys. Lett. B 94, 245-250 (1980).
  * [29] A. V. Efremov and A. V. Radyushkin, Asymptotical Behavior of Pion Electromagnetic Form-Factor in QCD, Theor. Math. Phys. 42, 97-110 (1980).
  * [30] F. M. Dittes and A. V. Radyushkin, Two loop contribution to the evolution of the pion wave function, Phys. Lett. B 134, 359-362 (1984).
  * [31] G. R. Katz, Two Loop Feynman Gauge Calculation of the Meson Nonsinglet Evolution Potential, Phys. Rev. D 31, 652 (1985).
  * [32] S. V. Mikhailov and A. V. Radyushkin, Evolution Kernels in QCD: Two Loop Calculation in Feynman Gauge, Nucl. Phys. B 254, 89-126 (1985).
  * [33] S. S. Agaev, V. M. Braun, N. Offen and F. A. Porkert, Light Cone Sum Rules for the pi0-gamma*-gamma Form Factor Revisited, Phys. Rev. D 83, 054020 (2011).
  * [34] G. T. Bodwin, H. S. Chung, J. H. Ee and J. Lee, New approach to the resummation of logarithms in Higgs-boson decays to a vector quarkonium plus a photon, Phys. Rev. D 95, 054018 (2017).
  * [35] G. T. Bodwin, X. Garcia Tormo, i and J. Lee, Factorization in exclusive quarkonium production, Phys. Rev. D 81, 114014 (2010).
  * [36] M. Beneke and V. A. Smirnov, Asymptotic expansion of Feynman integrals near threshold, Nucl. Phys. B 522, 321 (1998).
  * [37] S. A. Larin, The Renormalization of the axial anomaly in dimensional regularization, Phys. Lett. B 303, 113-118 (1993).
  * [38] T. Hahn, Generating Feynman diagrams and amplitudes with FeynArts 3, Comput. Phys. Commun 140, 418 (2001).
  * [39] R. Mertig, M. Bohm and A. Denner, FeynCalc - Computer-algebraic calculation of Feynman amplitudes, Comput. Phys. Commun 64, 345 (1991).
  * [40] V. Shtabovenko, R. Mertig and F. Orellana, New Developments in FeynCalc 9.0, Comput. Phys. Commun 207, 432 (2016).
  * [41] F. Feng, $Apart: A Generalized Mathematica Apart Function, Comput. Phys. Commun 183, 2158 (2012).
  * [42] A.V. Smirnov, Algorithm FIRE - Feynman Integral REduction, J. High Energy Phys. 0810, 107 (2008).
  * [43] T. Hahn and M. Perez-Victoria, Automatized one loop calculations in four-dimensions and D-dimensions, Comput. Phys. Commun 118, 153 (1999).
  * [44] G. T. Bodwin, H. S. Chung, D. Kang, J. Lee and C. Yu, Improved determination of color-singlet nonrelativistic QCD matrix elements for S-wave charmonium, Phys. Rev. D 77, 094017 (2008).
  * [45] R. L. Workman et al. [Particle Data Group], Review of Particle Physics, PTEP 2022, 083C01 (2022).
  * [46] G. T. Bodwin, J. Lee and C. Yu, Resummation of Relativistic Corrections to $e+e-\to J/psi+eta(c)$, Phys. Rev. D 77, 094018 (2008).
  * [47] W. M. Yao et al. [Particle Data Group], Review of Particle Physics, J. Phys. G 33, 1-1232 (2006).
  * [48] P. Marquard, J. H. Piclum, D. Seidel and M. Steinhauser, Three-loop matching of the vector current, Phys. Rev. D 89, 034027 (2014).
  * [49] M. Egner, M. Fael, F. Lange, K. Schönwald and M. Steinhauser, Three-loop nonsinglet matching coefficients for heavy quark currents, Phys. Rev. D 105, 114007 (2022).