Coupled-channel meson-meson scattering in the diabatic framework

We apply the diabatic framework, a QCD-based formalism for the unified study of quarkoniumlike systems in terms of heavy quark-antiquark and open-flavor meson-meson components, to the description of coupled-channel meson-meson scattering. For this purpose, we first introduce a numerical scheme to find the solutions of the diabatic Schr\"odinger equation for energies in the continuum, then we derive a general formula for calculating the meson-meson scattering amplitudes from these solutions. We thus obtain a completely nonperturbative procedure for the calculation of open-flavor meson-meson scattering cross sections from the diabatic potential, which is directly connected to lattice QCD calculations. A comprehensive analysis of various elastic cross sections for open-charm and open-bottom meson-meson pairs is performed in a wide range of the center-of-mass energies. The relevant structures are identified, showing a spectrum of quasi-conventional and unconventional quarkoniumlike states. In addition to the customary Breit-Wigner peaks we obtain nontrivial structures such as threshold cusps and minimums. Finally, our results are compared with existing data and with results from our previous bound-state--based analysis, finding full compatibility with both.


I. INTRODUCTION
Ever since the discovery of the χ c1 (3872) back in 2003 [1], many quarkoniumlike meson states whose properties defy the conventional description as heavy quarkantiquark bound states have been discovered [2]. The mass of these unconventional states generally lies close below or above the lowest open-flavor threshold with the same quantum numbers. This suggests a possible relevant role of open-flavor meson-meson components. As a matter of fact many theoretical descriptions incorporating meson-meson degrees of freedom have been attempted. These include (but are not limited to) effective field theories, molecular models, and "unquenched" quark models. We refer the reader to [3][4][5][6][7][8][9][10][11][12][13][14] and references therein for a more comprehensive review on the present theoretical and experimental landscape.
In recent years, new theoretical input from lattice QCD has sparked further progress. In particular, the static energy levels for a bottomonium bb configuration mixing with one or two open-bottom meson-meson configurations have been calculated in lattice QCD [15,16], and significant progress has been made in the understanding of the properties of some charmoniumlike states near open-charm meson-meson thresholds [17,18]. These results have been used to study the effect of up to two spin-averaged thresholds on the bottomonium spectrum [19,20]. A more general framework for the interaction of a QQ pair (where Q stands for a heavy quark, b or c) with an arbitrary number of meson-meson channels has been obtained through the implementation of the diabatic approach in QCD [21][22][23]. This type of approach, first developed in molecular physics [24], has been applied * roberto.bruschini@ific.uv.es † pedro.gonzalez@uv.es to calculate the mass spectrum and OZI-allowed strong decay widths of charmoniumlike [21,22] and bottomoniumlike [23] states. In the diabatic framework the dynamics, including the Q-Q interaction and the QQ-meson-meson mixing induced by string breaking, is completely described by a potential matrix whose elements are directly related to the static energy levels calculated in lattice QCD. This diabatic potential matrix is then plugged into a multichannel Schrödinger equation involving all the QQ and mesonmeson components for a given set of J P C quantum numbers. Then, the solutions of this diabatic Schrödinger equation allow for the description of quarkoniumlike meson states within the same J P C family [21]. Specifically, bound state solutions for energies below the lowest openflavor J P C threshold can be directly assigned to quarkoniumlike meson states. In contrast, solutions for energies above one or more thresholds, containing as many free-wave meson-meson components, cannot be directly assigned to physical mesons. Indeed, a dedicated formalism has to be developed in order to obtain a proper physical description.
As an alternative, one may try to avoid this difficulty following a bound state approximation even for energies above threshold. This is, quarkoniumlike meson states in this energy region can be assigned to bound state solutions of a reduced set of Schrödinger equations where the coupling with open thresholds is neglected. Then, decay widths and mass shifts induced from the coupling of these bound states with the continuum of free meson-meson states can be evaluated using a procedure well known in nuclear physics and hadron spectroscopy [25,26]. Actually, we have previously followed this approximation for a spectral analysis of charmoniumlike and bottomoniumlike mesons [21][22][23]. Notwithstanding the good results obtained, this approximation presents some shortcomings. On the one hand, there is an ambiguity regarding the number of bound states for a given J P C related to the possibility of choosing different sets of neglected thresholds when calculating some of them. In [21][22][23] this ambiguity was obviated ad hoc by assuming a oneto-one correspondence between calculated bound states and Cornell states. On the other hand, mixing between different QQ partial waves induced through their coupling to the same meson-meson thresholds is not properly taken into account in the case of open thresholds. Finally, there may be physical mesons escaping the bound state approximation.
In order to overcome these shortcomings we develop in this article a formalism to describe quarkoniumlike mesons from the continuum solutions of the diabatic Schrödinger equation above threshold. More precisely, we solve numerically the diabatic Schrödinger equation for energies above threshold. Then we decompose the asymptotic limit of each solution, which consists of one or more free waves in various meson-meson channels, as a superposition of stationary meson-meson scattering states. From this decomposition we obtain the on-shell S-matrix for the coupled-channel meson-meson scattering process. Finally, we identify quarkoniumlike meson states as resonances in the calculated scattering cross sections as a function of the energy. From this analysis we recover the quarkoniumlike meson states obtained in the bound state approximation as well as identify additional structures in the cross sections.
It should be pointed out that numerous treatments of the effects of coupled channels in the quarkoniumlike meson spectrum can be found in the literature, see for instance [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40]. A major difference of these treatments with respect to our diabatic approach is the use of a phenomenological mixing potential (mostly through a 3 P 0 model) instead of an interaction based on lattice QCD. Furthermore in some of these treatments a perturbative approach, whose validity is questionable (see Appendix in [22]), has been employed. Such differences make rather difficult a direct comparison between these studies and ours. Somewhat closer to our approach is a quite recent first attempt to study J P C = 0 ++ , 2 ++ charmoniumlike resonances in coupled DD-D sDs scattering on the lattice [41]. However, the several simplifying assumptions needed for this first lattice investigation of the coupledchannel system also prevent a direct comparison with our results.
The contents of this article are organized as follows. In Sec. II we briefly revisit the diabatic approach in QCD and write down the diabatic Schrödinger equation. Section III treats the numerical method used to solve the diabatic Schrödinger equation for continuum energies above the lowest open-flavor threshold. The theoretical framework used to describe coupled-channel mesonmeson scattering is detailed in Sec. IV, where we derive the formula for calculating the S-matrix. The calculated elastic cross sections are briefly discussed in Sec. V. Then in Sec. VI we summarize our main findings.

II. DIABATIC SCHRÖDINGER EQUATION
The diabatic approach in QCD is a recently proposed theoretical framework aiming at a unified description of both conventional and unconventional quarkoniumlike states based on lattice QCD results. Its construction, as well as its differences from other commonly used approaches in heavy-meson spectroscopy, has been explained extensively in [21]. In this section we briefly recap its main features.
The dynamics is governed by the diabatic Schrödinger equation where |Ψ is the state of the system containing one QQ and N meson-meson components M (i) 1M (i) 2 with i = 1, . . . , N , E the c.m. energy, K the kinetic energy operator, and V the diabatic potential operator. If we conveniently represent |Ψ as a column vector (1) . . .
where |ψ (0) and |ψ (i) are the states of the QQ and M components respectively, then we can represent the kinetic energy operator as a matrix (1) . . .
with µ (0) and µ (i) respectively the reduced mass of the QQ and M (i) 1M (i) 2 component and p 2 the squared relative momentum operator. Accordingly, in this notation the diabatic potential operator reads where vanishing matrix elements have been omitted. By projecting on the relative position space we have where the wave function r|Ψ is written as with ψ (0) (r) ≡ r|ψ (0) and ψ (i) (r) ≡ r|ψ (i) respectively the wave function of the QQ and M It is important to realize that, since we do not project explicitly on spin, each component of the wave function (6) is intended as a spin vector. In practice, since the diabatic potential must incorporate the symmetries of QCD, we impose J P C conservation. Therefore each solution to the diabatic Schrödinger equation is labeled by its J P C and m J quantum numbers, with J and m J the total angular momentum and its projection, P the parity, and C the charge-conjugation parity. Then the wave function components are more conveniently expressed in terms of the irreducible tensors of order J: where r = rr, C m l ,ms,m J l,s,J are Clebsch-Gordan coefficients, Y m l l (r) spherical harmonics, and ξ ms s spin vectors. Concretely, the wave function components are expanded as where t and k label the distinct (l (0) , s (0) ) and (l (i) , s (i) ) partial waves coupling to J P C in the QQ and M As for the potential energy the diagonal term V (0) , describing the Q-Q interaction, is directly connected to quenched lattice QCD results for the static light field energy [42]. It is represented by the kernel where V C (r) corresponds to the Cornell potential V C (r) = σr − χ r + m Q + mQ − β (13) with σ, χ, m Q (mQ), and β respectively the string tension, color Coulomb strength, heavy quark (antiquark) mass, and a constant. The other diagonal terms read with T (i) the threshold value where m M (i) 1 and mM(i) 2 are the masses of the corresponding mesons.
As for the nonvanishing offdiagonal elements V  (i) 2 mixing due to string breaking, they can be determined ab initio from the static energy levels calculated in lattice QCD. This is done by expanding their kernel in J P C -conserving tensors, where for simplicity we have assumed the radial mixing potentials V (i) mix (r) to be independent of the l, s, and J P C quantum numbers. Thus, using Eqs. (11)-(16), the diabatic Schrödinger equation (5) reduces to a multichannel radial Schrödinger equation with a radial diabatic potential matrix written in terms of V C (r), T (i) , and V (i) mix (r). It is then possible to determine V (i) mix (r) by virtue of the direct correspondence between the static energy levels and the eigenvalues of the diabatic potential matrix [24]. Our educated guess [21], where the factors ∆ (i) and ρ (i) respectively represent the strength and mixing length scale of the mixing, is based on the static energy levels recently calculated in unquenched lattice QCD [15,16]. The physical reasons behind this parametrization have been detailed in [21], which we refer to. It may be instructive to derive explicitly the multichannel radial equation in the simplest possible example of a QQ component interacting with a single meson-meson component M (1) 2 , each with only one partial wave contributing to J P C . In this case Eq. (5) reads . (18) If we now substitute Eqs. (9) and (10), multiply by Y J,m J (r) , and integrate over the solid angler, we get the multichannel radial equation  where we have introduced the reduced radial wave functions u (0) The generalization to any number of partial waves and an arbitrary number of meson-meson components is straightforward by considering each u It has to be remarked that by using Eqs. (4), (14), and (15) we are neglecting all meson-meson interactions from light meson exchange depicted in Fig. 1a. Instead, the diabatic potential incorporates a meson-meson interaction mediated by QQ depicted in Fig. 1b. Then our approximation could be regarded as the assumption that other interactions (like light meson exchange potentials) are subdominant with respect to the diabatic QQ-mediated ones. Although similar approximations for the potential matrix are commonly used in this context (compare for example with Eq. (7) of Ref. [16]), this assumption may be debatable. In fact there have been effective field theory studies, see for instance [43,44], where the dynamical energy levels have been described through a dominant contact meson-meson interaction. Therefore our treatment should be considered as an exploratory research on the effect of the mixing term. One should keep in mind though that additional contributions to the potential might be implicitly taken into account through the effective values of the parameters in the mixing potential.
For energies below the lowest J P C threshold the spectrum of solutions to the diabatic Schrödinger equation consists in a certain number of discrete energies with wave functions representing properly normalizable states.
Then, a discrete solution with energy E can be assigned to a quarkoniumlike meson with mass E. Concretely, solutions for energies far below the lowest threshold, containing a very dominant QQ component, differ very little from the quark model solutions obtained from the Cornell potential. They correspond to the lowest-lying conventional quarkonium (charmonium and bottomonium) states observed experimentally. In contrast, solutions with a relevant meson-meson component may appear close below the lowest threshold, as may be the case for the charmoniumlike meson χ c1 (3872) [21].
For energies above one or more thresholds the solutions, containing one or more free-wave meson-meson components, cannot be directly assigned to physical mesons. Instead, they have a more natural interpretation in terms of stationary meson-meson scattering states. More precisely, the diabatic Schrödinger equation for any energy E above threshold can be used to describe the J P C contribution to a meson-meson scattering process at c.m. energy E, where the meson-meson interaction is mediated by QQ (see Fig. 1b). From this interpretation, quarkoniumlike mesons observed experimentally can be properly assigned to meson-meson scattering resonances.

III. NUMERICAL SOLUTION ABOVE THRESHOLD
In this section we briefly revisit the finite difference method for the solution of the diabatic Schrödinger equation and show how the imposition of appropriate boundary conditions can be used to obtain a set of independent solutions for each energy above threshold.
In order to keep the discussion both simple and general, we start by examining the simplest case of a single radial Schrödinger equation where V (r) is some spherical potential, u(r) a reduced radial wave function, µ a reduced mass, and l a relative orbital angular momentum. The finite difference method consists in discretizing the radial configuration space in a lattice of equally spaced points r n . This gives rise to a matrix Schrödinger equation whose solution yields the numerical wave function u n ≡ u(r n ).
Discretization of the centrifugal and potential energy terms in (20) is straightforward, while that of the second derivative proceeds in two steps. The first step is to introduce the O(d 2 ) approximation for the second derivative where d is the discretization step. Notice however that this expression cannot be applied at the origin and extreme of the lattice, given respectively by r 0 = 0 and r extr ≡ r N+1 = (N + 1)d where N is the number of points in its interior. Thus the second step is to impose for these points Dirichlet boundary conditions u 0 = b 0 and u N+1 = b extr , so the second derivative for n = 1 and n = N reads while for n = 2, . . . , N − 1 it is given by Eq. (21). The discretized spherical Schrödinger equation is then written in matrix form as where u is the numerical reduced wave function (without the boundary values) H is the Hamiltonian matrix with K the tridiagonal kinetic energy matrix 2µr 2 , and B a constant numerical vector related to the boundary conditions, whose only nonzero components are in general the first and the last one: If the potential V (r) is either regular or diverges at most as r −1 for r → 0, as it is always in our case, the physical boundary condition at the origin is b 0 = 0. As for the boundary condition at r extr , assuming r extr to be very large, there are two distinct cases: These two scenarios lead to very different treatments of the numerical Schrödinger equation.
In the first case we have B = 0 and Eq. (24) reduces to a secular equation for the Hamiltonian matrix H, whose solution provides the bound state spectrum and the corresponding wave functions.
In the second case the energy spectrum is known a priori, being the continuum E ≥ lim r→∞ V (r), but there is no physical criterion to fix the boundary condition at the extreme. Then for each continuum energy E the wave functions are calculated solving the nonhomogeneous Eq. (24) with an arbitrary nonzero boundary condition b extr = 0. Indeed, as we have fixed the boundary condition at r 0 there is only one linearly independent solution for each value of E, therefore different choices of b extr will produce the same wave function up to a global multiplicative factor. Note that the assumption b extr = 0 is not restrictive in numerical applications. In fact, since the energies yielding a nontrivial solution that vanishes at r extr constitute a discrete subset of the continuum (i.e., the energy levels of a particle in a spherical box), the chance of picking accidentally the exact numerical value of one such energy is negligible.
The procedure we have detailed for a single-channel equation can be easily extended to the solution of the diabatic Schrödinger equation. We begin by observing that if instead of (20) one has a multichannel radial Schrödinger equation (like Eq. (19) for example), then the discretization procedure yields a system of coupled numerical equations that can be rearranged in the form (24), with the tridiagonal matrix H substituted by a more general Hermitian banded matrix.
Regarding the boundary condition at r extr , note that each diabatic channel falls into one of the two aforementioned cases i) or ii) depending on the energy. In particular, QQ and closed meson-meson channels fall into category i) whereas open meson-meson channels fall into category ii). So for energies below the lowest threshold, where all meson-meson channels are closed, all reduced wave function components vanish at "infinity", which corresponds to bound state solutions. Otherwise above threshold the boundary condition at r extr becomes a numerical vector whose only nonzero block is that corresponding to open meson-meson channels, whereñ is the total number of partial waves coupling to J P C in the open meson-meson channels.
Finally, notice that one can choose up toñ linearly independent boundary conditions b open , yielding as many independent solutions with the same energy. Therefore, for an arbitrary continuum energy one is able to build a basis set of solutions to the diabatic Schrödinger equation by solving (24) with boundary conditions (30) for n linearly independent numerical vectors b open . In practice, we build and solve the nonhomogeneous linear system (24) in Python using the NumPy and SciPy libraries [45,46] and taking the canonical basis of Cñ as our set ofñ linearly independent vectors b open .
For the sake of completeness, we report here that in this paper we have used d = 10 −3 fm and r extr = 200 fm.
1 ) coupling to J P C .
If there were no mixing, V mix (r) = 0, the spectrum of the diabatic Schrödinger equation would decompose in a discrete spectrum of pure quarkonium states (ψ (1) (r) = 0) and a continuum E ≥ T of pure mesonmeson (ψ (0) (r) = 0) free states given by a normalization factor introduced to facilitate the connection with the scattering states (see below). Then, the asymptotic behavior of a continuum solution would be written as where we have defined the symbol for the asymptotic equality relation, meaning that the two sides are equal in the limit r → ∞, and used If now we introduce the mixing potential (17), the continuum solutions cease being pure meson-meson states as they acquire a QQ component as well. Notice though that as the mixing potential (17) goes to zero exponentially for r → ∞, the QQ behaves as a bound component vanishing very quickly in the asymptotic limit. Hence, these mixed free-bound states correspond to a mesonmeson pair interacting at short distances through the mixing with QQ. The asymptotic behavior of these states is known from elastic scattering theory, see for example Eq. (11.17) in [47], as ψ (1) where the effect of the short-range interaction amounts to the introduction of a phase shift η J P C ;1 with respect to the solution without mixing with QQ. Notice that here and in the following we omit bound components and focus instead on the asymptotic behavior of the open meson-meson components.
Let us now generalize to an arbitrary number of partial waves (l k ) coupling to J P C in the open channel. Recalling from the previous section that whenever there is more than one (partial-wave) channel there are many independent solutions with the same energy, the asymptotic behavior of any continuum solution can be conveniently expressed as (compare for example with Chap. X, Eqs. (12), (14) in [48]) where we have introduced the additional label h to distinguish independent solutions with the same energy. Notice that, as the diabatic mixing couples the various meson-meson partial waves to each other, the coefficients a J P C ;k;h give account of the flow of probability from one partial-wave channel to the others through the coupledchannel interaction.
In the most general case of an arbitrary number of meson-meson components, one must realize that there is a different number of open channels depending on the energy. Then, if E is the energy and n the number of open thresholds at that energy, E ≥ T (j) for j = 1, . . . , n, the asymptotic behavior of the solutions is given by the open meson-meson components ψ (j) where h = 1, . . . ,ñ withñ the total number of partial waves coupling to J P C in all the open meson-meson channels. For more than one open threshold, these solutions correspond to a meson-meson scattering process where the coupled-channel interaction is provided by the mixing of QQ with all the meson-meson components.

B. Meson-meson S-matrix
Let us now see how to extract the on-shell mesonmeson S-matrix from the solutions of the diabatic Schrödinger equation above threshold. As shown in the Appendix, the meson-meson scattering states from an initial open channel j , with quantum numbers J P C , m J , l (j ) , s (j ) , to a final open channel j can be expressed as where k and k are used to label distinct partial waves (l (j ) , s (j ) ) and (l (j) , s (j) ) coupling to J P C in the initial and final channel respectively, and f j←j J P C ;k,k is the corresponding partial-wave scattering amplitude. Notice that the normalization factors have been chosen according to the state normalization by energy instead of momentum. This is motivated by the fact that the open channels in general differ in their threshold and reduced mass, which makes normalization by momentum unpractical as there is a different value of the relative meson-meson momentum for each open channel. It is important to realize that there are in totalñ independent scattering states (one for each partial wave k coupling to J P C ), that is as many as independent solutions to the diabatic Schrödinger equation at the same energy. Hence Eqs. (36) and (37) give two equivalent representations of the space of solutions at energy E, related to each other through a change of basis transformation. The change of basis transformation is written as where Γ We can further simplify Eqs. (40) and (41) using matrix notation. For the sake of simplicity, let us rearrange momentarily the indices (j, k), (j , k ) as, = 1, . . . ,ñ, so that each or specifies uniquely a partial-wave channel for the scattering process. In this way we can drop the subscripts k, k and rewrite Eqs. (40) and (41) as h a () and respectively. Let us now introduce theñ×ñ Jost matrices F ± J P C , where labels the rows and h the columns, and theñ ×ñ change of basis matrix Γ J P C , with h and the row and column index respectively. Then Eq. (42) can be rewritten in matrix notation as where 1 is theñ-dimensional identity matrix, meaning that the change of basis matrix Γ J P C is just the inverse of F − J P C . Plugging the change of basis matrix back into (43) yields the scattering amplitude from which, recalling that J P C conservation implies that the S-matrix is block-diagonal one can recognize the J P C block of the S-matrix as Equation (49) is a general formula known in multichannel scattering theory, see for example Eq. (20.18) in [47] (notice that we lack the momentum factors because we adopt a different normalization), that allows us to calculate the on-shell S-matrix numerically from the solution of the diabatic Schrödinger equation above threshold. Indeed, the numerical values of the amplitudes a (j) J P C ;k;h and phase shifts η (j) J P C ;k;h , defining the Jost matrices (44), can be obtained simply by fitting Eq. (36) to the long-distance numerical wave functions from Sec. III.
Finally, it is shown in the Appendix that the total unpolarized cross-section can be calculated as with each J P C cross section σ j←j J P C obtained as 2 , and we have restored the notation (j, k), (j , k ) for the partial-wave channels.

V. RESULTS
In this section we identify quarkoniumlike mesons from the structures in the calculated J P C cross sections as a function of the energy. More precisely, we do not use the total cross section (50) but rather a scaled J P C cross section defined as This scaled J P C cross section is more convenient for our theoretical analysis for a number of reasons: • it makes easier the distinction of resonances; • it is a dimensionless quantity; • it is not affected by the purely kinematical (p (j) ) −2 behavior of the cross section, what allows for a more detailed study for energies close above threshold; • it is symmetric under exchanges j ↔ j , as per the detailed balance principle (i.e., time reversal symmetry); • its values are bounded by the unitarity condition j,j σ j←j J P C ≤ 1, where the maximum value of 1 is expected at the mass of an isolated resonance with little nonresonant background.
We restrict our study to meson-meson pairs with pretty small widths and thresholds well separated in energy. In this manner we avoid the technical complications deriving from the treatment of the widths and the possible presence of overlaping thresholds. This constrains our formalism to be safely applicable only up to energies of 4.1 GeV in the charmoniumlike sector, see [21,22], and 10.8 GeV in the bottomoniumlike one, see [23]. For the sake of simplicity we shall discuss here only elastic processes, as they are the most convenient for the current theoretical analysis.
For practical purposes we shall compare the scattering resonances with existing quarkoniumlike states and with our former predictions from the bound state approximation. For this comparison to be meaningful we use the same parameters as in Refs. [21][22][23], that we list here for completeness.
For the Cornell potential we take the standard and for the heavy quark masses we use As for the mixing potential (17), in order not to spoil our predictive power we choose universal values ρ and ∆ for all thresholds. Moreover, we adopt the prescription of taking ∆ as the effective mixing strength of QQ with the isospin-singlet combination of two approximately degenerate thresholds containing up and down quarks. Then the effective mixing strength of QQ with a hidden-strange threshold is given by ∆/ √ 2 (see Appendix in Ref. [23] for more details). We use for charmoniumlike systems and for bottomoniumlike ones. ii) There appear additional resonances that do not correspond to the energy of any Cornell state. These resonances, with a significant meson-meson component, are located close to some meson-meson threshold coupling to J P C .
iii) Superposition of different resonances and their overlap with thresholds may result in the creation of complex structures in the cross section which deviate from the customary Breit-Wigner form.
The emerging physical picture is that of a quarkoniumlike spectrum consisting in quasi-conventional plus unconventional states from the meson-meson interaction provided by the diabatic mixing of Fig. 1b. The appearance of the unconventional states, as well as their eventual composition, depends on the strength of the QQ-meson-meson mixing and the position of the Cornell potential bound states with respect to the meson-meson thresholds coupling to J P C .
In what follows we simply identify the enhancements in the calculated open-charm and open-bottom cross sections which can be associated to resonances. A more thorough analysis, especially of the most complex structures such as threshold cusps and minimums, shall be the subject of a forthcoming paper.

A. Elastic cross section for open-charm mesons
The calculated scaled cross section for elastic opencharm meson-meson scattering processes with J P C = (0, 1, 2) ++ and J P C = 1 −− are shown in Fig. 2.
Visual comparison of these plots with the masses and widths reported in Table V of Ref. [22] provides us with an independent check of the results obtained in the bound state approximation. So the presence of a narrow resonance with J P C = 0 ++ close below the D sDs threshold with a mass of 3920.9 MeV, coming from the interaction of the D sDs threshold with the 2P Cornell state, is confirmed from the calculated peak in the DD cross section. Another peak, without correspondence in the bound state approximation, is visible in the D sDs cross section close above threshold, which may include some enhancement due to the presence of the resonance close below threshold. This second peak at about 3950 MeV can be partly attributed to a quasi-conventional 0 ++ resonance with a very dominant 2P cc component coming as well from the interaction of the 2P Cornell state at 3953.7 MeV with the D sDs threshold. Notice that this state was also present in the bound state approximation when neglecting the (open) D sDs threshold, but it was discarded because a one-to-one correspondence with Cornell states was assumed and the 2P cc core had already been assigned to the calculated bound state at 3920.9 MeV [22]. These two peaks may be related to the experimental candidates X(3915) (see below) and the more uncertain χ c0 (3860). Indeed, it is quite possible that due to the proximity in energy of the peaks their effect is present in both candidates. If this is the case, dedicated experiments would be required to unveil their presence. In this regard it is important to emphasize that these peaks are more easily distinguished in the scaled cross section (52) than in the physical one (51).
Apart from these peaks a tiny bump at about 4030 MeV can be seen in the D * D * cross section. It might correspond to a new resonance due to the interaction between the D * D * threshold and the distant 2P Cornell state, but its experimental observation in the physical cross section seems unlikely. (Even more difficult should be the observation of the other tiny structure at the DD threshold.) For J P C = 2 ++ the calculated DD elastic cross section shows two peaks, a broader one at about 3.9 GeV, close above the DD * threshold, and a narrower one at about 4 GeV, close below the D * D * threshold. They correspond to the two excited 2 ++ resonances at 3881.1 MeV and 4003.9 MeV calculated in the bound state approximation, the broader (narrower) one resulting mainly from the threshold interactions with the 2P (1F ) Cornell state. Notice that the narrower peak around 4 GeV is also visible in the DD * and, to a much lesser extent, in the D sDs cross sections. From the experimental point of view there is a 2 ++ candidate, the χ c2 (3930), which may be assigned to the broader peak.
It may be worth mentioning that there is a tiny enhancement of the DD * cross section on the background of the broader peak, near 3.95 GeV, that can be associated to a quasi-conventional resonance with a very dominant 2P cc component. Moreover, a quasi-conventional resonance with a very dominant 1F cc component can be inferred from the bump in the D * D * cross section just above threshold. Neither of these resonances find correspondence in the bound state approximation since the 2P and 1F cc cores were assigned to the bound states corresponding to the peaks near 3.9 GeV and 4 GeV respectively.
For J P C = 1 ++ two peaks are clearly visible in the DD * cross section, one very close to threshold, which can be correlated to the presence of the χ c1 (3872) state close below threshold, and another one at about 3.95 GeV, corresponding to a quasi-conventional 2P cc resonance, that has no correspondence in the bound state approximation. This last state was discarded in Ref. [22] because the 2P cc core had already been assigned to the calculated bound state at 3871.7 MeV. At present, there is no suitable experimental candidate in the PDG [2] for this 1 ++ resonance at about 3.95 GeV. We note however that recent studies of e + e − → γ ωJ/ψ from the BESIII Collaboration [49] show that a good fit to data is obtained by introducing either two or three resonant structures. In the three-resonance fitting scenario, along with the X(3872) and X(3915) there is an additional resonance, labeled X(3960), at 3963.7 ± 5.5 MeV (see Table I of that reference). Then our calculated 1 ++ peak at 3.95 GeV may be tentatively identified with X(3960). Moreover, notice that according to [49] the fitting including X(3872), X(3915), and X(3960) implies a width for X(3915) much smaller than the current PDG average value, opening an alternative interpretation of this state as the experimental counterpart of the very narrow 0 ++ calculated peak at 3920.9 MeV (see above). We strongly encourage further experimental efforts in establishing the possible existence of a 1 ++ charmoniumlike resonance near 3.95 GeV as well as the quantum numbers of the X(3915), what could provide a crucial test of our predictions.
For J P C = 1 −− a quasi-conventional resonance around 3.77 GeV is clearly visible in the calculated DD cross section, in perfect correspondence with the state of mass 3771.7 MeV (with a very dominant 1D cc component) obtained in the bound state approximation and with the experimental ψ(3770). On the other hand the D * D * cross section shows a peak at 4080 MeV, close below the D sD * s threshold, which may be assigned to the experimental ψ(4040). This is particularly relevant because, as explained in Ref. [22], the bound state approximation is not suited for the description of this resonance. In this regard it is worth mentioning that the calculated bump is not a standard Breit-Wigner curve. This complex structure seems to include the effect of a quasi-conventional 3S cc resonance located around 4097 MeV, and maybe also some effect from a third resonance above 4.1 GeV (as an extended numerical calculation beyond 4.1 GeV seems to point out). The minimum around 4.14 GeV in the calculated cross sections has to do with this complexity.
In addition to the aforementioned minimum, the calculated open-charm meson-meson scaled cross sections reveal the presence of other peculiar structures such as the threshold cusps in the 1 −− D sDs and DD * cross sections at the D * D * threshold. Threshold cusps have been studied extensively in the literature, see for example [50] and references therein, but it is not clear a priori which thresholds may produce a nontrivial structure. In this respect, the diabatic formalism may be an ideal framework to study their occurrence as it provides a unified description for energies across various thresholds. An analysis of the calculated threshold cusps in the charm as well as in the bottom sectors will be given in a separated paper.

B. Elastic cross section for open-bottom mesons
The calculated scaled cross section for elastic openbottom meson-meson scattering processes with J P C = (0, 1, 2) ++ and J P C = 1 −− are shown in Fig. 3.
As in the charmoniumlike case, the calculated cross section confirm the bottomoniumlike resonances calculated within the bound state approximation, see Table VIII of Ref. [23]. We see that for 0 ++ there is a peak at about the energy of the 4P Cornell state at 10782.2 MeV, corresponding to a quasi-conventional 4P bb resonance (recall that the nearby B sB * s threshold does not couple to 0 ++ ). This peak, much more visible in the B * B * than in the BB and B sBs cross sections, corresponds to the bound state predicted at 10778.1 MeV. Although no experimental candidate is known at present, it has to be remarked that this resonance shows up as a standard Breit-Wigner peak since it is quite isolated from any other structure. This may facilitate the experimental discovery of the corresponding bottomiumlike resonance, notwithstanding its relatively small width. Actually, apart from this peak there is only a tiny cusp in the B * B * cross section at the B sBs threshold and a small bump in BB just above threshold which may be due to the interaction of BB with the 3P Cornell state at 10536.6 MeV.
For 2 ++ there are two relevant structures: a Breit-Wigner peak in the BB cross section, somewhat below the BB * threshold, and a very broad and asymmetric bump (most visible in the B * B * cross section) reaching its maximum around the B sB * s threshold. The first peak corresponds to a quasi-conventional resonance, with a very dominant 2F bb component, predicted in the bound state approximation (the effect of the BB * threshold being just a small shift in the mass of the 2F bb Cornell state). On the other hand, the second more complex structure may be involving the effect of more than one resonance: the one at 10782.3 MeV from the bound state approximation and some other resonance above 10.8 GeV, as the presence of a minimum might be indicating (and an extended numerical calculation beyond 10.8 GeV seems to confirm). It may also include some effect from the interaction of the B * B * threshold with the 2F Cornell state at 10601.0 MeV. This complexity may make difficult the identification of the underlying resonance(s). In contrast, we expect the Breit-Wigner peak in BB close below the BB * threshold to be clearly identified when data are available.
For 1 ++ there is a peak close below the B sB * s threshold that is clearly visible in the B * B * and BB * cross sections, corresponding to the resonance at 10778.9 MeV from the bound state approach. In addition to this peak there appears an enhancement in the B sB * s cross section just above the corresponding threshold, at about the mass of the 4P Cornell state (10782.2 MeV). The presence of the quasi-conventional 4P bb resonance in this structure is somehow diluted within the enhancement of the B sB * s cross section due to the presence of the resonance close below threshold. In [23] we suggested to look experimentally for a narrow resonance at about 10780 MeV, which is in perfect agreement with the sharp peak calculated just below the B sB * s threshold. We realize now that the presence of the second peak just above the B sB * s threshold could hinder its experimental detection. In this regard, the expected decay of the narrow resonance to Υ(1S)φ could be determinant for its disentanglement from data.
For 1 −− the two prominent peaks observed in the figure are in perfect agreement with the bound state approximation prediction. The first peak, close below the BB * threshold, comes mostly from the interaction of the BB * threshold with the 4S bb state at 10615.0 MeV. Also related to this interaction there is a bump in the BB * cross section close above threshold which is associated to a quasi-conventional 4S bb resonance. As for the second prominent peak it can be associated to a quasi-conventional 3D bb resonance. Experimentally, the prominent peaks can be assigned to the well established Υ(10580) and the not well established Υ(10753) respectively, whereas the effect of the smaller bump is expected to be diluted within the data assigned to Υ(10580). Notice that both prominent peaks are slightly asymmetric, possibly indicating additional contributions from a threshold interaction with the 3D Cornell state in the first peak and with the 4S Cornell state in the second one.
Altogether, the results for charmoniumlike and bottomoniumlike states show, as anticipated, a spectrum of quasi-conventional resonances in one-to-one correspondence with the Cornell spectrum plus additional unconventional resonances located close to some meson-meson thresholds.
It is worth remarking that these results, which overcome the shortcomings of the bound state approximation, allow for a description of all existing candidates in the energy range under study as well as for definite prescriptions which could guide the experimental searches of yet unknown predicted resonances.

VI. SUMMARY
The diabatic framework, a QCD-based formalism for the description of quarkoniumlike systems in terms of QQ and open flavor meson-meson components, has been extended to the study of coupled-channel meson-meson scattering. This has allowed to overcome the shortcomings of the bound state approximation previously used for the analysis of quarkoniumlike states with masses above the lowest meson-meson threshold.
Our study of coupled-channel meson-meson scattering has gone through two separate steps. First, we have solved numerically the diabatic Schrödinger equation above threshold by using an extension of the well known grid method which allows for nonvanishing values of the wave function at the boundaries. Then, we have used the asymptotic behavior of these continuum solu- tions to extract the meson-meson scattering amplitudes and cross sections. This procedure has been applied to the calculation of the elastic cross section as a function of the energy for open-charm as well as open-bottom meson-meson scattering. Quarkoniumlike resonances have been identified from standard Breit-Wigner peaks and more general structures in the cross section. We have shown that quarkoniumlike resonances previously obtained from the bound state approximation constitute only a part (as a direct consequence of the constraints inherent to such approximation) of the whole spectrum. We find that this is formed by quasi-conventional resonances, with masses close to those of Cornell bound states, as well as by unconventional ones, with masses close to the energies of some meson-meson thresholds. It must be observed that sometimes the quasi-conventional resonances are not isolated from other enhancements in the cross section, which hinders their identification. All existing experimental candidates have been assigned to calculated states. As for predicted but yet undiscovered resonances, an initial analysis to guide experimental searches has been outlined. In this regard, the experimental confirmation of our predictions of a 1 ++ quasi-conventional charmoniumlike resonance with a mass about 3950 MeV, a narrow 0 ++ unconventional charmoniumlike resonance with a mass about 3920 MeV assigned to X(3915), and a 1 ++ unconventional bottomoniumlike resonance with a mass about 10780 MeV, would provide strong further support to our theoretical description.