The Dynamical Diquark Model: First Numerical Results

We produce the first numerical predictions of the dynamical diquark model of multiquark exotic hadrons. Using Born-Oppenheimer potentials calculated numerically on the lattice, we solve coupled and uncoupled systems of Schroedinger equations to obtain mass eigenvalues for multiplets of states that are, at this stage, degenerate in spin and isospin. Assuming reasonable values for these fine-structure splittings, we obtain a series of bands of exotic states with a common parity eigenvalue that agree well with the experimentally observed charmoniumlike states, and we predict a number of other unobserved states. In particular, the most suitable fit to known pentaquark states predicts states below the charmonium-plus-nucleon threshold. Finally, we examine the strictest form of Born-Oppenheimer decay selection rules for exotics and, finding them to fail badly, we propose a resolution by relaxing the constraint that exotics must occur as heavy-quark spin-symmetry eigenstates.


I. INTRODUCTION
The existence of multiquark exotic hadrons is now a well-established feature of strong-interaction physics. At least 35 such states in the heavy quark-antiquark (Q = c or b) sector have been experimentally established at various levels of statistical significance. Many have been observed beyond the 5σ level at either multiple facilities, or through multiple production or decay channels, or both. A number of reviews in the recent literature [1][2][3][4][5][6][7][8][9] summarize both experimental and theoretical developments.
Yet, even after many years of intensive study, no single theoretical model has emerged to provide a successful, unified picture for understanding the spectroscopy, decay patterns, and structure of these novel states. The most heavily studied alternatives in this regard, including hadronic molecules, diquarks, hadroquarkonium, hybrids, and kinematical threshold effects, both benefits and drawbacks, are amply discussed in the aforementioned reviews. The complete spectrum might turn out to rely upon a delicate interplay of several of these physical effects, meaning that each one must be fully understood before a global model of the exotics can be confidently constructed.
In this work we employ a model in which the exotics are constructed from quasi-bound heavy-light diquarks δ,δ, which are formed via the attractive channels 3 ⊗ 3 →3 [δ ≡ (Qq)3] and3⊗3 → 3 [δ ≡ (Qq )3] between the colortriplet quarks. The most influential early application of diquarks to the problem of heavy exotics [10] treats tetraquarks as bound (δδ) molecules in a Hamiltonian formalism, using as interaction operators the spin-spin couplings between the various component quarks. A later variant of this approach [11] restricted the interactions to spin-spin couplings between quarks within either the δ * jfgiron@asu.edu † Richard.Lebed@asu.edu ‡ curtistaylor.peterson@asu.edu or theδ. Reference [12] provides a detailed review of diquark phenomenology prior to the discovery of the heavy exotics.
Such an approach inspired the development of the dynamical diquark picture [13], in which some of the light quarks created in the production process of the QQ pair coalesce with these heavy quarks to form a δ-δ pair. Due to the large energies available in either b → c or collider processes in which exotics are produced, the δ-δ pair can achieve through recoil a large spatial separation (> 1 fm) before being forced by confinement either to form a single tetraquark state, or if the energy is sufficiently high, for the color flux tube between the δ-δ pair to fragment to create a baryon-antibaryon pair. The key feature of this picture is a mechanism for producing multiquark states that are spatially large yet strongly bound. Indeed, the successive accretion of additional quarks through the color-triplet binding mechanism [14] can be used to interpret pentaquark states as triquarkdiquark states,θδ ≡ [Q(q 1 q 2 )3] 3 (Qq 3 )3 [15]. The effectiveness of this mechanism is clearly limited by competition from the attractive QQ andQq channels, and a full theory of multiquark hadrons would allow for including all possible configurations simultaneously. One important step in this much larger project is to uncover the predictions of the δ-δ mechanism and see if the current data provides support for the existence of such states.
The means by which the dynamical diquark picture may be realized as a fully predictive model, including spectroscopy and decay selection rules, is the subject of Ref. [16]. In the original proposal of the picture [13], the estimated size of the Z − c (4430) resonance appearing in B 0 → (ψ(2S)π − )K + follows from taking a δ-δ pair (of known masses) produced in the decay to recoil against the K + . Since the diquarks, like quarks, are color triplets, a Cornell Coulomb-plus-linear potential [17,18] was assumed. The final separation of the δ-δ pair upon coming relatively to rest was calculated to be 1.16 fm.
The significance of this classical turning point in forming the exotic state from the δ-δ pair was tied in Refs. [13,16] to the Wentzel-Kramers-Brillouin (WKB) approximation. The WKB transition wave function scales as 1/ √ p, where p is the classical relative momentum of the constituents. One therefore expects the color flux tube between the δ-δ pair to stretch nearly to its classical limit. Such a state is spatially large but still exhibits unscreened strong interactions between all of its components. It also possesses two heavy, slow-moving sources (δ andδ) connected by a lighter (mostly gluonic) field susceptible to more rapid changes. Analogous comments apply toθ-δ states. These properties indicate that the system can be characterized well by use of the Born-Oppenheimer (BO) approximation [19]. Although more familiar from its applications to atomic and molecular systems, the BO approximation has also been implemented in particle physics. In fact, its use as a fundamental tool in lattice-QCD calculations was initiated decades ago [20]. The relevant physical observables are the energies of the light degrees of freedom (d.o.f.) that connect a heavy, and hence static, QQ pair; such energies as a function of QQ separation and orientation are called BO potentials. While multiple aspects of this problem in strong-interaction physics have been studied in the intervening years, the ones most relevant to the present work involve the calculation of the BO potentials and its eigenvalues, which in turn give the masses of heavy-quark hybrid mesons; short overviews of the key lattice papers in this regard appear in Refs. [21,22]. For many years, the most accurate lattice results of hybrid static potentials for substantial QQ separation have been those of Juge, Kuti, and Morningstar [23][24][25][26]. Very recently, however, a new collaboration [27] has begun to improve upon these results. In addition, high-quality simulations focusing upon small QQ separations have been performed [28]. Also of note are lattice simulations of two heavy quark-antiquark pairs, which study the crossover between hadron molecule and δ-δ configurations [29].
A prototype of the approach underlying the present work is provided by Ref. [21], which supposes the known (neutral) exotics to be hybrid mesons, and then develops an effective field-theory formalism for computing their spectrum, using the BO potentials as input to the Schrödinger equations. The first treatment of multiquark exotic hadrons using the BO formalism appeared in Refs. [30][31][32]. In these works, the valence light qq pair (for tetraquarks) is treated as belonging to the light d.o.f. In contrast, the light quarks in the dynamical diquark model belong to the diquarks, which in turn are treated as the heavy, pointlike sources, while the light d.o.f. are purely gluonic (or include sea quarks).
This paper carries out one of the central proposals of Ref. [16]: a numerical calculation of the spectrum of δδ andθ-δ hidden-charm states, under the assumption that their basic structure-at least in the last moments of their evolution prior to decay-consists of heavy, slowmoving compact diquarks interacting through the same hybrid BO potentials as those appearing in the QQ sector. We solve the resulting Schrödinger equations numer-ically and identify known exotic states with the eigenstates of the lowest-lying BO potentials. One can already identify a number of assumptions implicit in this strategy; these and several others of equal significance are discussed below. Nevertheless, the initial results are quite encouraging: Choosing to fix to either X(3872) or Z c (4430) as a reference δ-δ state, one obtains a spectrum broadly consistent with the pattern of the known tetraquark states, and for which the excited states above the ground-state band [Σ + g (1S)] naturally have substantial spatial extent. In the pentaquark case, fixing to P c (4380) and P c (4450) predicts the masses of numerous unseen states.
In carrying out calculations in this scheme, one must keep in mind the central difficulty with diquark models (as discussed in any of the reviews [1][2][3][4][5][6][7][8][9]): the proliferation of many more potential multiquark states than have been observed to date. For example, the J P C = 1 ++ X(3872) appears to lack an isospin partner [33] that would arise from replacing its u → d quarks. On the other hand, the 1 +− Z c (3900) lies so close in mass to X(3872) as to suggest some sort of multiplet structure. A truly complete diquark model must explain both the absence of the former and the presence of the latter. While we indicate how the details of this fine structure might be resolved in Sec. V, for the present work we seek only to establish the broader pattern of multiplets, in particular, as collections of states sorted in mass by parity eigenvalue.
The organization of this paper is as follows: In Sec. II we reprise the notation for dynamical diquark-model states developed in Ref. [16]; the reader unfamiliar with BO notation appropriate to the "diatomic" system is referred to Appendix A. Section III discusses the relevant Schrödinger equations, both uncoupled and coupled versions, as introduced in Ref. [21]. Details of our numerical approach to solving these equations appear in Appendix B. In Sec. IV we present our results, and outline the approximations used in obtaining them in Sec. V, describing how these simplifications can be lifted one by one. Finally, Sec. VI presents our conclusions and directions for subsequent development of the model.

A. Tetraquarks
The notation adopted in Ref. [16] for states in the dynamical diquark model begins with the notation introduced in Ref. [11] for diquark-antidiquark (δ-δ) states of good total J P C with zero orbital angular momentum: The number before each δ(δ) subscript is the diquark (antidiquark) spin, and the outer subscript on each ket is the total quark spin J. By straightforward use of 9j symbols, these states can also be expressed in terms of states of good heavy-quark (QQ) and light-quark (qq) spins [from which eigenvalues of the charge-conjugation parity C given in Eq. (1) are immediately determined]. The states in Eq. (1) represent the tetraquark analogues to heavy S-wave quark-model states such as η Q and ψ(Υ). One then allows for nonzero relative orbital angular momentum L between the δ-δ pair. Using the generic symbol Y for X 0 , X 1 , Z, Z , X 2 in Eq. (1), and J = L+S, where now S is the total quark spin, one obtains the states Y (J) L . In terms of C Y , the C-parity of the underlying S-wave state Y , these states have P = (−1) L and C = (−1) L C Y . Such states include the analogues of the P -wave quark-model states χ Q , h Q , as well as the Dwave, F -wave, etc. states. All of these states also possess radial excitations, labeled by the quantum number n.
Lastly, one appends the quantum numbers from the Born-Oppenheimer (BO) excitation of the gluon field with respect to the quark state. Given the BO potential with quantum numbers Γ ≡ Λ η as defined in Appendix A, states receive multiplicative factors to their P and C quantum numbers of ρ ≡ (−1) Λ and κ ≡ η (−1) Λ = ηρ, respectively. In total, and suppressing the radial quantum number n, the physical tetraquark eigenstates may be labeled Y The resulting states associated with the lowest BO potentials (as calculated on the lattice) are listed in Table I. If the light d.o.f. carry nonzero isospin I [e.g., (cu)(cd)], then C-parity eigenvalue of the state is replaced by the G-parity eigenvalue, G ≡ C(−1) I , where the C eigenvalue is that of the neutral member of the isospin multiplet.

B. Pentaquarks
Much of the same construction holds for theθ-δ pentaquarks. In that case, the states analogous to those in Eq. (1) are denoted by [16]: The number before eachθ(δ) subscript is the triquark (diquark) spin, and the outer subscript on each ket is the total quark spin J. In this list, the light diquark internal toθ is restricted to carry spin 0 (as well ud flavor content with isospin 0), since the known heavy pentaquark candidates P c (4380) and P c (4450) both appear in the decay of Λ b , whose light ud valence quarks carry these attributes. In general, 6 additional states, for which the light diquark inθ carries spin 1, can be defined. The relative orbital angular momentum L and the BO potential quantum numbers Γ ≡ Λ are then incorporated. Since theθ and δ components cannot form a charge-conjugate pair, the core states of Eq. (3) are not C eigenstates, and the BO potentials Λ are not η eigenstates (see Appendix A). Defining ρ ≡ (−1) Λ as before, one obtains Suppressing the radial quantum number n, the physical pentaquark eigenstates may be labeled P (J)ρ SL , where S now denotes the total quark spin. The resulting states associated with the lowest BO potentials (as calculated on the lattice) are listed in Table II.

III. SCHRÖDINGER EQUATIONS FOR THE BORN-OPPENHEIMER POTENTIALS
At its core, the calculation of the spectrum of hybrids in Ref. [21] amounts to the use of the BO potentials calculated on the lattice in Refs. [25,28] to find the energy eigenvalues of Schrödinger equations between a static QQ pair (cc, cb, and bb are all considered). The relevant Schrödinger equations actually arise directly from QCD through a systematic 1/m Q expansion by the application of effective field theories: first NRQCD [34,35] (in which the hard scale m Q is integrated out), and then pN-RQCD [36,37] (in which the softer scale of momentum transfer between the QQ pair is integrated out). The gluonic pNRQCD static energies between the QQ pair are then none other than the BO potentials, which are obtained numerically on the lattice.
Since the fundamental quark mass m Q appears directly in the analysis of Ref. [21], the authors take care to identify the details of their renormalization scheme for both m Q and the perturbative short-distance behavior of the potential between the fundamental QQ pair. In our case, the corresponding mass m δ is that of the diquark (or mθ for the triquark), which is of course not a fundamental Lagrangian parameter, and therefore in this analysis m δ , TABLE I. Quantum numbers of the lowest tetraquark states expected in the dynamical diquark picture. For each of the expected lowest Born-Oppenheimer potentials, the full multiplet for given nL eigenvalues is presented, using both the state notation developed in Ref. [16] and the corresponding J P C eigenvalues. States with J P C not allowed for conventional qq mesons are indicated in boldface.
More central to the current calculation is that the full Schrödinger equations for the "diatomic" system contain not one, but two special points, corresponding to the two heavy-constituent positions, separated by a distance r. One expects additional symmetries between the static energies to arise in the limit r → 0, where the cylindrical D ∞h symmetry for the "homonuclear" δ-δ case or the conical C ∞v symmetry for the "heteronuclear"θ-δ case (see Appendix A and Fig. 1) is supplanted by the higher spherical O(3) symmetry. These r → 0 static energy configurations, transforming as color adjoints in the light d.o.f. and corresponding to degenerate BO potentials in the limit r → 0, are called gluelumps. One finds, for instance, that the δ-δ Σ − u and Π + u BO potentials approach a single gluelump with J P C = 1 +− . Additionally, the loss of the single (spherical) symmetry center for r > 0 means that the angular part of the Laplacian in the Schrödinger equation is no longer solved by familiar spherical harmonics, but by slightly more complicated forms [38].
All δ-δ BO potentials arising from a particular gluelump appear together in a coupled system of Schrödinger equations. The ground-state BO potential Σ + g does not mix with others, and therefore appears in an uncoupled equation. While the degenerate BO potentials Π + u and Π − u carry opposite parities, only Π + u (nL) produces states with the same J P C quantum numbers as Σ − u (nL) for L > 0 [Π + u (nS) is forbidden by Eq. (A2)]. Indeed, Π + u is found in lattice simulations to approach Σ − u as r → 0, to form the 1 +− gluelump [37]. The Schrödinger equations for Π + u (nL) and Σ − u (nL) with L > 0 therefore must be solved as a coupled system, while those for Σ + g (nL), Π − u (nL) (with L > 0), or Σ − u (nS) remain uncoupled. One finds different mass eigenvalues emerging from Π ± u (1P ), a lifting of the parity symmetry called Λdoubling [38]. Higher-mass gluelumps have been found to split into even more BO potentials [37]. For example, the 2 −− gluelump supports the BO potentials Σ − g , Π + g , and ∆ − g ; its Schrödinger equations for D-wave solutions would couple all three of them. Analogous comments apply to theθ-δ BO potentials, once the η = g, u subscript is removed.
The radial Schrödinger equations for the uncoupled BO potentials V Γ assume the conventional form (with = 1): and for the coupled potentials Π + u , Σ − u , they read: TABLE II. Quantum numbers of the lowest pentaquark states expected in the dynamical triquark-diquark picture. For each of the expected lowest Born-Oppenheimer potentials, the full multiplet for given nL eigenvalues is presented, using both the state notation developed in Ref. [16] and the corresponding J P eigenvalues.
BO states State notation The details of the numerical methods for solving these coupled Schrödinger equations, using the state-of-the-art techniques of Refs. [39] and [40], are discussed in Appendix B.

IV. RESULTS
Using the techniques of Sec. III, one predicts the full spectrum of states in the dynamical diquark model with only two further inputs: the diquark masses m δ and mδ (or triquark mass mθ) and specific functional forms V (r) for the BO potentials Γ. The potentials are intrinsically nonperturbative in nature and can only be computed from first principles by using lattice QCD simulations. In this regard we apply the results of Refs. [23][24][25][26], especially the online summary of results in Ref. [26], to which we refer as JKM; and separately, we apply the results of the very recent calculations of Ref. [27], to which we refer as CPRRW. Furthermore, since the corresponding calculation for the ground-state (Σ + g ) BO potential for the cc system with both component masses given by m c would generate the conventional charmonium spectrum, we also include for Σ + g cases the phenomenological Cornell potential used in the fit of Ref. [41] (but suppressing spin-dependent couplings), to which we refer as BGS.
The numerical results are summarized in Tables III,  IV, and V for hidden-charm tetraquarks and Table VII for hidden-charm pentaquarks. We take mδ = m δ in each case. Fine-structure spin-dependent mass splittings among the states of each level Γ(nL) (as enumerated in Tables I and II) are neglected in the present calculation, and the approach to include them in future calculations is discussed in detail in Sec. V.
The tables are organized by identifying particular exotic states of known mass and J P C eigenvalues [42] as reference states for particular levels Γ(nL) that contain a state of the given J P C . For charged states, the C value used is that of its neutral isospin partner. Our reference states are: In Tables III, IV, and V, the fits in the left-hand columns correspond to choosing X(3872) to be the unique Σ + g (1S) 1 ++ state, one of the two 1 ++ states in Σ + g (1D), and one of the two 1 ++ states in Π + u (1P )-Σ − u (1P ), respectively (these being the only n = 1, J P C = 1 ++ states in Table I).
The fits in the right-hand columns correspond to choosing Z c (4430) to be one of the two 1 +− states in Σ + g (2S), one of the two 1 +− states in Σ + g (2D), and one of the , respectively (these being the only n = 2, J P C = 1 +− states in Table I). The value of m δ is obtained from these fits, and is used for all other states in the spectrum of Tables III, IV, and V. Also calculated in the tables are values of the typical length scales 1/r −1 and r for the states; as described in Ref. [40] and in Appendix B, expectation values can be calculated using the same procedure as one uses to compute eigenvalues, without the need to generate explicit eigenfunctions.
The choice of X(3872) as an n = 1 state and Z c (4430) as an n = 2 state is not logically necessary. However, if X(3872) is chosen as the lowest n = 2 state [Σ + g (2S)], then the mass of the Σ + g (1S) state is about 3270 MeV, which lies squarely in the region of conventional charmonium between J/ψ and χ c0 (1P ). Such a state with any of the Σ + g (1S) quantum numbers given in Table I would certainly have been discovered decades ago, thus rendering the n = 2 assignment untenable. Similar comments apply to the Σ + g (1S) states predicted in Tables IV and V; indeed, the fit of Table IV fixing X(3872) to Σ + g (1D) and the fit of Table V GeV m c , which means that these portions of the tables effectively reproduce the conventional charmonium spectrum plus its lowest hybrids, once δδ are replaced withcc. The other fits in Tables IV and V predict exotic Σ + g (1S) states that lie below the opencharm threshold at some distance from the conventional charmonium states that are known to be the only ones populating this region, and hence again produce conflicts with observation.  ) for hidden-charm dynamical diquark states that are eigenstates of the indicated BO potentials corresponding to quantum numbers nL, for given diquark masses m δ (in GeV). The particular form of the BO potential used is that given by lattice simulations JKM [26] or CPRRW [27], or (for the Σ + g potential) by the Cornell potential obtained by a fit BGS [41] to conventional charmonium. Also predicted are the corresponding expectation values for the length scales 1/r −1 and r (in fm). Fixing to the experimental mass of X(3872) or Z − c (4430) predicts m δ and the whole spectrum, either under the assumption that X(3872) is a Σ + g (1S) state or that Zc(4430) is a Σ + g (2S) state (as indicated by boldface).   As for the Z c (4430) as an n = 1 state, the assignments Σ + g (1D) or Π + u (1P )-Σ − u (1P ) (not tabulated) are logically possible, but they again lead to a mass for Σ + g (1S) that is too small (3815 and 3450 MeV, respectively) compared to data, and also a Σ + g (2S) mass (4390 and 4045 MeV, respectively) that, combined with the Σ + g (1S) states, would generate a total of at least four 1 +− states lying below the Z c (4430). Experimentally, only two candidates [Z c (3900) and Z c (4200)] have confirmed J P C = 1 +− , but even the existence of the Z c (4200) remains unconfirmed. In contrast, the existence of the Z c (4020) is confirmed, and it is widely expected to be 1 +− , but only its C = − quantum number has been confirmed. As discussed previously, the assignment of Z c (4430) to an n = 1 state is problematic due to the prediction of numerous unseen light exotic states, although the choice of Z c (4430) as Σ + g (1D) is not yet definitively excluded, particularly if a pair of 1 +− exotics near 4390 MeV [from Σ + g (2S)] is observed.
One unique level assignment appears to work particularly well with observation: In Table III, X(3872) and Z c (4430) are identified as states in the multiplets Σ + g (1S) and Σ + g (2S), respectively.  [13]. The lowest levels in order of increasing mass are Σ + g (1S), Σ + g (1P ), Σ + g (2S), Σ + g (1D), Σ + g (2P ), Σ + g (3S) (not tabulated, ≈ 4890 MeV), Π + u (1P )-Σ − u (1P ), and Σ + g (2D). In particular, only one of the BO potentials beyond Σ + g is represented among the lowest states, and even then, it is in the 6 th excited level. What one might call the "hybrid" exotic levels begin ∼ 1 GeV above the Σ + g (1S) states, just as for hybrid charmonium [21]. The lowest observed states in this assignment are therefore "quark-model" δ-δ states, in that the gluonic field does not contribute to the valence spin-parity quantum number. The states enumerated in Ref. [11] are included in this list, although their order and spacing as determined here depends intrinsically upon the calculated BO potential Σ + g (r). The particular mass eigenvalues calculated here apply to all states in each multiplet Γ(nL) listed in Table I, which are degenerate at this level of the calculation: Fine-structure spin-(and isospin-) dependent corrections are therefore neglected here. However, one may estimate the magnitude of these mass splittings by examining those for conventional charmonium. Note first that m J/ψ −m ηc = 113 MeV and m χ c2 (1P ) −m χ c0 (1P ) = 141 MeV, and the charmonium fine-structure splittings tend to decrease somewhat with higher excitation number. These splittings arise from spin-spin, spin-orbit, and tensor cc operators all proportional to 1/m 2 c . The corresponding coefficient in the dynamical diquark model is slightly smaller (1/m 2 δ ), but the diquark spins can be as large as 1 (compared with 1 2 for c orc), so that the spinoperator expectation values in the numerator can be substantially larger. Based upon this reasoning, we crudely (but conservatively) estimate the largest fine-structure splittings in each exotic multiplets to be ∼ 150 MeV; the true value might turn out to be even larger, but then one runs the risk of producing multiple overlapping bands of exotic states, resulting in diminished model predictivity. Since the X(3872) appears to be the lowest exotic candidate, Table I (2P )]. Since the heaviest charmoniumlike state currently observed is X(4700), we do not analyze the higher levels in any further detail.
In fact, X(4700) and several of the other exotic candidates [Y (4140), Y (4274), X(4350), X(4500)] have only been observed as resonances decaying to J/ψ φ, which makes them good candidates for ccss states [43], and if so, then they do not belong to the current analysis; instead, they would appear as part of an identical analysis using heavier (cs), rather than (cu) or (cd), diquarks.
Perhaps the most evident pattern among the spectrum of charmoniumlike exotic bosons [1], once the five states listed above are removed, is the appearance of fairly well-separated clusters: one between the X(3872) and at least as high as Z c (4020) [and possibly as high as the less well-characterized states Z c (4050) and Z c (4055)]; another from Z c (4200) to Y (4390) [and possibly as low as the less well-characterized X(4160)]; the Z c (4430) by itself; and the 1 −− states X(4630) and Y (4660). Even among the five ccss candidates, only Y (4140) appears to lie starkly outside this band structure. The recently observed [44] (at > 3σ) η c π − resonance Z c (4100) also appears to fall into this gap.
We now confront this hypothesis with the full data set. The charmoniumlike states for which the J P C quantum numbers are either unambiguously determined or "favored" experimentally [42] are listed in Table VI. In addition, some dispute remains that the X(3915) could be a 2 ++ state [45] [possibly the same as the conventional charmonium χ c2 (2P )], and the recently observed χ c0 (3860) [46] has properties consistent with being the conventional 0 ++ χ c0 (2P ), but the existence of this state has not yet been confirmed. Charged [47] and neutral [48] structures from 4032-4038 MeV are not included because of their uncertain nature. Lastly, several states (confirmed and unconfirmed) have unknown J P but known C eigenvalues [42]: Z 0 c (4020) (noted above) and Z 0 c (4055) are C = −, while Z 0 c (4050), Z 0 c (4250), X(4350) are C = +. Under our hypothesis, in Σ + g (1S) the X(3872) is the sole 1 ++ state, and X(3915) is one of two 0 ++ states (although it could instead be the ccss ground state [43]), with the other 0 ++ state possibly being χ c0 (3860) [instead of the cc χ c0 (2P ) assignment], or even possibly (if fine-structure splittings turn out to be large in this case and J P = 0 + is confirmed) Z 0 c (4100). The two 1 +− states are Z 0 c (3900) and Z 0 c (4020). As for 2 ++ , the state χ c2 (3930) has still not been confirmed as the cc χ c2 (2P ), making it a potential candidate to complete the Σ + g (1S) multiplet. The states X(3940), Z 0 c (4050), and Z 0 c (4055) are also potential members (noting the known C eigenvalues of the latter two). The existence of the only other state claimed in this range, the 1 −− Y (4008) [49], is being challenged by increasingly adverse evidence, and the state may disappear completely with newer data and analysis. In fact, Y (4008) has P = − and thus would not fit into Σ + g (1S), which represents a success of the model: All hidden-charm exotics below 4100 MeV are predicted to have positive parity.
The Σ + g (1P ) states all have P = −, and indeed, Y (4220), Y (4360), and (with a small stretch of the band mass range) Y (4390) fit into this multiplet. The multiplet actually contains a fourth 1 −− state, but note that the famous 1 −− Y (4260) may actually be a composite of the other 1 −− states [50]. The 0 −− Z 0 c (4240) is especially notable because, if confirmed, it is the lightest state with exotic J P C (i.e., not allowed for conventional qq mesons). The Σ + g (1P ) also allows for a pair of 1 −+ (J P C -exotic) states. The remaining unassigned states of Σ + g (1P ) are two each of 0 −+ , 2 −+ , and 2 −− , and one 3 −− . The unassigned observed states in this mass range are X(4160), Y (4274) and X(4350) (both identified with ccss above), Z 0 c (4200), and Z 0 c (4250). Of these, Z 0 c (4200) is problematic because it is a 1 +− (P = +) state, but again, its existence remains unconfirmed.
The only clear candidate in the Σ + g (2S) multiplet is Z 0 c (4430), although X(4500) [and possibly X(4700), if the allowed mass range for the band is stretched] are the two potential 0 ++ members, although again, they have been suggested as ccss states. X(4700) can also fit naturally into Σ + g (1D).
Finally, the 1 −− X(4630) and Y (4660) states fit into Σ + g (2P ), assuming the lower bound of the mass range (given above as 4750-4900 MeV) can be stretched downward slightly. In fact, the greatest difficulty of our full level assignment is the tension between the Y (4390)-X(4630) mass difference (∼ 240 MeV) and the multiplet Σ + g (2P )-Σ + g (1P ) mass difference (∼ 460 MeV). If one supposes that the Y (4390) lies at the top of the Σ + g (1P ) multiplet and X(4630) lies at the bottom of the Σ + g (2P ) multiplet, then the assignment remains sensible. Clearly a full analysis of fine-structure splittings will be necessary to assess the fate of this assumption.
We now turn to the hidden-charm pentaquarks, where to date only the two states P c (4380) and P c (4450) have been observed. Since they have opposite parities, these states must belong to distinct BO potential multiplets. While one expects the diquark mass m δ to assume the same value as that appearing in the best fits to the hidden-charm tetraquarks (Table III), the triquark mass mθ may be freely adjusted to fix one of the masses. 1 Such a fit, assuming that P c (4450) is either the J P = 5 2 + or 3 2 + state in Σ + (1P ), is presented in Table VII. However, then the P c (4380) must have J P = 3 2 − or 5 2 − . The latter assignment places it in the higher-mass Σ + (1D) level (therefore excluded), while the former assignment places it in the Σ + (1S) level, 370 MeV lower (in contrast with the observed mass splitting 4450−4380=70 MeV). Such a huge mass difference appears to be impossible to accommodate simply by using fine-structure effects of a natural size that place P c (4380) at the top of its band and P c (4450) at the bottom of its band.
Much more natural for matching with the known spectroscopy is to fix P c (4450) as the J P = 3 2 − state in Σ + (2S) (not tabulated) and to identify P c (4380) as the 5 2 + state in Σ + (1P ) , since the Σ + (2S)-Σ + (1P ) multiplet-average mass difference is then calculated to be only 200 MeV. Clearly, measuring P for either of the observed states would distinguish these scenarios. In fits with this level assignment, we notably find mθ ≈ 1.9 GeV, which is only slightly larger than m δ . We also predict Σ + (1S) hiddencharm pentaquark ground states in this fit to lie near 3880 MeV; such states would be stable against decay to J/ψ N and possibly even to η c N . They would decay through annihilation of the cc pair to light hadrons plus a nucleon, and would have narrow widths, comparable to that of η c (1S) [O(10) rather than O(100) MeV].
It should be noted that the fits in Table VII and those discussed in the previous paragraph use the full "homonuclear" BO potentials determined on the lattice, including the reflection quantum number η. However, this reflection symmetry disappears for the "heteronuclear" system. Nevertheless, the ground-state potential Σ + g in the "homonuclear" case is well separated from any other BO potential (in particular, from Σ + u ), so that nothing is lost by using it as the "heteronuclear" ground-state Σ + potential. For the higher potentials (which we did not use in the phenomenological analysis), the "homonuclear" BO potentials represent the interactions of a system with two equal masses 2µ but the same separation parameter r.
The final results involve the BO decay selection rules first discussed for exotics in Refs. [31,32] and obtained for this model in Ref. [16]. These rules assume not only that the light d.o.f. decouple from the heavy QQ pair, and that they adjust more rapidly in a physical process than the QQ, but also that the quantum numbers of the QQ and the light d.o.f. are separately conserved in the decay. As noted in Ref. [16], the strictest tests of the BO decay selection rules occur for single light-hadron decays (π, ρ, ω, φ) of exotics to conventional charmonium (Σ + g ) states. Since Z c (3900) has J P = 1 + , the observed decay Z + c (3900) → J/ψ π + requires π + to be in  ) for hidden-charm dynamical diquark-triquark states that are eigenstates of the indicated BO potentials corresponding to quantum numbers nL, for given diquark m δ and triquark mθ masses (in GeV). The particular form of the BO potential used is that given by lattice simulations JKM [26] or CPRRW [27], or (for the Σ + g potential) by the Cornell potential fit to conventional charmonium BGS [41]. Also predicted are the corresponding expectation values for the length scales 1/r −1 and r (in fm). Fixing m δ from the corresponding simulation in Table III  an even partial wave, and furthermore requiring conservation of heavy-quark spin [s QQ = 1 for J/ψ, and hence also for Z + c (3900)], Ref. [16] found Π + u (1P )-Σ − u (1P ) to be the most likely home for Z c (3900). However, as we have seen, the Π + u (1P )-Σ − u (1P ) level lies ∼ 1 GeV above the ground-state level Σ + g (1S), which is untenable for phenomenology.
The situation with the vector-meson decays is even worse, as the decay selection rules, when strictly applied as above, require the introduction of BO potentials beyond those listed in Table I. Lattice simulations predict these levels to lie even higher in mass (> 1 GeV) above Σ + g (1S). And in the pentaquark decays, either Π + (1D) or Σ − (1S), which are highly excited levels, is given as the favored home for the P = + candidate. The strictest application of the BO decay selection rules appears to conflict with the known spectroscopy.
The simplest way to resolve such issues is to note that the evidence for the conservation of QQ spin in exotics (as discussed in Ref. [16]) is imperfect, meaning that the requirement of separate conservation of QQ spin and light d.o.f. quantum numbers, which forced unacceptably high BO potentials to appear, may also be called into question. More precisely, in contrast to conventional quarkonium states, the exotics do not obviously occur in eigenstates of heavy-quark spin-symmetry. A better approach, as suggested by this work, appears to be that of obtaining the spectroscopy in a robust calculation, and then from the observed decays identifying the behavior of the states' internal quantum numbers-as is done for conventional quarkonium transitions.

V. APPROXIMATIONS OF THE MODEL
We have modeled mass eigenstates formed from a δ-δ (orθ-δ) pair of sufficient relative momentum to create between them a color flux tube of substantial spatial extent, but not so large as to induce immediate fragmentation of the flux tube. In this section, almost all the comments applied to δ-δ systems also apply toθ-δ systems.
The first and most obvious question is whether the quantized states of such a configuration of dynamical origin are best described in terms of the static configuration provided by lattice simulations. We have argued that the transition from the former to the latter paradigm is facilitated by the WKB approximation, specifically by the enhancement of the amplitude when the δ-δ system approaches its classical turning point.
The original result [13] r = 1.16 fm for the Z − c (4430) spatial extent, at which point the δ-δ pair in the process B 0 → (ψ(2S)π − )K + comes completely to rest, is only slightly smaller than the 1.224(15) fm lattice calculation of the string-breaking distance very recently presented in Ref. [51]. However, the result of Ref. [13] explicitly depends upon available phase space for the δ-δ pair, and hence upon m B 0 and m K + (in addition to m δ ). Of course, the mass eigenvalue of a ccdū bound state should depend only upon its internal dynamics, and not upon the details of the process through which it is produced; therefore, 1.16 fm should be viewed as a theoretical maximum for the possible size of the exotic state Z − c (4430), in contrast to the smaller values of r computed above. In particular, the constituents of the actual bound state should carry nonzero internal kinetic energy, which the naive calculation ignores. Since lattice-calculated static potentials provide the best available ab initio information on the nature of gluonic fields of finite spatial extent, they provide the most natural framework for modeling δ-δ bound states, even ones of dynamical origin.
The next obvious approximation is that this model assumes (effectively) structureless, pointlike δ quasiparticles. A true diquark quasiparticle-setting aside the fact that (like a quark) it carries nonzero color charge and therefore is a gauge-dependent object-should have a finite size, comparable to that of a heavy meson (a few tenths of a fm). But then, the states obtained above have a natural size only 2-4 times larger, meaning that the notion of a tetraquark state with well-separated components comes into question. Corrections that probe the robustness of the present results by including the finite size of diquarks as a perturbation are planned for subsequent work [52].
Another consequence of the assumption of structureless diquarks is the absence of both spin-and isospindependent effects. Each row of Tables I and II lists all the eigenstates of specific quantum numbers n and L for a particular BO potential Γ, which are degenerate at this stage of the calculation. Inclusion of the requisite finestructure corrections is necessary to lift the degeneracies and to produce a full spectrum of states. At this juncture, the present numerical results are identical to those one would obtain by using the methods of Ref. [21] for hybrid mesons, except with the heavy-quark mass replaced by the somewhat heavier diquark/triquark mass, and with the full quantum numbers for the states obtained only after including the light-quark spins and isospins.
The most important fine-structure corrections identified here fall into two categories: First are the spin-spin interactions within each of δ andδ (orθ); their importance in understanding the fine structure of the exotics spectrum was first emphasized in Ref. [11]. Second, since each of δ andδ contains a light quark (or two forθ), one expects in general long-distance spinand isospin-dependent corrections to modify the spectrum; without the latter, the δ-δ states would fall into degenerate uū, ud, dū, dd quartets rather than into the experimentally observed I = 0 singlets and I = 1 triplets. Were the exotic states instead composed of molecules of two isospin-carrying hadrons, a natural differentiation in the spectrum based upon isospin would arise [53], e.g., through distinct couplings of the hadrons via (longdistance) π exchange vs. η exchange.
In contrast, the dominant interaction in the dynamical diquark model between δ andδ (orθ) occurs through the color-adjoint flux tube. However, even in this case one can identify isospin-dependent interactions through the exchange of colored π-like quasiparticles, owing to a variant of the Nambu-Goldstone theorem of chiralsymmetry breaking originally discussed in the context of color-flavor locking [54]. Thus, one expects modifications to the spectrum arising from long-distance spin-and isospin-dependent interactions. As suggested in Ref. [16], lattice simulations in which the heavy, static sources are assigned not only color and spin, but also isospin, quantum numbers would have excellent investigative value for this scenario. Incorporation of both the diquarkinternal and long-distance δ-δ fine-structure corrections are planned for the next refinement of model calculations [52].
Absent in the discussion up to this point is perhaps the most consequential of all corrections for the meson sector: For most J P C quantum numbers, the physical heavy hidden-flavor meson spectrum likely contains not only possible δ-δ states, but also ordinary QQ quarkonium, as well asQQg hybrids, in addition to molecules of heavy-meson pairs. Coupled-channel mixing effects to include all of these states can have a profound effect on the observed spectrum. For example, a commonly held view in the field [1][2][3][4][5][6][7][8][9] is that the peculiar properties of the X(3872)-particularly, its extreme closeness to the D 0 -D * 0 threshold, its small width, and its substantial collider prompt-production rate-can be explained by X(3872) being an admixture of a D 0 -D * 0 molecule and the yet-unseen χ c1 (2P ) charmonium state. The addition of δ-δ states clearly makes the complete spectrum all the more rich and complex. At this stage of the dynamical diquark study, we do not attempt to address this intricate and deeply interesting problem.

VI. DISCUSSION AND CONCLUSIONS
We have produced the first numerical predictions of the dynamical diquark model, which is the application of the Born-Oppenheimer (BO) approximation to the dynamical diquark picture. In turn, this picture describes a multiquark exotic state as a system of a compact diquark δ and antidiquarkδ for a tetraquark (or triquark θ for a pentaquark) interacting through a gluonic field of finite extent. Using the results of lattice simulations for the lowest BO potentials, we have obtained the mass eigenvalues of the corresponding Schrödinger equations (both uncoupled and coupled), and found that all known hidden-charm multiquark exotic states can be accommodated by the lowest (Σ + g ) potential, for which the gluonic quantum numbers are J P = 0 + . In this sense, our explicit calculations support a type of "quark-model" classification of the lowest multiquark exotics, in which one obtains the tetraquark or pentaquark quantum numbers by combining δ andδ (orθ) quantum numbers and their relative orbital angular momentum, exactly as one does for qq mesons.
Each level Γ(nL) for each BO potential Γ produces a distinct mass eigenvalue, but differences due to the spin (and isospin) quantum numbers of the δ,δ(θ) are ignored in this calculation, meaning that each mass eigenvalue corresponds to a degenerate multiplet of states with various I G , J P C quantum numbers. We estimated the maximum size of the neglected fine-structure splittings (∼ 150 MeV), and found that the spectrum of tetraquarks should consist of a lowest [Σ + g (1S)] band, all members of which have P = +, followed by a gap of about 100 MeV, and then a Σ + g (1P ) band of P = − states, then another gap and overlapping Σ + g (2S) and Σ + g (1D) bands of P = + states, and finally a band of Σ + g (2P ), P = − states. The order for pentaquark states is the same, except that the reflection quantum number "g" is no longer present, and the P eigenvalues are opposite those for tetraquarks. Many higher levels are predicted, but are not yet needed to accommodate currently observed exotics.
As of the present, the computed band structure with alternating P is supported by the known states, such that ] multiplet. The exceptions are the 1 ++ Y (4140), which lies in the first band gap, but may be a ccss state and thus fall outside the current analysis; and the 1 +− Z 0 c (4200) (if its existence is confirmed), since it lies in the region of the Σ + g (1P ), P = − band. For the pentaquarks, the small P c (4450)-P c (4380) mass difference is most easily accommodated by assigning the heavier state to the 2S band and the lighter one to the 1P band, leading to the prediction of 1S-band hidden-charm pentaquarks that may be stable against decay into charmed particles.
However, the BO decay selection rules, based upon separate conservation of heavy quark-antiquark and light degree-of-freedom quantum numbers in observed decay processes, mandate that known, low-lying exotics must appear in highly excited BO potentials, and thus contradict observation. We propose that the selection rules fail badly because they are based in part upon assigning the exotics to heavy-quark spin eigenstates, for which the experimental evidence appears to be quite mixed.
To develop the model further, one must perform a detailed analysis of the fine-structure corrections (both spin and isospin dependence) to determine whether the specific level structure suggested by data is supported by experiment. One must also include effects arising from finite diquark (triquark) sizes. These refinements will be implemented in subsequent work to be carried out by this collaboration.
Additional improvements rely upon, first of all, a reassessment of lattice simulations: How much does the spin of the heavy sources ( 1 2 for c andc, 0 or 1 for δ,δ), or the light-flavor content in the diquark/triquark case, modify the BO potentials? Finally, this work assumes that every exotic state in the charmoniumlike system is a δ-δ or δ-θ state. Ignored completely in this analysis are the possibilities that some of these states are high conventional cc, or that some are genuinely hadronic molecules, or threshold effects, or even mixtures of these types. Only a global analysis including observables such as detailed branching fractions and production lineshapes can truly disentangle the full spectrum. Γ ≡ Λ η . Λ is an angular momentum projection, and and η are inversion parities. With reference to Fig. 1, one defines a plane containing the axisr and a unit normal n to the plane. Denoting the total angular momentum of the light d.o.f.-a conserved quantity, thanks to the decoupling in the BO approximation-as J light , and the orbital angular momentum of the heavy d.o.f. as L QQ , one obtains the total orbital angular momentum of the system, Sincer·L QQ = 0, the axial angular momentumr·J light = r·L of the light d.o.f. is a good quantum number for the whole system, and its eigenvalues are denoted by λ = 0, ±1, ±2, . . .. One further notes that Analogous to the use of labels S, P, D, . . . for the quantum numbers L = 0, 1, 2, . . ., one denotes potentials with the eigenvalues Λ = 0, 1, 2, . . . by Σ, Π, ∆, . . .. Reflection R light of the light d.o.f. through the plane with unit normaln (which is spatial inversion P light of the light d.o.f., combined with a rotation by π radians about n originating from the δ-δ midpoint) transforms λ → −λ.
Since a glance at Fig. 1 reveals that the physical system (and hence its energy) must be invariant under R light , one finds that the energy eigenvalues must be a function only of Λ ≡ |λ|.
The eigenvalues of R light itself are denoted = ±1. Strictly speaking, specifying is essential only for Σ potentials, for which the = ±1 states may have distinct energies; however, one may also form eigenstates for (degenerate) Λ > 0 BO potentials in the same way as one uses linear combinations of f (+x) and f (−x) to form both even and odd functions from an arbitrary function f . Following this procedure [16], one obtains BO potentials for all values of Λ , in which case one also finds P light = (−1) Λ in the "homonuclear" case.
Even for the "homonuclear" system, a complete inversion of coordinates of the light d.o.f. through the midpoint of the QQ pair does not by itself produce an equivalent state to the original; rather, one must also exchange the QQ pair, or instead, perform charge conjugation C light upon the light d.o.f. Thus, the full system in the "homonuclear" case is physically invariant under the combination (CP ) light , and its eigenvalues η = +1, −1 are labeled as g, u, respectively.
In total, the three eigenvalues Λ η completely specify the irreducible representation of D ∞h for any "homonuclear diatomic" system. 2 In the "heteronuclear diatomic" case (such as B c tetraquarks or the triquark-diquark pentaquark), the (CP ) light symmetry is lost (leaving the symmetry group C ∞v ), but the good quantum numbers Γ ≡ Λ remain.
innermost root of any of the eigenvalues of −Q(r), and the outer classical turning point is defined as the outermost root of any of the eigenvalues of −Q(r). This definition of the classical turning points ensures that no nodes are missed in the counting procedure due to having skipped over a portion of the classically allowed region of one of the potentials. Moreover, such a definition is established through multichannel generalizations of the WKB approximation [56].
To find the desired value energy eigenvalue E N for which the solution U N has N nodes, one first chooses a window of E values that is bounded below by E low , bounded above by E high , and that contains E N . One then counts the number of nodes N of det U(r) in the classically allowed region by numerically integrating Eq. (B3) at E = E mid = (E low + E high )/2. Once the outer bound has been reached, one then sets E low = E mid if N ≤ N , or sets E high = E mid if N > N . This procedure is repeated until E high − E low meets a pre-specified tolerance.

Renormalized Numerov Integration Procedure
We now derive the renormalized Numerov method of Ref. [39]. Discretize the radial coordinate r such that the initial value r 0 is located to the left of the inner classical turning point and is sufficiently close to the origin, while the final value r f is located sufficiently to the right of the outer classical turning point. Moreover, choose the grid to be uniformly spaced, with r i+1 − r i = h for all i ∈ {0, 1, 2, ..., f }, and h some arbitrarily small real number. Moreover, let U(r i ) ≡ U i .
Numerical integration of Eq. (B3) can be accomplished by a coupled generalization of the Numerov recurrence relation: The recurrence relation of the renormalized Numerov scheme results from three substitutions: The first two substitutions turn Eq. (B5) into and the last substitution renormalizes Eq. (B8) as Equation (B9) provides a stable and efficient method for propagating U i in regions where the entries of Q diverge to ±∞, since R i ∼ Q i+1 Q −1 i . In most situations, one can choose R −1 0 = 0; however, this choice can present complications [39]. We instead set R −1 0 to be If r 0 = h, then this choice is equivalent to supposing that R −1 −1 = 0 precisely at the origin, r = r −1 = 0. We now describe a method for counting the number of nodes along the integration without explicitly monitoring det U i at each integration step. First suppose that there is only one node of det U between r i+1 and r i . From the last of Eq. (B7), one has that Since det(1 − T i+1 ) > 0 for any r i in the classically allowed region, one has that det R i+1 < 0, since det U encountering one node between r i+1 and r i implies that det U i and det U i+1 have opposite signs. It follows that one may monitor det R i at each integration step to count the number of nodes. Such a method is effective if only one node exists between each pair of grid points, as the existence of even numbers of nodes between any two grid points (due to degeneracies, for example) implies that det U i+1 and det U i have the same sign, and hence det R i+1 > 0. One can avoid missing nodes by instead monitoring the individual eigenvalues of R i at each integration point. If det R i < 0, then there is an odd number of negative eigenvalues. However, if det R i > 0, there is an even number of negative eigenvalues, and an incorrect node count occurs. Therefore, if one instead monitors each time one of the eigenvalues of R i changes sign, there is no danger of missing a node. This procedure is equivalent to monitoring the signature of R i , and one can therefore equivalently count the number of times one of the diagonal entries of the U matrix in an LU -decomposition of R i changes sign.

Calculation of Expectation Values
From nondegenerate perturbation theory, one learns that if the Hamiltonian H splits into some reference Hamiltonian H (0) and a small perturbation H , then the energy eigenvalues can be calculated perturbatively. To first order in , this perturbative expansion is: The expectation value of H is then found as the limit This result is true for any Hermitian operator H . Therefore, Eq. (B13) provides a method [40] for numerically calculating the expectation values of arbitrary operators using just the energy eigenvalues of the original problem and the energy eigenvalues of Eq. (B12), so long as these energies are non-degenerate. This observation is very powerful, because it means that one does not need to numerically compute the wave functions of Eq. (B3) in the calculation of expectation values.