Diabatic description of charmonium-like mesons

We apply the diabatic formalism, first introduced in molecular physics, to the description of heavy-quark mesons. In this formalism the dynamics is completely described by a diabatic potential matrix whose elements can be derived from unquenched lattice QCD studies of string breaking. For energies far below the lowest open flavour meson-meson threshold, the resulting diabatic approach reduces to the well-known Born-Oppenheimer approximation where heavy-quark meson masses correspond to energy levels in an effective potential between quark and antiquark. For energies close below or above that threshold, where the Born-Oppenheimer approximation fails, this approach provides a set of coupled Schr\"{o}dinger equations incorporating meson-meson components nonperturbatively, i.e. beyond loop corrections. A spectral study of heavy mesons containing $c\bar{c}$ with masses below 4.1 GeV is carried out within this framework. From it a unified description of conventional as well as unconventional resonances comes out.


I. INTRODUCTION
The discovery of the χ c1 (3872) in 2003 [1] may be considered as the initio of a new era in heavy-quark meson spectroscopy. This resonance and a plethora of new states (ψ(4260), ψ(4360), X(3915), and many others, see [2]) discovered since then have masses and decay properties that do not correspond to the conventional heavy quark (Q) -heavy antiquark (Q) meson description, such as the one provided by nonrelativistic or semirelativistic quark models that has been so successful in the past [3][4][5]. A feature of any of these unconventional states is that its mass lies close below or above the lowest open flavour meson-meson threshold with the same quantum numbers. This suggests a possible relevant role of open flavour meson-meson thresholds in the explanation of the structure of the new states. As a matter of fact, the nonrelativistic Cornell quark model [3,4] incorporates some of these effects through meson loops where the interaction connecting QQ and open flavour meson-meson is derived from the QQ binding potential. Similar kind of loop contributions, with quark pair creation models like the 3 P 0 one providing the valence-continuum coupling, have been extensively studied in the literature (see for instance [6] and [7]). However, these perturbative loop contributions seem to be insufficient for a detailed description of the new structures. This has led to the building of phenomenological models involving implicit or explicit meson-meson components, for example in the forms of tetraquarks, meson molecules, and hadroquarkonium (see [8][9][10][11] and references therein).
Ab initio calculations from QCD have been also carried out. From lattice QCD, a Born-Oppenheimer (B-O) approximation for heavy-quark mesons has been developed [12] (for a connection with effective field theories see [13] * roberto.bruschini@ific.uv.es † pedro.gonzalez@uv.es and references therein). In this approximation, based on the large ratio of the heavy quark mass to the QCD energy scale associated with the gluon field, the heavyquark meson masses correspond to energy levels of a Schrödinger equation for QQ in an effective potential. This potential is defined by the energy of a stationary state of light-quark and gluon fields in the presence of static Q and Q sources, which is calculated in lattice QCD. Thus, conventional quarkonium masses are the energy levels in the ground state potential calculated in quenched (without light quarks) lattice QCD whose form is Cornell-like [14], whereas quarkonium hybrid (QQg bound state where g stands for a gluon) masses are energy levels in the quenched excited state potentials. Although no tetraquark potentials have been calculated yet from lattice QCD, some information on them has been also extracted [15]. The immediate question arising is whether these hybrid and tetraquark B-O potentials may correctly describe or not the new states. The answer to this question can be derived from [15], where an assignment of the masses of some of the new states to energy levels in these potentials has been pursued. In essence, quoting this reference, although the B-O approximation provides a starting point for a coherent description of the new states based firmly on QCD, a detailed description of them requires to go beyond quenched lattice calculations and beyond the B-O approximation.
An intermediate step in this direction was taken in [16] by identifying the unquenched lattice energy for static Q and Q sources, when the QQ configuration mixes with one or two open flavour meson-meson ones [17,18], with a QQ potential. This unquenched approximation allows for some physical understanding of threshold effects beyond hadron loops. However, the description in terms of effective QQ channels does not give detailed account of the configuration mixing.
In this article we take a step further to go beyond the B-O approximation. For this purpose we use the diabatic approach developed in molecular physics for tack-ling the configuration mixing problem (see for instance [19]). This allows us to establish a general framework for a unified description of conventional and unconventional heavy-quark meson states. This framework is applied to the calculation of J ++ and the low-lying 1 −− meson states with Q = c (charm quark) where there are sufficient data available to test its validity.
To our knowledge this is the first time that a complete treatment of heavy-quark meson states involving heavy quark-antiquark and meson-meson degrees of freedom incorporates the results from ab initio calculations in quenched and unquenched lattice QCD.
The contents of the paper are organized as follows. In Sec. II the mathematical formalism and the physical picture leading to the B-O approximation for heavy-quark mesons is revisited. In Sec. III we detail the diabatic approach and in Sec. IV we adapt it to the description of heavy-quark meson states. The application to meson states containing cc is detailed in Sec. V. For the sake of simplicity we consider states involving non-overlapping thresholds with small widths. The comparison of our results to existing data serves as a stringent test of our treatment. Finally, in Sec. VI our main conclusions are summarized.

II. BORN-OPPENHEIMER APPROXIMATION IN QCD
The Born-Oppenheimer (B-O) approximation was developed in 1927 for the description of molecules [20], and since then it has been a fundamental approximation in chemistry. More recently it has been employed for the description of heavy-quark meson bound states from QCD [12,15]. Next, we briefly recall the main steps in its construction for the description of a heavy-quark meson system containing a heavy quark-antiquark (QQ) interacting with light fields (gluons and light quarks), with Hamiltonian where K QQ is the QQ kinetic energy operator with µ QQ being the reduced QQ mass, p (P ) the QQ relative (total) three-momentum, and H lf QQ the part of the Hamiltonian containing the light field energy operator and the QQ -light-field interaction. Notice that H lf QQ depends on the Q and Q positions but does not contain any derivative with respect to the Q and Q coordinates.
A heavy-quark meson bound state |ψ is a solution of the characteristic equation where E is the energy of the state. Note that |ψ contains information on both the QQ and light fields.

A. Static limit
The first step in building the B-O approximation consists in solving the dynamics of the light fields by neglecting the QQ motion, i.e. setting the kinetic energy term K QQ equal to zero. This corresponds to the limit where Q and Q are infinitely massive, what can be justified because the Q and Q masses, m Q and m Q , are much bigger than the QCD scale Λ QCD , which is the energy scale associated with the light fields.
As we are interested in the internal structure of the system and this does not depend on the centre of mass motion (which coincides with the QQ centre of mass motion in the infinite mass limit) it is convenient to use the QQ relative position r = r Q − r Q , and work in the QQ centre of mass frame where P = 0.
In this static limit r is fixed, ceasing to be a dynamical variable. This is, the components of r can be considered as parameters, rather than operators, in the expression of H lf QQ that will depend operationally on the light fields only. We shall indicate this parametric dependence renaming H lf QQ as H lf static (r). It is then possible to solve the dynamics of the light fields for any value of r: where |ζ i (r) are the light field eigenstates, V i (r) the corresponding eigenvalues, and i stands for the set of quantum numbers labelling the eigenstates. Note that both the eigenvalues and the eigenstates depend parametrically on r, and that for every value of r the eigenstates {|ζ i (r) } form a complete orthonormal set for the light fields: As for the eigenvalues V i (r), they correspond to the energies of stationary states of the light fields in the presence of static Q and Q sources placed at a relative position r, and can be calculated ab initio in lattice QCD.
More precisely, in quenched (with gluon but not lightquark fields) lattice QCD [14] the ground state of the light fields is associated with a QQ configuration, and up to spin dependent terms that we shall not consider the static energy of this ground state mimics the form of the phenomenological Cornell potential with σ, χ and β standing for the string tension, the colour coulomb strength, and a constant fixing the origin of the potential respectively. On the other hand, unquenched (with gluon and lightquark fields) lattice QCD calculations [17,18] have shown that due to string breaking the association of the light field ground state with a QQ configuration holds only for small values of the relative QQ distance r ≡ |r|. When increasing r the QQ configuration mixes significantly with meson-meson configurations. More in detail: below (above) an open-flavour meson-meson threshold the energy of a stationary state of the light fields changes with r, from the one corresponding to the QQ (mesonmeson) configuration to the one of meson-meson (QQ) configuration, avoiding in this manner the crossing of the static light field energies corresponding to pure QQ and meson-meson configurations that would take place at the threshold mass in absence of string breaking. In Fig. 1 we have represented graphically this situation for QQ and one meson-meson threshold (the representation for two meson-meson thresholds can be seen in [17,18]).

B. Adiabatic expansion
Having solved the static problem for the light fields, the next step in the construction of the B-O approximation consists in reintroducing the QQ motion. This is done by solving the bound state equation where E denotes the mass of the bound state, making use of the so-called adiabatic expansion for |ψ : where |r ′ is a state indicating the QQ relative position and we have omitted spin degrees of freedom for simplicity. The qualifier "adiabatic" refers to the fact that each term in the expansion depends only on a single value of r ′ , what can be related to the physical situation where the light fields respond almost instantaneously to the motion of the quark and antiquark. However, as will be shown in what follows, this physical expansion is not mathematically convenient when configuration mixing takes place. Note that as the states |ζ i (r ′ ) depend on r ′ , so do the coefficients ψ i , one for each light field state.
Using (8) and multiplying on the left by r| the bound state equation can be rewritten as then multiplying on the left by ζ j (r)| yields The first term on the left hand side of (10) can be developed as with τ ji (r) ≡ ζ j (r)|∇ζ i (r) and τ being the so-called Non-Adiabatic Coupling Terms (NACTs) of the first and second order respectively. Furthermore, using ∇ ζ j (r)|ζ i (r) = ∇δ ji = 0 we have from which it follows so that (∇τ (r)) ji = ζ j (r) ∇ 2 ζ i (r) + ∇ζ j (r)|∇ζ i (r) = τ (2) ji (r) − (τ (r) 2 ) ji (15) and finally ζ j (r)|∇ 2 ψ i (r)|ζ i (r) = δ ji ∇ 2 ψ i (r) + 2τ ji (r) · ∇ψ i (r) + ((∇ · τ (r)) ji + (τ (r) 2 ) ji )ψ i (r) ≡ ((∇ + τ (r)) 2 ) ji ψ i (r). (16) The bound state equation (10) then reads This is a multichannel equation where ψ i (r) stands for the i-th component of the heavy-quark meson wave function, that is in general a mixing of QQ and mesonmeson components. Notice though that this is not the usual Schrödinger equation because of the presence of the NACTs τ inside the kinetic energy operator. These terms introduce a coupling between the wave function components and reflect the non-trivial interaction between the QQ motion and the light field states.

C. Single channel approximation
The last step in the construction of the B-O approximation consists in neglecting the NACTs inside the kinetic energy operator: This is called the single channel approximation because the bound state equation (17) then factorizes in a set of decoupled single channel Schrödinger equations where V j (r), corresponding to the energy of the stationary j-th state of the light fields in the presence of static Q and Q sources, plays the role of an effective potential. Eqs. (4), (8), (18) and (19) define the B-O approximation.
Notice that the single channel approximation can be deemed reasonable only up to QQ distances for which the NACTs can be neglected, i.e. for distances where the QQ and meson-meson configuration mixing associated with the light field eigenstates is negligible (for a specific calculation see Sec. IV C). This makes the B-O approximation to be justified only for bound state energies far below the lowest open flavour meson-meson threshold. In particular, conventional heavy-quark meson masses, far below the lowest open flavour meson-meson threshold, can be described as the energy levels in the potential corresponding to the quenched ground state of the light fields, i.e. the Cornell potential.

III. DIABATIC APPROACH
For energies close below or above an open flavour meson-meson threshold the mixing between the QQ and meson-meson configurations gives rise to nonvanishing NACTs, so that the single channel approximation (18) cannot be maintained. Instead, one has to deal with the set of coupled equations (17), which is not practicable for two reasons: i) There is no yet direct lattice QCD calculation of the NACTs τ .
ii) When τ = 0, the wave function components in the expansion (8) do not correspond to pure QQ or meson-meson but rather to a mixing of both, the amount of mixing depending on r.
These drawbacks can be overcome through the use of the diabatic approach, where one expands the bound state |ψ on a basis of light field eigenstates calculated at some fixed point r 0 . As the {|ζ i (r) } form a complete set for the light fields whatever the value of r, switching from a {|ζ i (r) } to {|ζ i (r 0 ) } is equivalent to a r-dependent change of basis in the light degrees of freedom.
The diabatic expansion of the bound state reads where the coefficients ψ i , one coefficient for each light field state, are functions of r ′ that depend parametrically on r 0 . A nice physical feature of this expansion is that the light field state |ζ i (r 0 ) corresponding to each component ψ i does not depend on the QQ relative position r ′ . This means that if one chooses the fixed point r 0 far from the avoided crossing, then the wave function components correspond to either pure QQ or meson-meson for any value of r ′ . In other words, in the diabatic approach one expands the bound states in terms of the more intuitive Fock components (pure QQ and pure meson-meson) instead of components which are a mixing of QQ and meson-meson.
Substituting (20) in the bound state equation (7) and projecting on r| yields where all the derivatives are taken with respect to r. If we now multiply on the left by ζ j (r 0 )|, as ∇ |ζ i (r 0 ) = 0 the equation reads is the so-called diabatic potential matrix.
The complete equivalence between Eqs. (17) and (22) has been shown elsewhere [19] and is reproduced, for the sake of completeness, in Appendix A. In short, the troublesome NACTs in (17) that break the single channel approximation when configuration mixing is present (thus invalidating the B-O framework) are taken into account in (22) through the diabatic potential matrix. This is utterly convenient since, as we shall see in Sec. IV B, the elements of this matrix are directly related to the static light field energy levels calculated in quenched and unquenched lattice QCD.
It is also easy to show that when the single channel approximation (18) holds the diabatic potential matrix (23) becomes a diagonal matrix containing the static light field energy levels calculated in quenched lattice QCD, and consequently Eq. (22) reproduces the set of single channel Schrödinger equations (19).
Therefore, the diabatic approach is a complete general framework appliable to conventional heavy-quark mesons lying far below the lowest open flavour mesonmeson threshold as well as to unconventional ones lying close below or above that threshold.

IV. HEAVY-QUARK MESONS IN THE DIABATIC FRAMEWORK
In order to apply the diabatic framework to the description of heavy-quark meson bound states we examine first the case of a single meson-meson threshold. Then we proceed to the generalization to an arbitrary number of thresholds.

A. Spectroscopic equations
Let us consider one meson-meson threshold. Let us fix a value for r 0 such that the ground state of the light fields is associated with the QQ configuration and the first excited state with the meson-meson one. To make this more clear we relabel the diabatic light field states as and the diabatic wave function components as Accordingly, we rename the diabatic potential matrix components (23) as Let us realize that having associated each component of the wave function with pure QQ or pure meson-meson, we can easily incorporate to the kinetic energy operator the fact that the reduced mass of the meson-meson component, µ M1M2 , is different from µ QQ . Hence, we shall use − Then, the bound state equations read or in matrix notation where K is the kinetic energy matrix and Ψ(r) is a column vector notation for the wave function: The multichannel Schrödinger equation (27), or equivalently (28), defines formally the diabatic approach for the description of the heavy-quark meson system.

B. Mixing potential
To solve (27) we need to know the diabatic potential matrix Eq. (26). Regarding the diagonal element V QQ (r), we see from (26a) that it corresponds to the expectation value of the static energy operator in the light field state associated with a pure QQ configuration. This can be identified with the ground state static energy calculated in quenched lattice QCD, see Fig. 1, given by the Cornell potential In the same way, from (26b) we identify the other diagonal term V M1M 2 (r) with the static energy associated with a pure meson-meson configuration, given by the threshold mass T M1M 2 (the sum of the meson masses) up to one pion exchange effects that we do not consider here.
As for the off-diagonal term, the mixing potential V mix (r), we can use the eigenvalues of the diabatic potential matrix to derive its form. As shown in Appendix A, these eigenvalues correspond to the static energy levels that are calculated in unquenched lattice QCD which have been pictorially represented in Fig. 1. More precisely, the eigenvalues of the diabatic potential matrix are the two solutions V ± (r) of the secular equation where I is the identity matrix. These solutions read from which we obtain where we have dropped the vector notation for r as the energy levels calculated in lattice QCD depend only on the modulus r = |r|. Eq. (36) tells us that a detailed calculation of the mixing potential |V mix (r)| from ab initio lattice data on V ± (r) is possible. As a matter of fact, an effective parametrization of V mix (r) from lattice data has been proposed [18]. However, a general parametrization at all distances, which is indispensable to solve the multichannel Schrödinger equation (27), is still lacking. While we encourage work along this direction, we must resort to general arguments to get the shape of |V mix (r)|. In this regard, the general form of the curves V + (r) and V − (r) near any threshold, reflecting the physical picture of the QQ -meson-meson mixing, is expected to be similar as it happens to be the case when two thresholds are incorporated into the lattice calculation [17,18]. Furthermore, the same form is expected for Q = b and Q = c since the underlying mixing mechanism (string breaking) is the same. Therefore, we shall proceed to a parametrization of |V mix (r)| according to this general form, and we shall rely on phenomenology to fix the values of the parameters.
Let us begin by observing that unquenched lattice QCD results show that for every value r, and that at the crossing radius r M1M 2 c , defined by |V mix (r)| gets approximately its maximum value with ∆ being the distance of the static energy levels at the crossing radius On the other hand we have far from the crossing radius r M1M2 c . Consequently, from (36) we obtain that V mix (r) vanishes in both asymptotic limits: To summarize, lattice QCD indicates that the mixing potential |V mix (r)| approaches a maximum value of ∆/2 at r ≈ r M1M2 c and vanishes asymptotically as the distance from the crossing radius increases. The simplest parametrization that takes into account these behaviours, thus providing a good fit to lattice QCD calculations of V ± (r), is a Gaussian shape: where Λ is a parameter with dimensions of energy. To better understand the physical meaning of Λ we write it in terms of the string tension σ as where ρ has now dimensions of length. Then at distances for which V C (r) ≈ σr+m Q +m Q −β the mixing potential can be also written as from which it is clear that ρ, the width of the Gaussian curve, fixes a radial scale for the mixing.

D. General case
The multichannel Schrödinger equation (27) defines the heavy quark meson system when only one threshold is considered, but in general it may be necessary to incorporate several meson-meson thresholds. In such a case one has to extend the formalism, what is is more easily done in the matrix notation (28).
The generalization of the kinetic energy matrix is straightforward: where µ (i) MM with i = 1, . . . , N is the reduced mass of the i-th meson-meson component, N is the number of meson-meson thresholds, and matrix elements equal to zero are not displayed.
As for the extension of the diabatic potential matrix (30), the presence of interaction terms between different meson-meson components would make not practicable our procedure to extract the mixing potentials. Following what it is usually done in molecular physics [19], we neglect some interactions between components. Namely, in line with lattice QCD studies of string breaking [18], we assume that different meson-meson components do not interact with each other. It seems reasonable to think that this naturally occurs when dealing with relatively narrow, well-separated thresholds. With this assumption the diabatic potential matrix with N thresholds reads where V C (r) stands for the Cornell potential, T MM for the mass of the i-th threshold and V (i) mix (r) for the mixing potential between the QQ and the i-th meson-meson components.
In Fig. 2 we draw the eigenvalues of this matrix for cc and the first three open flavour meson-meson thresholds.
The diabatic potential matrix (56) can be regarded as a generalization of the two threshold model of string breaking introduced in [18], the two main differences being that in our study each dynamical quark flavour can introduce more than one threshold and that we have parametrized the coupling between quark-antiquark and meson-meson components with a Gaussian instead of a constant.
Let us add that even tough there is presumably an infinite number of possible meson-meson components, in practice one needs to consider only a limited subset of them when searching for bound states. As a matter of fact, a meson-meson component hardly plays any role in the composition of a bound state whose mass lies far below the corresponding threshold.

E. Quantum numbers
Heavy-quark meson states are characterized by quantum numbers I G J P C where I, G, J, P , C stand for the isospin, G-parity, total angular momentum, parity, and charge conjugation quantum numbers respectively.
Let us focus on isoscalars I = 0 heavy-quark mesons, for which G = C. Since the diabatic potential matrix is spherically symmetric and spin-independent, the QQ component of the wave function can be characterized by the relative orbital angular momentum quantum number l QQ , the total spin s QQ , the total angular momentum J and its projection m J so that where Y m l l (r) is the spherical harmonic of degree l, ξ ms s is the eigenstate of the total QQ spin and Y l QQ (r)ξ s QQ mJ J is a shorthand notation for the sum where u (QQ) E,l QQ (r) is the QQ radial wave function. The same can be done for the meson-meson components of the wave function, considering the meson-meson relative orbital angular momentum l M1M 2 and the sum of their spins s M1M 2 . Therefore, with a straightforward extension of the above notation we write Note that for the spectroscopic state to have a definite value of J, the QQ and all the meson-meson components must have the same total angular momentum, hence the unified notation for J. A bound state made of QQ and meson-meson has definite parity and C-parity only if all the wave function components have the same parity under these transformations. This requirement translates into different conditions depending on whether the wave function component is associated with QQ or meson-meson. For the QQ component, P and C quantum numbers are given by P = (−1) l QQ +1 and C = (−1) l QQ +s QQ . (61) On the other hand, for each meson-meson component one has where P M is the parity of the meson. As for C-parity, one has to consider two distinct cases: if M 1 = M 2 the C-parity of the meson-meson component is given by if otherwise M 1 = M 2 one can build both positive and negative C-parity states taking the linear combinations with |M 1 M 2 0 being the isospin singlet state obtained from the combination of the M 1 and M 2 isomultiplets and

F. Bound state solutions
Given a spherically-symmetric and spin-independent diabatic potential matrix, each QQ configuration with a distinct value of (l QQ , s QQ ) can be treated as a channel per se, and the same can be said for each meson-meson configuration with a distinct value of (l M1M 2 , s M1M2 ). Then finding the spectrum of a given J P C family boils down to solving a multichannel, spherical Schrödinger equation involving only those channels with the corresponding J P C quantum numbers.
One should realize though that a complete numerical nonperturbative solution of the spectroscopic equations (28) is only possible for energies below the lowest J P C threshold. Above it the asymptotic behaviour of its meson-meson component as a free wave, against the confined QQ wave, prevents obtaining a physical solution. Nonetheless, an approximate physical solution for energies above threshold is still possible, under the assumption that the effect of an open threshold on the abovelying bound states can be treated perturbatively. More in detail, we proceed in the following way: i) We build the effective J P C diabatic potential matrix out of the Cornell QQ potential, the threshold masses, and the QQ -meson-meson mixing potentials.
ii) We solve the spectroscopic equations for energies up to the lowest J P C threshold mass, and we analyse the (n 2S+1 L J ) QQ and meson-meson content of the bound states.
iii) We build a new J P C diabatic potential matrix neglecting the QQ coupling to the lowest (first) threshold. We solve it for energies in between the lowest and the second thresholds and discard as spurious any solution containing a (n 2S+1 L J ) QQ state entering in the bound states calculated in ii). The rationale underlying this step is that a given spectral state in between the lowest and the second thresholds containing such a (n 2S+1 L J ) QQ component would become, when the lowest threshold were incorporated, the bound state below threshold containing it found in ii).
iv) We build a new J P C diabatic potential matrix by neglecting the coupling to the lowest threshold and to the second one. We solve it for energies in between the second and the third thresholds and discard as spurious any solution containing a (n 2S+1 L J ) QQ state entering in the bound states calculated in ii) and iii), and so on.
v) We assume that corrections to the physical states thus obtained due to the coupling with open thresholds can be implemented perturbatively.
The formulation of an appropriate perturbative scheme for the calculation of these corrections, giving rise to mass shifts as well as to decay widths to open flavour mesonmeson states, will be the subject of a forthcoming paper. Here we shall concentrate on the calculation of the "unperturbed" heavy-quark meson spectrum. The technical procedure followed to solve the spectroscopic equations is detailed in Appendices C and D.

V. CHARMONIUM-LIKE MESONS
The formalism we have developed in the previous sections can be tested in charmonium-like mesons (heavy mesons containing cc) where, unlike in the bottomoniumlike case, there are several well-established experimental candidates for unconventional isoscalar states, presumably containing significant meson-meson components. In particular, we centre on isoscalar states with masses up to about 4.1 GeV, for which the relevant thresholds have very small widths and do not overlap. A list of these thresholds is shown in Table I.
The possible values of the meson-meson relative orbital angular momentum contributing to any given set of quantum numbers J P C are shown in Table II. Note that   [2]. we use the common notation D (s) to refer to charmed as well to charmed strange mesons and the shorthand notation D (s) D * (s) for the meson-meson C-parity eigenstate defined by Eq. (65).
In order to calculate the heavy-quark meson bound states we have to fix the values of the parameters. For the Cornell potential (6)  in order to fit the 2s centre of gravity. Let us note that one could alternatively choose to fit the 1s or 1p centres of gravity, or to get a reasonable fit to the three of them. Our choice is based on the assumption that relativistic mass effects in the higher states, which are at least in part incorporated in β, are expected to deviate less from those in the 2s states. We should also mention that the value of the charm quark mass we use is completely consistent with the one needed to correctly describe cc electromagnetic decays within the Cornell potential model framework [22].
The low-lying spectrum from this Cornell potential for J ++ and 1 −− isoscalar states is shown in Table III. For the lowest J ++ states it is worth to remark, apart from the good average mass description, the excellent fit to the mass of the lowest 1 ++ state, χ c1 (1p) (3510.9 MeV versus the experimental mass 3510.7 MeV). However, an accurate fit of the lowest (0, 2) ++ masses, in particular  Table IV. Correlated values of the mixing potential parameters giving rise to a 0 + (1 ++ ) bound state with a mass close below the DD * threshold.
for χ c0 (1p), would require the incorporation of correction terms (e.g. spin-spin, spin-orbit, tensor) to the Cornell radial potential. As for the first excited J ++ states one could expect a similar situation (the 2s states lie in between the 1p and 2p ones) in the absence of threshold effects that we analyse in what follows.
As for the parameters of the mixing potential (45), we have to rely on phenomenology since the only lattice information available is for bb. We fix them by requiring that our diabatic treatment fits the mass of some unconventional experimental state lying close below threshold. In particular, we can use the mass of χ c1 (3872), a wellestablished experimental resonance lying just below the DD * threshold, to infer the possible of values for ∆ cc and ρ cc . As the crossing of the Cornell potential with the DD * threshold takes place around r DD * c = 1.76 fm, we conservatively vary ρ cc from 0.1 fm to 0.8 fm, this last value corresponding to almost half of r DD * c . Then, for every value of ρ cc we get the minimal value of ∆ cc to accurately fit the mass of χ c1 (3872). The calculated values are listed in Table IV. It should be pointed out that large values of ∆ cc would deform the shape of the avoided energy crossings as compared to the one calculated in lattice for bb, against our bbcc universality arguments for the shape of the mixing potential. On the other hand, large values of ρ cc would make the mixing angle between the cc and a single M 1 M 2 threshold, calculated from Eq. (52), to have an asymp- totic behaviour in conflict with the one observed in the lattice under the natural assumption that this behaviour is similar for bb and cc. More precisely, unquenched lattice QCD calculations of the mixing angle [17] show that θ approaches π/2 quite rapidly for r > r M1M2 c , thus ruling out a large radial scale for the mixing. Henceforth we use for this value gives the most accurate asymptotic behaviour of the mixing angle, see Fig. 3, and consequently The resulting mixing potential is drawn in Fig. 4 for For any other threshold the only difference comes from the substitution of the threshold mass. Notice that we have drawn |V mix (r)| with no sign prescription for V mix (r). This sign can be reabsorbed as a relative phase between the charmonium and mesonmeson components. For the calculations in this paper a positive sign has been taken. We have checked that for the observables considered in this article the same results are obtained with a negative sign. It should be realized though that this could not be the case for other observables.
The calculated spectrum of J ++ states, containing one cc state with l cc = 1 (1p or 2p), is shown in Table V.
It is illustrative to compare these results with the cc masses in Table III obtained with the Cornell potential. A glance at these tables makes clear that the presence of the thresholds gives rise to attraction in the sense that the resulting masses are reduced with respect the corresponding Cornell cc masses. For the lowest-lying 0 + (J ++ ) states (J = 0, 1, 2) there is a very small mass difference indicating an almost negligible attraction for these states.  This is understood for the thresholds are far above in energy (≥ 200 MeV) so that no significant mixing occurs (less than 1% meson-meson probability). The situation is completely altered for the first excited 0 + (J ++ ) states. Thus, the fitting of the first excited 1 ++ resonance, χ c1 (3872) with a measured mass of 3871.69 ± 0.17 MeV, requiring a mass reduction of 81 MeV with respect to the Cornell cc mass, implies a very strong mixing, 99% of DD * component, whereas for the 0 ++ and 2 ++ states the predicted mixing is about 40% (mainly from D s D s ) and 15% (shared by D s D s and D * D * ) respectively, with corresponding mass reductions of 33 MeV and 20 MeV. It is amazing that these (0, 2) ++ mass predictions are in complete agreement with data regarding their positions with respect to the D s D s threshold, both below it. Moreover, their calculated numerical values are pretty close to the measured ones. So, the calculated 2 ++ mass, 3933.5 MeV, is very close to that of the experimental resonance χ c2 (3930): 3927.2 ± 2.6 MeV. And the 0 ++ calculated mass, 3920.4 MeV, is consistent with the ones of the experimental candidates: χ c0 (3860), with a measured mass of 3862 +26+40 −32−13 MeV, and X(3915) with a measured mass of 3918.4 ± 1.9 MeV, although in this last case the assignment to a 2 ++ state cannot be completely ruled out, see [2] and references therein. This suggests that further mass corrections for these sates as the ones due to spin-dependent terms in the cc potential, or those taking into account the effect of the lower threshold DD, or the deviations from the assumption of the same values of the mixing potential parameters for all the thresholds, are small and might be implemented perturbatively.
It should also be emphasized that our nonperturbative formalism provides us with the meson wave functions in terms of their cc and meson-meson components.
For χ c1 (3872) the radial cc and DD * (l DD * = 0, 2) wave function components are plotted in Fig. 5.  For the calculated 0 ++ state the radial wave function is drawn in Fig. 6.
As can be checked, the wave function with a rms radius of 1.26 fm is made mainly of cc and D s D s with a 59% and 37% probability respectively. This indicates a dominant DD strong decay mode from cc as it is experimentally the case for χ c0 (3860). On the other hand, a J/ψ ω decay mode may get a significant contribution from D s D s since it is OZI allowed through the small ss content of ω. This could cause this mode to be also a dominant one as it is experimentally the case for X(3915). Hence, it could be that χ c0 (3860) and X(3915) are just the same resonance observed through two different decay modes. As for the calculated 2 ++ state, the wave function, with a rms radius of 1.06 fm, is plotted in Fig. 7. It is mostly that of the cc component. This is in accord with a very dominant DD strong decay mode as it is experimentally the case for χ c2 (3930).
Certainly these qualitative arguments on the dominant strong decay modes should be supported by trustable and predictive quantitative calculations. As mentioned before, the development of a consistent formalism for the calculation of the decay widths to open flavour mesonmeson states, which is out of the scope of this article, is in progress. One should keep in mind though that the dearth of current detailed quantitative decay data for comparison will be a serious drawback to test it. We strongly encourage experimental efforts along this line.
Regarding electromagnetic radiative transitions, although important progress for the accurate calculation of decays from the cc component has been reported [22], a reliable and consistent calculation incorporating the meson-meson contribution as well is lacking. We encourage a theoretical effort along this line.
One can do better, as we show next, for leptonic decays from the low-lying 1 −− states since the decay widths depend on the wave function at the origin and the contribution from meson-meson components is suppressed as they are not in s-wave, see Table II.
The calculated 1 −− spectrum of states is listed in Table VI.
Again, a comparison with the cc masses in Table III makes clear that the presence of the thresholds gives rise to attraction. As it was the case for J ++ , the lowest state, lying far below the lowest threshold, has no mixing at all being the 1s cc state. A pretty small mixing is present for the next two higher states that can be mostly assigned to the 2s (95%) and 1d (97%) cc states respectively. It is worth to mention that for the 1d state with a Cornell cc mass of 3795.8 MeV, the DD threshold lying 66 MeV below does not produce enough attraction to bring the state below threshold.
The first state with a significant mixing, 36% of D s D * s , is predicted at 4071 MeV and contains a 60% of 3s cc and a 4% of 2d cc as well. Its wave function is drawn in Fig. 8.
In this case the vicinity of the D s D * s threshold at 4080 MeV to the 3s cc Cornell mass at 4097 MeV produces sufficient attraction to bring the state below threshold, in agreement with data under its assignment to the ψ(4040) resonance with a measured mass of Hence, our results agree with data within the experimental intervals. The reason for this agreement has to do with the reduced probability of the 3s cc component, 60%, induced by the mixing with the D s D * s threshold. This mixing is also responsible for the 4% of 2d cc component. This small (big) 2d (3s) probability could be increased (decreased) if a tensor interaction were incorporated as a correction term to the Cornell potential. Maybe the bias we observe in our results, both agreeing with the maximum allowed experimental values, is an indication in this sense. In any case a modest additional probability reduction of the 3s cc component should be expected.
It is worth to mention that the explanation of the leptonic width for ψ(4040) has been linked in the literature to that of ψ(4160) through a very significant s-d mixing [23]. Our results do not support this idea. Instead the D s D * scc(3s) mixing appears to be the main physical mechanism underlying the ψ(4040) decay to e + e − .
Unfortunately, at the current stage of our diabatic development we cannot properly evaluate ψ(4160), the main reason being that the dominant Cornell 2d cc state lies only 100 MeV below the first s-wave 1 −− threshold, DD 1 , which is composed of two overlapping thresholds, DD 1 (2420) and DD 1 (2430), the last one with a large width. Quite presumably this double threshold gives a significant contribution by itself to the leptonic width of ψ(4160).
This current limitation applies as well to the description of unconventional states with masses above 4.1 GeV such as ψ(4260) lying close below the DD 1 double threshold, or ψ(4360) and ψ(4415) lying close below a multiple threshold at 4429 MeV. The same limitation applies for J ++ states. Work along this line is in progress.

VI. SUMMARY AND CONCLUSIONS
A general formalism for a unified description of conventional and unconventional heavy-quark meson states has been developed and successfully applied to isoscalar J ++ and 1 −− charmonium-like states with masses below 4.1 GeV.
The formalism adapts the diabatic approach, widely used in molecular physics to tackle the configuration mixing problem, to the study of heavy-quark meson states involving quark-antiquark as well as meson-meson components. A great advantage of using this approach, against the Born-Oppenheimer (B-O) approximation commonly used for heavy-quark mesons, is that the bound states are expanded in terms of QQ and meson-meson configurations instead of the mixed configurations that correspond to the ground and excited states of the light fields. Then instead of being forced to use a single channel approximation to solve the bound state problem as in B-O, what in practice is equivalent to neglect the configuration mixing, one can write a treatable multichannel Schrödinger equation where the interaction between configurations is incorporated through a diabatic potential matrix. Moreover, the diagonal and offdiagonal elements of this potential matrix can be directly related to the static energies obtained from ab initio quenched (only QQ or meson-meson configuration) and unquenched (QQ and meson-meson configurations) lattice calculations. This connection defines the diabatic approach in QCD.
It is worth to emphasize that this approach goes also beyond the incorporation of hadron loop corrections to the B-O scheme that have been used sometimes in the literature to deal with unconventional charmonium-like mesons. Indeed, the diabatic bound state wave functions, given in terms of quark-antiquark and meson-meson components, allow for a complete nonperturbative evaluation of observable properties.
This theoretical framework has been tested in the charmonium-like meson sector where there is compelling evidence of the existence of mixed-configuration states, in particular the very well-established 0 + (1 ++ ) resonance χ c1 (3872) that we use to fix our parametrization of the mixing potential.
Although a complete (at all energies) spectral description would require additional theoretical refinements, as for example the incorporation of threshold widths, the results obtained for states with mass below 4.1 GeV, for which the significant thresholds are very narrow, are encouraging. All the mass values are well reproduced and their locations with respect to the thresholds correctly predicted making clear the cc -threshold attraction. This points out to the diabatic approach as an appropriate framework for a unified and complete nonperturbative description of heavy-quark meson states. stationary point is precisely the corresponding energy eigenvalue.
Technically speaking, the results presented here are analytically valid only when using a complete, i.e. infinite, orthonormal set. Since in realistic applications one deploys a limited set, the correspondence drawn here is only approximate and so are the energies and eigenstates obtained with the variational method.
A shortcoming of the variational method is that the degree of approximation is not known a priori. To assure this not to be any problem we choose an appropriate orthonormal set of states reflecting some of the properties of the physical states and employ a very high number of states in the set.