Studies of $Z$ boson decay into double $\Upsilon$ mesons at the NLO QCD accuracy

Cong Li Ying-Zhao Jiang Zhan Sun sunzhan˙[email protected] Department of Physics, Guizhou Minzu University, Guiyang 550025, People’s Republic of China.

Abstract

In this paper, we employ the nonrelativistic QCD factorization to conduct a comprehensive examination of the $Z$ boson decay into a pair of $\Upsilon$ mesons, achieving accuracy at the next-to-leading-order (NLO) in $\alpha_{s}$ . Our calculations demonstrate that the QED diagrams are indispensable in comparison to the pure QCD diagrams, and the implementation of QCD corrections markedly enhance the QCD results, whereas it substantially diminish the QED results. To ensure consistency with the experimental methodology, we have taken into account the feed-down transitions originating from higher excited states, which exhibit significant relevance. Combining all the contributions, we arrive at the NLO prediction of $\mathcal{B}_{Z\to\Upsilon(nS)+\Upsilon(mS)}\sim 10^{-11}$ , which is notably lower than the upper limits set by CMS.

pacs:

12.38.Bx, 12.39.Jh, 14.40.Pq

I Introduction

Heavy-quarkonium production in $Z$ boson decay has garnered significant interest from both theoretical and experimental physicists over the past decades z decay 1 ; z decay 2 ; z decay 3 ; z decay 4 ; z decay 5 ; z decay 6 ; z decay 7 ; z decay 8 ; z decay 9 ; z decay 10 ; z decay 11 ; z decay 12 ; z decay 13 ; z decay 14 ; z decay 15 ; z decay 16 ; z decay 17 ; z decay 18 ; z decay 19 ; z decay 20 ; z decay 21 ; z decay 22 ; z decay 23 ; z decay 24 ; z decay 25 ; z decay 26 ; z decay 27 ; z decay 28 ; z decay 29 ; z decay 30 ; z decay 31 ; z decay 32 ; z decay 33 ; z decay 34 ; z decay 35 ; z decay 36 ; z decay 37 ; z decay 38 ; z decay 39 ; z decay 40 ; z decay 41 ; z decay 42 ; z decay 43 . In the various decay channels of the $Z$ boson into heavy quarkonium, the yield of double $J/\psi$ (or $\Upsilon$ ) mesons holds a notable advantage due to the unique experimental signature arising from their subsequent decay into four muons, which can distinctly be recognized and analyzed.

In 2019, the CMS Collaboration measured the branching ratio of the $Z$ boson decaying to double $J/\psi$ CMS:2019wch , denoted as $\mathcal{B}_{Z\to J/\psi+J/\psi}<2.2\times 10^{-6}$ . This upper limit was recently updated to be $1.4\times 10^{-6}$ by analyzing a larger sample of data CMS:2022fsq . At leading order (LO) in $\alpha_{s}$ , the standard model provides a prediction of $\mathcal{B}_{Z\to J/\psi+J/\psi}\sim 10^{-12}$ Likhoded:2017jmx , which was upgraded to $10^{-10}$ by further evaluating the QED contributions resulting from the virtual-photon effects Gao:2022mwa . Recently, Li $et~{}al.$ suggested that the QCD corrections play a crucial role in enhancing the contributions from QCD diagrams and simultaneously diminishing the contributions from QED diagrams, ultimately maintaining the next-to-leading-order (NLO) results at the order of $10^{-10}$ z decay 40 .

In addition to the double $J/\psi$ yield, the CMS group has also explored the Z boson decaying into a pair of $\Upsilon$ mesons, thereby establishing upper limits on the branching ratio CMS:2019wch ; CMS:2022fsq

	$\displaystyle\mathcal{B}_{Z\to\Upsilon(mS)+\Upsilon(nS)}\leq 3.9\times 10^{-7},$
	$\displaystyle\mathcal{B}_{Z\to\Upsilon(1S)+\Upsilon(1S)}\leq 1.8\times 10^{-6}.$		(1)

By taking into account both the QCD and QED diagrams, Gao $et~{}al.$ provided an estimation of $\mathcal{B}_{Z\to\Upsilon(1S)+\Upsilon(1S)}$ at the LO QCD accuracy Gao:2022mwa . In light of the considerable effect that NLO QCD corrections have on double-charmonia production in $e^{-}e^{+}$ annihilation Zhang:2005cha ; Gong:2007db ; Zhang:2008gp ; Brambilla:2010cs ; Dong:2011fb ; Sun:2018rgx ; Sun:2021tma , it is prudent to investigate whether high-order terms in $\alpha_{s}$ could similarly engender a significant boost in the yield of double $\Upsilon$ mesons in $Z$ boson decay. To achieve this, in the present work, we will examine the decay of $Z$ boson into $\Upsilon$ pair with NLO precision in $\alpha_{s}$ , incorporating both QCD and QED diagrams within the nonrelativistic QCD (NRQCD) formalism NRQCD1 . Note that, the measured branching ratio of $Z$ boson decay into double $\Upsilon$ mesons includes considerations for feed-down transitions. Consequently, this evaluation will also consider the impact of higher excited states.

It is noteworthy that, the considerable mass of the $b$ quark typically leads to a more rapidly converging perturbative series in the expansion of $\alpha_{s}$ and $v^{2}$ , compared to the case of the $c$ quark. Coupled with the significant production rates of Z bosons at the LHC, the proposed HL-LHC, or the future CEPC CEPC , the yields of double $\Upsilon$ mesons in $Z$ decay would offer an excellent laboratory for probing the mechanisms governing heavy quarkonium formation.

The remainder of this paper is structured as follows: Section II provides an outline of the calculation formalism. This is followed by the presentation of phenomenological results and discussions in Section III. Finally, Section IV is dedicated to a summary.

II Calculation formalism

II.1 Theoretical Framework

As previously discussed, the CMS measurements have incorporated feed-down transitions; consequently, our calculations will assess the contributions stemming from higher excited states, including those of the $\Upsilon(mS)$ and $\chi_{b}(nP)$ mesons. Within the NRQCD framework NRQCD1 ; NRQCD2 , the decay width of $Z\to H_{1}(b\bar{b})+H_{2}(b\bar{b})$ can be factorized as

\displaystyle\Gamma=\hat{\Gamma}_{Z\to b\bar{b}[n_{1}]+b\bar{b}[n_{2}]}\langle\mathcal{O}^{H_{1}}(n_{1})\rangle\langle\mathcal{O}^{H_{2}}(n_{2})\rangle,

(2)

where $\hat{\Gamma}_{Z\to b\bar{b}[n_{1}]+b\bar{b}[n_{2}]}$ represents the perturbative calculable short-distance coefficients (SDCs), denoting the production of the intermediate states consisting of $b\bar{b}[n_{1}]$ and $b\bar{b}[n_{2}]$ . With the restriction to color-singlet contributions, we have $n_{1}=^{3}S_{1}^{1}$ and $n_{2}=^{3}S_{1}^{1}(\textrm{or}~{}^{3}P_{J}^{1})$ ( $J=0,1,2$ ).¹¹1When addressing the feed-down transitions, we exclude the evaluation of $Z\to\chi_{b}(nP)+\chi_{b}(n^{\prime}P)$ . This process is anticipated to contribute $Z\to\Upsilon(mS)+\Upsilon(m^{\prime}S)$ through a multi-step decay, which involves at least two decay stages. The universal nonperturbative long distant matrix element (LDME) $\langle\mathcal{O}^{H_{1(2)}}(n_{1(2)})\rangle$ stands for the probabilities of transitions from $b\bar{b}[n_{1(2)}]$ into the $H_{1(2)}$ meson states.

The $\hat{\Gamma}_{Z\to b\bar{b}[n_{1}]+b\bar{b}[n_{2}]}$ can further be expressed as

\displaystyle\hat{\Gamma}_{Z\to b\bar{b}[n_{1}]+b\bar{b}[n_{2}]}=\frac{\kappa}{2m_{Z}}\frac{1}{N_{I}N_{s}}|\mathcal{M}|^{2},

(3)

where $|\mathcal{M}|^{2}$ is the squared matrix elements, $1/N_{s}$ is the spin average factor of the initial $Z$ boson multiplied by the identity factor $N_{I}$ of the two final $b\bar{b}$ states. For $b\bar{b}[n_{1}]$ and $b\bar{b}[n_{2}]$ representing the same states, $N_{I}$ is set to $1/2!$ ; however, if they are distinct states, $N_{I}$ equals 1. The symbol $\kappa$ denotes the factor originating from the standard two-body phase space.

Based on the framework used to handle $Z-J/\psi+J/\psi$ z decay 40 , we include terms up to the $\alpha^{3}$ order and achieve NLO accuracy in $\alpha_{s}$ . The squared matrix elements specified in Eq. (3) can subsequently be expressed as follows

	$\displaystyle\bigg{\|}\left(\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha_{s}}+\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}\right)+\left(\mathcal{M}_{\alpha^{\frac{3}{2}}}+\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}\right)\bigg{\|}^{2}$	(4)
$\displaystyle=$	$\displaystyle\big{\|}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha_{s}}\big{\|}^{2}+2\textrm{Re}\left(\mathcal{M}^{*}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}\right)$
	$\displaystyle+2\textrm{Re}\left(\mathcal{M}^{}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{3}{2}}}\right)+2\textrm{Re}\left(\mathcal{M}^{}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}+\mathcal{M}^{*}_{\alpha^{\frac{3}{2}}}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}\right)$
	$\displaystyle+\big{\|}\mathcal{M}_{\alpha^{\frac{3}{2}}}\big{\|}^{2}+2\textrm{Re}\left(\mathcal{M}^{*}_{\alpha^{\frac{3}{2}}}\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}\right)+\cdots.$

Accordingly, we decompose the SDCs into three distinct components,

\displaystyle\hat{\Gamma}=\hat{\Gamma}^{(0,1)}_{1}+\hat{\Gamma}^{(0,1)}_{2}+\hat{\Gamma}^{(0,1)}_{3},

(5)

where

$\displaystyle\hat{\Gamma}^{(0)}_{1}$	$\displaystyle\propto$	$\displaystyle\big{\|}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha_{s}}\big{\|}^{2},$
$\displaystyle\hat{\Gamma}^{(1)}_{1}$	$\displaystyle\propto$	$\displaystyle 2\textrm{Re}\left(\mathcal{M}^{*}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}\right),$
$\displaystyle\hat{\Gamma}^{(0)}_{2}$	$\displaystyle\propto$	$\displaystyle 2\textrm{Re}\left(\mathcal{M}^{*}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{3}{2}}}\right),$
$\displaystyle\hat{\Gamma}^{(1)}_{2}$	$\displaystyle\propto$	$\displaystyle 2\textrm{Re}\left(\mathcal{M}^{}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}+\mathcal{M}^{}_{\alpha^{\frac{3}{2}}}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}\right),$
$\displaystyle\hat{\Gamma}^{(0)}_{3}$	$\displaystyle\propto$	$\displaystyle\big{\|}\mathcal{M}_{\alpha^{\frac{3}{2}}}\big{\|}^{2},$
$\displaystyle\hat{\Gamma}^{(1)}_{3}$	$\displaystyle\propto$	$\displaystyle 2\textrm{Re}\left(\mathcal{M}^{*}_{\alpha^{\frac{3}{2}}}\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}\right).$	(6)

The subscripts 1, 2, and 3 represent the order in $\alpha$ , while the superscript $0(1)$ denotes the terms at the LO (or NLO) level in $\alpha_{s}$ .

Refer to caption — Figure 1: Representative QCD Feynman diagrams for $Z\to b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}S_{1}^{1},^{3}P_{J}^{1}]$ . Diagram (a) ( ${\alpha^{\frac{1}{2}}\alpha_{s}}$ order) is the QCD tree-level diagram. Diagrams (b)-(f) $({\alpha^{\frac{1}{2}}\alpha^{2}_{s}}~{}\textrm{order})$ depict the NLO QCD corrections to (a). Diagram (b) specifically represents the counter-term diagram.

The representative Feynman diagrams of $\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha_{s}}$ , $\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}$ , $\mathcal{M}_{\alpha^{\frac{3}{2}}}$ , and $\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}$ are displayed in Figs. 1 and 2. Fig. 1(a) ( $\alpha^{\frac{1}{2}}\alpha_{s}$ order) demonstrates the QCD tree-level diagram (4 diagrams), whereas Figs. 1(b)-(f) elucidate the subsequent QCD NLO corrections, which comprise 56 individual one-loop diagrams and 20 counter-term diagrams. In Figs. 2(a,b) ( $\alpha^{\frac{3}{2}}$ order) we observe the tree-level diagrams pertinent to QED, encompassing 8 diagrams in total. Furthermore, Figs. 2(c)-(k) depict the high-order corrections in $\alpha_{s}$ , including 52 one-loop diagrams and 32 counter-term diagrams.

II.2 The decay width

According to Eqs. (2) and (3), the decay width can generally be written as

\displaystyle\Gamma^{\textrm{(0,1)}}_{i}=\frac{\kappa}{2m_{Z}}\frac{1}{N_{I}N_{s}}\times\zeta_{i}\times\mathcal{C}^{(0,1)}_{i}\times\langle\mathcal{O}^{H_{1}}(^{3}S_{1}^{1})\rangle\langle\mathcal{O}^{H_{2}}(^{3}S_{1}^{1},^{3}P_{J}^{1})\rangle,~{}i=1,2,3,

(7)

where $\kappa$ is equal to $\frac{\sqrt{m^{2}_{Z}-16m^{2}_{b}}}{8\pi m_{Z}}$ and $N_{s}=3$ .

For the yield of $b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}S_{1}^{1}]$ , the $\zeta_{i}$ takes the following form ( $e_{b}=-\frac{1}{3}$ )

\displaystyle\zeta_{i}

\displaystyle=

\displaystyle\frac{\pi^{3}}{16m_{b}^{6}\sin^{2}\theta_{\textrm{w}}\cos^{2}\theta_{\textrm{w}}}e^{2(i-1)}_{b}\alpha^{i}\alpha^{3-i}_{s},

(8)

while for the scenario involving $b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{1}]$ ,

\displaystyle\zeta_{i}

\displaystyle=

\displaystyle\frac{\pi^{3}(4\sin^{2}\theta_{\textrm{w}}-3)^{2}}{144m_{b}^{6}\sin^{2}\theta_{\textrm{w}}\cos^{2}\theta_{\textrm{w}}}e^{2(i-1)}_{b}\alpha^{i}\alpha^{3-i}_{s}.

(9)

The LDMEs $\langle\mathcal{O}^{H}(^{3}S_{1}^{1},^{3}P_{J}^{1})\rangle$ can be expressed in terms of the wave functions at the origin by utilizing the following formulae

	$\displaystyle\langle\mathcal{O}^{\Upsilon(nS)}(^{3}S_{1}^{1})\rangle$	$\displaystyle=$	$\displaystyle\frac{1}{4\pi}\|R_{nS}(0)\|^{2},$
	$\displaystyle\langle\mathcal{O}^{\chi_{bJ}(mP)}(^{3}P_{J}^{1})\rangle$	$\displaystyle=$	$\displaystyle\frac{3}{4\pi}\|R^{\prime}_{mP}(0)\|^{2}.$		(10)

II.2.1 LO

The coefficients of $\mathcal{C}^{(0)}_{i}$ derived from the LO processes, as illustrated in Figs. 1(a) and 2(a,b) are free from divergences. In the following, we provide the expressions for various processes ( $r\equiv\frac{m^{2}_{Z}}{4m_{b}^{2}}$ ).

$b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}S_{1}^{1}]$

	$\displaystyle\mathcal{C}^{(0)}_{1}=\frac{65536m_{b}^{2}(r^{2}-10r+24)}{9r^{4}},$
	$\displaystyle\mathcal{C}^{(0)}_{2}=\frac{16384m_{b}^{2}(3r+2)(r^{2}-10r+24)}{3r^{4}},$
	$\displaystyle\mathcal{C}^{(0)}_{3}=\frac{1024m_{b}^{2}(3r+2)^{2}(r^{2}-10r+24)}{r^{4}}.$		(11)

$b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{0}]$

	$\displaystyle\mathcal{C}^{(0)}_{1}=\frac{16384(r^{4}+182r^{3}-428r^{2}+152r+144)}{27r^{5}},$
	$\displaystyle\mathcal{C}^{(0)}_{2}=\frac{32768(7r^{4}-7r^{3}-59r^{2}+56r+36)}{9r^{5}},$
	$\displaystyle\mathcal{C}^{(0)}_{3}=\frac{512(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}{3r^{5}}.$		(12)

$b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{1}]$

	$\displaystyle\mathcal{C}^{(0)}_{1}=\frac{131072(2r^{3}-12r^{2}+13r+18)}{9r^{5}},$
	$\displaystyle\mathcal{C}^{(0)}_{2}=-\frac{16384(19r^{3}-24r^{2}-88r-72)}{3r^{5}},$
	$\displaystyle\mathcal{C}^{(0)}_{3}=\frac{1024(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}{r^{5}}.$		(13)

$b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{2}]$

	$\displaystyle\mathcal{C}^{(0)}_{1}=\frac{32768(r^{4}+20r^{3}-188r^{2}+308r+360)}{27r^{5}},$
	$\displaystyle\mathcal{C}^{(0)}_{2}=\frac{16384(10r^{4}-37r^{3}-104r^{2}+488r+360)}{9r^{5}},$
	$\displaystyle\mathcal{C}^{(0)}_{3}=\frac{1024(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}{3r^{5}}.$		(14)

II.2.2 NLO

We utilize the dimensional regularization with $D=4-2\epsilon$ to isolate the ultraviolet (UV) and infrared (IR) divergences. The on-mass-shell (OS) scheme is employed to set the renormalization constants for the $b$ -quark mass ( $Z_{m}$ ) and heavy-quark filed ( $Z_{2}$ ); the minimal-subtraction ( $\overline{MS}$ ) scheme is adopted for the QCD-gauge coupling ( $Z_{g}$ ) and the gluon filed $Z_{3}$ . The renormalization constants are taken as

$\displaystyle\delta Z_{m}^{OS}$	$\displaystyle=$	$\displaystyle-3C_{F}\frac{\alpha_{s}}{4\pi}N_{\epsilon}\left[\frac{1}{\epsilon_{\textrm{UV}}}+\frac{4}{3}+2\textrm{ln}{2}\right],$
$\displaystyle\delta Z_{2}^{OS}$	$\displaystyle=$	$\displaystyle-C_{F}\frac{\alpha_{s}}{4\pi}N_{\epsilon}\left[\frac{1}{\epsilon_{\textrm{UV}}}+\frac{2}{\epsilon_{\textrm{IR}}}+4+6\textrm{ln}{2}\right],$
$\displaystyle\delta Z_{3}^{\overline{MS}}$	$\displaystyle=$	$\displaystyle\frac{\alpha_{s}}{4\pi}(\beta_{0}-2C_{A})N_{\epsilon}\left[\frac{1}{\epsilon_{\textrm{UV}}}+\textrm{ln}\frac{4m_{b}^{2}}{\mu_{r}^{2}}\right],$
$\displaystyle\delta Z_{g}^{\overline{MS}}$	$\displaystyle=$	$\displaystyle-\frac{\beta_{0}}{2}\frac{\alpha_{s}}{4\pi}N_{\epsilon}\left[\frac{1}{\epsilon_{\textrm{UV}}}+\textrm{ln}\frac{4m_{b}^{2}}{\mu_{r}^{2}}\right],$	(15)

where $N_{\epsilon}=\frac{1}{\Gamma[1-\epsilon]}\left(\frac{4\pi\mu_{r}^{2}}{4m_{b}^{2}}\right)^{\epsilon}$ is an overall factor, $\gamma_{E}$ is the Euler’s constant, and $\beta_{0}=\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{f}$ is the one-loop coefficient of the $\beta$ function. $n_{f}(=n_{L}+n_{H})$ represents the number of the active-quark flavors; $n_{L}(=4)$ and $n_{H}(=1)$ denote the number of the light- and heavy-quark flavors, respectively.²²2In our calculations, the $c$ quark is considered a light quark in comparison to the $b$ quark; therefore, we assign $n_{l}=4$ for the number of light quarks (encompassing $u,d,s$ , and $c$ ) and $n_{h}=1$ for the number of heavy quarks (specifically, the $b$ quark). In ${\rm SU}(3)$ , the color factors are given by $T_{F}=\frac{1}{2}$ , $C_{F}=\frac{4}{3}$ , and $C_{A}=3$ .

Incorporating the QCD corrections allows us to derive the coefficients for $\mathcal{C}^{(1)}_{i}$ , which can be formulated in a generic manner

\displaystyle\mathcal{C}^{(0)}_{i}+\mathcal{C}^{(1)}_{i}=\mathcal{C}^{(0)}_{i}\left[1+\frac{\alpha_{s}}{\pi}\left(\xi_{i}\beta_{0}\textrm{ln}\frac{\mu_{r}^{2}}{4m_{b}^{2}}+a_{i}n_{L}+b_{i}n_{H}+c_{i}\right)\right],

(16)

where $\xi_{1}=\frac{1}{2}$ , $\xi_{2}=\frac{1}{4}$ , and $\xi_{3}=0$ . The coefficients $a_{i}$ , $b_{i}$ , and $c_{i}$ are dependent solely on the variables of $r$ and $m_{b}$ . Their fully analytical expressions can be found in Appendices 33-56.

To generate all the necessary Feynman diagrams and corresponding analytical amplitudes, we use the FeynArts package Hahn:2000kx . We then employ the FeynCalc package Mertig:1990an to handle the traces of the $\gamma$ and color matrices, which transforms the hard scattering amplitudes into expressions with loop integrals. When calculating the $D$ -dimensional $\gamma$ traces that incorporate a single $\gamma_{5}$ matrix and involve UV and/or IR divergences, we follow the scheme outlined in Refs. Korner:1991sx ; z decay 4 ; z decay 22 and choose the same starting point ( $Z$ -vertex) to write down the amplitudes without implementation of cyclicity. Subsequently, we utilize our self-written Mathematica codes that include implementations of Apart Feng:2012iq and FIRE Smirnov:2008iw to reduce these loop integrals to a set of irreducible Master Integrals (MIs). The fully-analytical expressions for these MIs can be found in Appendix 21-32. As a cross check, we simultaneously adopt the LoopTools Hahn:1998yk package to numerically evaluate these MIs, obtaining the same numerical results.

III Phenomenological results

The parameters incorporated into our calculations are defined as follows: $m_{b}=4.7\pm 0.1$ GeV, $m_{Z}=91.1876$ GeV, and $\alpha=1/128$ . Additionally, we employ the two-loop $\alpha_{s}$ running coupling constant in our analysis.

The wave functions in Eq. (10) read Eichten:1995ch

$\displaystyle\|R_{1S}(0)\|^{2}$	$\displaystyle=$	$\displaystyle 6.477~{}\textrm{GeV}^{3},~{}~{}~{}~{}~{}~{}\|R_{2S}(0)\|^{2}=3.234~{}\textrm{GeV}^{3},$
$\displaystyle\|R_{3S}(0)\|^{2}$	$\displaystyle=$	$\displaystyle 2.474~{}\textrm{GeV}^{3},$
$\displaystyle\|R^{{}^{\prime}}_{1P}(0)\|^{2}$	$\displaystyle=$	$\displaystyle 1.417~{}\textrm{GeV}^{5},~{}~{}~{}~{}~{}~{}\|R^{{}^{\prime}}_{2P}(0)\|^{2}=1.653~{}\textrm{GeV}^{5},$
$\displaystyle\|R^{{}^{\prime}}_{3P}(0)\|^{2}$	$\displaystyle=$	$\displaystyle 1.794~{}\textrm{GeV}^{5}.$	(17)

The branching ratios of diverse excited bottomonium states transitioning to lower energy states are documented in Refs. Br1 ; Br2 ; Br3 ; PDG .

Table 1: SDCs (in unit:

10^{-12}

\textrm{GeV}^{-5}

) corresponding to the process

Z\to b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}S_{1}^{1}]

with

m_{b}=4.7

GeV and

\mu_{r}=\frac{m_{Z}}{2}

\hat{\Gamma}^{\textrm{LO}}_{\textrm{total}}

incorporates the individual contributions of

\hat{\Gamma}_{1}^{(0)}

\hat{\Gamma}_{2}^{(0)}

, and

\hat{\Gamma}_{3}^{(0)}

, while

\hat{\Gamma}^{\textrm{NLO}}_{\textrm{total}}

refers to the combined sum of the contributions of

\hat{\Gamma}_{1,2,3}^{(0,1)}

$\hat{\Gamma}^{(0)}_{1}$	$\hat{\Gamma}^{(1)}_{1}$	$\hat{\Gamma}^{(0)}_{2}$	$\hat{\Gamma}^{(1)}_{2}$	$\hat{\Gamma}^{(0)}_{3}$	$\hat{\Gamma}^{(1)}_{3}$	$\hat{\Gamma}^{\textrm{LO}}_{\textrm{total}}$	$\hat{\Gamma}^{\textrm{NLO}}_{\textrm{total}}$
$165.2$	$244.9$	$229.4$	$134.4$	$79.67$	$-24.79$	$474.3$	$828.8$

Table 2: SDCs (in unit:

10^{-12}

\textrm{GeV}^{-7}

) corresponding to the process

Z\to b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{1}]

with

m_{b}=4.7

GeV and

\mu_{r}=\frac{m_{Z}}{2}

\hat{\Gamma}^{\textrm{LO}}_{\textrm{total}}

incorporates the individual contributions of

\hat{\Gamma}_{1}^{(0)}

\hat{\Gamma}_{2}^{(0)}

, and

\hat{\Gamma}_{3}^{(0)}

, while

\hat{\Gamma}^{\textrm{NLO}}_{\textrm{total}}

refers to the combined sum of the contributions of

\hat{\Gamma}_{1,2,3}^{(0,1)}

	$\hat{\Gamma}^{(0)}_{1}$	$\hat{\Gamma}^{(1)}_{1}$	$\hat{\Gamma}^{(0)}_{2}$	$\hat{\Gamma}^{(1)}_{2}$	$\hat{\Gamma}^{(0)}_{3}$	$\hat{\Gamma}^{(1)}_{3}$	$\hat{\Gamma}^{\textrm{LO}}_{\textrm{total}}$	$\hat{\Gamma}^{\textrm{NLO}}_{\textrm{total}}$
$J=0$	$92.12$	$81.42$	$8.631$	$3.358$	$0.323$	$-0.070$	$101.1$	$185.8$
$J=1$	$15.26$	$10.84$	$-0.373$	$2.282$	$1.916$	$-0.499$	$16.80$	$29.42$
$J=2$	$76.08$	$15.61$	$5.985$	$-0.700$	$0.652$	$-0.413$	$82.72$	$97.21$

Before proceeding further, let us first scrutinize the SDCs. Inspecting the data in Tables 1 and 2, we find

1)

In the production of double $b\bar{b}[^{3}S_{1}^{1}]$ states, considering the LO level in $\alpha_{s}$ , the inclusion of $\hat{\Gamma}_{2,3}^{(0)}$ enhances the pure QCD results, denoted as $\hat{\Gamma}_{1}^{(0)}$ , by approximately a factor of 2. When QCD corrections are incorporated, the QCD results are further augmented by about 1.5 times, as indicated by the ratio $\hat{\Gamma}_{1}^{(1)}/\hat{\Gamma}_{1}^{(0)}$ . However, high-order terms in $\alpha_{s}$ significantly reduce the QED results, with $(\hat{\Gamma}_{3}^{(1)}+\hat{\Gamma}_{3}^{(0)})/\hat{\Gamma}_{3}^{(0)}$ being approximately $70\%$ . Due to the interference between QCD and QED results, $\hat{\Gamma}_{2}^{(0)}$ experience a $40\%$ enhancement when QCD corrections are considered, as reflected in the ratio $\hat{\Gamma}_{2}^{(1)}/\hat{\Gamma}_{2}^{(0)}$ .
2)

Regarding the $b\bar{b}[^{3}S_{1}^{1}]$ yield associated with $b\bar{b}[^{3}P_{J}^{1}]$ , the LO SDCs $\hat{\Gamma}_{2,3}^{(0)}$ amplify $\hat{\Gamma}_{1}^{(0)}$ by approximately ten percents. With the inclusion of QCD corrections, $\hat{\Gamma}_{1}^{(1)}$ significantly enhances the $\hat{\Gamma}_{1}^{(0)}$ for the $b\bar{b}[^{3}P_{0,1}^{1}]$ process; however, the enhancement is milder for $b\bar{b}[^{3}P_{2}^{1}]$ . In terms of QED processes, $\hat{\Gamma}_{3}^{(1)}$ slightly diminishes the LO SDC $\hat{\Gamma}_{3}^{(0)}$ for $b\bar{b}[^{3}P_{0,1}^{1}]$ , while it substantially reduces that for $b\bar{b}[^{3}P_{2}^{1}]$ . Consequently, due to the combined influence of enhancement and reduction effects, the QCD corrections ( $\hat{\Gamma}_{2}^{(1)}$ ) have a substantial impact on $\hat{\Gamma}_{2}^{(0)}$ corresponding to $b\bar{b}[^{3}P_{0,1}^{1}]$ , whereas the effect is more moderate for $b\bar{b}[^{3}P_{2}^{1}]$ .

The ratio $\hat{\Gamma}^{\textrm{NLO}}_{\textrm{total}}/\hat{\Gamma}^{\textrm{LO}}_{\textrm{total}}$ in Tables 1 and 2 suggests that QCD corrections can substantially enhance the LO results, underscoring the importance of our newly-calculated higher order terms in $\alpha_{s}$ .

To illustrate the relative importance of the QCD and QED contributions, along with their interference effects, we plot $\hat{\Gamma}_{0,1,2}$ as a function of the renormalization scale in Fig. 3.

Table 3: Decay widths of

Z\to\Upsilon(mS)+\Upsilon(nS)

(in unit:

10^{-12}

GeV) with

m_{b}=4.7

GeV and

\mu_{r}=m_{Z}/2

. The subscripts “dir” and “fd” denote direct-production processes and feed-down effects, respectively. All predictions are composed of contributions from

\Gamma_{1,2,3}^{(0,1)}

	$\Gamma_{\textrm{dir}}$	$\Gamma_{\textrm{fd}}$	$\Gamma_{\textrm{total}}$	${\mathcal{B}}_{\textrm{theo}}$	$\mathcal{B}_{\textrm{exp}}$ CMS:2022fsq
$1S+1S$	$110.1$	$43.58$	$153.7$	$6.158\times 10^{-11}$	$<1.8\times 10^{-6}$
$1S+2S$	$109.9$	$29.16$	$139.1$	$5.575\times 10^{-11}$
$1S+3S$	$84.10$	$17.78$	$101.9$	$4.083\times 10^{-11}$
$2S+2S$	$27.45$	$7.329$	$34.78$	$1.394\times 10^{-11}$	$<3.9\times 10^{-7}$
$2S+3S$	$41.99$	$5.000$	$46.99$	$1.883\times 10^{-11}$
$3S+3S$	$16.06$	$0.838$	$16.90$	$6.773\times 10^{-12}$

In Table 3, we confront our predictions with the CMS measurements. The subscripts “dir” and “fd” signify direct-production processes and feed-down effects, respectively. The ratio $\Gamma_{\textrm{fd}}/\Gamma_{\textrm{dir}}$ reveals that the feed-down contributions from the transitions of excited bottomonium states are crucial for $Z$ boson decay into a pair of $\Upsilon$ mesons. For example, in the case of $\Upsilon(1S)+\Upsilon(1S)$ , feed-down contributions constitute approximately $30\%$ of the total results. Comparisons with experimental data indicate that the NLO predictions derived from the NRQCD framework are markedly below the upper threshold set by CMS.

Finally, we examine the uncertainties in our predictions resulting from variations in the $b$ -quark mass $m_{b}$ and the renormalization scale $\mu_{r}$ .

$\displaystyle\mathcal{B}_{Z\to\Upsilon(1S)+\Upsilon(1S)}$	$\displaystyle=$	$\displaystyle(6.158^{+0.254+2.868}_{-0.236-0.767})\times 10^{-11},$
$\displaystyle\mathcal{B}_{Z\to\Upsilon(1S)+\Upsilon(2S)}$	$\displaystyle=$	$\displaystyle(5.575^{+0.233+2.586}_{-0.216-0.695})\times 10^{-11},$
$\displaystyle\mathcal{B}_{Z\to\Upsilon(1S)+\Upsilon(3S)}$	$\displaystyle=$	$\displaystyle(4.083^{+0.168+1.905}_{-0.156-0.509})\times 10^{-11},$
$\displaystyle\mathcal{B}_{Z\to\Upsilon(2S)+\Upsilon(2S)}$	$\displaystyle=$	$\displaystyle(1.394^{+0.059+0.647}_{-0.054-0.174})\times 10^{-11},$
$\displaystyle\mathcal{B}_{Z\to\Upsilon(2S)+\Upsilon(3S)}$	$\displaystyle=$	$\displaystyle(1.883^{+0.078+0.877}_{-0.073-0.235})\times 10^{-11},$
$\displaystyle\mathcal{B}_{Z\to\Upsilon(3S)+\Upsilon(3S)}$	$\displaystyle=$	$\displaystyle(6.773^{+0.277+3.173}_{-0.258-0.844})\times 10^{-12},$	(18)

where the uncertainties listed in the first column stem from fluctuations in the value of $m_{b}$ between $4.6$ to $4.8$ GeV, centered around 4.7 GeV, whereas those in the second column result from adjustments to $\mu_{r}$ ranging from $2m_{b}$ to $m_{Z}$ around $m_{Z}/2$ . The above results suggest that a 0.1 GeV variation in $m_{b}$ around 4.7 GeV lead to a modest change in predictions; however, adjusting $\mu_{r}$ from $2m_{b}$ to $m_{Z}$ , particularly around $m_{Z}/2$ , significantly affects the predictions.

Note that, in $Z$ boson decay into double $\Upsilon$ , besides the color-singlet processes of interest, color-octet (CO) processes are also possible. The primary CO contributions are derived from the process $Z\to b\bar{b}[^{3}S_{1}^{8}]+g^{*}$ , followed by $g^{*}\to b\bar{b}[^{3}S_{1}^{8}]$ . Given the ratio of $\frac{\langle{\mathcal{O}}^{H}(^{3}S_{1}^{8})\rangle}{\langle{\mathcal{O}}^{H}(^{3}S_{1}^{1})\rangle}\sim 10^{-3}$ Br3 , the branching ratio $\mathcal{B}_{Z\to\Upsilon(1S)+\Upsilon(1S)}$ due to CO contributions are estimated to be $10^{-14}$ . Moreover, considering that experiments account for feed down from higher excited states, the semi-inclusive production of double $\Upsilon$ mesons in $Z$ boson decay, such as contributions from $b$ (or $\bar{b}$ ) quark fragmentation into $\Upsilon$ , should be considered. Utilizing the fragmentation function from Ref. z decay 4 , the branching ratio $\mathcal{B}_{Z\to\Upsilon(1S)+\Upsilon(1S)+X}$ is predicted to be $\sim 10^{-11}$ , which is comparable with the results in Eq. (18).

IV Summary

In this paper, we carry out a comprehensive study of $Z\to\Upsilon(mS)+\Upsilon(nS)$ utilizing the NRQCD framework at the QCD NLO accuracy. The LO results indicate that the contributions from QED diagrams significantly enhance those from the QCD processes. The QCD corrections substantially amplify the QCD results, while concurrently reducing the QED results. Additionally, we discover that the transitions involving higher excited states of bottomonium mesons exert an indispensable influence. Taking into account both the QCD and QED contributions, we estimate the branching ratio for $\mathcal{B}_{Z\to\Upsilon(mS)+\Upsilon(nS)}$ to be on the order of $10^{-11}$ . This value is significantly lower than the upper limit established by CMS and aligns with experimental observations.

Appendix A

In this section, we present the expressions for the NLO coefficients $a_{i}$ , $b_{i}$ , and $c_{i}$ in Eq. (16) which are formulated as a superposition of the MIs $\mathcal{I}$ . The coefficients related to the $Z\to b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}S_{1}^{1}]$ process can be obtained by substituting $m_{c}$ in Eqs. (A1)-(A3) of Ref. z decay 40 with $m_{b}$ . In the following, we only provide the expressions of $a_{i}$ , $b_{i}$ , and $c_{i}$ pertaining to the $Z\to b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{1}]$ process.

A.1 Master Integrals

To start, we introduce the following definitions and showcase the finite ( $\epsilon^{0}$ -order) terms of the MIs that are involved ( $r\equiv\frac{m^{2}_{Z}}{4m_{b}^{2}}$ ):

	$\displaystyle a=\sqrt{r},~{}b=\sqrt{r-1},~{}c=\sqrt{r-4},~{}d=2r+1,$
	$\displaystyle f=r+ac-4,~{}g=r-ac-4,~{}h=rg+2ac,$
	$\displaystyle j=rf-2ac,~{}j_{1}=(r-4)ab,~{}j_{2}=(2-r)bc.$		(19)

There is only one 1-point scalar integral,

\displaystyle\mathcal{I}_{1}

\displaystyle=

\displaystyle\frac{\lambda}{\mu_{r}^{4-D}}\int{\frac{d^{D}k}{k^{2}-m_{b}^{2}}}=m_{b}^{2}\left[1-2\ln(m_{b})\right],

(20)

where $k$ denotes the loop momentum, $\lambda=\frac{\mu_{r}^{4-D}}{i\pi^{\frac{D}{2}}\gamma_{\Gamma}}$ with $\gamma_{\Gamma}=\frac{\Gamma^{2}(1-\epsilon)\Gamma(1+\epsilon)}{\Gamma(1-2\epsilon)}$ .
There are five 2-point scalar integrals,

\displaystyle\mathcal{I}^{(1)}_{2}

\displaystyle=

\displaystyle\frac{\lambda}{\mu_{r}^{4-D}}\int{\frac{d^{D}k}{k^{2}[(k+\frac{2p_{1}+p_{2}}{2})^{2}-m_{b}^{2}]}}=2-2\ln{(m_{b})}-\frac{2r}{d}\left[\ln{(2r)}-i\pi\right],

(21)

\displaystyle\mathcal{I}^{(2)}_{2}

\displaystyle=

\displaystyle\frac{\lambda}{\mu_{r}^{4-D}}\int{\frac{d^{D}k}{k^{2}(k-\frac{p_{1}+p_{2}}{2})^{2}}}=2-\ln{({m_{b}}^{2}r)}+i\pi,

(22)

\displaystyle\mathcal{I}^{(3)}_{2}

\displaystyle=

\displaystyle\frac{\lambda}{\mu_{r}^{4-D}}\int\frac{d^{D}k}{(k^{2}-m_{b}^{2})[(k+\frac{p_{1}+p_{2}}{2})^{2}-m_{b}^{2}]}=\frac{c}{a}\left[\ln{\left(\frac{4r}{(f+4)^{2}}\right)}+i\pi\right]+2\left[1-\ln{(m_{b})}\right],

\displaystyle\mathcal{I}^{(4)}_{2}

\displaystyle=

\displaystyle\frac{\lambda}{\mu_{r}^{4-D}}\int\frac{d^{D}k}{(k^{2}-m_{b}^{2})[(k+p_{1}+p_{2})^{2}-m_{b}^{2}]}=\frac{b}{a}\left[\ln{(-2ab+d-2)}+i\pi\right]+2\left[1-\ln{(m_{b})}\right],

\displaystyle\mathcal{I}^{(5)}_{2}

\displaystyle=

\displaystyle\frac{\lambda}{\mu^{4-D}}\int\frac{d^{D}k}{k^{2}(k-p_{1})^{2}}=-\ln{\left(4m_{b}^{2}\right)}+i\pi+2,

(25)

There are seven 3-point scalar integrals,

	$\displaystyle\mathcal{I}^{(1)}_{3}$
	$\displaystyle=\frac{\lambda}{\mu_{r}^{4-D}}\int\frac{d^{D}k}{k^{2}[(k+\frac{p_{2}}{2})^{2}-m_{b}^{2}][(k+\frac{2p_{1}+p_{2}}{2})^{2}-m_{b}^{2}]}$
	$\displaystyle=\frac{1}{2acm_{b}^{2}}\left\{{\ln(2)\ln\left(\frac{(3(g+4)+h)^{2}}{d(g+h+4)^{2}}\right)+\left[\ln(r)-i\pi\right]\ln\left(\frac{d}{(g+3)^{2}}\right)+4\textrm{Li}_{2}\left(\frac{c}{a}\right)+\textrm{Li}_{2}\left(-\frac{2(ac+j)}{d}\right)}\right.$
	$\displaystyle\left.{-\textrm{Li}_{2}\left(\frac{c^{2}}{ac-h}\right)-\textrm{Li}_{2}\left(-\frac{c+ag}{a}\right)-2\textrm{Li}_{2}\left(-\frac{c}{a}\right)}\right\},$		(26)

	$\displaystyle\mathcal{I}^{(2)}_{3}$
	$\displaystyle=\frac{\lambda}{\mu_{r}^{4-D}}\int\frac{d^{D}k}{k^{2}[(k-\frac{p_{2}}{2})^{2}-m_{b}^{2}](k+\frac{p_{1}+p_{2}}{2})^{2}}$
	$\displaystyle=\frac{1}{acm_{b}^{2}}\left\{{\ln(2)\ln\left(\frac{d(f+4)^{2}}{(3r-ac)^{2}}\right)+\left[\ln(r)-i\pi)\right]\ln\left(\frac{4dr}{(3r-ac)^{2}}\right)+\textrm{Li}_{2}\left(-\frac{2(ac+j)}{d}\right)+2\textrm{Li}_{2}\left(\frac{c}{a}\right)}\right.$
	$\displaystyle\left.{+\textrm{Li}_{2}\left(-\frac{g}{2}\right)-\textrm{Li}_{2}\left(-\frac{f}{2}\right)-\textrm{Li}_{2}\left(\frac{c^{2}}{ac-h}\right)-\textrm{Li}_{2}\left(-\frac{c+ag}{a}\right)}\right\},$		(27)

	$\displaystyle\mathcal{I}^{(3)}_{3}$	$\displaystyle=$	$\displaystyle\frac{\lambda}{\mu_{r}^{4-D}}\int\frac{d^{D}k}{k^{2}[(k+\frac{p_{1}}{2})^{2}-m_{b}^{2}](k+\frac{p_{1}+p_{2}}{2})^{2}}$
		$\displaystyle=$	$\displaystyle\frac{1}{acm_{b}^{2}}\left\{\left[\ln(r)-i\pi\right]\ln\left(\frac{(f+4)^{2}}{4r}\right)+2\textrm{Li}_{2}\left(\frac{c}{a}\right)-2\textrm{Li}_{2}\left(-\frac{c}{a}\right)-\textrm{Li}_{2}\left(-\frac{g}{2}\right)+\textrm{Li}_{2}\left(-\frac{f}{2}\right)\right\},$

	$\displaystyle\mathcal{I}^{(4)}_{3}$
	$\displaystyle=\frac{\lambda}{\mu_{r}^{4-D}}\int\frac{d^{D}k}{k^{2}[(k+\frac{2p_{1}+p_{2}}{2})^{2}-m_{b}^{2}](k+\frac{p_{1}+p_{2}}{2})^{2}}$
	$\displaystyle=\frac{1}{acm_{b}^{2}}\left\{{\ln(2)\ln\left(\frac{d}{(g+3)^{2}}\right)+\left[\ln(r)-i\pi)\right]\ln\left(\frac{g+6}{2g+6}\right)+\textrm{Li}_{2}\left(-\frac{g}{2}\right)-\textrm{Li}_{2}\left(-\frac{f}{2}\right)-2\textrm{Li}_{2}\left(-\frac{c}{a}\right)}\right.$
	$\displaystyle\left.{-\textrm{Li}_{2}\left(\frac{2ac-2h}{d}\right)+\textrm{Li}_{2}\left(\frac{h-ac}{dr}\right)+\textrm{Li}_{2}\left(\frac{c}{a}-f\right)}\right\},$		(29)

			$\displaystyle\mathcal{I}^{(5)}_{3}$
			$\displaystyle=\frac{\lambda}{\mu_{r}^{4-D}}\int\frac{d^{D}k}{k^{2}[(k-\frac{p_{1}}{2})^{2}-m_{b}^{2}][(k+\frac{p_{1}+2p_{2}}{2})^{2}-m_{b}^{2}]}$
			$\displaystyle=\frac{1}{2acm_{b}^{2}}\left\{{\left[\ln\left(\frac{(r-ab)^{2}}{r}\right)+i\pi\right]\ln\left(\frac{br(2b-3c)+j-j_{1}}{br(2b+3c)+j+j_{1}}\right)+\left[\ln(2r)-i\pi\right]\ln\left(\frac{(g+3)^{2}}{d}\right)}\right.$
			$\displaystyle\left.{{-\textrm{Li}_{2}\left(-\frac{2(ac+j)}{d}\right)-\textrm{Li}_{2}\left(\frac{c^{2}-h+j_{1}+j_{2}}{r}\right)-\textrm{Li}_{2}\left(-\frac{-c^{2}+h+j_{1}+j_{2}}{r}\right)+\textrm{Li}_{2}\left(\frac{c^{2}-j+j_{1}-j_{2}}{r}\right)}}\right.$
			$\displaystyle\left.{-2\textrm{Li}_{2}\left(-\frac{c}{a}\right)+\textrm{Li}_{2}\left(\frac{c^{2}-j-j_{1}+j_{2}}{r}\right)+{\textrm{Li}_{2}\left(\frac{c^{2}}{ac-h}\right)+\textrm{Li}_{2}\left(-\frac{c+ag}{a}\right)}}\right\},$
					(30)

	$\displaystyle\mathcal{I}^{(6)}_{3}$
	$\displaystyle=\frac{\lambda}{\mu_{r}^{4-D}}\int\frac{d^{D}k}{k^{2}[(k+\frac{p_{2}}{2})^{2}-m_{b}^{2}][(k+\frac{p_{1}+2p_{2}}{2})^{2}-m_{b}^{2}]}$
	$\displaystyle=\frac{1}{acm_{b}^{2}}\left\{{\ln(2)\ln\left(\frac{(f+4)^{4}\left(3(2r-ac)+h\right)}{4(g+6)\left(4(r-ac)+h\right)}\right)+\ln(r)\ln\left(\frac{dr}{(d+f+3)^{2}}\right)-i\pi\ln\left(\frac{16d}{(g+6)^{2}r}\right)}\right.$
	$\displaystyle\left.{+\textrm{Li}_{2}\left(-\frac{g}{8}\right)+\textrm{Li}_{2}\left(-\frac{4ac+j}{2d}\right)-\textrm{Li}_{2}\left(\frac{4ac-h}{4r}\right)+2\textrm{Li}_{2}\left(\frac{c}{2a}\right)+\textrm{Li}_{2}\left(\frac{f}{2r}\right)-\textrm{Li}_{2}\left(-\frac{f}{8}\right)-\textrm{Li}_{2}\left(\frac{g}{2r}\right)}\right.$
	$\displaystyle\left.{-\textrm{Li}_{2}\left(\frac{c^{2}}{4ac-h}\right)}\right\},$		(31)

	$\displaystyle\mathcal{I}^{(7)}_{3}=\frac{\lambda}{\mu^{4-D}}\int\frac{d^{D}k}{k^{2}[(k+\frac{p_{1}}{2})^{2}-m_{b}^{2}](k-p2)^{2}}$
	$\displaystyle=\frac{1}{2acm_{b}^{2}}\left\{{\ln(2)\ln\left(\frac{64d}{(f+4)^{4}(g+6)^{2}}\right)+\ln(r)\ln\left(\frac{d(f+4)^{4}}{[(g+6)r]^{2}}\right)+i\pi\ln\left(\frac{(g+6)^{2}}{4d}\right)+\textrm{Li}_{2}\left(-\frac{h}{2}\right)}\right.$
	$\displaystyle\left.{+2\textrm{Li}_{2}\left(\frac{ac}{2}\right)-\textrm{Li}_{2}\left(\frac{jr-2ac}{4d}\right)+\textrm{Li}_{2}\left(\frac{2acr-jr^{2}}{2d}\right)-\textrm{Li}_{2}\left(\frac{2ac-fr}{2}\right)-\textrm{Li}_{2}\left(-\frac{2ac+hr}{4}\right)}\right\}.$
			(32)

A.2 NLO coefficients

A.2.1 $b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{0}]$

\displaystyle a_{1}

\displaystyle=

\displaystyle-\frac{1}{3}{\mathcal{I}^{(2)}_{2}}-\frac{2(r^{4}-25r^{3}-146r^{2}+176r+72)}{9(r^{4}+182r^{3}-428r^{2}+152r+144)}-\frac{2}{3}\ln(2m_{b}),

(33)

	$\displaystyle b_{1}$	$\displaystyle=$	$\displaystyle\frac{2(r^{4}+176r^{3}-716r^{2}+440r+288)}{3m_{b}^{2}r(r^{4}+182r^{3}-428r^{2}+152r+144)}{\mathcal{I}_{1}}-\frac{r^{5}+184r^{4}-76r^{3}-1280r^{2}+1024r+576}{3r(r^{4}+182r^{3}-428r^{2}+152r+144)}{\mathcal{I}^{(3)}_{2}}$		(34)
			$\displaystyle-\frac{2(r^{5}-13r^{4}+814r^{3}-2404r^{2}+1176r+864)}{9r(r^{4}+182r^{3}-428r^{2}+152r+144)}-\frac{2\ln(2m_{b})}{3},$		(34)

$\displaystyle c_{1}$	$\displaystyle=$	$\displaystyle\frac{120r^{7}+19848r^{6}-103058r^{5}+54543r^{4}+96834r^{3}-16648r^{2}-32352r-5760}{9m_{b}^{2}r(r-4)(2r+1)(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}_{1}$
		$\displaystyle+\frac{5796r^{6}-14278r^{5}-21720r^{4}+16167r^{3}+32542r^{2}+14136r+2016}{9r(2r+1)(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(1)}_{2}$
		$\displaystyle+\frac{-1407r^{4}+16354r^{3}-25036r^{2}+2136r+6264}{9(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(2)}_{2}$
		$\displaystyle+\frac{7r^{5}+1408r^{4}-8952r^{3}+27008r^{2}-20608r-4992}{18(r-4)(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(3)}_{2}$
		$\displaystyle-\frac{1469r^{4}+2814r^{3}-16728r^{2}+11440r+4704}{9(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(4)}_{2}$
		$\displaystyle+\frac{2m_{b}^{2}r(r^{4}+60r^{3}-684r^{2}+640r+240)}{3(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(1)}_{3}$
		$\displaystyle+\frac{m_{b}^{2}(-21r^{5}+208r^{4}-2137r^{3}+2010r^{2}+2444r-624)}{12(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(2)}_{3}$
		$\displaystyle-\frac{m_{b}^{2}(r^{5}-658r^{4}+708r^{3}+924r^{2}-816r-288)}{3(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(3)}_{3}$
		$\displaystyle+\frac{m_{b}^{2}(r^{5}-480r^{4}-3239r^{3}+6838r^{2}-1868r-528)}{12(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(4)}_{3}$
		$\displaystyle-\frac{2m_{b}^{2}(4r^{6}+1208r^{5}-6399r^{4}+10840r^{3}-6148r^{2}-720r+1344)}{3(r-2)(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(5)}_{3}$
		$\displaystyle+\frac{2m_{b}^{2}(3r^{5}+96r^{4}-406r^{3}+860r^{2}-528r-288)}{3(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(6)}_{3}$
		$\displaystyle-\frac{24m_{b}^{2}r^{2}(5r^{2}-2r-6)}{(r-2)(r^{4}+182r^{3}-428r^{2}+152r+144)}\mathcal{I}^{(7)}_{3}$
		$\displaystyle+\frac{(19r^{5}+3512r^{4}-7204r^{3}+200r^{2}+3968r+768)}{r^{5}+182r^{4}-428r^{3}+152r^{2}+144r}\ln(m_{b})$
		$\displaystyle+\frac{-25r^{6}-12935r^{5}+80893r^{4}-118558r^{3}-26064r^{2}+85824r+21888}{9(r-4)r(r^{4}+182r^{3}-428r^{2}+152r+144)}+11\ln(2),$

\displaystyle a_{2}

\displaystyle=

\displaystyle-\frac{1}{6}{\mathcal{I}^{(2)}_{2}}+\frac{7r^{4}+38r^{3}+16r^{2}-424r-144}{72(7r^{4}-7r^{3}-59r^{2}+56r+36)}-\frac{\ln(2m_{b})}{3},

(36)

	$\displaystyle b_{2}$	$\displaystyle=$	$\displaystyle\frac{14r^{4}-35r^{3}-136r^{2}+292r+144}{6m_{b}^{2}r(7r^{4}-7r^{3}-59r^{2}+56r+36)}{\mathcal{I}_{1}}-\frac{7r^{5}+7r^{4}-94r^{3}-80r^{2}+328r+144}{6r(7r^{4}-7r^{3}-59r^{2}+56r+36)}{\mathcal{I}^{(3)}_{2}}$		(37)
			$\displaystyle+\frac{7r^{5}-256r^{4}+328r^{3}+2288r^{2}-3216r-1728}{72r(7r^{4}-7r^{3}-59r^{2}+56r+36)}-\frac{\ln(2m_{b})}{3},$		(37)

	$\displaystyle c_{2}=$
	$\displaystyle\frac{21672r^{8}-69156r^{7}-219454r^{6}+398155r^{5}+899784r^{4}-46124r^{3}-616672r^{2}-309120r-43776}{144(r-4)r(2rm_{b}+m_{b})^{2}(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}_{1}$
	$\displaystyle+\frac{3840r^{8}+34740r^{7}-143528r^{6}-28915r^{5}+327531r^{4}+350456r^{3}+158500r^{2}+33648r+2880}{144r(2r+1)^{2}(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(1)}_{2}$
	$\displaystyle+\frac{-240r^{5}+714r^{4}+13853r^{3}-44450r^{2}+13812r+13176}{72(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(2)}_{2}$
	$\displaystyle+\frac{56r^{5}-1915r^{4}+13254r^{3}-30968r^{2}+15472r+8736}{36(r-4)(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(3)}_{2}$
	$\displaystyle-\frac{120r^{5}+1991r^{4}-4392r^{3}-1320r^{2}+6016r+2544}{36(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(4)}_{2}$
	$\displaystyle+\frac{r(2r^{4}-207r^{3}+213r^{2}+242r+12)m_{b}^{2}}{3(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(1)}_{3}$
	$\displaystyle+\frac{(21r^{5}-7333r^{4}+22129r^{3}-13620r^{2}-140r+4368)m_{b}^{2}}{96(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(2)}_{3}$
	$\displaystyle+\frac{(100r^{5}+965r^{4}-5556r^{3}+7692r^{2}-1680r-2016)m_{b}^{2}}{24(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(3)}_{3}$
	$\displaystyle-\frac{(13r^{5}-3993r^{4}+1633r^{3}+12124r^{2}-11084r-3696)m_{b}^{2}}{96(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(4)}_{3}$
	$\displaystyle-\frac{(628r^{6}-1309r^{5}-4890r^{4}+14188r^{3}-5944r^{2}-4704r+1920)m_{b}^{2}}{24(r-2)(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(5)}_{3}$
	$\displaystyle+\frac{(15r^{5}-171r^{4}+2044r^{3}-6236r^{2}+3552r+2016)m_{b}^{2}}{12(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(6)}_{3}$
	$\displaystyle-\frac{3r^{2}(3r^{3}+2r^{2}-20r-24)m_{b}^{2}}{4(r-2)(7r^{4}-7r^{3}-59r^{2}+56r+36)}\mathcal{I}^{(7)}_{3}$
	$\displaystyle+\frac{759r^{5}-708r^{4}-6736r^{3}+3312r^{2}+7312r+1536}{8r(7r^{4}-7r^{3}-59r^{2}+56r+36)}\ln(m_{b})$
	$\displaystyle+\frac{-16562r^{7}+78949r^{6}+154533r^{5}-887080r^{4}-10228r^{3}+1040208r^{2}+535104r+66816}{144(r-4)r(2r+1)(7r^{4}-7r^{3}-59r^{2}+56r+36)}$
	$\displaystyle+\frac{11}{2}\ln(2),$		(38)

\displaystyle a_{3}

\displaystyle=

\displaystyle b_{3}=0,

(39)

	$\displaystyle c_{3}$
	$\displaystyle=\frac{\mathcal{I}_{1}}{9(r-4)r(2rm_{b}+m_{b})^{2}(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}$
	$\displaystyle\times(1512r^{9}+3216r^{8}-31716r^{7}+29233r^{6}+208412r^{5}+192816r^{4}-374944r^{3}$
	$\displaystyle-536912r^{2}-231744r-32256)$
	$\displaystyle+\frac{4(1446r^{8}-1323r^{7}-10786r^{6}+27874r^{5}+53997r^{4}+26482r^{3}+2900r^{2}-1320r-288)}{9r(2r+1)^{2}(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(1)}_{2}$
	$\displaystyle-\frac{8(60r^{5}-429r^{4}-70r^{3}+3868r^{2}-3000r-1728)}{9(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(2)}_{2}$
	$\displaystyle-\frac{64(17r^{5}-103r^{4}+51r^{3}+850r^{2}-1640r-624)}{9(r-4)(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(3)}_{2}$
	$\displaystyle-\frac{4(255r^{5}-8r^{4}-1272r^{3}+1392r^{2}-784r+192)}{9(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(4)}_{2}$
	$\displaystyle-\frac{8r(55r^{4}-24r^{3}-480r^{2}-56r+96)m_{b}^{2}}{3(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(1)}_{3}$
	$\displaystyle-\frac{(273r^{5}-484r^{4}-2138r^{3}+3408r^{2}-56r-1248)m_{b}^{2}}{3(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(2)}_{3}$
	$\displaystyle+\frac{2(155r^{5}-530r^{4}-528r^{3}+4728r^{2}-3792r-2304)m_{b}^{2}}{3(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(3)}_{3}$
	$\displaystyle+\frac{4(86r^{5}+117r^{4}-805r^{3}-892r^{2}+1916r+528)m_{b}^{2}}{3(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(4)}_{3}$
	$\displaystyle-\frac{4(18r^{6}+85r^{5}-32r^{4}-772r^{3}+272r^{2}+1344r+192)m_{b}^{2}}{24(r-2)(3(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(5)}_{3}$
	$\displaystyle-\frac{32(3r^{5}-42r^{4}+20r^{3}+320r^{2}-672r-288)m_{b}^{2}}{3(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\mathcal{I}^{(6)}_{3}$
	$\displaystyle+\frac{2(36r^{6}-7r^{5}-668r^{4}-168r^{3}+2608r^{2}+3344r+768)}{r(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}\ln(m_{b})$
	$\displaystyle+\frac{-1404r^{8}+5174r^{7}+34583r^{6}-121380r^{5}-207200r^{4}+497824r^{3}+727440r^{2}+265536r+23040}{9(r-4)r(2r+1)(9r^{5}+2r^{4}-152r^{3}-16r^{2}+592r+288)}.$

A.2.2 $b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{1}^{1}]$

\displaystyle a_{1}

\displaystyle=

\displaystyle-\frac{1}{3}{\mathcal{I}^{(2)}_{2}}-\frac{2r^{3}-15r^{2}+16r+36}{18(2r^{3}-12r^{2}+13r+18)}-\frac{2\ln(2m_{b})}{3},

(41)

	$\displaystyle b_{1}$	$\displaystyle=$	$\displaystyle\frac{4(r^{3}-9r^{2}+14r+18)}{3m_{b}^{2}r(2r^{3}-12r^{2}+13r+18)}{\mathcal{I}_{1}}-\frac{2r^{4}-8r^{3}-23r^{2}+74r+72}{3r(2r^{3}-12r^{2}+13r+18)}{\mathcal{I}^{(3)}_{2}}$		(42)
			$\displaystyle-\frac{2r^{4}+45r^{3}-290r^{2}+264r+432}{18r(2r^{3}-12r^{2}+13r+18)}-\frac{2\ln(2m_{b})}{3},$		(42)

$\displaystyle c_{1}$	$\displaystyle=$	$\displaystyle-\frac{(-784r^{7}+4412r^{6}-3080r^{5}-9905r^{4}+5302r^{3}+13190r^{2}+6456r+960)}{12(r-4)r(2r^{3}-12r^{2}+13r+18)(2m_{b}r+m_{b})^{2}}\mathcal{I}_{1}$
		$\displaystyle+\frac{64r^{7}-2592r^{6}+396r^{5}+12116r^{4}+16867r^{3}+10206r^{2}+2930r+336}{12r(2r+1)^{2}(2r^{3}-12r^{2}+13r+18)}\mathcal{I}^{(1)}_{2}$
		$\displaystyle+\frac{-4r^{4}+323r^{3}-877r^{2}-134r+522}{12r^{3}-72r^{2}+78r+108}\mathcal{I}^{(2)}_{2}$
		$\displaystyle+\frac{19r^{4}-291r^{3}+1227r^{2}-1180r-624}{18(r-4)(2r^{3}-12r^{2}+13r+18)}\mathcal{I}^{(3)}_{2}$
		$\displaystyle+\frac{-12r^{4}+74r^{3}+473r^{2}-2020r-1176}{18(2r^{3}-12r^{2}+13r+18)}\mathcal{I}^{(4)}_{2}$
		$\displaystyle+\frac{2m_{b}^{2}r(4r^{3}-23r^{2}+41r+30)}{6r^{3}-36r^{2}+39r+54}\mathcal{I}^{(1)}_{3}$
		$\displaystyle-\frac{m_{b}^{2}(280r^{4}+450r^{3}-449r^{2}-1040r+312)}{48(2r^{3}-12r^{2}+13r+18)}\mathcal{I}^{(2)}_{3}$
		$\displaystyle+\frac{m_{b}^{2}(34r^{4}-135r^{3}-129r^{2}+276r+144)}{12(2r^{3}-12r^{2}+13r+18)}\mathcal{I}^{(3)}_{3}$
		$\displaystyle-\frac{m_{b}^{2}(408r^{4}-1010r^{3}-2523r^{2}+512r+264)}{48(2r^{3}-12r^{2}+13r+18)}\mathcal{I}^{(4)}_{3}$
		$\displaystyle-\frac{m_{b}^{2}(22r^{4}-247r^{3}+635r^{2}+146r-336)}{6(2r^{3}-12r^{2}+13r+18)}\mathcal{I}^{(5)}_{3}$
		$\displaystyle-\frac{m_{b}^{2}(10r^{4}-25r^{3}-162r^{2}+168r+144)}{6(2r^{3}-12r^{2}+13r+18)}\mathcal{I}^{(6)}_{3}$
		$\displaystyle+\frac{40r^{4}-223r^{3}+171r^{2}+444r+96}{2r^{4}-12r^{3}+13r^{2}+18r}\ln(m_{b})$
		$\displaystyle+\frac{-598r^{6}+6783r^{5}-19185r^{4}-5999r^{3}+40917r^{2}+29940r+5472}{18(r-4)r(2r+1)(2r^{3}-12r^{2}+13r+18)}+11\ln(2),$

\displaystyle a_{2}

\displaystyle=

\displaystyle-\frac{1}{6}{\mathcal{I}^{(2)}_{2}}+\frac{-5r^{3}+24r^{2}+68r+72}{18(19r^{3}-24r^{2}-88r-72)}-\frac{\ln(2m_{b})}{3},

(44)

	$\displaystyle b_{2}$	$\displaystyle=$	$\displaystyle\frac{r(r(19r-72)-184)-144}{3m_{b}^{2}r(r(r(19r-24)-88)-72)}{\mathcal{I}_{1}}-\frac{19r^{4}+14r^{3}-232r^{2}-440r-288}{6r(19r^{3}-24r^{2}-88r-72)}{\mathcal{I}^{(3)}_{2}}$		(45)
			$\displaystyle+\frac{-5r^{4}-234r^{3}+212r^{2}+960r+864}{18r(19r^{3}-24r^{2}-88r-72)}-\frac{\ln(2m_{b})}{3},$		(45)

	$\displaystyle c_{2}$
	$\displaystyle=\frac{116r^{8}+2088r^{7}-12033r^{6}-6098r^{5}+10386r^{4}+72762r^{3}+94260r^{2}+45360r+7296}{12(r-4)r(19r^{3}-24r^{2}-88r-72)(2m_{b}r+m_{b})^{2}}\mathcal{I}_{1}$
	$\displaystyle-\frac{128r^{8}+332r^{7}-3272r^{6}+31989r^{5}+76948r^{4}+67272r^{3}+28738r^{2}+5900r+480}{12r(2r+1)^{2}(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(1)}_{2}$
	$\displaystyle+\frac{-16r^{5}+113r^{4}-838r^{3}-5202r^{2}+1820r+4392}{12(-19r^{3}+24r^{2}+88r+72)}\mathcal{I}^{(2)}_{2}$
	$\displaystyle-\frac{9r^{5}-607r^{4}+4602r^{3}-17004r^{2}+15664r+17472}{36(r-4)(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(3)}_{2}$
	$\displaystyle+\frac{24r^{5}+207r^{4}-1346r^{3}+6268r^{2}+6448r+2544}{18(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(4)}_{2}$
	$\displaystyle-\frac{2m_{b}^{2}r(3r^{4}-62r^{3}+121r^{2}+131r+12)}{3(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(1)}_{3}$
	$\displaystyle+\frac{m_{b}(492r^{5}+2098r^{4}-4215r^{3}+2314r^{2}+184r-4368)}{48(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(2)}_{3}$
	$\displaystyle-\frac{m_{b}(18r^{5}-119r^{4}-1203r^{3}+2922r^{2}-480r-2016)}{12(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(3)}_{3}$
	$\displaystyle-\frac{636r^{5}+1578r^{4}+6869r^{3}-2382r^{2}+5176r+3696}{48(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(4)}_{3}$
	$\displaystyle+\frac{m_{b}(54r^{5}-458r^{4}+1745r^{3}+728r^{2}+788r-480)}{6(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(5)}_{3}$
	$\displaystyle-\frac{m_{b}(6r^{5}+41r^{4}+460r^{3}-1764r^{2}+2208r+2016)}{6(19r^{3}-24r^{2}-88r-72)}\mathcal{I}^{(6)}_{3}$
	$\displaystyle+\frac{-539r^{4}+442r^{3}+2656r^{2}+3240r+768}{-38r^{4}+48r^{3}+176r^{2}+144r}\ln(m_{b})$
	$\displaystyle-\frac{-468r^{7}+20074r^{6}-70839r^{5}-100506r^{4}+248630r^{3}+519516r^{2}+256080r+33408}{36(r-4)r(2r+1)(19r^{3}-24r^{2}-88r-72)}+\frac{11}{2}\ln(2),$
			(46)

\displaystyle a_{3}

\displaystyle=

\displaystyle b_{3}=0,

(47)

	$\displaystyle c_{3}$
	$\displaystyle=-\frac{648r^{9}-1628r^{8}-200r^{7}-5270r^{6}-29735r^{5}-58006r^{4}-89032r^{3}-79760r^{2}-34176r-5376}{3(r-4)r(2rm_{b}+m_{b})^{2}(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}_{1}$
	$\displaystyle+\frac{2(28r^{8}-1626r^{7}+1334r^{6}+14055r^{5}+19662r^{4}+11099r^{3}+2300r^{2}-196r-96)}{3r(2r+1)^{2}(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(1)}_{2}$
	$\displaystyle-\frac{4(4r^{5}-93r^{4}+276r^{3}+140r^{2}-616r-576)}{3(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(2)}_{2}$
	$\displaystyle+\frac{8(-9r^{5}+269r^{4}-1230r^{3}+1068r^{2}+3952r+2496)}{9(r-4)(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(3)}_{2}$
	$\displaystyle-\frac{4(39r^{5}-219r^{4}+1102r^{3}+1960r^{2}+664r+96)}{9(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(4)}_{2}$
	$\displaystyle-\frac{4r(57r^{4}-49r^{3}-127r^{2}+10r+96)m_{b}^{2}}{3(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(1)}_{3}$
	$\displaystyle+\frac{2(6r^{5}+71r^{4}-174r^{3}-226r^{2}+68r+624)m_{b}^{2}}{3(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(2)}_{3}$
	$\displaystyle+\frac{2(36r^{5}-161r^{4}+546r^{3}+294r^{2}-1164r-1152)m_{b}^{2}}{3(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(3)}_{3}$
	$\displaystyle-\frac{4(21r^{5}-72r^{4}+80r^{3}+12r^{2}-536r-264)m_{b}^{2}}{3(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(4)}_{3}$
	$\displaystyle-\frac{4(18r^{6}-45r^{5}+271r^{4}+476r^{3}+188r^{2}+464r+96)m_{b}^{2}}{3(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(5)}_{3}$
	$\displaystyle-\frac{16(6r^{5}-11r^{4}+50r^{3}-144r^{2}-480r-288)m_{b}^{2}}{3(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}\mathcal{I}^{(6)}_{3}$
	$\displaystyle+\frac{72r^{6}-72r^{5}+222r^{4}+2200r^{3}+3440r^{2}+2928r+768}{9r^{6}-9r^{5}+34r^{4}+228r^{3}+248r^{2}+144r}\ln(m_{b})$
	$\displaystyle+\frac{-1944r^{8}+10224r^{7}-9718r^{6}-62979r^{5}+122574r^{4}+343096r^{3}+350880r^{2}+131136r+11520}{9(r-4)r(2r+1)(9r^{5}-9r^{4}+34r^{3}+228r^{2}+248r+144)}.$
			(48)

A.2.3 $b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{2}^{1}]$

\displaystyle a_{1}

\displaystyle=

\displaystyle-\frac{1}{3}{\mathcal{I}^{(2)}_{2}}-\frac{2(r^{4}+11r^{3}-131r^{2}+224r+180)}{9(r^{4}+20r^{3}-188r^{2}+308r+360)}-\frac{2\ln(2m_{b})}{3},

(49)

	$\displaystyle b_{1}$	$\displaystyle=$	$\displaystyle\frac{2(r^{4}+14r^{3}-296r^{2}+776r+720)}{3m_{b}^{2}r(r^{4}+20r^{3}-188r^{2}+308r+360)}{\mathcal{I}_{1}}-\frac{r^{5}+22r^{4}-160r^{3}-284r^{2}+1912r+1440}{3r(r^{4}+20r^{3}-188r^{2}+308r+360)}{\mathcal{I}^{(3)}_{2}}$		(50)
			$\displaystyle-\frac{2(r^{5}+23r^{4}+73r^{3}-1366r^{2}+1968r+2160)}{9r(r^{4}+20r^{3}-188r^{2}+308r+360)}-\frac{2\ln(2m_{b})}{3},$		(50)

	$\displaystyle c_{1}$
	$\displaystyle=\frac{204r^{8}+1812r^{7}-57877r^{6}+146122r^{5}+217023r^{4}-115334r^{3}-239518r^{2}-104568r-14400}{9(r-4)r(2rm_{b}+m_{b})^{2}(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}_{1}$
	$\displaystyle+\frac{1248r^{7}+6628r^{6}+127916r^{5}+289545r^{4}+290345r^{3}+156466r^{2}+43890r+5040}{9r(2r+1)^{2}(r^{4}+20r^{3}-188r^{2}+308r+360}\mathcal{I}^{(1)}_{2}$
	$\displaystyle+\frac{-123r^{4}+868r^{3}-28744r^{2}-2868r+15660}{9(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(2)}_{2}$
	$\displaystyle-\frac{101r^{5}+146r^{4}-2496r^{3}-260r^{2}+1936r+12480}{18(r-4)(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(3)}_{2}$
	$\displaystyle-\frac{113r^{4}+798r^{3}+5286r^{2}+17944r+11760}{9(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(4)}_{2}$
	$\displaystyle+\frac{r(r^{3}+14r^{2}-78r-72)}{r^{4}+20r^{3}-188r^{2}+308r+360}\mathcal{I}^{(5)}_{2}$
	$\displaystyle+\frac{2r(r^{4}+24r^{3}-306r^{2}+340r+600)m_{b}^{2}}{3(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(1)}_{3}$
	$\displaystyle-\frac{(21r^{5}+1115r^{4}+1003r^{3}-3945r^{2}-2324r+1560)m_{b}^{2}}{12(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(2)}_{3}$
	$\displaystyle-\frac{(r^{5}-109r^{4}+498r^{3}+369r^{2}-996r-720)m_{b}^{2}}{3(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(3)}_{3}$
	$\displaystyle+\frac{(r^{5}+1251r^{4}+14599r^{3}+9475r^{2}+700r-1320)m_{b}^{2}}{12(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(4)}_{3}$
	$\displaystyle-\frac{2(4r^{6}+209r^{5}-1014r^{4}+4477r^{3}-5416r^{2}-3468r+3360)m_{b}^{2}}{3(r-2)(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(5)}_{3}$
	$\displaystyle+\frac{2(3r^{5}+45r^{4}-25r^{3}+230r^{2}-600r-720)m_{b}^{2}}{3(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(6)}_{3}$
	$\displaystyle-\frac{72(r-4)r(r+1)m_{b}^{2}}{(r-2)(r^{4}+20r^{3}-188r^{2}+308r+360)}\mathcal{I}^{(7)}_{3}$
	$\displaystyle+\frac{19r^{5}+416r^{4}-3664r^{3}+4940r^{2}+9296r+1920}{r^{5}+20r^{4}-188r^{3}+308r^{2}+360r}\ln(m_{b})$
	$\displaystyle+\frac{-104r^{7}-3176r^{6}+29871r^{5}-103006r^{4}+32774r^{3}+441486r^{2}+306648r+54720}{9(r-4)r(2r+1)(r^{4}+20r^{3}-188r^{2}+308r+360)}$
	$\displaystyle+11\ln(2).$		(51)

\displaystyle a_{2}

\displaystyle=

\displaystyle-\frac{1}{6}{\mathcal{I}^{(2)}_{2}}+\frac{-11r^{4}+65r^{3}+148r^{2}-628r-360}{18(10r^{4}-37r^{3}-104r^{2}+488r+360)}-\frac{\ln(2m_{b})}{3},

(52)

	$\displaystyle b_{2}$	$\displaystyle=$	$\displaystyle\frac{10r^{4}-79r^{3}-68r^{2}+1136r+720}{3m_{b}^{2}r(10r^{4}-37r^{3}-104r^{2}+488r+360)}{\mathcal{I}_{1}}-\frac{10r^{5}-17r^{4}-262r^{3}+352r^{2}+2632r+1440}{6r(10r^{4}-37r^{3}-104r^{2}+488r+360)}{\mathcal{I}^{(3)}_{2}}$		(53)
			$\displaystyle-\frac{11r^{5}+121r^{4}-730r^{3}-1724r^{2}+6096r+4320}{18r(10r^{4}-37r^{3}-104r^{2}+488r+360)}-\frac{\ln(2m_{b})}{3},$		(53)

	$\displaystyle c_{2}$
	$\displaystyle=\frac{3096r^{8}-24996r^{7}-30164r^{6}+323237r^{5}+146433r^{4}-521269r^{3}-663026r^{2}-317208r-54720}{18(r-4)r(2rm_{b}+m_{b})^{2}(10r^{4}-37r^{3}-104r^{2}+488r+360}\mathcal{I}_{1}$
	$\displaystyle+\frac{1536r^{8}+17268r^{7}+111220r^{6}+635537r^{5}+1298688r^{4}+1176860r^{3}+526618r^{2}+107484r+7200}{36r(2r+1)^{2}(10r^{4}-37r^{3}-104r^{2}+488r+360}$
	$\displaystyle\times\mathcal{I}^{(1)}_{2}-\frac{192r^{5}+1041r^{4}+12656r^{3}+93670r^{2}-18084r-65880}{36(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(2)}_{2}$
	$\displaystyle+\frac{-289r^{5}-193r^{4}+9762r^{3}-36644r^{2}+65680r+87360}{36(r-4)(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(3)}_{2}$
	$\displaystyle-\frac{48r^{5}+401r^{4}+2466r^{3}+11496r^{2}+23680r+6360}{9(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(4)}_{2}$
	$\displaystyle+\frac{r(2r^{3}-8r^{2}-57r-36)}{10r^{4}-37r^{3}-104r^{2}+488r+360}\mathcal{I}^{(5)}_{2}$
	$\displaystyle+\frac{2r(4r^{4}-189r^{3}+48r^{2}+247r+60)m_{b}^{2}}{3(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(1)}_{3}$
	$\displaystyle-\frac{(822r^{5}+4154r^{4}-12011r^{3}+14334r^{2}-13424r-21840)m_{b}^{2}}{48(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(2)}_{3}$
	$\displaystyle+\frac{(56r^{5}+427r^{4}-4071r^{3}+8370r^{2}+384r-10080)m_{b}^{2}}{12(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(3)}_{3}$
	$\displaystyle+\frac{(1414r^{5}+12714r^{4}+50113r^{3}+13030r^{2}-4208r+18480)m_{b}^{2}}{48(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(4)}_{3}$
	$\displaystyle-\frac{(340r^{6}-796r^{5}+3933r^{4}+6358r^{3}-28252r^{2}+1224r+4800)m_{b}^{2}}{6(r-2)(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(5)}_{3}$
	$\displaystyle+\frac{(30r^{5}+21r^{4}+1760r^{3}-3004r^{2}+7680r+10080)m_{b}^{2}}{6(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(6)}_{3}$
	$\displaystyle-\frac{18r(r^{3}-10r-8)m_{b}^{2}}{(r-2)(10r^{4}-37r^{3}-104r^{2}+488r+360)}\mathcal{I}^{(7)}_{3}$
	$\displaystyle+\frac{291r^{5}-915r^{4}-4006r^{3}+13632r^{2}+17032r+3840}{20r^{5}-74r^{4}-208r^{3}+976r^{2}+720r}\ln(m_{b})$
	$\displaystyle+\frac{-7688r^{7}+41692r^{6}+62841r^{5}-824998r^{4}+799850r^{3}+2723892r^{2}+1386288r+167040}{36(r-4)r(2r+1)(10r^{4}-37r^{3}-104r^{2}+488r+360)}$
	$\displaystyle+\frac{11}{2}\ln(2),$		(54)

\displaystyle a_{3}

\displaystyle=

\displaystyle b_{3}=0,

(55)

	$\displaystyle c_{3}$
	$\displaystyle=-\frac{1}{9(r-4)r(2rm_{b}+m_{b})^{2}(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}$
	$\displaystyle\times(-1512r^{9}-1032r^{8}+32016r^{7}-11119r^{6}-117845r^{5}+317166r^{4}+933424r^{3}+943568r^{2}$
	$\displaystyle+451200r+80640)\mathcal{I}_{1}$
	$\displaystyle+\frac{2(444r^{8}+4470r^{7}+24088r^{6}+165809r^{5}+332502r^{4}+265415r^{3}+89104r^{2}+6612r-1440)}{9r(2r+1)^{2}(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(1)}_{2}$
	$\displaystyle-\frac{4(48r^{5}-105r^{4}+952r^{3}+7196r^{2}-6360r-8640)}{9(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(2)}_{2}$
	$\displaystyle-\frac{8(55r^{5}-281r^{4}+462r^{3}+2516r^{2}-15472r-12480)}{9(r-4)(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(3)}_{2}$
	$\displaystyle-\frac{4(21r^{5}+619r^{4}+2706r^{3}+10896r^{2}+12008r+480)}{9(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(4)}_{2}$
	$\displaystyle-\frac{4r(110r^{4}+249r^{3}-339r^{2}+386r+480)m_{b}^{2}}{3(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(1)}_{3}$
	$\displaystyle-\frac{2(57r^{5}-367r^{4}-608r^{3}+3210r^{2}-2852r-3120)m_{b}^{2}}{3(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(2)}_{3}$
	$\displaystyle+\frac{2(83r^{5}-125r^{4}-336r^{3}+4554r^{2}-4044r-5760)m_{b}^{2}}{3(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(3)}_{3}$
	$\displaystyle+\frac{8(7r^{5}+93r^{4}+376r^{3}-347r^{2}+196r+660)m_{b}^{2}}{3(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(4)}_{3}$
	$\displaystyle-\frac{4(18r^{6}+103r^{5}+247r^{4}+2660r^{3}+2348r^{2}-720r+480)m_{b}^{2}}{3(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(5)}_{3}$
	$\displaystyle-\frac{16(6r^{5}-51r^{4}-80r^{3}+28r^{2}-1920r-1440)m_{b}^{2}}{3(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\mathcal{I}^{(6)}_{3}$
	$\displaystyle+\frac{2(36r^{6}+101r^{5}-731r^{4}-60r^{3}+7432r^{2}+7736r+1920)}{r(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}\ln(m_{b})$
	$\displaystyle+\frac{-2052r^{8}+4274r^{7}+57269r^{6}-136869r^{5}-528734r^{4}+1033984r^{3}+1897728r^{2}+747072r+57600}{9(r-4)r(2r+1)(9r^{5}+11r^{4}-134r^{3}+140r^{2}+1336r+720)}.$
			(56)

Acknowledgements.

This work is supported by the Natural Science Foundation of China under the Grant No. 12065006.

References

(1) B. Guberina, J. H. Kuhn, R. D. Peccei and R. Ruckl, Rare Decays of the Z0, Nucl. Phys. B 174 (1980), 317-334.
(2) W. Y. Keung, Off Resonance Production of Heavy Vector Quarkonium States in $e^{+}e^{-}$ Annihilation, Phys. Rev. D 23 (1981), 2072.
(3) V. D. Barger, K. m. Cheung and W. Y. Keung, Z BOSON DECAYS TO HEAVY QUARKONIUM, Phys. Rev. D 41 (1990), 1541.
(4) E. Braaten, K. m. Cheung and T. C. Yuan, Z0 decay into charmonium via charm quark fragmentation, Phys. Rev. D 48 (1993), 4230-4235.
(5) S. Fleming, Electromagnetic production of quarkonium in Z0 decay, Phys. Rev. D 48, R1914-R1916 (1993).
(6) G. Alexander et al. [OPAL], Prompt J $\psi$ production in hadronic Z0 decays, Phys. Lett. B 384 (1996), 343-352.
(7) P. Abreu et al. [DELPHI], Search for promptly produced heavy quarkonium states in hadronic Z decays, Z. Phys. C 69 (1996), 575-584.
(8) K. m. Cheung, W. Y. Keung and T. C. Yuan, Color octet quarkonium production at the $Z$ pole, Phys. Rev. Lett. 76 (1996), 877-880.
(9) P. L. Cho, Prompt upsilon and psi production at LEP, Phys. Lett. B 368 (1996), 171-178.
(10) S. Baek, P. Ko, J. Lee and H. S. Song, Color octet heavy quarkonium productions in Z0 decays at LEP, Phys. Lett. B 389 (1996), 609-615.
(11) P. Ernstrom, L. Lonnblad and M. Vanttinen, Evolution effects in $Z^{0}$ fragmentation into charmonium, Z. Phys. C 76 (1997), 515-521.
(12) E. M. Gregores, F. Halzen and O. J. P. Eboli, Prompt charmonium production in Z decays, Phys. Lett. B 395 (1997), 113-117.
(13) C. f. Qiao, F. Yuan and K. T. Chao, A Crucial test for color octet production mechanism in $Z^{0}$ decays, Phys. Rev. D 55 (1997), 4001-4004.
(14) S. Baek, P. Ko, J. Lee and H. S. Song, Color octet mechanism and $J/\psi$ polarization at LEP, Phys. Rev. D 55 (1997), 6839-6843.
(15) M. Acciarri et al. [L3], Inclusive $J/\psi$ , $\psi^{\prime}$ and chi( $c$ ) production in hadronic $Z$ decays, Phys. Lett. B 407 (1997), 351-360.
(16) M. Acciarri et al. [L3], Heavy quarkonium production in $Z$ decays, Phys. Lett. B 453 (1999), 94-106.
(17) C. G. Boyd, A. K. Leibovich and I. Z. Rothstein, J / psi production at LEP: Revisited and resummed, Phys. Rev. D 59 (1999), 054016.
(18) M. Acciarri et al. [L3], Heavy quarkonium production in $Z$ decays, Phys. Lett. B 453 (1999), 94-106.
(19) R. Li and J. X. Wang, The next-to-leading-order QCD correction to inclusive $J/\psi(\Upsilon)$ production in $Z^{0}$ decay, Phys. Rev. D 82 (2010), 054006.
(20) L. C. Deng, X. G. Wu, Z. Yang, Z. Y. Fang and Q. L. Liao, $Z_{0}$ Boson Decays to $B^{(*)}_{c}$ Meson and Its Uncertainties, Eur. Phys. J. C 70 (2010), 113-124.
(21) Z. Yang, X. G. Wu, L. C. Deng, J. W. Zhang and G. Chen, Production of the $P$ -Wave Excited $B_{c}$ -States through the $Z^{0}$ Boson Decays, Eur. Phys. J. C 71 (2011), 1563.
(22) C. F. Qiao, L. P. Sun and R. L. Zhu, The NLO QCD Corrections to $B_{c}$ Meson Production in $Z^{0}$ Decays, JHEP 08 (2011), 131.
(23) T. C. Huang and F. Petriello, Rare exclusive decays of the Z-boson revisited, Phys. Rev. D 92 (2015) no.1, 014007.
(24) J. Jiang, L. B. Chen and C. F. Qiao, QCD NLO corrections to inclusive $B_{c}^{*}$ production in $Z^{0}$ decays, Phys. Rev. D 91 (2015) no.3, 034033.
(25) Q. L. Liao, Y. Yu, Y. Deng, G. Y. Xie and G. C. Wang, Excited heavy quarkonium production via Z⁰ decays at a high luminosity collider, Phys. Rev. D 91 (2015) no.11, 114030.
(26) G. Aad et al. [ATLAS], Search for Higgs and $Z$ Boson Decays to $J/\psi\gamma$ and $\Upsilon(nS)\gamma$ with the ATLAS Detector, Phys. Rev. Lett. 114 (2015) no.12, 121801.
(27) Y. Grossman, M. König and M. Neubert, Exclusive Radiative Decays of W and Z Bosons in QCD Factorization, JHEP 04 (2015), 101.
(28) 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 (2018) no.1, 016009.
(29) A. K. Likhoded and A. V. Luchinsky, Double Charmonia Production in Exclusive $Z$ Boson Decays, Mod. Phys. Lett. A 33 (2018) no.14, 1850078.
(30) M. Aaboud et al. [ATLAS], Searches for exclusive Higgs and $Z$ boson decays into $J/\psi\gamma$ , $\psi(2S)\gamma$ , and $\Upsilon(nS)\gamma$ at $\sqrt{s}=13$ TeV with the ATLAS detector, Phys. Lett. B 786 (2018), 134-155.
(31) Z. Sun and H. F. Zhang, Next-to-leading-order QCD corrections to the decay of $Z$ boson into $\chi_{c}(\chi_{b})$ , Phys. Rev. D 99 (2019) no.9, 094009.
(32) J. P. Lansberg, New Observables in Inclusive Production of Quarkonia, Phys. Rept. 889 (2020), 1-106.
(33) A. M. Sirunyan et al. [CMS], Search for rare decays of Z and Higgs bosons to J $/\psi$ and a photon in proton-proton collisions at $\sqrt{s}=$ 13 TeV, Eur. Phys. J. C 79 (2019) no.2, 94.
(34) Z. Sun, The studies on $Z\rightarrow\Upsilon(1S)+g+g$ at the next-to-leading-order QCD accuracy, Eur. Phys. J. C 80 (2020) no.4, 311.
(35) Z. Sun and H. F. Zhang, Comprehensive studies of $\Upsilon(1S)$ inclusive production in Z boson decay, JHEP 06 (2021), 152.
(36) Z. Sun, X. Luo and Y. Z. Jiang, Impact of $Z\to\eta_{c,b}+g+g$ on the inclusive $\eta_{c,b}$ meson yield in Z-boson decay, Phys. Rev. D 106 (2022) no.3, 034001.
(37) X. C. Zheng, X. G. Wu, X. J. Zhan, H. Zhou and H. T. Li, Next-to-leading order QCD corrections to $Z\to\eta_{Q}+Q+\bar{Q}$ , Phys. Rev. D 106 (2022) no.9, 094008.
(38) W. L. Sang, D. S. Yang and Y. D. Zhang, Z boson radiative decays to a P-wave quarkonium at NNLO and LL accuracy, Phys. Rev. D 106 (2022) no.9, 094023.
(39) 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 (2023) no.1, 014021.
(40) 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 (2023), 120.
(41) 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]].
(42) G. Y. Wang, X. G. Wu, X. C. Zheng, J. Yan and J. W. Zhang, Improved analysis of double $J/\psi$ production in Z-boson decay, Eur. Phys. J. C 84 (2024) no.5, 544.
(43) 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 $\mathcal{O}(\alpha_{s}v^{2})$ corrections, Phys. Rev. D 109 (2024) no.7, 074004.
(44) 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 (2019), 134811.
(45) 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 (2023), 137534.
(46) A. K. Likhoded and A. V. Luchinsky, Double Charmonia Production in Exclusive $Z$ Boson Decays, Mod. Phys. Lett. A 33 (2018) no.14, 1850078.
(47) D. N. Gao and X. Gong, Note on rare Z-boson decays to double heavy quarkonia*, Chin. Phys. C 47 (2023) no.4, 043106.
(48) Y. J. Zhang, Y. j. Gao and K. T. Chao, Next-to-leading order QCD correction to e+ e- — $>$ J / psi + eta(c) at s**(1/2) = 10.6-GeV, Phys. Rev. Lett. 96 (2006), 092001.
(49) B. Gong and J. X. Wang, QCD corrections to $J/\psi$ plus $\eta_{c}$ production in $e^{+}e^{-}$ annihilation at $S^{(1/2)}$ = 10.6-GeV, Phys. Rev. D 77 (2008), 054028.
(50) Y. J. Zhang, Y. Q. Ma and K. T. Chao, Factorization and NLO QCD correction in $e^{+}e^{-}\to J/\psi(\psi(2S))+\chi_{c0}$ at B Factories, Phys. Rev. D 78 (2008), 054006.
(51) N. Brambilla, S. Eidelman, B. K. Heltsley, R. Vogt, G. T. Bodwin, E. Eichten, A. D. Frawley, A. B. Meyer, R. E. Mitchell and V. Papadimitriou, et al. Heavy Quarkonium: Progress, Puzzles, and Opportunities, Eur. Phys. J. C 71 (2011), 1534.
(52) H. R. Dong, F. Feng and Y. Jia, $O(\alpha_{s})$ corrections to $J/\psi+\chi_{cJ}$ production at $B$ factories, JHEP 10 (2011), 141, [erratum: JHEP 02 (2013), 089].
(53) Z. Sun, X. G. Wu, Y. Ma and S. J. Brodsky, Exclusive production of $J/\psi+\eta_{c}$ at the $B$ factories Belle and Babar using the principle of maximum conformality, Phys. Rev. D 98 (2018) no.9, 094001.
(54) Z. Sun, Next-to-leading-order study of $J/\psi$ angular distributions in $e^{+}e^{-}\to J/\psi+\eta_{c},\chi_{cJ}$ at $\sqrt{s}\approx 10.6$ GeV, JHEP 09 (2021), 073.
(55) G. T. Bodwin, E. Braaten and G. P. Lepage, Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium, Phys. Rev. D 51 (1995), 1125-1171, [erratum: Phys. Rev. D 55 (1997), 5853].
(56) J. B. Guimarães da Costa et al. [CEPC Study Group], CEPC Conceptual Design Report: Volume 2 - Physics & Detector, [arXiv:1811.10545 [hep-ex]].
(57) A. Petrelli, M. Cacciari, M. Greco, F. Maltoni and M. L. Mangano, NLO production and decay of quarkonium, Nucl. Phys. B 514 (1998), 245-309.
(58) T. Hahn, Generating Feynman diagrams and amplitudes with FeynArts 3, Comput. Phys. Commun. 140 (2001), 418-431 .
(59) R. Mertig, M. Bohm and A. Denner, FEYN CALC: Computer algebraic calculation of Feynman amplitudes, Comput. Phys. Commun. 64 (1991), 345-359.
(60) J. G. Korner, D. Kreimer and K. Schilcher, A Practicable gamma(5) scheme in dimensional regularization, Z. Phys. C 54 (1992), 503-512.
(61) F. Feng, $\tt{Apart}$ : A Generalized Mathematica Apart Function, Comput. Phys. Commun. 183 (2012), 2158-2164.
(62) A. V. Smirnov, Algorithm FIRE – Feynman Integral REduction, JHEP 10 (2008), 107.
(63) T. Hahn and M. Perez-Victoria, Automatized one loop calculations in four-dimensions and D-dimensions, Comput. Phys. Commun. 118 (1999), 153-165.
(64) E. J. Eichten and C. Quigg, Quarkonium wave functions at the origin, Phys. Rev. D 52 (1995), 1726-1728.
(65) B. Gong, L. P. Wan, J. X. Wang and H. F. Zhang, Complete next-to-leading-order study on the yield and polarization of $\Upsilon(1S,2S,3S)$ at the Tevatron and LHC, Phys. Rev. Lett. 112 (2014) no.3, 032001.
(66) H. Han, Y. Q. Ma, C. Meng, H. S. Shao, Y. J. Zhang and K. T. Chao, $\Upsilon(nS)$ and $\chi_{b}(nP)$ production at hadron colliders in nonrelativistic QCD, Phys. Rev. D 94 (2016) no.1, 014028.
(67) Y. Feng, B. Gong, L. P. Wan and J. X. Wang, An updated study of $\Upsilon$ production and polarization at the Tevatron and LHC, Chin. Phys. C 39 (2015) no.12, 123102.
(68) P. A. Zyla et al. [Particle Data Group], Review of Particle Physics, PTEP 2020 (2020) no.8, 083C01.

	$\displaystyle\bigg{\|}\left(\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha_{s}}+\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}\right)+\left(\mathcal{M}_{\alpha^{\frac{3}{2}}}+\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}\right)\bigg{\|}^{2}$	(4)
$\displaystyle=$	$\displaystyle\big{\|}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha_{s}}\big{\|}^{2}+2\textrm{Re}\left(\mathcal{M}^{*}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}\right)$
	$\displaystyle+2\textrm{Re}\left(\mathcal{M}^{}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{3}{2}}}\right)+2\textrm{Re}\left(\mathcal{M}^{}_{\alpha^{\frac{1}{2}}\alpha_{s}}\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}+\mathcal{M}^{*}_{\alpha^{\frac{3}{2}}}\mathcal{M}_{\alpha^{\frac{1}{2}}\alpha^{2}_{s}}\right)$
	$\displaystyle+\big{\|}\mathcal{M}_{\alpha^{\frac{3}{2}}}\big{\|}^{2}+2\textrm{Re}\left(\mathcal{M}^{*}_{\alpha^{\frac{3}{2}}}\mathcal{M}_{\alpha^{\frac{3}{2}}\alpha_{s}}\right)+\cdots.$

Studies of ZZ boson decay into double Υ\Upsilon mesons at the NLO QCD accuracy

Abstract

pacs:

I Introduction

II Calculation formalism

II.1 Theoretical Framework

II.2 The decay width

II.2.1 LO

II.2.2 NLO

III Phenomenological results

IV Summary

Appendix A

A.1 Master Integrals

A.2 NLO coefficients

A.2.1 bb¯[3S11]+bb¯[3PJ0]b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{0}]

A.2.2 bb¯[3S11]+bb¯[3P11]b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{1}^{1}]

A.2.3 bb¯[3S11]+bb¯[3P21]b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{2}^{1}]

Acknowledgements.

References

Studies of $Z$ boson decay into double $\Upsilon$ mesons at the NLO QCD accuracy

A.2.1 $b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{J}^{0}]$

A.2.2 $b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{1}^{1}]$

A.2.3 $b\bar{b}[^{3}S_{1}^{1}]+b\bar{b}[^{3}P_{2}^{1}]$