Exclusive $eta_b$ decay to double $J/\psi$ at next-to-leading order in $alpha_s$

Within the nonrelativistic QCD (NRQCD) factorization framework, we calculate the exclusive decay process \eta_b \to J/\psi J/\psi to next-to-leading order in the strong coupling constant, while at leading order in charm quark relative velocity. It is found that this new contribution to the amplitude is comparable in magnitude with the previously calculated relativistic correction piece, but differs by a phase about 90^\circ. Including this new contribution will increase the previous prediction to B(\eta_b \to J/\psi J/\psi) substantially, thus brightening the discovery potential of this clean hadronic decay channel of \eta_b in the forthcoming LHC experiment.

The quest for the η b meson, the lowest-lying pseudoscalar bottomonium state, started shortly after the discovery of Υ in 1977.A thorough understanding of its properties, such as its mass and decay pattern, will be great benefit to strengthening our knowledge towards heavy quarkonium physics [1].
Tremendous efforts have been invested in various high energy collision experiments during the past thirty years to seek this elusive particle, but convincing evidence about its existence has not emerged until very recently.After many futile searches, there finally comes out an exciting news from Babar collaboration last month [2], which reported the first unambiguous sighting of η b with a 10 σ significance, through the hindered magnetic dipole transition process Υ(3S) → η b γ.It is interesting to note the rather precisely measured hyperfine Υ(1S)-η b mass splitting, 71.4 +2. 3  −3.1 (stat) ± 2.7(syst) MeV, can play a decisive role in discriminating different quarkonium models [1].
Aside from its mass, almost nothing is known about the dynamical properties of the η b , e.g., its decay pattern.Owning to the weaker QCD coupling at bottom mass scale, together with the copious phase space opened up for innumerable decay channels, one in general expects that each η b exclusive decay mode may only be allotted a rather small branching ratio.Indeed, a rough estimate of some 2-body and 3-body hadronic decay channels of η b supports this expectation, and probably η b disintegrates primarily into final states of high multiplicities [3].Some "golden" modes have been proposed for hunting η b , such as η b → J/ψJ/ψ [4] and η b → J/ψγ [5,6].Despite very clean signature thanks to the presence of J/ψ, these decay modes unfortunately have rather suppressed branching ratios.Obviously, the e + e − collision experiments like B factory are not well suited to pin down these decay modes due to the relatively low η b yield.By contrast, high energy hadron collision experiments in general possess a much larger η b production rate, which in turn allows for triggering these rare but very clean decay modes.For the decay η b → J/ψJ/ψ, the original hope is that it should be discoverable in Tevatron Run 2 [4].Subsequent investigation reveals the situation may not be so optimistic [3,7] 1 .In particular, an explicit NRQCD calculation predicts B[η b → J/ψJ/ψ] ∼ 10 −8 , which implies that notwithstanding the pessimistic prospect of observing this decay channel in Tevatron Run 2, the chance may still be bright in the forthcoming LHC program [3].
The exclusive decay η b → J/ψJ/ψ, apart from its phenomenological interest, is worthy of theoretical exploration in its own right.Since this hard exclusive process comprises entirely heavy quarkonium, it can serve as an ideal laboratory to test the applicability of NRQCD and perturbative QCD [9].Note this process is quite similar to the double charmonium production processes at B factory, e.g. e + e − → J/ψη c , J/ψJ/ψ, which have stirred quite a few attention since the first Belle measurements were released in 2002 [10] (for an incomplete list about NRQCD-based investigation, see [11,12,13,14,15,16,17,18,19,20]).No doubt a systematic study of this hadronic decay process will offer complementary insight towards understanding the decay and production mechanisms of quarkonium.
A peculiar feature of the decay η b → J/ψJ/ψ is that2 , at the lowest order in α s and charm quark relative velocity in J/ψ, v c , the amplitude vanishes.Therefore, to obtain a non-vanishing result, one must proceed to the higher orders either in α s or v c .In Ref. [3], the leading relativistic correction to the tree-level process was explored.In this work, we make progress along the alternative route, i.e., we compute the O(α s ) correction but stay at the zeroth order in v c .Since α s ∼ v 2 c , the latter contribution is expected to be as important as the former one, so it will be definitely desirable to work out this missing piece.
Before presenting our main result, it is instructive to recall some dynamical properties about this process.Let λ and λ stand for the helicities of two J/ψ mesons in the η b rest frame.
According to the helicity selection rule and angular momentum conservation, the decay configuration (λ, λ)=(0, 0) exhibits the slowest asymptotic decrease with m b , B ∼ 1/m 4 b [21].However, it is impossible for η b to decay into two longitudinally polarized J/ψ because of conflict between parity and Lorentz invariance.The hadron helicity conservation, |λ+ λ| = 0, can be violated either by the nonzero charm mass m c or by the transverse momentum of c inside J/ψ, q ⊥ .For every unit of deviation of this rule, there is a further suppression factor of m 2 c /m 2 b or q 2 ⊥ /m 2 b .For the only physically allowed helicity configurations (λ, λ)=(±1, ±1), the selection rule is violated by two units, so one expects B ∼ 1/m 8 b no matter the cause of violation is due to nonzero q ⊥ or nonzero m c .
A noteworthy finding is that, at the lowest order in α s , the violation of the rule should be ascribable to the nonzero transverse momentum of c in J/ψ, instead of the nonzero m c [3].
The corresponding branching ratio exhibits the following asymptotic scaling where the subscript emphasizes the corresponding contribution comes from the relativistic correction.The factor v 6 c affiliated with the would-be "leading-power" scaling, stems from the squares of the wave functions at the origin of two J/ψ, since ψ J/ψ (0) ∝ (m c v c ) 3/2 .In the last equation q ⊥ ∼ m c v c has been assumed.
One may naturally wonder how the situation becomes when one goes to the next-toleading order (NLO) in α s .Some typical NLO Feynman diagrams contributing to η b → J/ψJ/ψ are shown in Fig. 1.Should it also be obligatory to retain finite relative momentum of c in order to obtain a non-vanishing result?To answer this question, it is helpful to draw some clue from an analogous hadronic decay process, Υ → J/ψη c , which also violates the helicity selection rule.In this Υ decay process, the lowest-order diagrams start already at one-loop, i.e., at O(α 6 s ), which in fact share the same topology as those represented by Fig. 1e).There it is found that keeping a nonzero m c but neglecting the relative momentum of c is sufficient to give rise to a non-vanishing result [22].This may be reckoned as a persuasive hint that a similar pattern may occur in our case.A more direct evidence comes from inspecting the cut structure of Fig. 1.It is easy to observe that all the diagrams in Fig. 1 [except Fig. 1e)] permit a specific cut that has a clear physical meaning, i.e., a cut through two intermediate gluon lines that divides the full process into η b → gg and gg → J/ψJ/ψ.
Since both sub-processes are finite at the zeroth order in v c , so does the corresponding portion of the imaginary part of the η b → J/ψJ/ψ amplitude which is obtained by stitching both on-shell sub-amplitudes together in accordance with the recipe of Cutkosky rule.Since a fraction of the imaginary part of the total amplitude is nonzero, there is no reason to expect the full imaginary part should vanish, thus we are reassured that at this order, the nonzero m c can adequately play the role of violating the helicity selection rule, and it suffices to set the relative velocity v c to zero in the calculation.As a result, we expect the NLO perturbative contribution scales as follows where the subscript stresses the corresponding contribution comes from the O(α s ) perturbative correction.As will be confirmed by the explicit calculation, this power-law scaling indeed holds besides logarithmical scaling violation effect.
Having said this much, let us now proceed to the actual calculation.We first set up some notations.Let P, ε ( P , ε) signify the momentum and polarization vector of each J/ψ.
Owning to parity and Lorentz invariance, the decay amplitude M is constrained of the following structure: Apparently, the only allowed helicity configurations indeed are (λ, λ) = (±1, ±1).All the nontrivial dynamics is encoded in the reduced amplitude A, and our task is then to unravel its explicit form.
In total 80 diagrams can be drawn for this process at O(α 6 s ).A full treatment of these diagrams, at first sight, may seem to be a daunting task.Luckily, the calculation turns out to be much simpler than one would imagine, thanks to some gratifying traits of this process.First, renormalization is not needed in this work.This originates from the fact that the amplitude at leading order in α s and v c vanishes.Since the contribution of the counterterm diagrams are proportional to the LO diagram, they simply vanish too.As a consequence, all the diagrams which contain primitive ultraviolet (UV) divergent subgraph must instead be free of the UV divergences, otherwise these unremoved divergences would signal the failure of the renormalization theory.A remarkable property of this process is that, these superficially UV-divergent diagrams are not only finite, in fact all of them vanish.Similarly, among the remaining manifestly UV-finite diagrams, a simplification further occurs, namely many diagrams with a symmetric topology also vanish, e.g., diagrams involving one 4-gluon vertex, diagrams involving two 3-gluon vertices, and diagrams containing a gluon ladder between the b and b lines.This can be mainly attributed to the specific Lorentz structure as indicated in (3), that is, in these diagrams, the Levi-Civita tensor that arises from the b b annihilation amplitude, is either contracted with some symmetric Lorentz indices, or there are no sufficient number of independent 4-vectors to contract with it.
There finally survive 40 diagrams with non-vanishing contribution, some of which are shown in Fig. 1a)-e).In fact, because of the interchanging symmetry between two J/ψ, diagrams of a given topology usually give rise to an identical answer, and practically only much fewer diagrams need to be calculated.We employ the automated Feynman Diagram Calculation package (FDC) to fulfill the analytic evaluation of the required one-loop diagrams.
FDC is a powerful program based on the symbolic language Reduce, which is designed to automate the perturbative quantum field theory calculation in computer.FDC was initially developed by one of us long ago [23], and the function of automatic one-loop calculation has recently been realized by two of us [24].Recently, FDC has been successfully applied to several important quarkonium production processes to compute the NLO perturbative corrections [19,20].
In passing, we would like to comment on some specific issues in the calculation.One technical complication is the occurrences of six-point one-loop integrals in Fig. 1a) and b).
FDC has implemented some systematic recursive algorithm to reduce higher-point integrals down to lower-point ones.In this work, it turns out that all the encountered 6-point integrals can be reduced to the 4-point integrals and lower ones.In addition, Passarino-Veltman algorithm [25] has also been built into FDC to expedite the reduction of tensor integrals to the scalar ones.
It turns out that each individual diagram is infrared finite, at least in the Feynman gauge as we worked with.At a glance, one may feel that some care needs to be paid in handling Fig. 1b), since the Coulomb divergences may potentially arise there.Hearteningly, a close examination shows that this diagram is in fact free of any singular 1/v c poles.The absence of Coulomb singularity can be most transparently understood in the NRQCD factorization ansatz [9].What we can calculate within perturbation theory is the so-called hard (short-distance) coefficient.All the long-distance singularities encountered in computing the one-loop QCD diagrams, must be exactly reproduced in the respective one-loop NRQCD diagrams with a potential gluon attached between the outgoing c-c lines that forms a J/ψ, so they ultimately cancel out in the intended hard coefficient via the matching procedure.
However, the corresponding diagrams in the NRQCD sector trivially vanish in our case, again because the η b → J/ψJ/ψ process vanishes at LO in α s and v c , therefore no any net long-distance divergence, such as Coulomb divergence, should be expected to appear in the diagrams in the QCD sector.
Finally it may be worthwhile to compare our process with Υ → J/ψη c .As stated earlier, the latter process is described by a small subset of Fig. 1, the class of diagrams exemplified by Fig. 1e), in which three "Abelian" gluons connect between b quark line and the charm quark lines.Since Υ has negative C parity, charge conjugation invariance demands that three gluons must form a totally symmetric color-singlet state, therefore only the fully symmetric part of the color factor d abc d abc can contribute in this Υ decay process [22].On the contrary, due to the positive C parity of η b , in our case the three gluons must instead be totally antisymmetric in color space, so the f abc f abc piece finally survive in the color factor.We have explicitly checked this is indeed the case.
After analytic manipulation of all the diagrams by FDC, together with some simplifying by hand, we can express the reduced amplitude as following: where the relativistic normalization is used for the quarkonium states, and is the magnitude of J/ψ momentum.The complex-valued f function, which has encompassed the loop contribution, is deliberately chosen as dimensionless, therefore it can depend on m b and m c only through their dimensionless ratio m 2 c /m 2 b .Furthermore, the f function is normalized such that it varies with m 2 c /m 2 b only logarithmically, which is slower than any power law scaling.After some quick algebra, one can verify that Eq. ( 4) is in conformity to B αs ∼ 1/m 8 b up to logarithmical modification.We display the shape of the f function in Figure 2. The analytic expression for the real part of the f function is too lengthy to be reproduced here.However, the analytical result for the imaginary part is quite compact, where It is instructive to deduce the asymptotic behavior of f in the small ξ limit.After some As can be clearly visualized in Fig. 2, as ξ → 0, the real and imaginary parts of the f function exhibit ln 2 ξ and ln ξ dependence, respectively.This seems to be a generic trait of exclusive bottomonium decays to double charmonium, which is also seen in the analogous f function in the Υ → J/ψη c process [22].It is interesting to note that, up to the physically relevant point ξ = m 2 c /m 2 b = 0.102, which is obtained by substituting m c = 1.5 GeV and m b = 4.7 GeV, the asymptotic expressions in ( 6) and ( 7) still constitute decent approximations to the true expressions, where the agreement for the real part is especially satisfactory, and the agreement for the imaginary part is at 70% level.This may suggest that, if one is willing to carry out this calculation by hand yet without loss of accuracy, the most efficient way is to first expand the integrand in m 2 c /m 2 b prior to performing the loop integrals.It is interesting to note that, provided that ξ is not overly small, say, ξ > 3.8 × 10 −4 , then −Im f is always bigger than |Re f |.In particular, at the phenomenologically relevant point ξ = 0.102, we find the corresponding f = 0.216 − 2.30i = 2.31 e −i 84.6 • , and the real part seems to be practically negligible.
Since the NLO perturbative amplitude is dominated by its imaginary part, one may like to speculate which type of cuts might constitute the most important contribution.It might sound attractive to inquire whether that specific cut as mentioned earlier, the one acting upon two intermediate gluon lines that partitions the full process into η b → gg and gg → J/ψJ/ψ, would dominate over all other possible cuts.If this were the case, we would be happy because this process then could be endowed with a "plausible" physical picture -that it proceeds first through η b annihilation into two gluons, and these two gluons subsequently rescatter into double J/ψ.Let us examine this conjecture concretely.With the aid of the Cutkosky rule, it is straightforward to deduce this particular fraction of the imaginary part of the f function, which we call Im f 2g : Note the occurrence of the 1/ξ term, which can be traced to the cut contribution from Fig. 1a), renders this function scale with quark masses as m 2 b /m 2 c , so that the true scaling B αs ∼ 1/m 8 b is violated.This is unacceptable and clearly indicates approximating the full imaginary part by this specific fraction is unjustified.Numerically for ξ = 0.1, we obtain Im f 2g = 1.15, which is far off the true value Im f = −2.30.A careful investigation reveals that each diagram in Fig. 1 [except Fig. 1e)] possesses a rich cut structure, for example, Fig. 1a), b), d) allow for 6, 5, and 4 different kind of cuts, respectively.Among all the possible cuts, those which pass through one gluon line and another nonadjacent quark line generally lead to more involved expressions in (5).For a given diagram, it seems not profitable to single out one specific cut to mimic its full imaginary part.For example, only after the remaining cuts are included, the imaginary part of Fig. 1a) then has a sensible answer, free of the dangerous 1/ξ singularity.
We are now very close to the intended formula.For completeness, we still need the expression of the reduced amplitude that arises from the tree-level relativistic correction [3] Roughly speaking, v 2 J/ψ is a quantity related to the second derivative of the wave function at the origin for J/ψ, which governs the size of relativistic corrections.In fact, this factor admits a rigorous definition as a ratio of NRQCD matrix elements [9,11]: where ǫ denotes the polarization vector of J/ψ at its rest frame, ψ and χ represent Pauli spinor fields in NRQCD, D is the spatial covariant derivative.The matrix elements appearing in the above ratio are understood to be the renormalized ones.It is worth stressing that the sign of v 2 J/ψ needs not necessarily to be positive.Plugging equations ( 4) and (9) in the formula we then obtain the desired partial width.Note the cubic power of momentum reflects that two J/ψ are in relative P -wave state.This formula has already taken into account the sum over transverse polarizations of two J/ψ as well as their indistinguishability.It may be more convenient to present an explicit expression for the branching ratio, where the wave function at the origin of η b drops out, In deriving this, we have approximated the total width of η b by its hadronic width: where the factor K gg encodes the magnitude of NLO perturbative correction to η b → gg [26] n f denotes the number of active light flavors.One usually takes n f = 4, by treating charm quark also as a massless parton.For α s (2m b ) = 0.18, we get K gg = 1.58.
Equation ( 12) constitutes the key formula of this work.From this equation one can readily confirm the asymptotic behaviors anticipated in (1) and ( 2).
Let us now explore the numerical outcome of (12).The input parameters are m b , m c , α s , ψ J/ψ (0) and v 2 J/ψ , respectively.We take m b = M η b /2 ≈ 4.7 GeV and m c = 1.5 GeV.Obviously, our predictions are quite sensitive to the value assigned for the strong coupling constant.To account for the affiliated ambiguity, we slide the renormalization scale µ from 2m b to 2m c , corresponding to varying α s between 0.18 and 0.26.The wave function at the origin of J/ψ can be extracted from its electric width: where the first-order perturbative correction has been included.Using the measured electric width 5.55 keV [27], we obtain ψ J/ψ (0) = 0.263 GeV 3/2 .
Among various input parameters, the least known about is v 2 J/ψ .Nevertheless, there is a useful relation, first derived by Gremm and Kapustin by employing the equation of motion of NRQCD [28], which expresses this quantity in terms of the charm quark pole mass and the J/ψ mass.To our purpose this relation reads [11] 1.5×10 -7 2.0×10 Owing to the intrinsic ambiguity associated with the charm quark pole mass, a wide spectrum of estimates for m c pole is scattered in literature.Ref. [11] takes m c pole = 1.4 GeV in the analysis of double charmonium production at B factory, which leads to v 2 J/ψ = 0.22 from (15).Note this value is compatible with v 2 J/ψ = 0.225 +0.106 −0.088 , which is obtained from a recent Cornell potential-model based analysis [29,30].A positive v 2 J/ψ is helpful to alleviate the discrepancy between the LO NRQCD prediction and the measured J/ψ + η c production rate at Belle.Notwithstanding the phenomenological inclination, however, one may be alert that a rather different value might not be excluded on theoretical ground.For example, if m c pole = 1.75 ± 0.15 GeV is assumed, which stems from a recent QCD moment sum rule analysis [31], one would instead obtain v 2 J/ψ = −0.22 ± 0.15, whose sign has can be entirely attributed to the NLO perturbative contribution.In contrast, the prediction including only the relativistic correction simply vanishes as v 2 J/ψ approaches zero.Finally we present an updated estimate about the discovery potential of this decay mode in the LHC experiments.In the hadron collision experiment J/ψ can be cleanly reconstructed via its muonic decay mode.Multiplying ( 16) by the branching ratios of 6% for each of the decay J/ψ → µ + µ − , we obtain B[η b → J/ψJ/ψ → 4µ] ≈ (0.7 − 6.7) × 10 −10 .
Assuming the η b production cross section at LHC to be about 15 µb [3], one then finds the cross section for the 4 µ events to be about 1-10 fb.For a 300 fb −1 data, which is expected to be gleaned in one year run at LHC design luminosity, the number of produced events may reach 300-3000.The product of acceptance and efficiency for detecting J/ψ decay to µ + µ − is estimated to be ǫ ≈ 0.1 [4], which is perhaps a conservative estimate for LHC.Multiplying the number of the produced events by ǫ 2 , we expect between 3 and 30 observed events per year.From this rough estimate, we are tempted to conclude that, the prospect of observing η b at LHC through the 4µ mode looks promising.However, we should be fully aware that it will be a highly challenging task for experimentalists to single out the relatively few signal events from the presumedly abundant continuum background.
To summarize, in this work we have studied the exclusive decay process η b → J/ψJ/ψ to NLO in α s while at LO in v c .We found that this new contribution to the amplitude is comparable in magnitude with the previously calculated tree-level relativistic correction piece, but differs by a phase about 90 • .Including this new piece of contribution will substantially enhance the previous estimate of B(η b → J/ψJ/ψ), thus the observation prospect of this clean η b hadronic decay mode in the forthcoming LHC experiment becomes brighter.

FIG. 1 :
FIG. 1: Lowest order diagram, together with some representative next-to-leading order diagrams [from a) to e)] that contribute to η b → J/ψ J/ψ.

FIG. 2 :
FIG.2:The solid lines represent the real and imaginary parts of f (ξ), the dashed lines represent their respective asymptotic expressions as given in (6) and(7).The vertical line marks the place ξ = m 2 c /m 2 b = 0.102, which is adopted in the phenomenological analysis.

FIG. 3 :
FIG. 3: B[η b → J/ψ J/ψ] as a function of v 2 J/ψ (left panel) and α s (right panel).The solid lines represent the full prediction from (12) and the dashed lines are obtained by keeping only the relativistic correction piece.In the left panel, three different pairs of curves (each pair tinted by the same color) are obtained by fixing the renormalization scale µ to be 2m b , m b and 2m c respectively, so that the corresponding strong coupling constant α s is 0.18, 0.22 and 0.26.