Prospects of discovering stable double-heavy tetraquarks at a Tera-$Z$ factory

Motivated by a number of theoretical considerations, predicting the deeply bound double-heavy tetraquarks $T^{\{bb\}}_{[\bar u \bar d]}$, $T^{\{bb\}}_{[\bar u \bar s]}$ and $T^{\{bb\}}_{[\bar d \bar s]}$, we explore the potential of their discovery at Tera-$Z$ factories. Using the process $Z \to b \bar b b \bar b$, we calculate, employing the Monte Carlo generators MadGraph5$\_$aMC@NLO and Pythia6, the phase space configuration in which the~$b b$ pair is likely to fragment as a diquark. In a jet-cone, defined by an invariant mass interval $m_{bb}<M_{T^{\{bb\}}_{[\bar q \bar q']}} + \Delta M$, the sought-after tetraquarks $T^{\{bb\}}_{[\bar q \bar q^\prime]}$ as well as the double-bottom baryons,~$\Xi_{bb}^{0,-}$, and $\Omega_{bb}^-$, can be produced. Using the heavy quark--diquark symmetry, we estimate $\mathcal{B} (Z \to T^{\{bb\}}_{[\bar u \bar d]} + \; \bar b \bar b) = (1.4^{+1.1}_{-0.5}) \times 10^{-6}$, and about a half of this for the $T^{\{bb\}}_{[\bar{u}\bar{s}]}$ and $T^{\{bb\}}_{[\bar d \bar s]}$. We also present an estimate of their lifetimes using the heavy quark expansion, yielding $\tau(T^{\{bb\}}_{[\bar q \bar q^\prime]}) \simeq 800$~fs. Measuring the tetraquark masses would require decays, such as $T^{\{bb\} -}_{[\bar u \bar d]} \to B^- D^- \pi^+$, $T^{\{bb\} -}_{[\bar u \bar d]} \to J/\psi \overline K^0 B^-$, $T^{\{bb\} -}_{[\bar u \bar d]} \to J/\psi K^- \overline B^0$, $T^{\{bb\} -}_{[\bar u \bar s]} \to \Xi_{bc}^0 \Sigma^-$, and $T^{\{bb\} 0}_{[\bar d \bar s]} \to \Xi_{bc}^0 \bar\Sigma^0$, with subsequent decay chains in exclusive non-leptonic final states. We estimate a couple of the decay widths and find that the product branching ratios do not exceed~$10^{-5}$. Hence, a good fraction of these modes will be required for a discovery of $T^{\{bb\}}_{[\bar q \bar q']}$ at a Tera-$Z$ factory.

, we explore the potential of their discovery at Tera-Z factories. Using the process Z → bbbb, we calculate, employing the Monte Carlo generators MadGraph5 aMC@NLO and Pythia6, the phase space configuration in which the bb pair is likely to fragment as a diquark. In a jet-cone, defined by an invariant mass interval m bb < M T {bb} → Ξ 0 bcΣ 0 , with subsequent decay chains in exclusive non-leptonic final states. We estimate a couple of the decay widths and find that the product branching ratios do not exceed 10 −5 . Hence, a good fraction of these modes will be required for a discovery of T {bb} [qq ′ ] at a Tera-Z factory.

I. INTRODUCTION
The experimental discovery of exotic, hidden-charm and hidden-beauty states, has opened a new field in hadron spectroscopy. The exotic states, called X, Y , Z and P c , have been analysed in a number of theoretical models. They have been claimed to be hybrid quarkonia, hadron molecules, coupled-channel or threshold effects, and multiquark states, see [1][2][3][4][5] for recent reviews and extensive references therein. Their dynamics is very much an open question with the models dividing themselves approximately in two classes, those reflecting residual QCD long-distance effects, dominated by meson exchanges, and those reflecting genuine short-distance forces, induced by gluon exchanges. A particular realisation of the latter class of models assumes that heavy baryons and tetraquarks may be viewed as diquark-quark and diquark-antidiquark objects, respectively, with the diquarks having well-defined color, spin and flavor quantum numbers. Indeed, if tetraquarks, which are stable against strong and radiative decays, could be found in experiments, this would provide an irrefutable evidence of compact diquarks as building blocks of hadronic matter.
The objects of our interest in this paper are the doublebottom J P = 1 + tetraquarks T  Fig. 1, which have evoked lately a lot of theoretical and phenomenological interest [6][7][8][9][10][11], though the possibility of stable multiquark states was already pointed out a long time ago [12,13]. Likewise, estimates of the tetraquark masses with two heavy quarks carried out in quark models well over a decade ago also predicted stable tetraquarks [14].
Concentrating on the T {bb} [ūd] state, which is a J P = 1 + , I = 0 tetraquark, consisting of the S-wave bound axialvector {bb} diquark and the scalar light [ūd] antidiquark, its mass is pitched at (10389 ± 21) MeV [6], lying about 215 MeV below the BB * threshold. Other esti- mates are somewhat higher in mass, with a Q-value of −121 MeV [7], −189 ± 10 MeV [8], and −60 +30 −38 MeV [9]. In all likelihood, T likewise are below their corresponding mesonic thresholds, and hence they are also expected to decay weakly. Since no weakly decaying multiquark state has so far been observed, their observation would herald a new era in genuine multiquark physics.
At a Z-factory, the underlying partonic process for the searches of such tetraquarks is which has been measured at LEP, having the branching ratio [15]: The production of a double-bottom tetraquark, T volves the formation of the bb-diquark and its fragmentation, producing a bb-diquark jet. The conducive configuration for the diquark formation is the one where the two b-quarks (or two b-antiquarks) are almost collinear, satisfying θ bb ≤ δ, where θ bb is the angle between the bquark momenta and δ is the cone apex angle, and have a small relative velocity, i. e., their energies E b1 and E b2 differ by a small amount, (E b1 − E b2 ) / (E b1 + E b2 ) ≤ ǫ, with (ǫ, δ) ≪ 1, with an example shown in Fig. 2. The bb-diquark jet can be defined, like a quark jet, by a jet resolution parameter, such as the invariant mass, or a Sterman-Weinberg jet cone [16]. A judicious choice of these cut-off parameters forces the two b-quarks to overlap in the phase space, which then fragment as a diquark. A tetraquark of the commensurate quark flavor is formed by picking up a light antidiquark pairūd, us, ords from the debris of the jet, which consists of mostly soft pions or kaons. Likewise, the bb-diquark will also fragment into double-bottom baryons, Ξ 0 bb (bbu), Ξ − bb (bbd), and Ω − bb (bbs), shown in the right-hand frame in Fig. 3. Outside the bb jet-cone, no double-bottom hadrons (tetraquarks or double-bottom baryons) will be produced. Hence, the hadronic texture of the bb-diquark jet is anticipated to be different from the fragmentation of two b-quark jets, whose fragmentation products are B-mesons or b-baryons.
We calculate the branching fraction by defining the (bb) jet with a cut-off on the bb-pair invariant mass, ∆M . An estimate of ∆M can be obtained from the inclusive B c -meson production in Z-boson decays. Using the inclusive cross section σ(e + e − → Z → bbcc), obtained via MadGraph [17] and Pythia6 [18], and the inclusive B c -meson cross section σ(e + e − → Z → B c +b+c), obtained by using the NRQCD-based calculations [19], one can evaluate the fraction f (cb → B c ) for the cb-pair fusion into the B c mesons 2 . The phase space of the fragmentation process to produce the B c -meson is limited by m cb < M cb cut , with M cb cut being the maximum value of the invariant mass in which the cb-fusion takes place. Beyond this cut, the b-andc-quarks fragment independently. We use M cb cut to estimate ∆M , which yields us the partial branching ratio in Eq. (3). The details of this calculation are given in Section II.
The final step in the calculation of the branching ratio is to estimate the probability of the diquark-jet (bb) jet (∆M ) to fragment into T into the double-bottom baryons, shown in the right panel in Fig. 3. This involves estimating the relative probability of emitting a light (anti)-diquark [ūd], [ūs], or [ds] from the vacuum in the presence of a bb-diquark color source to that of picking up a light quark q (q = u, d, s) from the qq-pair produced in the similar situation. The two probabilities are related by heavy quark-heavy diquark symmetry.
This probability is similar to the relative probability that a b-quark fragments into a b-baryon, such as Λ b = bud, which involves picking up the diquark [ud] from the vacuum, to that of a b-quark fragmenting into a heavylight meson B − = bū (orB 0 = bd). Denoting these fractions by f Λ b and f Bu (f B d ), respectively, we need to know the ratio f Λ b /(f Bu +f B d ). As is well-known, this ratio has been measured differently at the Fermilab-Tevatron [20]: −0.056−0.087 , at the LHC by LHCb Collab. [21], which finds the significant dependence of this ratio on p T  [15]. The dynamical reason behind these variations is not clear, but as we are estimating the Z-boson branching ratios, we shall use the average from the LEP experiments.
We also present an estimate of the total lifetimes of the tetraquarks T {bb} [qq ′ ] based on the heavy quark expansion. As each of the two b-quarks in the tetraquark will decay independently, the decay widths of T These baryons also remain to be discovered, and their discovery channels have been recently presented in 3 The LHCb analysis has been updated in Ref. [22]. Refs. [26,27]. This underscores the huge potential of the Tera-Z colliders in mapping out the landscape of the double-heavy baryons and double-heavy tetraquarks.
The rest of this work is organised as follows. In Sec. II, we present an estimate of the branching ratios of doublebottom tetraquarks in Z-boson decays. Lifetimes of the tetraquarks are discussed in Sec. III. In Secs. IV and V, we discuss the discovery modes for the T are anticipated to be produced in the partonic process Z → bbbb, with the two b-quarks (or two b-antiquarks) moving collinearly and having a small relative velocity. A reasonable cut on the invariant mass of the two b-quarks ensures these constraints.
To calculate the branching ratio in (4), we invoke the decay Z → B c + X, which has been measured at LEP. The idea is to calculate the fraction of the cb quarks which hadronise into the B c -mesons, f (cb → B c ), in the underlying process Z → bbcc. Theoretically, the B c production cross section has been calculated at the leading order (LO) in the NRQCD framework, yielding: [19]. The central (upper, lower) value corresponds to the input values of the quark masses m b = 4.9 GeV and m c = 1.5 GeV (m b = 5.3 GeV and m c = 1.2 GeV; m b = 4.5 GeV and m c = 1.8 GeV). With the same input, we have generated three sets of 10000 LO parton showered e + e − → bbcc events at the Z-boson mass using MadGraph [17] and Pythia6 [18], varying the bottom and charmed quark mass values. (The mass parameters used by MadGraph are the pole masses. The difference between the 1S mass value used in [19] and the pole mass of a heavy quark is of the α 2 s order (see (1) of [28]), which can be neglected safely.) The cross section σ(e + e − → Z → bbcc), evaluated using MadGraph  Table I. We then generate 10000 showered e + e − → bbbb events at the Z-boson mass as the centre-of-mass energy with MadGraph [17] and Pythia6 [18] at the NLO of QCD, with the b-quark pair invariant mass distribution displayed in Fig. 5. With the cross sections by Mad-Graph, the branching ratio B(Z → bbbb) is obtained to be 3.23×10 −4 , which is consistent with the world average (3.6 ± 1.3) × 10 −4 [15].
Next, we employ the same kinematics cuts, namely  +2.2 GeV, and find that 4.9% (5.3%), 8.9% (9.2%) 4 The NLO MadGraph result for the cross section is enhanced by more than a factor of 2 compared with the LO result, but as the corresponding NRQCD calculation is done only in the LO approximation, we use the LO simulation to be consistent with Ref. [19].  3.8% (3.9%) of the e + e − → Z → bbbb events fragment into a hadron with the bb (bb) quark pair. The numerical differences between the charged conjugated processes reflect statistical fluctuations, which will disappear in simulations with higher statistics. We average them and list the results of the ratios f (bb → H {bb} ) in Table II. Thus, our simulation yields The uncertainty in the fragmentation fraction is sizeable, and, within ±1σ, it lies in the range (4 − 9)%.
We show the dependence of the branching ratio B(Z → T

III. LIFETIMES
Before discussing the weak decays, we present an estimate of the lifetimes, the inverse of total decay rates. The decay width of T where m T and p µ T are the mass and four-momentum of the tetraquark, respectively. The factor 1/3 results after averaging over the polarization states λ of the tetraquark with the spin-parity J P = 1 + . The effective electroweak Hamiltonian H ew eff is governing the weak decays to be given in the following discussions. Using the optical theorem, one can rewrite the total decay width as with the transition operator T as the time-ordered nonlocal double insertion of the effective Hamiltonian: The heavy-quark expansion, the operator product expansion in essence, greatly simplifies the decay widths of inclusive decays. Integrating over the off-shell intermediate states in Eq. (10) results in a set of local operators with increasing dimensions. Higher dimensional operators are suppressed by inverse powers of the heavy-quark mass. (For a recent review see Ref. [29].) Up to dimension 6, we have the transition operator expanded as [29]: where G F is the Fermi constant and V CKM are the elements of the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix. The coefficients c i,b are the corresponding short-distance coefficients from the semileptonic and nonleptonic heavy quark decays. These short-distance coefficients have been precisely calculated using QCD perturbation theory [29]. The total decay width is then determined by matrix elements of these operators, and at the leading order in 1/m b , only thebb operator contributes: The matrix element corresponds to the bottom-quark number in the T state and is twice the matrix element for B-meson and Λ b -baryon. Accordingly, we expect that the lifetime of the double-bottom tetraquark is one half of that of the B-meson [15]: This is approximately twice the corresponding lifetime τ (T {bb} [qq ′ ] ) ≃ 367 fs estimated in [6].

IV. WEAK DECAYS OF T {bb} [ūd]
If a state is below the lowest hadronic threshold, its decay width is determined by weak interactions. The decay modes are dominated by flavor-changing charged currents in the effective Hamiltonian [30]: Here, α and β are the quark color indices, and P L = (1 − γ 5 ) /2 is the left-handed projector. The coefficients C 1 and C 2 are determined after matching the effective theory and the Standard Model, being scaledependent quantities.
We start by considering the two-body decays of T

{bb} [ūd]
and continue by taking up the three-body ones which can be divided into the purely mesonic and baryonic decays. decays due to the b-quark decay b → c + d +ū (left panel) and b → c + s +c (right panel).

A. Two-Body Baryonic Decays
As the doubly bottom tetraquark T {bb} [ūd] (10482) 5 consists of two b-quarks and two antiquarks,ū andd, the weak decay of one of the b-quarks results into two possible Cabibbo-allowed decay channels at the quark level: b → c+d+ū or b → c+s+c. From three quarks and three antiquarks one can easily construct baryon and antibaryon, as shown in Fig. 7. Baryonic states produced in these modes have in addition to the well-known antibaryons, also the so-far undiscovered bottom-charmed baryons shown in Fig. 4. They are being searched for at the LHC and are also a part of the experimental programme at future electron-positron colliders. Amplitudes of these decays are non-factorisable as a quark and antiquark produced in the weak transition hadronise into baryon and anti-baryon, respectively (see Fig. 7). Taking into account the axial-vector nature of the tetraquark, with J P = 1 + , the general form of the decay amplitude induced by the b → c + d +ū quark channel is as follows: where u(p Ξ bc ) and v(p p ) are the wave-functions of the Ξ 0 bc -baryon and antiproton with the four-momenta p Ξ bc and p p , respectively, q = p Ξ bc −p p , ε µ T (p T ) is the polarisation vector of the axial-vector tetraquark with the fourmomentum p T = p Ξ bc + p p and mass m T . The f Ξ bcp → Ω 0 bcΛ − c decay described by the right diagram in Fig. 7. 5 The mass assignment is according to the predictions of Ref. [7]. decays due to the b-quark decay b → c + d +ū.
One needs to further reconstruct the bottom-charmed baryons Ξ 0 bc and Ω 0 bc . Including the decay chain Ξ 0 bc → Λ b K − π + , we find that the two-body baryonic decay modes of the double-bottom tetraquarks are expected to have branching fractions of order of 10 −6 . decays can be identified. They are shown in Fig. 8. The factorisable amplitudes of these decays can be written as follows: where a eff 1 = C 1 + C 2 /N c , with N c = 3 being the number of quark colors, and the standard definition of the πmeson leptonic decay constant is used: Next, we need to parametrise the transition matrix elements from the tetraquark state to the double-meson one, having the total angular momentum J BD and zero electric charge (specified by the superscript index), in terms of the form factors. For the case of the D-and B-mesons, which are pseudoscalar particles, the total momentum J BD is completely determined by the angular momentum L BD of the system. For the BD-system in the S-wave, the transition matrix elements: [ūd] (p T ) , (22) define the S-wave generalised form factors.
One can also expect a production of (B * D) 0 JBD or (BD * ) 0 JBD pairs in the S-wave and the total angular momentum (J B * D = 1 or J BD * = 1) is determined by the spin of the vector meson (S B * = 1 or S D * = 1). This requires another set of S-wave generalised form factors.
Currently, we lack a reliable dynamical approach to calculate the generalised form factors, which prevents us to reliably predict the branching fractions for T →B 0 D 0 π − decays. However, as after the b-quark weak decay, the emitted charmed quark and three spectators have an invariant mass larger than the B − D + andB 0 D 0 thresholds, we expect the split into two hadrons will not cause any dynamical suppression. In this case, the transition T Accordingly, their branching fractions can reach the value of 10 −3 as well.
From Fig. 8, one may speculate that the final B − D + andB 0 D 0 mesons arise from an intermediate bottomcharmed tetraquark state T bc ud . In the case of a scalar (J P = 0 + ) tetraquark, the transition amplitude is governed by the four form factors: where q µ = p µ T − p µ T ′ , ε µ T is the polarization vector of the axial-vector T with the decay amplitude decays due to the b-quark decay b → c + s +c.
The factor 2 arises since there are two b-quarks in the initial state, while the factor 1/3 denotes the spin average. Neglecting the π-meson mass, the decay width can be written as follows: where G F = 1.166 × 10 −5 GeV −2 , f π ≃ 130 MeV, V ud = 0.974, and |V cb | = 40.5 × 10 −3 [15]. The effective coefficient a eff 1 is a scale-dependent quantity but in a wide range of the energy scale its value is close to unity, so we take a eff 1 = 1 in estimates. Using m T = 10.5 GeV, the Ξ bb → Ξ bc transition form factor, A 0 (q 2 = 0) = 0.44 [31], and m T ′ = m T + m c − m b ∼ 7.2 GeV, we obtain the estimate of the partial decay width: [ūd] → T bc ud π − ) ≃ 7.9 × 10 −16 GeV ≃ 1.2 ns −1 , (27) and the branching fraction with account of the lifetime (14): If one uses the B c → J/ψ transition form factor instead, A 0 (q 2 = 0) = 0.53 [32], the results can be enhanced by approximately 50%. As the intermediate state, T bc ud , will subsequently turn into the B − D + andB 0 D 0 states with equal probabilities, and thus we expect

Hidden-Charm Final States
There are some Feynman diagrams in which the hidden-charm mesons, such as J/ψ and ψ ′ , can be produced. The corresponding Feynman diagrams are shown in Fig. 9. Decay channels include decays due to the b-quark decay b → c + s +c. whose CKM factors are V * cb V cs . Their decay branching ratios can be comparable with the B → J/ψK decays [15]: So, it reasonable also to expect branching fractions of the T The other type of diagrams with hidden-charm have a light-and a b-baryon in the final state. A representative example of these decays T

{bb}− [ūd]
→ Ξ 0 bp J/ψ, which can be measured at the colliders, is presented in Fig. 10. As the two-body baryonic decays discussed earlier, the corresponding amplitudes are non-factorisable as the quark pair produced in the weak decay is devided and hadronised into a baryon and meson. An account of a quark pair picked up from the vacuum makes the theoretical analysis even more complicated and we postpone such a discussion for future. [ds] (10643) 0 , their masses are closer to the corresponding hadronic threshold,B * B s , but are still estimated to lie below it. In particular, these tetraquark states are predicted to be below the threshold (10691 MeV) by 48 MeV in Ref. [7]. If this is true, their decay modes are dominated by flavor-changing charged 6 The mass assignment is according to the predictions of Ref. [7]. currents in the effective Hamiltonian (15). Like the decays of the tetraquark T {bb} [ūd] , they have corresponding two-body and three-body decays, which can be divided into the purely mesonic and baryonic modes. They are briefly discussed below.

A. Two-Body Baryonic Decays
As the double-bottom tetraquark T → B 0 s D + π − for the neutral one.

Hidden-Charm Final State
There are channels in which the hidden-charm mesons like J/ψ, ψ ′ , and etc. can be produced: The final states mentioned here are well reconstructed at both electron-positron and hadron colliders.

C. Three-body baryonic decay modes
The other type of diagrams correspond to three-body decays with a light-and bottom baryon in the final state. The most interesting processes could be T

VI. W -EXCHANGE DIAGRAMS
The W -exchange diagrams result into two-body mesonic decays, some of which can be of interest at the LHC. The use of W -exchange decay modes of doubleheavy baryons has been advocated in Ref. [33]. The decays with the J/ψ-meson production are the most likely. So, the processes which can be searched for are T [ds] (10643) 0 →B 0 J/ψ. The corresponding diagrams for these processes are presented in Fig. 11.
The decay amplitude to leading order can be factorized. For the T {bb} [ūs] (10643) − → B − J/ψ decay, as an example, one can write it as follows: where a eff 2 = C 2 +C 1 /N c is the effective Wilson coefficient and the matrix element between the vacuum and J/ψmeson is parametrised in terms of the decay constant: The general decomposition of the transition matrix element can be written in the form similar to the B → A transition, where A is an axial-vector meson [34]. An advantage of these decay modes is that there are only two mesons in the final state. The J/ψ final state can be easily reconstructed in the µ + µ − channel, and the bottom meson can be studied in its decays into twobody final states. For the T [ds] (10643) 0 →B 0 J/ψ, the CKM factors in the transition are V cb V * cs , and thus they can have sizeable branching fractions. Without any additional power suppression, we expect the branching ratios for these two modes to be of order of 10 −3 . For the T {bb} [ūd] (10643) − → B − J/ψ decay, the decay width is suppressed due the CKM factor V cd , leading to a smaller branching fraction by approximately the factor 25.

VII. CONCLUSIONS
We have presented the estimate of the double-bottom hadron production in Z-boson decays. These include the tetraquarks T The lifetimes of the tetraquarks T {bb} [qq ′ ] are estimated, and for the SU (3) F -multiplet they are expected to be approximately similar, about one half of the B-meson lifetimes, τ ≃ 0.8 ps. We also discussed some signature decay modes of these tetraquarks. These include the two-body baryonic decays, such as T → Ω 0 bcΛ − c , whose branching ratios are estimated to be of order of 10 −3 . Likewise, threebody mesonic decays of T are also discussed, and the order of magnitude of their signature decay modes are also presented. Finally, we emphasise the so-called W -exchange diagrams, which give rise to potentially interesting two-body decays [ds] (10643) 0 →B 0 J/ψ. We estimate that the product branching ratios in many of the decay modes discussed here are not expected to exceed 10 −5 . Hence, several of these decay modes will be required to measure the masses of the tetraquarks to establish them firmly. Thus, we recommend to aim for an integrated luminosity of 10 12 Z-bosons at the Tera-Z factories being contemplated now. They will greatly advance the physics of the double-heavy hadrons, both in the mesonic and baryonic sectors.