A Free-Quark Amplitude

# Modeling Charmonium-η Decays of Jpc=1−− Higher Charmonia

## Abstract

We propose a new model to create a light meson in the heavy quarkonium transition, which is inspired by the NambuJona-Lasinio (NJL) model. Hadronic transitions of higher charmonia with the emission of an meson are studied in the framework of the proposed model. The model shows its potential to reproduce the observed decay widths and make predictions for the unobserved channels. We present our predictions for the decay width of and , where are higher and wave vector charmonia, which provide useful references to search for higher charmonia and determine their properties in forthcoming experiments. The predicted branching fraction is one order of magnitude smaller than the channel. Estimates of partial decay width are given for , and by assuming them as bound states with quantum numbers , and , respectively. Our results are in favor of these assignments for and . The corresponding experimental data for these states has large statistical errors which do not provide any constraint on the mixing angle if we introduce mixing. To identify , precise measurements on its hadronic branching fraction are required which are eagerly awaited from BESIII.

## 1 Introduction

Quantum chromodynamics (QCD), the gauge theory of strong interaction, received huge devolvement during the last few decades. However, it is still a subject of intensive research of various theoretical constructs (for instance, see recent review [1]). Study of heavy quarkonium decays is a good probe to understand the nonperturbative nature of QCD at different energy regimes. Thanks to the wealth of experimental facilities like CLEO, Belle, BABAR, CDF, D, and BESIII, now we have intensive experimental data in the charmonium () energy regime. This provides us with great opportunities to test and explore the nature of strong interactions in the heavy quark sector. Currently BESIII is taking data in the energy regime and it is easy to produce higher charmonia through annihilations. Figure 1 is the sketch of the intermediate production of vector charmonia, which further decay into . Since the center-of-mass energy of BESIII can go up to GeV [2], which is around the mass region of and , it is a good opportunity to study the production and decay mechanism of higher vector charmonia. In the future ANDA also plans to collect data [3] in the energy regime which colliders are not capable of producing directly. These experimental facilities will surely help us to deepen our understanding of heavy quarkonium physics and hence the nonperturbative aspects of strong interaction.

Heavy hadron spectroscopy has celebrated almost four decades since the discovery of in 1974. During this era, many theoretical studies have been carried out in the quark model framework to produce the spectrum of heavy quarkonium systems [4, 5]. An important manifestation of studying heavy is that the spectrum of these states can be explained by using nonrelativistic formalism. For instance, the Cornell potential model [6, 7], which incorporates a spin-independent color Coulomb plus scalar linear confined potential, was hugely successful at describing the spectrum of the charmonium systems. The Cornell spin-independent potential is an approximate heavy quark spin symmetry (HQSS) within systems. Deviations with experiments can be observed in such potential models due to HQSS breaking effects [8, 9]. One possible source of breaking HQSS is the spin-dependent potential which introduces relativistic corrections to the Cornell potential model. The widely used relativized quark potential model, sometimes also referred as the Godfrey-Isgur (GI) model [10], is so far considered the best available option to reproduce the spectra of heavy quarkonium systems.

In heavy systems, hadronic transitions serve as a crucial probe of their internal structures and help to establish the understanding of light quark coupling with a heavy degree of freedom. In QCD, the well-established formalism for hadronic transitions is multipole expansion (ME) [11, 12, 13, 14], which assumes that these transitions take place due to the intermediate process of gluon emission. These gluons are supposed to be soft, having wavelengths much larger than the size of a heavy quarkonium system. These soft gluons further couple to (s) and to complete such kinds of hadronic transitions.

Development of heavy meson chiral Lagrangians (HMCL) [15] is the foremost simplification to QCDME. HMCL serve as an effective field theory (EFT) to QCDME in a soft exchange where the gluonic exchanges are predominantly of limited momenta. With the assumptions that (i) the heavy involved in the process is well separated to consider it in a stringlike picture and (ii) the momentum of the emitted light meson is not too large, the HMCL are successful at reproducing the hadronic transitions among lower charmonia [16, 17]. The experimental status of the spectrum of higher vector charmonium(like) states is very rich now and several precise measurements have been recorded for their hadronic transitions [18]. To describe the observed transitions of higher systems there is a potential need for a theoretical model which can predict the transitions in the high momentum regime and help to identify the missing higher states through their hidden-flavor decays. We try to fulfill this need by modeling hadronic transitions of higher vector charmonia. Our proposed model is away from all the assumptions [(i) and (ii)] of HMCL and useful to predict the transitions involving much large momenta.

Another possibility is that the transition between two waves, to or to wave charmonia with the emission of () might occur through intermediate open-charm contributions. Heavy quarkonium states can couple to intermediate heavy mesons through the creation of a light quark-antiquark pair. The formalism which incorporate intermediate heavy mesons within hadrons is sometimes referred to as coupled-channel effects. For instance, using the quark pair creation mechanism [19], intermediate meson loop contributions are found to be essential to explain the suppression of dielectric decay widths of higher bottomonium [20]. Coupled-channel effects have also been taken into account in the QCDME framework to study the hadronic transitions with two-pion emission for the charmonium system and found a good agreement with experimental measurements [21]. In this paper, we neglect the coupled-channel effects for simplicity, which can be included in the future in the unquenched quark model.

To investigate the intermediate charmed meson loop effects on decay, nonrelativistic effective field theory (NREFT) formalism was constructed [9, 22, 23]. It is noted that if we go to much higher waves e.g., or with , the decay momentum is not so small, as it lies in the relativistic regime; hence, the NREFT formalism is not very suitable for studying hadronic transitions of higher charmonia.

These indications bring out the fact that there is a need for a comprehensive model which is capable of producing hadronic transitions with the emission of light meson(s) for higher mass charmonium systems. We attempt to fill this gap by modeling the transitions and , where refers to and vector charmonia with (). We present our predictions for hadronic transitions of higher vector charmonia into and , which provide useful references to determine their properties in ongoing and forthcoming experiments.

In hadron physics, the most widely used model to study open-flavor strong decays is or the quark pair creation (QPC) model. Within the framework of the model, quark pair creation is induced from QCD vacuum. Hence, the generated quark pair shares the quantum numbers of vacuum (); therefore, it is referred to as the pair creation mechanism. The traditional model has been widely used in hadron spectroscopy and decays [20, 24, 25, 26]. In the model, the probability to generate pairs is independent of the distance of the generation point from the valence quarks. In this work, we introduce the pair creation triggered from the NambuJona-Lasinio (NJL) four-point-like effective interaction (). Since is a mixture of scalar and pseudoscalar interactions, it raises the quantum numbers of the created pair as a mixture of and . It should be noted that the dynamics of the creation of a vertex in the NJL framework is totally different from the conventional mechanism.

Recently, the interaction of has been used to study the mixing of baryons with its pentaquark partner states and it is found that the leads to strong mixing between three-quark and five-quark (where ), with . It was reported that this expected mixing results in the decrease of the energy of the lowest state [27]. There is no hint of such kind of mixing within the conventional model (which only involves the scalar interaction). Hence these are charming motivations to consider this interaction to study the hadronic decays of higher quarkonia.

The paper is organized as follows. In Sec. 2, we review the development of the NJL model and its application in hadron spectroscopy and decays, where we deduce the effective Lagrangian for hadronic transitions with the emission of light meson(s). We conclude this section with an overview of well-established mixing formalism. Section 3 is devoted to discussing the results for , , and estimates of the branching fraction for , , and . Finally, we give a short summary in Sec. 4.

## 2 Theoretical Framework

### 2.1 NJL Motivated Effective Lagrangian

Effective field theories are very useful when the dynamics of the system involves only a few relevant degrees of freedom instead of all. The NJL model is one of the best examples of such kinds of effective theories which have the capability to recover almost all of the features of the exact leading theory. Historically, Nambu and Jona-Lasinio modeled a scheme to explain the pions as nucleon-antinucleon bound states [28, 29]; afterwards, the scheme gained more appreciation for being used at a more microscopic level by changing the nucleon field into a quark field . NJL model only involves the quark degree of freedom, while the gluon degree of freedom is frozen in its point-like interaction vertex. NJL four-point-like interaction between quarks can be described by

 LNJL=12gsNc∑a=0[(¯ψλaψ)2+(¯ψλaiγ5ψ)2], (1)

where indicates the color degree of freedom; are Gell-Mann matrices in flavor space with flavor singlet , where is the unit matrix in the three-dimensional flavor space, represents the quark field and is the coupling strength. This four-point-like color interaction with only one free parameter () has the capability to produce quite good results for the spectrum of low-lying and excited light mesons [30, 31]. Using the NJL motivated SU(2)SU(2) chiral Lagrangian for the excited pions, and mesons, the strong decay widths for the , , transitions ( and are the excited vector and pseudoscalar meson decaying into the vector and pseudoscalar meson, respectively) have been computed and found to be in good agreement with the experimental data [32, 33]. During the s, attempts were made to extend this approach to study the radiative transitions and strong decays of the charmed mesons. It is noted that the qualitative estimates of the strong decay widths using this approach agree well with the experimental data [34].

In light of these phenomenological studies we propose another quark pair creation mechanism which is inspired by the NJL model. We model the coupling of the light scalar and pseudoscalar meson with the charm quark. The effective Lagrangian of our model contains both the scalar and pseudoscalar interactions as present in the NJL model. In the case when we link the light production with an (anti)quark line, the effective Lagrangian of our proposed model can be written as

 LI=g(¯ψψ<σ>+¯ψiγ5ψ<η>), (2)

where is the overall coupling strength, is the heavy quark field, and and are singlet scalar and pseudoscalar meson, respectively. Since is the singlet, the light sector should also be in a singlet. That is why the above Lagrangian does not have flavor matrices as present in . The color index can also be suppressed. The above Lagrangian allows the coupling of the (anti)quark line only to a scalar or isospin singlet pseudoscalar. The possible Feynman diagram for the process is shown in Fig. 2.

In principle, can also couple to the antiquark line, so one needs to calculate two diagrams. But going through the details, it becomes clear that both diagrams are equally contributing.

Since the experimental data are available for , it provides us with an opportunity to test this model. No experimental data are available for the emission of a light scalar meson, so we mainly focus on emission transitions to check the validity of the model. It is quite possible that our proposed model can be extended to produce di-pion transitions. It would be through the intermediate production of a meson, which further decays into . But it will involve final state interactions (FSI) between and as investigated in Ref. [46]. This would be interesting but leads to intensive work which we will consider in our future studies.

To calculate the hadronic matrix elements, we prefer the mock hadron prescription [35] to express the initial and final meson wave functions. The mock hadron is defined as a collection of free quarks with the wave function of the bound quarks in a physical hadron normalized to the physical mass. The advantage of using this prescription is that it lets us calculate hadronic amplitudes as integrals over free-quark amplitudes. The initial mock meson wave function can be expressed as

 |A⟩=√2EA∑LS⟨Lm,SSz|JAmJA⟩∫d3p1ϕA(p1)χ12A|q1(p1),¯q2(−p1)⟩, (3)

where is the total energy of the meson; , , and are the orbital angular momentum between the quark and antiquark, the total spin of the quark-antiquark pair, and the total angular momentum of the meson, respectively; is the Clebsch-Gordan coefficient; and and are the spatial and spin wave functions of the initial meson , respectively. The relative momentum between the quark and antiquark, , is integrated over all values. We adopt the relativistic normalization

 ⟨A(p′)|A(p)⟩=2EAδ3(→p−→p′). (4)

For two-body decay, we define the transition amplitude as

 ⟨BC|HI|A⟩=2π√8EAEBECδ4(pi−pf)MmJAmJBmJC. (5)

Considering the standard relativistic phase space, we define the decay width in the center-of-mass (CM) frame as

 ΓA→BC=2πkEBECmA∑mJB,mJC∫dΩB|MmJAmJBmJC|2, (6)

for any fixed . Since the decay width is independent of the polarization of the initial state, we set in the following calculations. Here expresses the momentum of the outgoing mesons or , which is given as

 k=√[m2A−(mB−mC)2][m2A−(mB+mC)2]2mA, (7)

with and . The overlap of the wave functions of the initial meson and the final mesons can be expressed as

 MmJAmJBmJC=g∫d3p1ϕA(→p1)ϕ∗B(→p1−xB→PB)M0, (8)

where , with , with the charm quark mass. is the free-quark amplitude which is discussed in detail in the Appendix A. To find the overlap of the wave functions, we use simple harmonic oscillator (SHO) wave functions, which can be written in momentum space as

 ψnrlm(→p)=Rnrl(p)Yml(p,θ,φ), (9)

where , , and represent the radial, orbital, and magnetic quantum numbers, respectively. is the solid harmonic defined as a function of spherical harmonic. The radial wave function is given as

 Missing or unrecognized delimiter for \big (10)

where is an oscillatory parameter and is the associated Laguerre polynomial. Indeed, SHO wave functions serve as a coarse approximation to the true wave functions. However, qualitatively, SHO wave functions are similar to the realistic wave functions and useful for producing analytical results.

### 2.2 S−D Mixing

It is predicted that the charmonia near or above the open-charm threshold are an admixture of and waves [36, 37, 38]. An wave dominant state has a wave component in its wave function and vice versa. A well-established formalism of this mixing is based on reproducing the dielectric decay widths to deduce the mixing angle. If we neglect the open-charm contributions due to coupling to corresponding decay channels, under the assumption that the state only mixes with , and wave dominant states can be expressed as

 ψphys=cosθ|n3S1⟩+sinθ|(n−1)3D1⟩ (11)
 ψ′phys=−sinθ|n3S1⟩+cosθ|(n−1)3D1⟩ (12)

Here and represent the wave and wave dominant state, respectively. The relative sign between and is just a matter of convention. One may follow the other convention as in Ref. [37], but the effect of relative sign can be compensated by swapping . A rough estimate of the mixing angle can be made by computing the ratio of the dielectric decay widths [36, 37]. This has been done in the quark model framework by computing the wave functions using a potential model and then tuning the mixing angle to reproduce the dielectric decay widths. Considering as the dominant state with small component, there exist two sets of possible ranges for the values of the mixing angle: and [36, 37, 39].

There also exist a couple of quark-model-based phenomenological studies in favor of large mixing [38, 40]. A large mixing angle such as is used [38] to produce almost the same dielectric decay widths of and , which is consistent with experimental measurements [18]. The experimental fact that the dominant states have relatively small dielectric decay widths while dominant states have rather large widths can also be well described by considering the large mixing as the underlying mechanism [38].

## 3 Results and Discussions

### 3.1 Γ(Ψ→J/ψη)

To compute the decay widths for the process , it would be better to analyze the dependence of the wave function on the oscillatory parameter in the case of SHO wave functions. For light systems, the best value of from spectroscopy and decays is 4 [24, 41]. But for heavy states it is bit larger and also has a range in the literature [25, 42, 43, 44, 45]. The parameter relates to the size of quark-antiquark bound state. Since the size of heavy is smaller than the light system, for and for , where is the mass and is the velocity of the heavy quark. To choose the same for initial and final mesons is quite reasonable [46], although there exist some predictions that each meson has its own effective  [47]. Quark model studies show that the effective for higher multiplets is smaller than the corresponding lower ones. It is due to the fact that the excited states have large spatial extensions. For instance, for multiplets, and for multiplets, it is [42]. This is an indication that for higher charmonium states, the favorable value of should be around . The value has also been used to incorporate the spin counting predictions of open-flavor strong decays of higher wave states under the framework [44, 45]. We tune the parameter along with the coupling strength of the model. The coupling constant is fitted by choosing for all involved mesons, which agree with recent similar studies [43]. In principle, one can choose different values of for initial and final mesons to be more accurate, but for simplicity we choose the same for all charmonium states.

We explicitly show the dependence of the decay width in Fig. 3 for , , , , and , to clarify the possible acceptable range of the . It is worthy noting that the decay width does not change drastically around . It is clear from Fig. 3 that one can choose any value of within the safe region, i.e., . Our preferred value , lies in the safe region and, hence, is perfectly adequate.

The parameters used in our calculations are listed in Table 1. Table 2 shows the fitted results with the choice of best-fit values of the parameters of Table 1. We get quite impressive agreement with the experimental data. Although there exist only upper limits for the and decay processes, our computed decay widths for the former decay process lie within this limit, while for the latter process our predicted width is slighter larger than the central value. It is worthy noting that the experimental value of ) has large statistical errors. Considering this error range, our prediction in this case still lies within the upper limit.

Error estimation in the theoretical model is still an open question. It became an important debate among theoretical constructs in the last few years. In general, there are no surefire prescriptions for assigning error bars to theoretical models. In the case of parameter dependence, by doing numerical analysis one can confine the model space to a physically reasonable domain. Within this domain, there is a range of reasonable parametrizations that can be considered as delivering a decent fit [48]. Using this prescription we scan the parameters of our model and find the physically reasonable range of GeV and . For giving an idea of the uncertainties arising from model parameters, we explicitly give the decay widths in Table 3 by varying and in the described physical range.

The dependence of the decay widths on the mixing angle is very crucial to understand the behavior of this hadronic transition. We show the dependence of in Fig. 4 along with a band gap, which actually represents the decay width range when we consider the corresponding state as pure and wave states. Due to sinusoidal behavior, small mixing angles are adequate, while large mixing angles may ruin the predictions. Although Fig. 4 contains the predictions with a specific value of , still it can be seen that the decay widths of the process at and are exactly the same, which causes an exact overlap of the two curves in Fig. 3.

We also give the estimates of the decay widths of , , , , , and states decaying to in Fig. 5, which provide useful information to search and understand the missing higher vector states. Due to the fact that these higher states are poorly understood experimentally, we are not able to predict the widths exactly. To plot the decay width as a function of the mass of the corresponding higher vector state, we consider the mass range based on serval quark model predictions of the mass spectra listed in Table 4.

We include estimates of pure higher and wave states along with the predictions with small mixing. In all considered cases, the decay width of the wave state is smaller than that of the corresponding wave one. Given that the wave interferes destructively with the wave, its decay width is suppressed significantly with the use of a small mixing angle. After a particular mass value, the decay width of both and wave states becomes very sensitive to the initial mass. For example in the case, the decay width rises exceptionally after GeV. This critical mass value for other higher states can be observed from Fig. 5. With reference to the observed order of the emission rate, these critical mass values might provide the upper limits on the mass of the corresponding higher vector charmonium.

### 3.2 Predictions for ψ(nS/(n−1)D)→hc(1P)η with (n=3,4,5)

The hidden-charm decay of higher charmonia into the lowest wave meson, i.e., , is also possible. The threshold for this decay process is MeV. So, the first vector state which can decay into is . There exists an experimental evidence for around MeV mass by the CLEO Collaboration [53]. Their reported measurement on the branching fraction is with a confidence level. The HQSS violating transition requires the spin flip to be significantly suppressed relative to the corresponding heavy quark spin-conserving transitions like  [54]. The observed ratio is fully consistent with the earlier theoretical predictions [55]. It has been argued in [23, 56] that the coupled-channel effects due to intermediate charmed mesons for these transitions are quite small.

Among well-known higher vector charmonia, only and have enough phase space to decay into . Table 5 contains our predictions for these states. We predict the decay width of with the same order of magnitude as to the similar final state:

 Γ(ψ(4160)→hc(1P)η)Γ(ψ(4160)→J/ψη)=7.887×10−2, (13)
 Γ(ψ(4415)→hc(1P)η)Γ(ψ(4415)→J/ψη)=6.736×10−2. (14)

It is not easy to give the estimates of the decay width for , , and or higher ones because these states have not been experimentally well established up to now and hence, their masses are unknown. We give the initial mass dependence of the decay width of the transition of these higher vector states in Fig. 6, both for the pure and wave and for the standard mixing case.

For , the decay width of the pure wave is larger than that of the corresponding pure wave. It is due to that the overlap of the wave function of the wave with is larger than the corresponding wave. The wave interferes destructively with the wave; the decay width of wave, in this case, is not suppressed significantly as we have observed for the channel. For the case, the reason why the decay width of the pure wave is smaller is that the overlap of the wave with is smaller than the corresponding wave state. The swapping of the curves can be seen in Fig. 6 with the increase of the mixing angle. The constructive interference of the wave with the wave causes a significant increase in the decay width of the wave dominate state.

### 3.3 Y(4360), Y(4390) and Y(4660) Assignments

Just after the experimental observation of and in the initial state radiation (ISR) process at Belle [57], many theoretical studies were carried out to incorporate these states into conventional and exotic spectra (for an overview, see the discussion in Sec. 4.8 of the recent review [58]). Among these configurations, there exist a couple of interpretations by considering and as canonical  [50] and charmonium [50, 52], respectively. Predictions were also made for the dielectric widths, E1 and M1 transitions, and open flavor strong decays. In the screened potential model, the state is also interpreted as , while was considered as  [49]. Assuming as and as dominate, their dielectric decay widths can be reproduced to get agreement with experimental data by introducing the large mixing [38].

The predicted mass of in the BGS and RS potential model (as shown in Table 4) is somewhat larger than the experimental mass of . However, we notice that the mass predictions of various potential models for the higher states may differ by MeV [4, 59]. Therefore, it is not irrational to treat as dominant with a small component of , exactly in the same way as for the lower waves in Sec. 2.2. Mass predictions for of various potential models do not differ much from except the one produced by the screened potential model. Hence, it is also not irrational to treat as dominant having a small component of .

Very recently, the BESIII Collaboration announced the observation of a new resonant structure, in the process [60]. Its measured mass and total decay width are MeV and , respectively. The width of seems to be broader. It might be possible that this higher mass region could be described either by a resonance, or by a phase-space background as has already been noticed in the case of  [61]. Since the reported quantum numbers of this state are and it lies in the mass region of charmonium, hence, this new resonance may also be considered as a candidate for . As an estimate, we give our predictions of its transition branching fraction in Table 7 by assigning .

One may argue that the vector state can also be assigned as a or charmonium state. As listed in Table 4, in various potential models, and are predicted between GeV. Therefore, with reference to the observed mass of , it lies below the mass region. The value of its is two to three orders of magnitude smaller than those of the well-established conventional vector charmonia [62]. It is nearly impossible to accommodate as a conventional charmonium state. Many theoretical constructs consider as an exotic charmonium [61], for example, as a molecule with a binding energy of 29 MeV [63], or as a heavy hybrid meson with a gluonic excitation of about 1 GeV higher than the lowest and  [58, 62].

Despite the fact that these states do not decay into open-charm channels, it would be interesting to study their hidden-charm strong decays. By assuming and as and dominant states, respectively, we give our predictions for and , which might be helpful to understand the properties of these vector states. Because only experimental upper limits [64] exists for the product of the branching fraction and for and , we need to know the dielectric decay width of these states. There exist few theoretical predictions for for these states as summarized in Table 6.

As shown in Table 7, the only available experimental information is for these states. Therefore, we need to compute . Table 6 contains different theoretical predictions for of these states both for pure and waves, and with the large mixing (). For , we take the average value of keV from Table 6 by considering it as pure . The experimental upper limit for is given in Table 7. To give a comparison with Ref. [38], we also compute the upper limit of with large mixing, i.e., . In all three cases, pure , small, and large mixing, our predictions are in agreement with the experimental measurements. We conclude that could be considered as a potential candidate for dominant charmonium state.

For , we also include the predictions for pure and the mixed case. For the pure state, we take the average value of keV from Table 6 and list the experimental upper limit along with our prediction for in Table 7. In the large mixing case (), the dielectric decay width keV [38] allows us to give an upper limit on . Our predicted value, in this case, is within this upper limit as shown in Table 7. In all cases, our predictions agree with the experimental data. Hence, our results are consistent with the experimental data and the state can be considered as dominant with a small component.

In the case of , for the sake of completeness, we give our predictions of its hidden charm decay, with and without mixing. To identify this state, measurements on its hadronic branching fraction are required. We think that these estimates might be useful to clarify the picture of these vector states and give some references to search for the missing higher and wave vector charmonia.

## 4 Summary

A model to create a light meson for heavy quarkonium transition is proposed. This model is used to study the decays of higher vector charmonia into and . Computed decay widths are in excellent agreement with experimental data. The ratio is predicted for and . The initial state’s mass dependence of for higher vector charmonium is given. We suggest that the ongoing (Belle and BESIII) and forthcoming (ANDA and BelleII) experiments look for suggested unobserved decay channels. We also give the estimates of transition branching fractions for , , and by assuming them as bound states with quantum numbers , , and , respectively. Our predictions reflect that the state can be considered as a potential candidate for the charmonium state. Assuming to be , the predictions are consistent with the experimental upper limit. For a broader state, the update on its hadronic branching fraction from BESIII is eagerly awaited. We hope that our predictions might provide useful references to determine the properties of higher charmonium states in ongoing and forthcoming experiments.

## Acknowledgements

The authors are grateful to Feng-Kun Guo, Qiang Zhao, Martin Cleven, Eric Swanson and Ahmed Ali for useful discussions and suggestions. We are indebted to Shan-Gui Zhou for providing us Ref. [48] during the revision of this manuscript. M. N. A. gratefully acknowledges the hospitality at the DESY theory division, where part of this work was carried out. This work is supported by the National Natural Science Foundation of China under Grant No. 11261130311 (CRC110 by DFG and NSFC). M. N. A. received support from the CAS-TWAS President’s Fellowship for International Ph.D. Students.

After the revision of this manuscript, BESIII published an evidence of at center-of-mass energy GeV [65]. Along with their earlier measurement [66], this evidence will help to delve into the through its HQSS violating hadronic transitions.

## Appendix A Free-Quark Amplitude

To evaluate the matrix element of decay, we need spin matrix elements. At the quark level, these matrices involve the matrix elements of the Dirac bilinear (with and in our case) and Pauli matrix elements. Matrix elements of the Dirac bilinear in the nonrelativistic limit can be expressed as

 limq→0¯uq′s′Γuqs={δss′Γ=Ii2mq⟨s′|→σ|s⟩⋅(→q−→q′)Γ=iγ5 (15)
 limq→0¯v¯q¯sΓv¯q′¯s′={−δ¯s¯s′Γ=Ii2mq⟨¯s|→σ|¯s′⟩⋅(→¯q−→¯q′)Γ=iγ5. (16)

We have already shown that it is easy to handle the wave function overlap integration in a Cartesian basis; therefore, we express the elements of Pauli spinors in terms of Cartesian basis vectors as

 ⎧⎨⎩⟨↓|→σ|↑⟩=(^x+i^y)⟨↑|→σ|↓⟩=(^x−i^y)⟨↑|→σ|↑⟩=−⟨↓|→σ|↓⟩=^z (17)

For antiquark case, these relations are

 ⎧⎪⎨⎪⎩⟨¯↓|→σ|¯↑⟩=−(^x−i^y)⟨¯↑|→σ|¯↓⟩=−(^x+i^y)⟨¯↑|→σ|¯↑⟩=−⟨¯↓|→σ|¯↓⟩=−^z (18)

From Fig. 2, one can get the following relation of momentum conservation by considering the initial state in the center-of-mass reference frame:

 ⎧⎨⎩p2=−p1p′1=p1−kPB=−k. (19)

The meson’s space wave function can be written as

 ⎧⎨⎩ϕA{μA(p1m1−p2m2)}=ϕA(p1)ϕB{μB(p′1m1−p2m2)}=ϕB(p1−m1m1+m2PB) (20)

where the are the reduced masses of the constituent quarks of the mesons and , respectively. It should be noted that all the momenta should and can be expressed by the integration variable and the momentum of the meson, viz. and in this work. The free-quark amplitude for Fig. 2 is

 iM0=[¯uq′s′(iγ5)uqs][¯v¯q¯sIv¯q′¯s′], (21)
 iM0=i2mc⟨1′|→σ|1⟩⋅(→p1−→p′1)⋅⟨2|δss′|2′⟩, (22)

where is the mass of the charm quark. Collecting all the pieces together, the full amplitude becomes

 MmJAmJBmJC=gi2mc∫d3p1ϕA(→p</