Hunting for states in the recent LHCb di-$J/\psi$ invariant mass spectrum

Partial wave analysis is performed, with effective potentials as dynamical inputs, to scrutinize the recent LHCb data on the di-$J/\psi$ invariant mass spectrum. Coupled-channel effects are incorporated in the production amplitude via final state interactions. The LHCb data can be well described. A dynamical generated pole structure, which can be identified as the $X(6900)$ state, is found. Our analysis also provides hints for the existence of three other possible states: a bound state $X(6200)$, a broad resonant state $X(6680)$ and a narrow resonant state $X(7200)$. The $J^{PC}$ quantum numbers of the $X(6680)$ and $X(6900)$ states should be $2^{++}$, while the $X(6200)$ and $X(7200)$ states prefer $0^{++}$. To determine the above states more precisely, more experimental data for the channels, such as $J/\psi\psi(2S)$, $J/\psi\psi(3770)$, di-$\psi(2S)$, are required.


I. INTRODUCTION
In 2003, the first tetraquark candidate state χ c1 (3872) was discovered by Belle collaboration [1], and later confirmed by BABAR [2], CDF [3] and D0 [4] groups. The discovery of this state has inaugurated a new era of studying multiquark states that provide us a unique platform to gain more insights into the low-energy quantum chromodynamics (QCD). Since then, many charmonium-and bottomonium-like states have been observed by various experiments and, meanwhile, intrigued intensive theoretical investigations, see e.g. Refs [5][6][7][8][9][10] for reviews on the studies of XY Z particles.
Attempts have also been made to search for states containing only four heavy quarks by, e.g., the LHCb [11,12] and CMS [13] collaborations. Recently, the LHCb collaboration declared the observation of a narrow structure around 6.9 GeV and a broad structure located in the energy range [6.2, 6.8] GeV in the di-J/ψ invariant mass event distribution, using the proton-proton collision data at centre-of-mass (c.m.) energies of √ s = 7, 8 and 13 TeV [12]. The data of the di-J/ψ mass spectrum also hint a possible structure in the vicinity of 7.2 GeV. The narrow peak can be described by employing the Breit-Wigner parametrization and an associated X(6900) state is established with significance larger than 5σ. The mass and width of this state are determined to be M[X(6900)] = 6886 ± 11 ± 11 MeV , (1) Γ[X(6900)] = 168 ± 33 ± 69 MeV , in the scenario where the interference between the resonant contribution and the nonresonant single parton scattering is implemented. In the case without interference, the resultant determinations are M[X(6900)] = 6905 ± 11 ± 7 MeV, The observation of X(6900) is the first experimental evidence of the fully charmed tetraquark T cccc , in another word, a new member in the XY Z particle zoo is gained.
Turning back to the newly observed X(6900) state, plenty of theoretical works have been accumulated aiming at studying its properties, see e.g. Refs. [17,19,20,27,30,31,33,[37][38][39][40][41][42][43][44][45]. For instance, the mass of X(6900) obtained by QCD sum rules [30,31] agrees well with the experimental values. The X(6900) state is interpreted as an excited state both in the constituent quark model [27] and in the QCD sum rule approach [33]. As for its inner structure, it is pointed out in Ref. [19] that the X(6900) resonance is most likely to be a compact fourcharm quark state with little component of molecular nature. Nevertheless, based on spectral density function sum rule and pole counting rule, Ref. [17] suggests that more experimental data are needed in order to distinguish the nature of X(6900). Interestingly, it is claimed in Ref. [37] that the X(6900) state can be described as a Higgs-like boson, which conveys some signal beyond standard model. It is worth noting that another fully charmed state, denoted as X(7200), is found in e.g.

arXiv:2104.08589v2 [hep-ph] 31 Aug 2021
Refs. [17,30,38], which actually is also hinted by the LHCb data as mentioned above. In this work, we intend to hunt for all possible states in the di-J/ψ spectrum below 7.6 GeV, and try to determine their J P C quantum numbers.
To that end, we construct the amplitude for the di-J/ψ production on the pp collision, in which the coupledchannel effects are taken into account via final state interactions (FSI). The {J/ψJ/ψ, J/ψψ(2S), J/ψψ(3700), ψ(2S)ψ(2S)} channels are included in the T matrix with the aim of reproducing the di-J/ψ invariant mass spectrum, reported by the LHCb collaboration [12], in the range from 6.2 to 7.6 GeV. Of course, there are other double-charmonium channels with thresholds below 7.6 GeV, which can couple to the di-J/ψ system, such as η c η c , h c h c and χ cJ χ cJ . However, their contributions are either suppressed by heavy quark spin symmetry or subordinate due to the smallness of the relevant couplings estimated by meson exchanges, as pointed out by Ref. [18], and the readers are referred to Appendix B for detailed discussions. We further perform partial wave projection of the T matrix within helicity formalism. Explicit expressions for the S-wave amplitudes with J P C = 0 ++ , 2 ++ are obtained. We assume that the invariant mass spectrum is dominated by the S waves, and hence those partial waves beyond S-wave are neglected.
Various fits are performed to the LHCb data with the obtained production amplitudes in S wave. Due to the closeness of the X(6900) peak to the J/ψψ(3770) threshold, we first take only three channels, i.e., {J/ψJ/ψ, J/ψψ(2S), J/ψψ(3770)}, into account, and do fits with partial-wave amplitudes of 0 ++ and 2 ++ , respectively. It is found that both fits can well reproduce the experimental data in the energy range below 7.2 GeV. A di-J/ψ subthreshold state, named X(6200), is discovered both in the 0 ++ and 2 ++ amplitudes. However, there is no dynamically generated state responsible for the peak around 6.9 GeV. Thus, we make further fits, where the ψ(2S)ψ(2S) channel is now switched on and the fitting range is extended up to 7.6 GeV to cover more data. Following the same procedure as the three-channel fits, we find that the 0 ++ and 2 ++ amplitudes behave differently now. Furthermore, four dynamically generated states are established. We also carry out a combined fit with both 0 ++ and 2 ++ contributions, which indicates that the two partial-wave amplitudes are simultaneously sizable.
This manuscript is organized as follows. Section II serves to illustrate the basic ingredients of our theoretical framework, including coupled-channel potentials, partialwave projection, unitarized amplitudes and production amplitude of di-J/ψ. Numerical results are shown in Sec. III. A brief summary is given in Sec. IV. The explicit expressions of the helicity amplitudes are collected in Appendix A.

II. THEORETICAL FRAMEWORK
In this section, the S-wave production amplitude, with 0 ++ and 2 ++ quantum numbers, is constructed in order to reproduce the recent di-J/ψ invariant mass spectrum.
For a given angular momentum J, the partial wave projection of the helicity amplitudes can be obtained via with λ = λ 1 −λ 2 and λ = λ 3 −λ 4 . Here s = (p 1 +p 2 ) 2 and t = (p 1 − p 3 ) 2 are Mandelstam variables. Furthermore, d J λλ (z s ) are standard Wigner functions, z s ≡ cos θ and θ is the scattering angle. Hereafter, we focus on the S wave, therefore, L = 0 and J = S with S denoting the total spin of the initial or final states. The generalized Bose symmetry for identical particles and the conservation of J P C quantum numbers imply that L + S should be even, which means that J should take values of 0 or 2.
For S wave, the partial wave amplitudes in the LJS basis can be obtained through where the transformation matrix regarding the initial states is given by Here S 1 and S 2 are spins of the incoming vector particles V 1 and V 2 , respectively. Clesch-Gordon coefficients are represented by S 1 λ 1 S 2 − λ 2 |Sλ . The transformation matrix concerning the final states can be obtained exactly in the same way. Finally, the partial wave with J P C = 0 ++ reads Likewise, the partial wave with J P C = 2 ++ reads

C. Unitarized amplitudes
The conservation of probability implies a realistic amplitude should be unitary. This can be achieved by applying the on-shell-approximation version of the Bethe-Salpeter equation [46,47] to the partial wave amplitudes derived in the above subsection; see Refs. [48][49][50] for recent reviews on various unitarization approaches. In this manner, the unitarized amplitude is given by with where the subscripts of the entries of the matrix are channel indices. G(s) is a diagonal matrix defined by with its elements given by [51] g i (s) = 1 16π 2 a(µ) + ln Here , and a(µ) is a subtraction constant defined at the renormalization scale of µ.
Bound, virtual and resonant states correspond to pole singularities in different Riemann sheets (RSs) of the unitary T J matrix. RSs can be defined by performing analytical continuation of the g i (s) function through with ξ i = 0, 1 and p i (s) being the c.m. momentum in the ith channel. Furthermore, each RS can be denoted by a number where i is channel index and N is the total number of the involved channels. In practice, the nth RS is often marked by Roman numerals. For example, the fourth RS is equivalent to RS-IV.

D. Production amplitude
The di-J/ψ production amplitude, with given quantum numbers J P C , can be written as Here G is given by Eq. (15); A i 's represent direct productions of the di-J/ψ, J/ψψ(2S) and J/ψψ(3700) states, respectively. The second term on the right-hand side of the above equation stands for FSI where coupled channels rescattering effects of {J/ψJ/ψ, J/ψψ(2S), J/ψψ(3700), ψ(2S)ψ(2S)} are implemented. A diagrammatic representation of such a production mechanism is displayed in Fig. 1. The production amplitude M 1 (s) can be rewritten as The parameters r i 's are defined by r i = A i /A 1 and r 1 = 1. In principle, r 2 and r 3 are unknown complex parameters. Eventually, the invariant mass of the double-J/ψ spectrum is proportional to the amplitude squared multiplied by a phase space factor ρ(s), namely where the explicit expression of ρ(s) is Here λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz is the so-called Källén function.

III. NUMERICAL RESULTS AND DISCUSSIONS
The LHCb collaboration has recently reported the invariant mass spectrum of J/ψ pairs, coming from protonproton collision data at 7, 8 and 13 TeV [12]. In general, the di-J/ψ mass spectrum is dominated by the S-wave interaction, however, there are two S-wave candidates with different J P C quantum numbers either being 0 ++ or 2 ++ . Since we have both the 0 ++ and 2 ++ amplitudes, as given in the previous section, this enables us to perform a partial wave analysis of the di-J/ψ spectrum data and reveal possible underlying states with definite J P C . Therefore, in what follows, we are going to confront our theoretical model, given in Eq. (21), with the recent LHCb data.

A. Fits with three-coupled channels
We first consider the case with the three channels of {J/ψJ/ψ, J/ψψ(2S), J/ψψ(3770)}, which, in principle, are sufficient to reproduce the data in the energy region below 7.2 GeV where the structure around 6.9 GeV is covered. With Eqs. (20) and (21), the invariant mass spectrum of di-J/ψ is described by It should be emphasized that the unit in Eq. (20) is replaced by a constant γ in order to simulate the coherent background. Furthermore, r 2 and r 3 are set as 1, such that the number of free parameters is reduced. Actually, the effects of r 2 and r 3 may be partly absorbed by the γ constant and other free parameters involved in the rescattering amplitudes T i1 (s). In Eq. (23), the amplitude |A 1 (s)| 2 stands for direct production of J/ψ pairs, which encodes the information on short-distance interactions. Following Ref. [18], its modulus squared can be parametrized as where α and β are unknown overall constant and slope parameter, respectively. The distribution of events corresponding to the double-parton scattering (DPS) [52][53][54], shown in the LHCb paper [12], can be described by |A 1 (s)| 2 multiplied by the phase-space factor ρ(s), which determines α = 134 and β = 0.0123. Nevertheless, the overall constant is released as a free fitting parameter in our fits to be discussed below.
Coupled channel effects are implemented through FSI by the summation term in Eq. (23). We employ a uniform subtraction constant a(µ) for all the two-point G ii (s) functions. 2 Besides, the subtraction constant is fixed, i.e., a(µ = 1 GeV) = −3.0, where the renormalization scale is set to 1 GeV. We have checked that a(µ) tends to take the value of −3 even if it is set as a free fitting   parameter. In fact, this value also corresponds closely to the natural value of a subtraction constant [19,47], which can be directly calculated by using the formula derived in Ref. [56]. Note that we are using the same value for the subtraction constant for all the fits performed in this work. As for the scattering amplitudes T i1 (s), there are a few coupling constants stemming from the effective Lagrangian. To reduce the number of fitting parameters and to obtain natural values for those couplings, we employ where m j are masses of the vector particles involved in the scatterings, as specified in Eq. (A1). The values of the masses are taken from PDG [57]. In this way, theh i 's are chosen to be fitting parameters instead of h i 's. of Fig. 2. Fit-A is performed up to 7.2 GeV. Excellent agreement between our model and the LHCb data is achieved with a χ 2 /d.o.f. 0.96. The obtained values of the coupling constantsh i turn out to be of natural size as expected. It can also be seen from the figure that, the di-J/ψ spectrum is well described up to 9 GeV, even though the fitting range of Fit-A is from 6.2 to 7.2 GeV. Fit-B corresponds to the case that the production amplitude consists of the partial wave with J P C = 2 ++ only. The fitting range of Fit-B is the same as that of Fit-A. Results of Fit-B are compiled in the third column of Table II, and plots are displayed in the lower panel of Fig. 2. It can be found that the fit quality of Fit-B is as good  as Fit-A below 7.2 GeV, and the data beyond the fitting range are again well described within 1-σ uncertainty.
Since the coupling constants involved in the potentials are determined by the fits, it is now ready for us to investigate the pole structures of the unitarized amplitude in the RS. The pole positions and their residues are obtained using the fitted values of parameters of Fit-A and Fit-B, which are collected in Tables III and IV, respectively. In each case, a near-threshold bound state is found: √ s pole = 6173.9 +19.6 −41.5 MeV for the fit with partial wave of 0 ++ and √ s pole = 6169.3 +23.9 −44.2 MeV for the fit with partial wave of 2 ++ . This state is referred to as X(6200) in Ref. [18]. Compared with the three-coupled potentials used in Ref. [18], we systematically incorporate energy dependent terms via an effective Lagrangian approach. We conclude that the inclusion of those terms does not affect the existence of the X(6200) bound state. However, the LHCb data of the di-J/ψ invariant mass distribution in the range [6.2, 7.2 GeV] are described almost equally well by the 0 ++ and 2 ++ fits. Namely, the difference between the 0 ++ and 2 ++ partial wave amplitudes is not sensitive to the events distributed in the fitting range at all, preventing us from disentangling the J P C quantum numbers of this state.
The peak around 6.9 GeV is well described both in Fit-A and Fit-B, however, to which the dominate contribution is from the effect of the J/ψψ(3770) threshold. No resonant poles, responsible for this peak, were found in the nearby energy region. In Tables III and IV, the RS-II poles are presented for easy comparison with the results in Ref. [18].

B. Fits with four-coupled channels
As discussed in the previous subsection, the LHCb data below 7.2 GeV cannot distinguish between the par- tial waves of J P C = 0 ++ and 2 ++ in the three-coupled cases. To tackle thus issue, one may extend the fitting range to include more data at higher energies, such that the partial wave amplitude involved in the invariant mass spectrum formula Eq. (23) will be much more constrained. On the other hand, there exists a bump of events distribution just above 7.2 GeV, as can be seen in Figure 2, and the prediction, based on the coupled J/ψJ/ψ-J/ψψ(2S)-J/ψψ(3770) model, actually fails to describe it. Thus, we improve our formulation by incorporating one more channel, the ψ(2S)ψ(2S) channel, at the price of introducing three extra unknown parameters h 7 ,h 8 andh 9 .   6 GeV]. The above value for M cut J/ψJ/ψ is chosen such that a plateau-like behavior for Fit-D and a downtrend for Fit-C end simultaneously, as can be seen from Fig. 3. Compared to the three-coupled-channel fits, the fitting quality of Fit-C is improved in view of the obtained χ 2 /d.o.f. value, due to the inclusion of the fourth channel with more free coupling constants. However, we find that in the current cases the fitting parameters are now much more correlated, and consequently the resultant statistical errors become larger than those of the above threecoupled-channel fits. A comparison between the LHCb data and our predictions is shown in Fig. 4. It can be seen from the figure that the theoretical results of invariant mass spectrum, contributed either by the 0 ++ partial wave or by the 2 ++ one, behave quite differently now, as expected.
For Fit-C where the 0 ++ partial wave is used, the agreement is excellent in the sense that the theoretical line shape coincides with all the details of the LHCb data in the entire fitting range. Based on Fit-C, a pole is found in the RS-VIII and its location reads This resonant state is called X(7200) in literature, e.g., Ref. [17]. The coupling strengths of the X(7200) state with the four channels, i.e., the square roots of the magnitudes of the residues, are specified in Table VI. As can be seen from the Table, this state is most strongly coupled to ψ(2S)ψ(2S). The emergence of the X(7200) state accounts for the bump around 7.3 GeV in the events distribution. In addition to the discovery of the X(7200) state, a pole of √ s pole = 6124.8 +23.9 −121.8 MeV is obtained as well. It can be regarded as the X(6200) state discussed in the preceding three-coupled-channel analysis, but the pole location is shifted by an amount of about 50 MeV due to the incorporation of the di-ψ(2S) channel. It is worth noting that the above two states observed in the 0 ++ partial wave amplitudes are absent in the 2 ++ ones, providing us evidence that their J P C quantum numbers should be 0 ++ . Detailed pole information on the two states is given in Table VI. For Fit-D in which the 2 ++ partial wave is implemented, our prediction is in good agreement below 7.2 GeV, however, deviates from the events data of the bump around 7.3 GeV. Nevertheless, the procedure of pole hunting reveals a pole structure located at √ s pole = (6919.8 +17.2 −23.7 − i58.8 +10.9 −12.2 ) MeV . (27) which can be associated with the X(6900) state reported by the LHCb collaboration [12]. This finding indicates that the J P C quantum numbers are much more likely to be 2 ++ . It can also be found from Table VI that the coupling strength of X(6900) to the di-ψ(2S) pair is much larger than the states in the other three channels. Hence, the di-ψ(2S) channel is of particular importance to dynamically generate the X(6900) resonant state. In other words, the X(6900) state cannot be dynamically generated unless the heavier channel ψ(2S)ψ(2S) is included. In the absence of the ψ(2S)ψ(2S) channel, an alternative way of obtaining the X(6900) state is to introduce a Castillejo-Dalitz-Dyson (CDD) pole in the amplitude, as done in Ref. [19]. It seems compatible between Ref. [19] and this study that the former requires a CDD pole and the latter finds that the resonance is mostly due to the dynamical contributions from the ψ(2S)ψ(2S) channel, which is not included in Ref. [19] because of its rather large threshold compared to the X(6900) mass. Namely, a CDD pole could mimic dynamically generated states by heavier channels, which are not included explicitly in the coupled-channel study. Interestingly, a RS-II pole with a large imaginary part is also achieved in the 2 ++ partial wave amplitude, which reads The broad structure, ranging from 6.2 to 6.8 GeV in the spectrum and dubbed threshold enhancement by the LHCb collaboration [12], probably contains also an inherent ingredient, here represented by the obtained dynamically generated pole of J P C = 2 ++ , in addition to the kinematical threshold effects. For more information on the poles we obtained, see Table VI. In Fig. 5, we plot the pole locations in the complex √ s plane. For comparison, the results of X(6900) reported by the LHCb collaboration are also shown. Our determination of the X(6900) state is consistent with the experimental values within 1-σ uncertainties.

C. Combined fit
It is believed that the invariant mass spectrum should be dominated by S-wave contributions. However, there is no a priori criterion for us to judge if the 0 ++ partial wave or the 2 ++ one is more important than the other, since both of them have the same orbital angular momentum L = 0. It might be an appropriate way to treat them as equally important ingredients and, meanwhile, let the events data make the judgement. To that end, the invariant mass spectrum is recast as where we have used J = 0, 2 to denote the 0 ++ and 2 ++ amplitudes in S wave. Here, the incoherent backgrounds are assumed to be uniform. That is A J=0 1 (s) = A J=2 1 (s), which are parametrized in the same form as Eq. (24). On the other hand, we adopt different coherent backgrounds for J = 0 and J = 2 cases, represented by γ 0 and γ 2 in Eq. (29), respectively.  With the above preparations, fits incorporating both the 0 ++ and 2 ++ waves can be performed. It is found that, when we make a combined fit with four-coupled channels up to 7.6 GeV, the data around 6.9 GeV seem not sufficient enough to distinguish the nearby threshold effect of the J/ψψ(3770) channel from the contribution of a X(6900) state. As a result, a stable solution cannot be obtained in the four-channel fit. There-TABLE VIII. Poles and their residues based on the combined fit. The RSs, on which the poles are located, are given in the first column. The J P C quantum numbers of the obtained states are specified in the brackets.  Fig. 6. It can seen from Fig. 6 that both the 0 ++ and 2 ++ waves contribute sizably to the invariant mass spectrum. The former yields a bound state, located at √ s pole = 5979.6 +3.6 −9.1 MeV. It should correspond to the X(6200) state, if considering the influence caused by the closure of the J/ψψ(3770) channel. The latter allows the existence of a broad resonance and a narrow pole for the X(6900) state, as shown in Table VIII.
However, unlike the results of poles in Fit-C, there is no pole structure that accounts for the event enhancement around 7.3 GeV now. We owe the absence of such pole to the fact that the J/ψψ(3770) channel is excluded, which could have considerable impact on the formation of this pole. More data from experiments are indispensable for performing a comprehensive combined fit with four channels of {J/ψJ/ψ, J/ψψ(2S), J/ψψ(3770), ψ(2S)ψ(2S)}, in order to draw a solid conclusion on the existence of the X(7200) state discovered in Fit-C.

IV. SUMMARY AND CONCLUSIONS
In this work, we have presented a partial wave analysis of the recent di-J/ψ invariant mass spectrum. Coupledchannel effects, due to the rescatterings among the {J/ψJ/ψ, J/ψψ(2S), J/ψψ(3770), ψ(2S)ψ(2S)} states, are included in the production amplitude of pp → J/ψJ/ψ + anything. The unitary S-wave amplitudes with J P C = 0 ++ and 2 ++ are obtained with the help of the Bethe-Slapeter equation under on-shell approximation. Various fits are performed and the di-J/ψ invariant mass spectrum can be well reproduced. Nevertheless, it is found that, when fits are performed up to 7.6 GeV, the 0 ++ and 2 ++ amplitudes behave differently, enabling one to determine the J P C quantum numbers of the dynamically generated poles. Our final results are based on the four-coupled channel fits. Four states are found. In the case of 0 ++ wave, a bound state X(6200) is located below the di-J/ψ threshold and a narrow resonant state X(7200) resides at one can derive explicitly the 25 independent helicity amplitudes, which are listed below.
where p c.m. = | p 1 | = | p 2 | andp c.m. = | p 3 | = | p 4 | are the modulus of the momenta of the initial and final states in the center-of-mass (c.m.) frame, respectively. Furthermore, ω i = m 2 i + | p i | 2 (i = 1, . . . , 4) are the energies of the involved states. To obtain the above expressions, we have used the following polarization vectors for particles 1 and 3. It is straightforward to obtain the ones for the states V 2 and V 4 in the c.m. frame. Note that the scattering plane has been chosen to be the plane spanned by the x-axis and z-axis, such that the azimuthal angle φ = 0.
Feynman diagrams of the two types of scattering processes in meson-exchange picture: (a) ψψ → χcJ χ cJ , (b) ψψ → ψψ. The quantum numbers I G of the involved particles are shown in the brackets, where I and G are isospin and G parity, respectively.
can be found in Ref. [58]. In the S-wave doublet, ψ µ and η c are the vector mesons and pseudoscalar mesons, respectively. For the P -wave charmonium bound states, the multiplet J µ is composed of four states, χ c2 , χ c1 , χ c0 and h c . The conjugations of the fields are defined bȳ J = γ 0Ĵ † γ 0 andJ µ = γ 0 J µ † γ 0 . Note that the radial numbers n of the charmonia are suppressed for brevity. Now we expand the HQSS Lagrangian by inserting Eqs. (B2) and (B3) into Eq. (B1). The g 1 term gives In the heavy quark formalism, ψ µ (η c ) annihilates a state, while its conjugation ψ † µ (η † c ) creates a state. Therefore, one can use the g 1 term to describe the scattering processes of ψψ → ψψ, ψη c → ψη c , and η c η c → η c η c . However, the process such as ψψ → η c η c does not show up in the Lagrangian, as expected.
The g 2 term contains the χ cJ and h c states, which can be expanded as g 2 J µĴ J µĴ + H.c. = 4g 2 χ †µα c2 ψ α χ †ρ c2µ ψ ρ − 2g 2 (χ † c1γ ψ β χ †β c1 ψ γ − χ † c1γ ψ β χ †γ c1 ψ β ) + where "· · · " denotes the pieces with the number of ψ fields less than 2. It can be seen that the channels of ψψ → χ cJ χ cJ (J, J = 0, 1, 2) appear explicitly in the HQSS Lagrangian. On the contrary, the reaction ψψ → h c h c is absent, since it violates HQSS by flipping the charm-quark spin. We can further estimate the strength of ψψ → ψψ and ψψ → χ cJ χ cJ in the meson-exchange picture. Both of the two types of interactions are OZI allowed. Taking into the SU(3) flavor and isospin symmetries together with G-parity conservation into account, the lowest meson can be exchanged for the ψψ → ψψ is the f (500) meson (i.e. σ), while for ψψ → χ cJ χ cJ it is the ω meson, as illustrated in Fig. 7. 3 By integrating out the exchanged meson, one can expect that where M σ and M ω are masses of f (500) and ω, respectively. Since M ω > M σ , and hence g 2 < g 1 . That is, the coupling of ψψ → ψψ is larger than the one of ψψ → χ cJ χ cJ . Therefore, it might be a good approximation to take only the type of process of ψψ → ψψ into consideration, when studying the di-J/ψ invariant mass spectrum with limited experimental data. For convenience, we denote the 4ψ interaction in the HQSS Lagrangian by and briefly discuss its relationship with the effective Lagrangian given in Eq. (5) in Sec. II. To that end, the radial quantum numbers of the charmonium states have to be invoked and the L 4ψ Lagrangian is recast into L 4ψ = 4g 1 (mn; m n ) × ψ † µ (m S)ψ µ (mS)ψ † ν (n S)ψ ν (nS) . (B8) Here g 1 (mn; m n ) denotes the coupling constant corresponding to the process of ψ(mS)ψ(nS) → ψ(m S)ψ(n S). Note that ψ(1S) should be identified as the J/ψ state. Correspondence between the Lagrangian in Eq. (B8) and the one in Eq. (5) can be found. For instance, 4g 1 (11; 11) → h 1 , 4g 1 (11; 12) → h 2 , · · · (B9) However, it should be pointed out that the Lagrangian in Eq. (5) is constructed in explicitly relativistic formalism, and all kinds of contraction of Lorentz indices are taken into account. Therefore, the relativistic Lagrangian has more terms than the HQSS Lagrangian. Taking the interaction of J/ψJ/ψ → ψ(2S)ψ(2S) for example, there are two terms accompanied by h 4 and h 4 in the relativistic Lagrangian but only one term in the HQSS one indicated by 4g 1 (11; 22). The merit of the use of relativistic formalism is that the original analytical properties of the obtained scattering amplitudes are respected.
Finally, it is also worth noting that the incorporation of the ψ(3700) in the HQSS Lagrangian demands the introduction of the ψ(1D) charmonium state. The ψ(3770) state is usually regarded as a mixture of ψ(1D) and ψ(2S) states, and the 1 3 D 1 charmonium component is generally considered to be predominant [61]. We refer the readers to Ref. [58] for how to introduce the charmonium states with L = 2. Nevertheless, in our relativistic effective Lagrangian given in Eq. (5), the ψ(3700) state is directly included as an explicit degree of freedom and its relevant interactions are constrained by symmetries.