Coupled-channel interpretation of the LHCb double-$J/\psi$ spectrum and hints of a new state near the $J/\psi J/\psi$ threshold

Recently, the LHCb Collaboration reported pronounced structures in the invariant mass spectrum of $J/\psi$-pairs produced in proton-proton collisions at the Large Hadron Collider. In this Letter, we argue that the data can be very well described within two variants of a coupled-channel approach employing $T$-matrices consistent with unitarity: (i) with just two channels, $J/\psi J/\psi$ and $\psi(2S)J/\psi$, as long as energy-dependent interactions in these channels are allowed, or (ii) with three channels $J/\psi J/\psi$, $\psi(2S)J/\psi$ and $\psi(3770)J/\psi$ with just constant contact interactions. Both formulations hint at the existence of a near-threshold state in the $J/\psi J/\psi$ system with the quantum numbers $J^{PC}=0^{++}$ or $2^{++}$, which we refer to as $X(6200)$. We suggest experimental tests to check the existence of this state and discuss what additional channels need to be studied experimentally to allow for distinctive tests between the two mechanisms proposed. If the molecular nature of the $X(6200)$, as hinted by the three-channel approach, is confirmed, many other double-quarkonium states should exist driven by the same binding mechanism. In particular, there should be an $\eta_c\eta_c$ molecule with a similar binding energy.

Introduction.-Quantum chromodynamics (QCD) is highly nonperturbative at low energies. As a result, how hadrons emerge from QCD and how the hadron spectrum is organized are still challenging open questions. The quest of exotic hadrons beyond the conventional quark model classification of quark-antiquark mesons and three-quark baryons has been one of the central issues in the study of nonperturbative QCD. In the past decades, dozens of new resonant structures with exotic properties were reported by various experiments in particular in the spectrum of hadrons containing at least one heavy-flavor (charm or bottom) quark. However, these observations brought even more challenges as they seem not to fit into a single uniform classification scheme, and various interpretations were proposed for each of them, see Refs. [1][2][3][4][5][6][7][8][9][10] for recent reviews of such exotic states. Recently, the LHCb Collaboration reported resonant structures in the double-J/ψ invariant mass distribution using data for pp collisions at the c.m. energies 7, 8, and 13 TeV [11]. The form of the signal is nontrivial, departing significantly from the expected phase space distribution as well as single and double-parton scattering: An enhancement in the near-double-J/ψ threshold region from 6.2 to 6.8 GeV is seen, which is followed by a narrow peak around 6.9 GeV. Between the broad bump and the narrow peak, there is a dip at around 6.8 GeV. The narrow peak is now dubbed X(6900), and has spurred a flood of model explanations [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. Naturally, a fully-charmed compact tetraquark resonance is the most straightforward candidate. How-ever, most of the theoretical studies indicate that the cccc ground state should have a mass lower than 6.9 GeV [29][30][31][32][33][34][35][36][37][38]. Furthermore, the 700 MeV energy gap between the double-J/ψ threshold and 6.9 GeV is larger than a typical energy gap between the ground and radially/orbitally excited states. Thus, lower states should exist, if there is a cccc resonance with a mass around 6.9 GeV. Due to a smaller phase space, such lighter states are expected to have smaller widths. However, there are no obvious narrower peaks in the reported double-J/ψ spectrum.
It is well-known that threshold effects play an important, sometimes crucial, role for the properties of hadrons residing above open-flavor thresholds. For example, there is always a cusp at an S-wave threshold due to the analytic structure of the two-body Green's function (for a review, see Ref. [8]). It may lead to either a peak or a dip, depending on the interference with other contributions. The visibility of the corresponding structure in the line shape depends on whether or not it is enhanced by a nearby pole in the amplitude [39]. Thus, in order to properly interpret the new observations it is important to understand the role played by various thresholds located nearby. There are quite a few double-charmonium channels with the thresholds below 7.2 GeV which can couple to the double-J/ψ system, such as η c η c , h c h c , χ cJ χ cJ (J, J = 0, 1, 2), η c η c (2S), ψ(2S)J/ψ and ψ(3770)J/ψ. In this work we assume that the interaction between the quarkonia is dominated by the exchange of light modes (soft gluons or, e.g., pion pairs). Then the coupling of the double-J/ψ to the η c η c or h c h c flips the charmquark spin, and is expected to be suppressed due to the heavy quark spin symmetry (HQSS). Indeed, HQSS implies that the interactions involving the spin of a heavy quark are suppressed as O(Λ QCD /m Q ) [40], where Λ QCD denotes the momentum scale where QCD gets nonperturbative and m Q is the heavy-quark mass. From the point of view of the meson-exchange picture, the lowest meson that can be exchanged for the coupling of the χ cJ χ cJ to the double-J/ψ, keeping the SU(3) flavor and isospin symmetries, is the ω. It is heavier than the f 0 (500) (or, effectively, two pions) that can be exchanged for the transitions J/ψJ/ψ → ψ(2S)J/ψ or ψ(3770)J/ψ to happen. In this regard, it is interesting to notice that indeed the dip prior to the X(6900) peak appears around the ψ(2S)J/ψ threshold at 6783 MeV. Therefore, from this phenomenological point of view, among the doublecharmonium channels, the ψ(2S)J/ψ and ψ(3770)J/ψ ones are expected to play the most crucial role in describing the double-J/ψ spectrum up to the energies covering the X(6900) peak.
In this Letter, we aim at constructing minimal coupledchannel models able to describe the LHCb data on the double-J/ψ invariant mass distribution in the energy interval from the double-J/ψ threshold to 7.2 GeV and studying their predictions for pole locations and line shapes in the other double-charmonium channels. In particular, we consider a two-(J/ψJ/ψ and ψ(2S)J/ψ) and three-channel (J/ψJ/ψ, ψ(2S)J/ψ, and ψ(3770)J/ψ) models and find that (i) both models provide a remarkably good description of the data which, therefore, do not allow one to distinguish between them, (ii) both models predict the existence of a near-threshold pole around 6.2 GeV (we call it the X(6200)) corresponding to a shallow bound or virtual J/ψJ/ψ state, (iii) the structure of the other, above-threshold poles appears to be very different for the two models considered and so are the predicted line shapes in the ψ(2S)J/ψ channel. We conclude, therefore, that the existence of the X(6200) is a robust consequence of the proposed coupled-channel approach, while additional measurements of the other double-charmonium channels are necessary in order to better understand the nature of the higher poles.
Coupled-channel model.-Contrary to earlier attempts to understand the role played by the relevant double-charmonium thresholds for the double-J/ψ spectrum [24], the key idea of our approach is to present a minimal model consistent with the data and able to extract the poles responsible for the structures in the data. We, therefore, confine ourselves to those double-charmonium channels which are consistent with HQSS, and constrain the T -matrix with unitarity and causality.
As detailed below we work with a separable potential V . Then the T -matrix of the coupled-channel system can be written as where E is the double-J/ψ center-of-mass energy and G is a diagonal matrix for the intermediate two-body propagators. We use the dimensionally regularized twopoint scalar loop function [41], where s = E 2 , m i1 and m i2 are the particle masses in the i-th channel, is the corresponding threemomentum with λ(x, y, z) = x 2 +y 2 +z 2 −2xy−2yz−2xz for the Källén triangle function. Here µ denotes the dimensional regularization scale, and a(µ) is a subtraction constant. The T -matrix in Eq. (1) respects the constraints of unitarity.
Since the narrow dip and peak near 6.9 GeV are around the ψ(2S)J/ψ and ψ(3770)J/ψ thresholds, respectively, we consider only S waves which are able to produce nontrivial threshold structures. For the two-channel model, the potential V is parameterized as where a 1,2 , b 1,2 , and c are real free parameters. The momentum dependence of the potential is necessary here to produce nontrivial structures above the higher threshold, since purely constant contact-term potential can only produce bound or virtual state poles below threshold. For the three-channel model, the potential V is a 3 × 3 matrix, where a ij 's are real parameters of the model, and, aiming at the simplest possible formulation of the model consistent with the data, we do not consider an explicit momentum dependence in this case. The production amplitude in the J/ψJ/ψ channel (labelled as channel 1) can be constructed as where T ij are the elements of the T -matrix in Eq. (1), the ratios r i mimic potentially different production mechanisms for different channels (r 3 = 0 for the 2-channel fit), and the parameter b = 1 accounts for violation of unitarity in the production mechanism which should be present in a 2-body treatment given the complexity of the inclusive reaction from which the data were extracted. To describe the details of the short-distance production encoded in the function P (E) above, we take it in an exponential form, P (E) = αe −βE 2 , and fix the slope parameter β = 0.0123 GeV −2 from fitting to the doubleparton scattering (DPS) distribution quoted in the LHCb paper [11]. The energy dependence of the production operator accounts for the fact that the double-J/ψ and ψ(2S)J/ψ two-particle systems can be produced at the parton level and interact before the final double-J/ψ particles are detected. The overall strength parameter α is treated as a free parameter of the model.
Fit results.-Before fitting the data we get rid of the parameters which weakly affect the distribution or can be recast into other constants. In particular, we set µ = 1 GeV and the subtraction constant in the loop function is fixed as a(µ = 1 GeV) = −3; its variance within Eq. (1) can be absorbed into the redefinition of the contact interactions in the potential. Also, we choose the r i -parameters in the amplitude (5) equal to 1 since the fit does not call for their different values. Two-channel model: The two-channel parameterisation has 7 parameters. These are {a 1 , a 2 , b 1 , b 2 , c, b, α}. The fit was performed with randomly chosen 2 × 10 4 sets of initial values of the parameters and constrained by causality to ensure that there are no pole on the first Riemann sheet of the complex energy except on the real axis below threshold (see, e.g., Ref. [42]). The χ 2 function is then minimised using the MINUIT algorithm [43][44][45]. The best fit describes the data remarkably well with χ 2 /dof = 0.99-see Fig. 1. Interestingly, although the fit was only performed up to 7.2 GeV, a good description of the data is achieved in the entire energy interval up to 9 GeV. In this model, the dip in the line shape is produced due to a destructive interference of the ψ(2S)J/ψ threshold cusp which emerges from a coupled-channel dynamics with the background, see Eq. (5). The abovethreshold narrow hump is due to the energy dependence of the two-channel potential (3) which leads to a nearby resonance pole. A detailed analysis of the poles is given below.
Three-channel model: The three-channel model has 8 real parameters, {a ij (i j), b, α}. Two fits of similar quality are found with χ 2 /dof = 0.97 (Fit 1) and χ 2 /dof = 1.05 (Fit 2). All parameters of these fits coincide within their 1σ uncertainty except a 22 . A comparison of both fits with the data is given in Fig. 2  in the two-channel model, the description of the data is remarkably good, including the (not fitted) large-energy tail up to 9 GeV. In this model, the nontrivial structures in the line shape at approximately 6.8 and 6.9 GeV are due to the effect from the ψ(2S)J/ψ and ψ(3770)J/ψ thresholds, amplified by a nearby pole.
Pole analysis.-To study the pole structure of the fitted T -matrix, we generate more than 300 parameter sets within the 1σ contours in the parameter space for all combinations of the fit parameters and find all poles of the amplitude from the near-threshold region up to 7.  We focus first on the mass region of the pronounced structures in the data and study the T -matrix poles in the complex energy plane. For each pole we quote RS ±±... in parentheses to indicate the Riemann sheet where the pole is located, with the subscript composed of the signs of Im k i in all coupled channels involved, from the lowest in energy to the highest [46]. We find that the pole locations are quite different for the different models employed. In particular, there exist two such poles for the two-channel model (hereinafter the pole positions are given in MeV), while, for the 3-channel fits, there is only one (badly determined) remote pole on RS −++ -see Fig. 4. Meanwhile, both models confidently predict a pole very near the J/ψJ/ψ threshold. The two-channel fit allows Obviously, further modifications of the model to extend the coupled-channel set or include higher-order terms in the potential cannot destroy this pole simply because such modifications would only affect the highenergy tail of the distribution, far away from the J/ψJ/ψ near-threshold region. We conclude, therefore, that the existence of a pole near the J/ψJ/ψ threshold is a robust consequence of the coupled-channel dynamics within the suggested approach. For definiteness, we name this state X(6200). Its quantum numbers are either 0 ++ ( 2S+1 L J = 1 S 0 ) or 2 ++ ( 5 S 2 ), as required to have an Swave threshold composed of two identical vector bosons. As these two partial waves cannot interfere within our coupled-channel approach, we do not consider both amplitudes simultaneously. The latter scenario would require additional parameters not called for by the data.
Further predictions and tests.-As one can see from Figs. 1 and 2, although the models used to analyse the data in the J/ψJ/ψ channel are based on a different dynamical content, they can provide a description of the data of a comparable quality. However, further predictions of these models differ substantially, allowing for a direct experimental discrimination (or falsification of the whole approach). As one of such tests we propose measurements of the line shapes in other double-charmonium channels. As a representative example, in Fig. 5 we show the predictions of the two models employed in this work for the invariant mass spectrum in the ψ(2S)J/ψ final state. Indeed, the models predict quite different spectra above the ψ(3770)J/ψ threshold, so that experimental data for this channel as well as for the ψ(3770)J/ψ one should help to better understand the physical origin of the structures reported by LHCb.
Also, a direct experimental confirmation or refutation of the existence of the X(6200) state is very important for a better understanding of the double-charmonium spectrum. In particular, a distinct signal from this state could be seen in the final states like J/ψµ + µ − , µ + µ + µ − µ − , and η c η c which can be studied at energies below the nominal J/ψJ/ψ threshold.
To better understand the nature of the X(6200), we estimate its compositeness,X A (X A = 1 for molecules andX A = 0 for compact states), which was introduced in Ref. [47] to characterize near-threshold bound states, to extract the S-wave scattering length a 0 and the effective range r 0 , and then usē The results presented in Table I imply that while the two-channel model supports the X(6200) as a compact state, the three-channel approach is better compatible with its molecular interpretation. If the latter is true, the same mechanisms (for example, the two-pion exchange) which drive the X(6200) can provide sufficient binding also in other double-charmonium channels, so that many more double-charmonium molecular states can exist near relevant thresholds. In particular, if the X(6200) is a molecule, HQSS predicts a near-threshold scalar doubleη c state bound by a similar mechanism [48,49], to be searched for in, e.g., the 2[KK]π final state.
Conclusions.-In this Letter, we demonstrated that the recent LHCb data on the double-J/ψ invariant mass spectrum are consistent with a coupled-channel description. The best fits to the data imply the existence of a state near the J/ψJ/ψ threshold which we called the X(6200). This state can have the quantum numbers of a scalar or a tensor. Further experimental tests are outlined to verify the hypothesis of the existence of this state and shed light on its nature. If confirmed, this discovery may start a new era in the spectroscopy of doublecharmonium and double-bottomonium states. It would also be valuable to simulate the double-J/ψ (or doubleη c ) scattering on the lattice.

Supplemental material
Here we quote the values of the parameters found from the fits discussed in the main text and the corresponding correlation matrices. The parameters (without bars) used in the matrices V 2ch and V 3ch in the main text are obtained from the parameters (with bars) given here by multiplying the latter by the factor 2m i , where m i 's are the involved charmonium masses for which we use [10] m J/ψ = 3.0969 GeV, m ψ(2S) = 3.6861 GeV. (12) In the convention used throughout the paper the Swave T -and S-matrix are related as where k i is the absolute value of the c.m. 3-momentum in the i-th channel.

Two-channel model
The fitted values of the parameters of the two-channel model are listed in Table II. Their correlation matrix reads (see Table II The moduli of the T -matrix elements from the 2channel fit are shown in Fig. 6.

Three-channel model
The fitted values of the parameters of the three-channel model are listed in Table III. Their correlation matrices read (see Table III for the parameters order) for Fit 1, and for Fit 2.
The moduli of the T -matrix elements from the 3channel fits are shown in Fig. 7