application/xmlImplications for new physics from [formula omitted] and [formula omitted]Jian-Feng ChengYuan-Ning GaoChao-Shang HuangXiao-Hong WuPhysics Letters B 637 (2006) 260-265. doi:10.1016/j.physletb.2006.04.052journalPhysics Letters BCopyright © 2006 Elsevier B.V. All rights reserved.Elsevier B.V.0370-26936374-515 June 20062006-06-15260-26526026510.1016/j.physletb.2006.04.052http://dx.doi.org/10.1016/j.physletb.2006.04.052doi:10.1016/j.physletb.2006.04.052http://vtw.elsevier.com/data/voc/oa/OpenAccessStatus#Full2014-01-01T00:14:32ZSCOAP3 - Sponsoring Consortium for Open Access Publishing in Particle Physicshttp://vtw.elsevier.com/data/voc/oa/SponsorType#FundingBodyhttp://creativecommons.org/licenses/by/3.0/JournalsS300.2PLB23004S0370-2693(06)00531-410.1016/j.physletb.2006.04.052Elsevier B.V.PhenomenologyFig. 1The constraints on z from the branching ratios of B¯0→π0π0 and B¯0→K¯0K0. The big lattice denotes the constraint from B¯0→π0π0 and the small lattice for the constraint from B¯0→K¯0K0.Fig. 2The correlation between |C8g(mW)−C8g′(mW)| and |C13(mW)−C13′(mW)|.Fig. 3The correlation between B¯0→π0π0 and B¯0→K¯K in MSSM with the LL insertion.Implications for new physics from B¯0→π0π0 and B¯0→K¯0K0Jian-FengChenga⁎chengjf@mail.tsinghua.edu.cnYuan-NingGaoaChao-ShangHuangbXiao-HongWucaCenter for High Energy Physics, Tsinghua University, Beijing 100084, ChinabInstitute of Theoretical Physics, Academia Sinica, P.O. Box 2735, Beijing 100080, ChinacSchool of Physics, Korea Institute for Advanced Study, Seoul 130-722, South Korea⁎Corresponding author.Editor: M. CvetičAbstractWe have analyzed the B¯0→π0π0 puzzle in three kinds of models beyond the standard model (SM). It is shown that the minimal flavor violation (MFV) models, the minimal supersymmetric standard model (MSSM), and the two Higgs doublet models (2HDM) I and II cannot give an explanation of the B¯0→π0π0 puzzle within 1σ experimental bounds and the model III 2HDM can explain the puzzle without a conflict with other experimental measurements. If the constraint on C8g from b→sg is not imposed, for all kinds of insertions considered there are regions of parameter space, where the scalar quark mass is larger (much larger) than the gluino mass in the case of LR or RL (LL or RR), in which the puzzle can be resolved within 1σ experimental bounds.The branching ratios of B¯0 decays into two pions have been recently observed [1]:(1)Br(B¯0→π0π0)=(1.45±0.29)×10−6,Br(B¯0→π+π−)=(5.0±0.4)×10−6. The large branching ratio of the B decay into neutral pion final states is unexpected. The decay amplitudes of B¯0→ππ can be generally parameterized as(2)2A(B¯0→π0π0)=T[(PT−PEWT)eiα−CT],A(B¯0→π+π−)=−T[1+PTeiα], where T, C, P, and PEW are the tree, color-suppressed tree, penguin, and electro-weak penguin (EWP) amplitudes respectively, and α=arg(−λuλt) is the weak phase, where λp=VpbVpd∗(p=u,c,t). In SM one has the counting rules: the color-suppressed tree and penguin amplitudes are suppressed by a factor of λ (λ∼0.22 is the Wolfenstein parameter) and the EW penguin is suppressed by a factor of λ2, with respect to the tree amplitude [2]. So one should expect by the naive counting rules that the branching ratio of the B decay into neutral pion final states is O(λ2) of that for charged pion final states. However, the data (Eq. (1)) indict that the former is O(λ) of the later. The observed branching ratio of B¯0→π0π0=(1.45±0.29)×10−6 is much larger than the theoretical prediction, about 0.3×10−6, up to the αs order in the BBNS approach (QCDF) [3,4] in SM. In Li et al.'s approach (PQCD) the leading order (LO) prediction ∼10−7[5] is the same order as that of the QCDF prediction. In the recent paper [6] the next leading order (NLO) PQCD calculations have been carried out and the results are that the πK puzzle, the expected relation ACP(B±→π0K±)≈ACP(B0→π±K∓) disagreed significantly with the data, can be resolved but the predicted branching ratio of B¯0→π0π0 is about 0.3×10−6 which is still much smaller than the data, i.e., the π0π0 puzzle remains. If the large branching ratio persists it could indicate new physics.Though B¯0→π0π0 is not a pure-penguin process and has the contributions from tree operators, the tree contributions are of the order same as the penguin contributions because of the almost completely cancellation between the two terms in C2+C1/Nc where C1,2 are Wilson coefficients of tree operators, so B¯0→π0π0 is sensitive to new physics. Therefore, it seems that a lot of new models beyond SM could enhance the branching ratio and consequently resolve the puzzle [7]. However, any new model must simultaneously give an explanation for the branching ratio of B¯0→K¯0K0 since the two processes are closely related at quark level: the flavor changing neutral current b→d transition controls B¯0→K¯0K0 and the same transition gives significant contributions to B¯0→π0π0 which are of the order same as the tree contributions in SM. Recently the branching ratio of B¯0→K¯0K0 has been measured as (0.96±0.25)×10−6[1], which is consistent with the prediction from both the QCDF [4] and PQCD approaches [8]. Therefore, new physics (NP) contributions must satisfy the condition that they make Br of B¯0→π0π0 enhanced but keep Br of B¯0→K¯0K0 basically unchanged, compared with those in SM respectively, which will impose the significant limit on NP models.In the Letter we search for new models beyond the SM which can account for the data of branching ratios for both the B¯0→π0π0 and B¯0→K¯0K0 processes. To be specific, we concentrate on the well-known three kinds of models: the minimal flavor violation (MFV) models, the two Higgs doublet models (2HDM) and the minimal supersymmetric standard model (MSSM).The effective Hamiltonian relevant for the two processes in the SM can be expressed as [3](3)HeffSM=GF2∑p=u,cλp(C1Q1p+C2Q2p+∑i=3,…,10CiQi+C7γQ7γ+C8gQ8g)+h.c., where λp=VpbVpd∗, Q1,2 and Qi(i=3,…,10) are the tree and penguin operators respectively. Explicit forms for C1, C2, Ci, C7γ, C8g and Q1p, Q2p, Qi, Q7γ, Q8g can be found, e.g., in Ref. [3].In the QCD factorization approach, the dominant contributions to the decay amplitudes are given by:(4)M(B¯0→π0π0)=GF2fπFB→πmB2×12∑p=u,c[a2λu−(a4p+rπa6p)λp],M(B¯0→K¯0K0)=GF2fKFB→KmB2×∑p=u,c[(a4p+rKa6p)λp], where the definitions of the parameters ai and the chiral enhancement factors rπ, rK can be found in Ref. [3]. We take the values of running masses in the MS¯ scheme for light quarks such that rπ=rK≡r hereafter. The electro-weak penguin and annihilation contributions are neglected in above formula, which leads to the 10% theoretical uncertainty.First we consider the MFV models. The MFV models beyond the SM discussed in the Letter mean a class of models in which the general structure of flavor changing neutral current (FCNC) processes present in the SM is preserved. In particular, all flavor violating and CP-violating transitions are governed by the CKM matrix and the only relevant local operators are the ones that are relevant in the SM [9]. New parameters in the MFV models, e.g., the masses of charginos, squarks, Higgs particles in the MFV scenario of the MSSM, enter into Wilson coefficients of relevant local operators. Therefore, in the MFV models the amplitudes of the two decays are given same as Eq. (4) (with values of ai generally different from those in SM).We can model-independently determine |z|≡|∑p=u,c[(a4p+ra6p)λp]| from the measured branching ratios of B¯0→π0π0 and B¯0→K¯0K0 and a2=(0.0502−0.0689i)±(0.0025+0.0035i) which comes from the known tree contributions. The result is given in Fig. 1 where the 60% theoretical uncertainty (coming mainly from non-perturbative parameters such as form factors, distribution amplitudes and CKM matrix elements) has been taken into account. There is a narrow region which can simultaneously fit the data. However, assuming Wilson coefficients of relevant operators, except the chromo-magnetic dipole operator, change a little compared with SM, which is the case in the MFV models, z in the region corresponds to(5)|C8g(mW)|⩾2.6, which cannot be reached in the MFV models [10,11]. That is, the MFV models are excluded within 1σ experimental bounds.Next we consider models in which there are new operators in addition to those in the SM, e.g., the 2HDM and MSSM. The effective Hamiltonian in the 2HDM and MSSM can be written as [12,13](6)Heff=HeffSM+Heffnew,(7)Heffnew=GF2∑p=u,cVpbVpd∗(∑i=11,…,16[CiQi+Ci′Qi′]+∑i=3,…,10Ci′Qi′+C7γ′Q7γ′+C8g′Q8g′)+h.c., where Qi(′), i=11,…,16, are the neutral Higgs penguin operators and their explicit forms can be found in Refs. [12,13] with the substitution s→d. The primed operators, the counterpart of the unprimed operators, are obtained by replacing the chirality in the corresponding unprimed operators with opposite ones.From the effective Hamiltonian, Eq. (6), it follows that the main contributions to the decay amplitudes from the SM and new physics are given by:(8)M(B¯0→π0π0)=GF2fπFB→πmB2×12∑p=u,c[a2λu−(a4p+ra6p)λp],M(B¯0→K¯0K0)=GF2fKFB→KmB2×∑p=u,c{(a4p+ra6p)λp+msmb[h1(C11−C11′)+h2(C13(μ)−C13′(μ))]}, where we have set md=0. Due to the renormalization group equation (RGE) running, Wilson coefficients Ci, i=14,15,16, are related to C13 and the known constants h1,2 represent the running effects. The largest contributions to the hadronic elements of the neutral Higgs penguin operators at the αs order arise from penguin contractions with b quark in the loop, which are the same for B¯0→π0π0 and B¯0→K¯0K0 and have been included in a4 (see, for example, Ref. [14]). Therefore, although they can enhance the branching ratio of B¯0→π0π0, they alone cannot resolve the puzzle because the branching ratio of B¯0→K¯0K0 will also be enhanced by them, which will not agree with the data.The new physics contribution, the terms proportional to C11 and C13, respectively, which contributes to one mode but not to the another (precisely speaking, the contribution to the another mode is md/ms suppressed) comes from the hadronic matrix elements of Higgs penguin operators at the leading order in αs. It is the contribution that gives a possibility to make the data be account for without a conflict with all relevant experimental measurements. The key point is if one can have a sizable C13(′)(μ) and/or C11(′)(μ) in the 2HDM and MSSM.Let us analyze how large C13(′) and/or C11(′) are needed to fit the data. Letz=∑p=u,cλpλt(a4p+ra6p),z1=λuλta2,z2=msmb[h1(C11−C11′)+h2(C13(μ)−C13′(μ))], we have(9)r11≡232πmBBr(π0π0)min(GFfπFB→πλt)2τB⩽|z−z1|⩽r12≡232πmBBr(π0π0)max(GFfπFB→πλt)2τB,r21≡32πmBBr(K¯0K0)min(GFfKFB→Kλt)2τB⩽|z+z2|⩽r22≡32πmBBr(K¯0K0)max(GFfKFB→Kλt)2τB. From the data, 1.16×10−6⩽Br(π0π0)⩽1.74×10−6 and 0.71⩽Br(K¯0K0)⩽1.21×10−6[15], we have r12>r11>r22>r21. To satisfy the above two relations, it is necessary to have(10)|z2+z1|⩾r11−r22.In the model I and II 2HDMs and MSSM the Wilson coefficients of QCD penguin operators are not changed significantly, compared with those in SM, and the Wilson coefficient of chromo-magnetic operator can have a significant change [16]. Taking the SM values of Wilson coefficients of relevant operators but the chromo-magnetic operator and using RQE running, we can obtain the correlation between |C8g(mW)−C8g′(mW)| and |C13(mW)−C13′(mW)| from |z−z1|⩾r11 and Eq. (10), which is shown in Fig. 2 where C11=C13 has been assumed for simplicity, without losing the generality of discussions.11In the figure |C8g(mW)−C8g′(mW)| has no upper bound because we do not impose |z−z1|⩽r12. In all models beyond SM known so far, |C8g(mW)−C8g′(mW)| never reach very large value (say, 5) when all relevant experimental constraints are imposed. We do not need to know the upper bound for the analysis in the Letter. It follows from the figure that |C8g(mW)−C8g′(mW)|min=2.6 when C13(mW)−C13′(mW)=0, which reduce to Eq. (5) in the MFV models, as it should be.It is well known that the experimental upper bound of branching ratio for Bs→μ+μ− constrains severely parameters in the MSSM and model I and II 2HDMs [17]. Similarly, we show that the corresponding bound for Bd→μ+μ− implies that the Wilson coefficients of new operators in the MSSM and model I and II 2HDMs cannot be large. The branching ratio Bd→μ+μ− in the 2HDM and MSSM is given as(11)Br(Bd→μ+μ−)=GF2αem264π3mBd3τBdfBd2|λt|21−4mˆ2[(1−4mˆ2)|CQ1(mb)−CQ1′(mb)|2+|CQ2(mb)−CQ2′(mb)+2mˆ(C10(mb)−C10′(mb))|2], where mˆ=mμ/mBd. In the moderate and large tanβ cases the term proportional to (C10−C10′) in Eq. (11) can be neglected. The new CDF and D0 combined experimental upper bound of Br(Bd→μ+μ−) is 3.2×10−8[18] at 90% confidence level. We have the constraint(12)|CQ1(mW)−CQ1′(mW)|2+|CQ2(mW)−CQ2′(mW)|2≲2.2, where CQ1,2(′) are the Wilson coefficients of the operators Q1,2(′) which are Higgs penguin induced in leptonic and semileptonic B decays and their definition can be found in Refs. [19,20]. By substituting the quark-Higgs vertex for the lepton-Higgs vertex, it is straightforward to obtain Wilson coefficients relevant to hadronic B decays in the MSSM and model I and II 2HDMs. To translate CQ1,2 into CQ11,13, we have CQ11,13(′)(mW)∼0.037. Then it follows from Fig. 2 that |C8g(mW)−C8g′(mW)| must be larger than 2.4 in order to resolve the puzzle.The Wilson coefficients C8g(′) in the b→d transition are constrained by Br(B→Xdg). Because there is no data for Br of the B→Xdg decay and the difference between the B→Xdg and B→Xsg decays in the SM comes from CKM matrix elements, we assume the constraint on C8g(′) same as that from b→sg. In the presence of new physics a model-independent analysis gives that |C8g(mW)−C8g′(mW)|<2.01 when Br(b→sg)<9%[21]. That is, |C8g(mW)−C8g′(mW)| cannot satisfy the condition, larger than 2.4, because of the b→sg constraint. Therefore, we come to the conclusion that the puzzle of B¯0→π0π0 cannot get resolved within 1σ experimental bounds in the MSSM and model I and II 2HDMs.If one does not impose the b→sg constraint it is possible to resolve the puzzle in the MSSM because the Wilson coefficient C8g(∗) can reach values larger than 2.6 in some regions of parameter space. We have carried out detailed numerical calculations in the MSSM, imposing the constraints from the B¯d0–B¯d0 system, the mass difference ΔMd=(0.509±0.004)ps−1, mixing induced CP-violation phase angle β measured in charmonium B decays, sin2β=0.687±0.032[1], and B¯0→Xdγ, in addition to the constraint from Bd→μ+μ−, however, without imposing the b→sg constraint. δ13LL,RR and δ13LR,RL are constrained to be order of 10−1 and 10−2 respectively with moderate sparticle masses (say, 500 GeV) [22,23]. In particular, Br(B¯0→Xdγ)⩽1×10−5 extracted from exclusive B→ρ(ω)γ, as advocated in Ref. [23], gives a more stringent constraint. Br(B¯0→Xdγ) directly constrains |C7γ(mb)|2+|C7γ′(mb)|2 at the leading order. Due to the strong enhancement factor mg˜/mb associated with single δ13LR(RL) insertion term in C7γ(′)(mb), δ13LR(RL) (∼10−2) are more severely constrained than δ13LL(RR). However, if the left–right mixing of scalar bottom quark δ33LR is large (∼0.5), δ13LL(RR) is constrained to be order of 10−2 since the double insertion term δ13LL(RR)δ33LR(LR∗) is also enhanced by mg˜/mb. In the case of LL (RR) insertion the precisely measured mass difference ΔMd imposes a more severe constraint on δ13LL(RR).In numerical analysis we fix tanβ=10, vary mg˜ and mq˜ in the region between 300 GeV and 2 TeV, and the NHB masses in the ranges of 91GeV⩽mh⩽135GeV, 91GeV⩽mH⩽200GeV with mh<mH and 200GeV⩽mA⩽240GeV for the fixed mixing angle π/2 of the CP even NHBs, and scan δ13dAB in the range |δ13dAB|⩽0.1 for A=B and 0.05 for A≠B (A=L,R). The numerical result for the correlation between branching ratios of B¯0→π0π0 and B¯0→K¯0K0 is shown in Fig. 3 for the case of LL insertion. Due to the combined constraints mentioned above, in most of parameter space C8g(mW) is not large enough to accommodate the data in 1σ region. However, there are small regions of parameter space with x≫1 (say x∼40, x is the square of the ratio between scalar quark and gluino masses) where C8g(mW) is large enough to resolve the puzzle. In the case of both LL and RR insertion, the result is similar. In the cases of LR, both LR and RL insertions, we also have similar results. For x>1 we have the regions with large enough C8g(mW) and the regions are larger than those in the cases of LL, both LL and RR insertions. In short, numerical results of Br show that the MSSM can explain the puzzle within 1σ experimental bounds under all relevant experimental constraints except that from Br(b→sg).Finally we consider the model III 2HDM [24]. In the model III 2HDM there are tree-level flavor changing neutral currents (FCNC). After diagonalizing the mass matrix of quark fields, the flavor changing (FC) part of the Yukawa Lagrangian is [24](13)LY,FC(III)=ξijUQ¯i,LH˜2Uj,R+ξijDQ¯i,LH2Dj,R+h.c. In order to obtain naturally small FCNC one assumes the Cheng–Sher parameterization [25](14)ξijD=λijmimjv.Phenomenological constraints on parameters of the models have been extensively discussed [26]. For b→dss¯, the couplings λbd,db,ss are involved. λbd can reach 0.4 without a conflict with the measured mass difference ΔMBd if the mass of pseudo-scalar Higgs boson MA is large (say, ∼1TeV) [27] and it is also allowed by the recent data for B¯0→ρ(ω)γ. The constraint on λss from the analysis on B¯s0–B¯s0 shows that λss can reach O(100)[28] which means that the coupling of Higgs to s quark is O(10−2). It should be emphasized that the constraint from Bd→μ+μ− is irrelevant in the model III 2HDM because the decay involves λμμ besides λbd and λμμ has no relation to λss, which is different from the MSSM and 2HDMs I and II.In numerical calculations, we use mh=120GeV, md=6MeV, λbd=0.3, λss=150 and consequently obtain C13(mW)=0.41. Corresponding to this value, |C8g(mW)|⩾0.6, which can be satisfied under all relevant constraints. The numerical result for the values of parameters given above shows thatBr(B¯0→π0π0)=1.3×10−6,Br(B¯0→K¯0K0)=0.9×10−6. That is, the data of branching ratios of B¯0→π0π0 and B¯0→K¯0K0 are accounted for, as expected. At the same time, we have checked that the NP contribution to Br(B¯0→π+π−) is very small and negligible.In conclusion, we have analyzed the B¯0→π0π0 puzzle in three kinds of models beyond the SM. In the analysis 1σ experimental bounds and 60% theoretical uncertainty which mainly comes from the input of non-perturbative parameters have been taken into account. It is shown that the minimal flavor violation models, the minimal supersymmetric standard model, and the two Higgs doublet models I and II cannot give an explanation of the B¯0→π0π0 puzzle within 1σ experimental bounds when all relevant experimental constraints are imposed and the model III 2HDM can explain the puzzle without a conflict with other experimental measurements. Therefore, if the data of Br for B¯0→π0π0 and B¯0→K¯0K0 persist in the future, the MFV models, MSSM and model I, II 2HDMs will be excluded within 1σ experimental bounds and the model III 2HDM will be survived to resolve the puzzle. As it is obvious, the analysis depends on the estimate of theoretical uncertainty. Our estimate comes from the uncertainties of input parameters (form factors, distribution amplitudes, CKM matrix elements, etc.) as well as the error estimate of neglecting electro-weak penguin and annihilation contributions. If the uncertainty were 70%, |C8g(mW)−C8g′(mW)|min, which should be reached in a model in order to account for the data, would be 2.2 in the MFV models and 2.0 in the MSSM respectively so that the MFV models could not give an explanation of the data and the MSSM could. We have also analyzed the puzzle in the case without imposing the b→sg constraint in the MSSM. For all kinds of insertions there are regions of parameter space where the puzzle can be resolved within 1σ experimental bounds. We expect that similar effects appear in decays with PV final states, B→π0ρ0 and B→K¯0K0∗.AcknowledgementsThe work was supported in part by the Natural Science Foundation of China (NSFC), grant 10435040, grant 90503002, and grant 10225522.References[1]The Heavy Flavor Averaging Grouphttp://www.slac.stanford.edu/xorg/hfag/[2]M.GronauO.F.HernándezD.LondonJ.L.RosnerPhys. Rev. D5219956356[3]M.BenekeG.BuchallaM.NeubertC.T.SachrajdaNucl. Phys. B6062001245[4]M.BenekeM.NeubertNucl. Phys. B6752003333hep-ph/0308039[5]C.D.LüK.UkaiM.Z.YangPhys. Rev. D632001074009Y.Y.KeumA.I.SandaPhys. Rev. D672003054009[6]H.-n.LiS.MishimaA.I.Sandahep-ph/0508041[7]Y.-D.YangR.-M.WangG.R.LuPhys. Rev. D732006015003hep-ph/0509273S.BaekF.J.BotellaD.LondonJ.P.SilvaPhys. Rev. D722005114007[8]C.H.ChenH.N.LiPhys. Rev. D632001014003[9]A.J.BurasP.GambinoM.GorbahnS.JagerL.SilvestriniPhys. Lett. B5002001161hep-ph/0007085A.J.BurasActa Phys. Pol. B3420035615hep-ph/0310208[10]A.AliE.LunghiC.GreubG.HillerPhys. Rev. D662002034002[11]G.DegrassiP.GambinoP.Slavichhep-ph/0601135[12]J.-F.ChengC.-S.HuangPhys. Lett. B5542003155[13]J.-F.ChengC.-S.HuangX.-H.WuPhys. Lett. B5852004287hep-ph/0306086[14]J.-F.ChengC.-S.HuangX.-H.WuNucl. Phys. B701200454[15]C.W.BauerI.Z.RothsteinI.W.Stewarthep-ph/0510241[16]See, e.g.C.-S.HuangP.KoX.-H.WuY.-D.Yanghep-ph/0511129and references therein[17]C.-S.HuangQ.-S.YanPhys. Lett. B4421998209C.-S.HuangW.LiaoQ.-S.YanPhys. Rev. D591999011701S.R.ChoudhuryN.GaurPhys. Lett. B451199986K.S.BabuC.KoldaPhys. Rev. Lett.842000228C.-S.HuangPhys. Rev. D632001114021C.-S.HuangPhys. Rev. D642001059902ErratumP.H.ChankowskiL.SlawianowskaPhys. Rev. D632001054012C.BobethPhys. Rev. D642001074014hep-ph/0204225For a recent review, see, e.g.C.Koldain: Proceedings of the 12th International Conference on Supersymmetry and Unification of Fundamental Interactions, Tsukuba, Japan, June 2004hep-ph/0409205[18]R.BernhardCDF CollaborationD0 Collaborationhep-ex/0508058[19]Y.-B.DaiC.-S.HuangH.-W.HuangPhys. Lett. B3901997257[20]C.-S.HuangX.-H.WuNucl. Phys. B6572003304[21]C.GreubP.LinigerPhys. Lett. B4942000237C.GreubP.LinigerPhys. Rev. D632001054025G.HillerF.KrügerPhys. Rev. D692004074020[22]D.BecirevicNucl. Phys. B6342002105hep-ph/0112303[23]P.KoJ.-h.ParkG.KramerEur. Phys. J. C252002615hep-ph/0206297[24]W.S.HouPhys. Lett. B2961992179A.AntaramianL.J.HallA.RasinPhys. Rev. Lett.6919921871L.J.HallS.WeinbergPhys. Rev. D481993R979M.J.SavagePhys. Lett. B2661991135L.WolfensteinY.L.WuPhys. Rev. Lett.7319942809[25]T.P.ChengM.SherPhys. Rev. D3519873484T.P.ChengM.SherPhys. Rev. D4419911461[26]D.B.ChaoK.CheungW.-Y.KeungPhys. Rev. D591999115006E.O.IltanG.TuranI.TuranJ. Phys. G282002307Z.XiaoL.GuoPhys. Rev. D692004014002Y.-F.ZhouJ. Phys. G302004783R.A.DiazR.MartinezC.E.SandovalEur. Phys. J. C412005305[27]D.AtwoodL.ReinaA.SoniPhys. Rev. D5419963296[28]C.-S.HuangT.J.LiInt. J. Mod. Phys. A202005161and references therein