Heavy hadronic molecules with pion exchange and quark core couplings: a guide for practitioners

We discuss selected and important features of hadronic molecules as one of promising forms of exotic hadrons near thresholds. Using examples of $D \bar D^*$ systems such as $X(3872)$ and $Z_c$, emphasis is put on the roles of the one pion exchange interaction between them and their coupling to intrinsic quark states. Thus hadronic molecules emerge as admixtures of the dominant long-range hadron structure and short-range quark structure. For the pion exchange interaction, properties of the tensor force are analyzed in detail. More coupled channels supply more attractions, and heavier constituents suppress kinetic energies, providing more chances to form hadronic molecules of heavy hadrons. Throughout this article, we show details of basic ideas and methods.


Exotic phenomena
Since the discovery of X(3872) in 2003 at Belle/KEK and BaBar/SLAC, many candidates of new hadrons have been observed [1,2]. ‡ Their observed properties such as masses and life times are not easily explained by conventional methods and models of QCD. Thus they have been called exotic hadrons or simply exotics. Historically, exotic hadrons of multiquarks are already predicted by Gell-Mann in his original work of the quark model [3]. The states with quantum numbers which are not accessed by the standard quark model, mesons as quark and antiquark (qq) and baryons as three quarks (qqq), are often referred to as manifest or genuine exotics. In this regard, X(3872) is not manifestly exotic, but it shows up with many unusual properties. By now X(3872) has been observed by many experimental facilities, and is well established with its quantum numbers determined by LHCb, J P C = 1 ++ [4]. In the latest PDG data base, more than thirty particles are listed as candidates of exotic hadrons. Many of them are considered to contain charm quarks as constituents, while some of them only light quarks [5].
Those exotic candidates are observed near thresholds. For charmonia (cc pairs), the thresholds are the energies of D ( * )D( * ) above which an excited charmonium may decay into a D ( * )D( * ) pair §. Therefore, near the threshold region systems may contain an extra light qq (q = u, d) in addition to the heavy quark-antiquark pair (cc or bb). The natures of hadrons near thresholds and of those well below thresholds are qualitatively different each other. Quarkonium like states of cc or bb well below the threshold are essentially non-relativistic systems of a slowly moving heavy quark pair [6]. In contrast, exotics containing both heavy and light quarks may show up with various configurations such as compact multiquarks [7][8][9][10], hadronic molecules [11,12] and hybrids or complicated structure of quarks and gluons [13][14][15]. The question of how and where these different structures show up is an important issue in hadron physics and has been discussed in references [16,17].
In multi-quark systems, the quarks may rearrange into a set of colorless clusters. For instance, a hidden charm four quarks rearrange as ccqq → (cc)(qq) ∼ J/ψπ, or (cq)(qc) ∼ DD * . J/ψπ dominantly appears in decays because the pion is light and unlikely to be a constituent of hadrons. In the chiral limit massless pions play just as chiral radiations. In contrast, DD * may form quasi-stable states if suitable interactions are provided via light meson exchanges, in particular pion exchanges between light quarks. This is the crude but basic idea of how hadronic molecules are formed. The idea of hadronic molecule is dated back to the discussion of Λ(1405) as aKN molecule [11], and more were conjectured in the context of cc productions after the discovery of J/ψ [18].

Clusterization
The rearrangement of multi-quarks shares a general feature of clustering phenomena by neutralizing the original strong force among the constituents. Then among the neutralized clusters only relatively weak force acts. In the present case the color force is strong, while the meson exchange force is weak. In this clustering process hierarchies of matter, or separation of the energy scale occurs. Strong color force is of order hundred MeV while the weak meson exchange force is of order ten MeV. This qualitatively explains how hadronic molecules are bound with a binding energy of order ten MeV. In table 1, several candidates of hadronic molecules are shown. From these small binding energies one can verify that the spatial sizes of these systems are of order one fm or larger. With this inter-distance, the constituent hadrons in molecules can maintain their identity. Table 1. Candidates of exotic hadrons near thresholds, where ∆E = Mass − Threshold mass. Data are taken from PDG [5] for Λ, X and Z c , and from [19] for P c 's. The lower raw of Λ(1405) is for the higher pole of the two-pole scenario.

State
Mass ( Analogous phenomena are found in nuclear excited states with alpha clusters developed. A well known example is the Hoyle state, the first 0 + excited state of 12 C [20]. The formation of alpha clusters near the threshold of alpha decays is known as the Ikeda rule which predicts the dominance of alpha cluster components in nuclear structure in the threshold region of 4N (N = integer) nuclei [21]. Threshold phenomena are now regarded as universal phenomena and are discussed in the context of universality which covers various systems from quarks to atoms and molecules [22,23].
By now there are many articles which discuss hadronic molecules including comprehensive reviews [24]. Here in this article we do not intend to list all of the previous works, but rather focus on limited subjects that we believe important for the discussions of hadronic molecules. To elucidate the points we discuss D ( * )D( * ) systems, especially for X(3872) and some related states. We do not discuss baryons; for Λ(1405) there are many discussions including the summary one in PDG [5], for P c 's discussions has just started and we need more studies to make conclusive statements. In this way, this article is not inclusive. However, we try to emphasize general features by using a few specific examples. We also try to show some details of how basic ingredients are derived. Sometimes, we discuss items which are by now taken for granted. We think that this strategy is important because many current discussions seem to be based on ad hoc assumptions, and many explanations and predictions depend very much on them.

Pions and interactions
Now the most important ingredient is the interaction which is provided by light meson exchanges at long and medium distances. Among them best established is the one pion exchange potential. The pion is the Nambu-Goldstone boson of the spontaneous symmetry breaking (SSB) of chiral symmetry [25,26]. Its interaction with hadrons is dictated by the low-energy theorems. The leading term is the Yukawa term of πhh (h, h : hadrons). By repeating this twice, the one pion exchange potential emerges in the t-channel as shown in (t) of figure 1, where general structure of two-body amplitudes is shown. Hence, in hadronic molecules, pions play a role of the mediator of the force between constituent hadrons. Microscopically, the pion couples to the constituent quarks which are dynamically generated by SSB. Combined with the quark model wave functions of hadrons, the coupling strengths as well as form factors are estimated, schematically by where the sum is taken over the light quarks (i) in the hadrons as shown in figure 2. This method works qualitatively well for nuclear interactions and is now extended to other hadrons for the study of hadronic molecules. Other meson interactions of such as σ, ω and ρ mesons are also employed but with more parameters are needed. In fact, the masses of these mesons are of the same order of the inverse of hadron size, and their contributions may be masked by the form factors. Thus, the pion interaction is the best known and under control. In most part of this paper, we test models of hadronic molecules with the pion interaction. Another feature of the pion interaction is in its tensor structure. This is the consequence of SSB of chiral symmetry which leads to pseudoscalar nature of the pion with spin-parity J P = 0 − . Therefore, the coupling structure of the pion to hadrons is of σ · r type. This leads to the tensor force causing mixing of orbital motions of different angular momenta by two units. This provides extra attraction which contributes significantly to the formation of molecules. Although the importance of the tensor force has been long recognized in nuclear physics [27,28], quantitative understanding is being progressed by the developments in the microscopic treatment of many-body systems and in computer power [29][30][31][32].
In addition to the pion exchange interaction at long distances, we also discuss s-channel interactions at short distances where the incoming hadrons merge into a single hadron (one-particle) as an intermediate state (see figure 1). Hence this process leads to the mixing of configurations, an extended molecular structure of two particles and a compact one-particle state. The problem is also related to the question of the so-called compositeness [33][34][35][36][37]. We emphasize an importance of such a mixing for X(3872); a molecular component of DD * at long distances and a cc component at short distances [38].

Contents of this paper
This paper is organized as follows. In section 2, we show how the Yukawa vertex of a pion to heavy hadrons are derived. Coupling constants in different schemes are discussed in some detail. Estimation in the quark model is also discussed. In section 3 one pion exchange potential is derived with emphasis on general features of the potential. Special attention is payed to the tensor structure and form factors. A non-static feature is also discussed when the mass of the interacting hadrons changes, which is taken into account by an effective mass of the pion. In section 4, we discuss the structure of X(3872). After briefly reviewing experimental status, we discuss the molecular nature made by the one pion exchange. An important role of the short distance dynamics is also discussed, and consider a mixing structure of hadronic molecule coupled by a compact quark cc component. Lastly, discussion is given for Z c . In section 5 we summarize and discuss prospects.

Heavy hadron interactions
Hadronic molecules are composite systems of hadrons which are loosely bound or resonate. "Loosely" means that the binding or resonant energies are small as compared to the QCD scale of Λ QCD ∼ several hundred MeV, which is relevant to intrinsic structure of hadrons by quarks. In such a situation, the constituent hadrons can retain their intrinsic structure in the molecules. The interaction among the hadrons is colorless and its dominant part is expected to be dictated by meson exchanges. Among them, pion exchange interaction is the best under control. The pion couples to the light u, d quarks, and their dynamics is determined by the nature of Nambu-Goldstone bosons of spontaneous breaking of chiral symmetry. This is the case if hadrons contain light u, d quarks as constituents such as protons, neutrons and also heavy open flavor hadrons such as D mesons (cū) and Σ c baryons (cuu).
In this section, we discuss basic interactions for heavy mesons, that is the Yukawa vertices for P and P * with the pion, where P stands for D orB meson, and P * for D * or B * . We also employ the notation P ( * ) for either P or P * . In addition to chiral symmetry features associated with light quarks, heavy quark spin symmetry also applies in the presence of heavy quarks (either charm c or bottom b quark) [39,40]. In particular, heavy quark spin symmetry relates the mesons with different spins under spin transformations. For example, P -meson of spin-parity J P = 0 − is a spin partner of P * meson of J P = 1 − ; they are the same particles under heavy quark spin symmetry.
To implement the aspects of heavy quark spin and chiral symmetries in the effective Lagrangian, we shall quickly overview several issues such as a convention for heavy quark normalization, representations of heavy fields for D and D * mesons, and their properties under the heavy quark spin and chiral symmetry transformations. We also discuss how the relevant coupling constants are determined. We see that the constituent picture of the light quark coupled by the pion consistently describes the decay properties of the D ( * ) mesons as well as axial properties of the nucleon.

Heavy fields
When considering quantum fields of heavy particles of mass m H , it is convenient to redefine the effective heavy fields in which rapidly oscillating component in time, exp(−im H t), is factored out. For QCD the "heavy" means that m H is sufficiently larger than the QCD scale, Λ QCD m H , and the heavy quarks almost stay on mass-shell with quantum fluctuations being suppressed. In accordance with the redefinition of the field, the normalization of the effective heavy fields are naturally modified from the familiar one of quantum fields by the factor √ m H .
To show this point let us consider the standard Lagrangian for a complex scalar meson field of heavy mass The factor 1/2 is recovered when using the real components φ 1,2 . From this Lagrangian, the current is given by The field expansion may be done as in the standard conventions.
In the heavy mass limit m H → ∞, the particle is almost on-mass shell, and it is convenient to define the velocity v µ , v 2 = 1 which defines the on-shell momentum m H v µ .
Thus the momentum fluctuation k µ around it is considered to be small, k µ m H Moreover the Hilbert space of different heavy particle numbers decouple because particle-antiparticle creation is suppressed in the considering energy scale. Hence we define the heavy field by This means that the energy ofφ(x) is measured from m H . Inserting the relation into (2), we find and for the current, where in both equations we have only shown the leading term of order O(m H ). Note that the mass term in (2) disappears in the Lagrangian as expected because the energy is measured from m H . By absorbing the factor m H into the field as then we have In this convention, the heavy (boson) field φ H carries dimension 3/2 in units of mass, unlike 1 in the standard boson theory. In this paper, as in many references, we follow this convention, while we also come back to the ordinary convention of dimension 1. In terms of one-particle states, these two conventions correspond to different normalizations [40] v|v Moreover the one-particle to vacuum matrix element of the field is given by

Interaction Lagrangian
Let us consider the pseudoscalar P and vector P * mesons as a pair of heavy quark Q and light antiquarkq in the lowest S-wave orbit. Having overviewed the features of heavy particles in the previous subsection, the heavy meson field is defined in the frame of a fixed velocity v µ and contains multiplet of spin 0 pseudoscalar and spin 1 vector mesons. They are (D, D * ) in the charm sector, and (B,B * ) for the bottom sector. For convenience, namings and quark contents of various mesons are given in table 2. In this article through out, we place symbols without bar on the left of those with bar. This is a convention that is consistent with quark model calculations. Mesons P ( * ) K ( * )K ( * ) D ( * )D( * ) B ( * )B( * ) Quark contents qs sq cq qc qb bq A convenient way to express such heavy meson fields (including antiparticles) is where P a and P * a µ carry an index of isospin 1/2, a = 1, 2 ∼ u, d. The factor (1 + v /)/2 is a projector to constrain the heavy quark velocity at v. We employ the convention for γ-matrices and For later convenience, we express the meson fields explicitly in terms of the quark fields, Let us consider D + =diγ 5 c, D * + µ =dγ µ c, where d, c express the Dirac fields for the down and charm quarks. Under charge conjugation transformations, we can verify that In this convention, again using the notation P 's, the operators for the charge conjugated anti-particles are P † for pseudoscalars and −P * † µ for vectors. The corresponding states are defined by These relations will be used when forming eigenstates of charge conjugation of molecules formed by P ( * ) andP ( * ) mesons. The spin multiplet nature of P and P * is verified by writing (14) in the rest frame v µ = (1, 0, 0, 0), where only the spatial components remain for the vector meson as its degrees of freedom, Here the isospin label a is suppressed, since it is irrelevant under spin transformations. The combination P + σ · P * indicates that P and P * are the spin multiplet of (1/2, 1/2) representation of the heavy and light spin group SU (2) Q × SU (2) q . They transform under the heavy and light spin transformations as where Σ are the four component spin matrices defined by Σ k = i 2 ijk [γ i , γ j ], and θ Q and θ q are the rotation angles of heavy and light quark spins, respectively. The diagonal part of θ Q = θ q corresponds to the total spin rotation.
Chiral symmetry property of the heavy meson field (14) is inferred by the constituent nature of the light quark q. It is subject to non-linear transformations of chiral symmetry [41][42][43][44]. In this article, we consider two light flavors and therefore SU (2) L ×SU (2) R is the relevant chiral symmetry group, where the left (L) and right (R) transformations act on the two isospin groups. Explicitly, the quark field q of isospin 1/2 are transformed as where the isospin SU (2) matrix function h(g R , g L , π(x)) characterizes non-linear chiral transformations determined by global chiral transformations of g L,R ∈ SU (2) L,R at the pion field π(x). Therefore, the heavy meson fields of isospinor transforms under chiral symmetry transformations as The isovector pion field parametrizes unitary matrices as which linearly transforms as The non-linear transformation for the pion field is then conveniently expressed in terms of the square root which is subject to Here f π is the pion decay constant for which our convention is f π ∼ 93 MeV. The ξ-field defines the vector and axial-vector currents which are transformed as Note that the currents (28) are anti-Hermitian. Moreover, the vector and axial-vector currents are of even and odd power with respect to the pion field (see (31)), with properly satisfying the parity of the currents, 1 − and 1 + , respectively. Once established the heavy quark spin and chiral transformation properties we can write down the invariant Lagrangian. To the leading order of derivative expansion, we find By expanding the vector and axial-vector currents, V µ and A µ with respect to the pion field, the vector current leads to the pion-hadron interaction of the so-called Weinberg-Tomozawa interaction, while the axial vector current to the Yukawa coupling. The former strength is determined by the pion decay constant while the latter contains one unfixed parameter, the axial coupling constant g A . In the present scheme it corresponds to the one of the constituent quark as discussed in subsection 2.5. By inserting the expansion (31), we find the P P * π and P * P * π interaction Lagrangians as As anticipated, the strengths of these interactions are given by one coupling constant g A which is a consequence of heavy quark spin symmetry.

Meson decays I, D * → Dπ
To see the use of (32) together with the heavy quark normalization, let us consider the simplest and important example for meson decays, D * + (λ, p) → D + (p )π 0 (q), where λ labels the polarization of D * + . The relevant matrix element for these charged states is where ν (λ) is the polarization vector of the D * meson, and we have used the relation for the matrix elements; Here the masses of D and D * mesons are regarded sufficiently heavy and set equal to m H . For later convenience, we summarize the other matrix elements which are needed for the computations of the transition amplitudes P P * → P * P andPP * →P * P in table 3. By using the relations of (19), one can verify these relative signs. Table 3. Relative signs of various P ( * ) P ( * ) π couplings including antiparticles.
Now the decay width is computed by Note that in the heavy mass limit, m D * ∼ E D ∼ m H , heavy meson mass m H dependence in E D (p ) in the denominator and that in |V | 2 in the numerator cancel. The heavy mass independence is reasonable because the decay D * + → D + π 0 occurs through the spin flip of the light quark which should not depend on the heavy mass of the spectator heavy quark. Fixing the polarization (λ) ν of initial D * meson and performing the phase space integral together with the angle average of |q ν ν | 2 → |q · | 2 → q 2 /3, we find Using the experimental data we find This value is obtained in the limit m D , m D * → ∞ using the formula (37). If we take into account their finite values, we find g A ∼ 0.53. This estimates uncertainties of few percent at minimum in the discussions based on the leading terms of the heavy quark symmetry.

Meson decays II, K * → Kπ
Now it is instructive to demonstrate another textbook like calculation for the decay K * → Kπ, which is the analogue of D * → Dπ by replacing the charm quark by the (anti)strange quark. A commonly used Lagrangian in a flavor SU(3) symmetric form is written as where g is the coupling constant of a vector meson with two pseudoscalar mesons. In the SU(3) limit it is the ρ → ππ coupling constant and is given to be g ∼ 6 [44][45][46]. In this convention, the normalization is such that there are 2E particles in a unit volume. The matrix element of the above Lagrangian is then (again for the neutral pion decay) Therefore, we find the formula Here if we break SU (3) symmetry and take heavy mass limit for the strange quark, m K * ∼ E K → m H , we find the total decay width By using the experimental data we find which is close to the coupling constant of the ρ meson decay, g(ρ → ππ) ∼ 6. Comparing Eqs. (37) and (44) we find This is nothing but the generalized Goldberger-Treiman relation, implying that the coupling constant g scales as the meson mass m H , when g A is independent of m H as we shall discuss in the next subsection. In other words, flavor symmetry breaking for the coupling constant g defined by the Lagrangian (40) scales as that of the corresponding meson masses. As a matter of fact, estimation of g for the decay of D * implies g ∼ 12, obviously different from g ∼ 6 for the decay of K * . Indeed this is explained by the mass difference of K * and D * mesons, m K * ∼ 890 MeV and m D * ∼ 2000 within about 10 % accuracy.

Quark model estimate
In this subsection we show that the coupling constant g A in the Lagrangian (32) is nothing but the quark axial-vector coupling constant in the non-relativistic quark model. Let us start with the πqq Lagrangian in the axial vector type where we have denoted the quark axial-vector coupling by g q A , and ignored isospin structure for simplicity. This Lagrangian operates to the light quark-pion vertex at position x 1 as shown in figure. 3, where assignments of various variables are also shown. For example, the the center of mass and relative coordinates are defined by with m, M being the masses of the light and heavy quarks, In the non-relativistic limit, the matrix element of (48) forD * (at rest) → D(−q)π(+q) is given as [47] where χ i,f are the quark wave functions of the initial D * and final D meson, (ω π , q) the energy and momentum of the pion, p i the momentum of the quark in the initialD * meson. For a notational reason in the definition of the quark model wave function as explained below, here we consider the decay ofD * rather than D * . The wave functions χ i,f are written as a product of the plain wave for the center of mass and internal part including spin, φ i,f (r) where ω i,f are the energies of the initialD * and finalD mesons. Expressing p i by the relative momentum p r as we can perform the t and X-integral leading to the total energy-momentum conservation, leaving the r integral as where effective momentum transfer is defined bỹ and the relative momentum p r is replaced by −i∇ r . Using the harmonic oscillator wave function for both φ i and φ f , with the size parameter α, after some computation, we find For small |q|, we may set the form factor F (q 2 ) → 1. For the spin matrix element σ ·q , we evaluate the transition Having the spin wave functions and with the understanding that the spin operator acts on the first (left) spin state for the light quark, we find S = 0|σ ·q|S = 1 = 1 .
These results are compared with the matrix element (33), where we may set q µ = (0, 0, 0, |q|) and q µ µ = −|q|. Suppressing the second term in the heavy quark limit M → ∞, we find the relation in the limit |q| → 0.
Usually in the quark model g A is assumed to be unity, g q A = 1. However, it is know that this overestimates the axial couplings of various hadrons. For the nucleon it is known that the quark model predicts g N A = 5/3 [44]. In the quark model, the nucleon g A is defined to be the matrix element of spin and isospin operator n=1,2,3 Therefore, effectively the reduction of g q A ∼ 0.7 is needed to reproduce the data. Similarly heavy baryon transitions such as Σ ( * ) c → Λ c consistently implies small g q A [47]. How baryon g A is computed is found in Refs [44,47]

Meson exchange potential
In this section, we derive meson exchange potentials for the study of hadronic molecules. Starting from the classic method for the derivation, we revisit the one pion exchange potential (OPEP) for the nucleon (N ). We find it is useful to recognize important and universal features of meson exchange potentials.

Simple exercise
Let us illustrate a simple example for a scalar field φ interacting with the exchange of a scalar π meson of mass m. The extension to the case of physical pion will be done later in a rather straightforward manner. The model Lagrangian is from which the equation of motion and a special solution for π are obtained as We note that in this example, the coupling constant g carries dimension of unit mass. The potential energy for φ is given by the energy shift ∆E due to the interaction; where in the second line we have used the solution of (63). Inserting the complete set and representing the propagator in the momentum representation, d 3 p/(2π) 3 |p p|, we find the expression Now we regard φ as a field operator and expand in momentum space. Then consider a scattering process of p 1 , p 2 → p 1 , p 2 as shown in figure 4. Taking the matrix element p 1 , p 2 |∆E|p 1 , p 2 and performing x, y, and q integrals, we find Remarks are in order.
• The energy shift ∆E (65) is for the entire volume V ∼ (2π) 3 (p 1 + p 2 − p 1 − p 2 ) → δ 3 (0), and also for the normalization of 2E φ-particles per unit volume. In the center-of-mass frame, energies of the two particles are the same and conserved, E = p 2 + M 2 , where M is the mass of φ. Therefore, The energy shift per unit volume and per particle is given by • In the non-relativistic limit for the φ particle, we can take the static limit where the energy transfer q 0 is neglected such that −q 2 → +q 2 , and E ∼ M , This defines the potential in momentum space, and in turn in coordinate space.
• Though obvious but not very often emphasized is the fact that the potential appears always attractive when the coupling square g 2 is positive. If the coupling structure has spin dependence this is no longer the case, otherwise always so. This is understood by the formula of second order perturbation theory where the intermediate state of φφπ in the pion-exchange process is in higher energy state than the initial φφ (see figure. 4) . This fact is in contrast with what we know for the Coulomb force. The reason is that the latter is given by the unphysical component of the photon, which is manifest in the sign of the metric. Figure 4. One meson exchange potential. The vertex structure of σ · q is needed for the one pion exchange potential (OPEP) of the nucleon. For φφ or N N the line widths of bold or normal are irrelevant. It will become relevant when discussing the potential for P ( * ) P ( * ) .

OPEP for the nucleon-nucleon N N
Now the most familiar and important example is the one pion exchange potential (OPEP). In this section we will discuss OPEP for the nucleon, because the nucleon system is the best established, and can share common features with heavy hadrons. Since the pion is the pseudoscalar particle, the pion nucleon coupling is given either by the pseudoscalar or axial vector (pseudovector) form, When the nucleons are on mass-shell, it is shown that the matrix elements for N (p) → N (p )π(q) in the two schemes are equivalent by using the equation of motion for the nucleon. The equivalence of the two expressions leads to the familiar Goldberger-Treiman relation In the non-relativistic limit, the equivalent matrix elements reduce to where the two-component nucleon spinors are implicit and a on π and τ -matrix is an isospin index. The extra factor 2m N on the right hand side appears due to the normalization of the nucleon (fermion) field when the state is normalized such that there are 2E ∼ 2m N nucleons in a unit volume. The positivity of the coupling square as discussed in the previous subsection is ensured by ±i in (72). Inserting these coupling structures into the general form of (68), we find the OPEP for the nucleon in the momentum space Sometimes, this form is called the bare potential because the Lagrangian (70) does not consider the finite size structure of the nucleons and pions. The OPEP depends on q, a feature consistent with the low energy theorems of chiral symmetry; interactions of the Nambu-Goldstone bosons contain their momenta. At low energies the N N interaction (73) is of order O(q 2 ). In particular at zero momentum the interaction vanishes. In contrast, when q → ∞, the interaction approaches a constant. This requires a careful treatment for the large momentum or short range behavior of the interaction.
To see this point in more detail, let us decompose the spin factor σ 1 · qσ 2 · q into the central and tensor parts, where the tensor operator is defined by The first term of (74) is the spin and isospin dependent central force, which has been further decomposed into the constant and the Yukawa terms. The constant term takes on the form of the δ-function in the coordinate space. This singularity appears because the nucleon is treated as a point-like particle. In reality, nucleons have finite structure and the delta function is smeared out. In the chiral perturbation scheme starting from the bare interaction of (73) the constant (δ-function) term is kept and higher order terms are systematically computed by perturbation. In this case, low energy constants are introduced order by order together with a form factor with a cutoff to limit the work region of the perturbation series [48,49]. To determine the parameters experimental data are needed. This is possible for the N N force but not for hadrons in general. Alternatively, in nuclear physics the constant term is often subtracted. One of reasons is that the hard-core in the nuclear force suppresses the wave function at short distances and the δ-function term is practically ineffective. Then form factors are introduced to incorporate the structure of the nucleon, and the cutoff parameters there are determined by experimental data.
In this paper, we employ the latter prescription, namely subtract the constant term and multiply the form factor. As in (74), the constant and Yukawa terms in the central force have opposite signs, and hence part of their strengths are cancelled. The inclusion of the form factor in the Yukawa term is to weaken its strength, which is partially consistent with the role of the constant term. In practice, the central interaction of the OPEP is not very important for low energy properties. Rather the dominant role is played by the tensor force. We will see the important role of the tensor force in the subsequent sections.
So far we have discussed only the OPEP. In the so-called realistic nuclear force, to reproduce experimental data such as phase shifts and deuteron properties, more boson exchanges are included such as σ, ρ and ω mesons [50][51][52]. Their masses are fixed at experimental data except for less established σ. Coupling constants and cutoff masses in the form factors are determined by experimental data of N N phase shifts and deuteron properties. The resulting potentials work well for N N scatterings up to several hundred MeV, and several angular momentum (higher partial waves). However, if we restrict discussions to low energy properties, which is the case for the present aim for exotic hadronic molecules, meson exchanges other than the pion exchange are effectively taken into account by the form factors. As discussed in the next section, we will see this for the deuteron.
Having said so much, let us summarize various formulae for the OPEP for N N . Subtracting the constant term with the form factor included we find For the form factor we employ the one of dipole type The potential in the r-space is obtained by performing the Fourier transformation as where C(r; m, Λ) and T (r; m, Λ) are given by by

Deuteron
It is instructive to discuss how the OPEP alone explains basic properties of the deuteron by adjusting the cutoff parameter in the form factor. It implies the importance of OPEP especially for low energy hadron dynamics. Furthermore, we will see the characteristic role of the tensor force which couples partial waves of different orbital angular momenta by two units. Inclusion of more coupled channels gains more attraction, and hence more chances to generate hadronic molecules. The importance of the OPEP for the N N interaction is discussed nicely in the classic textbook [53]. The deuteron is the simplest composite system of the proton and neutron. It has the spin 1 and isospin 0. The main partial wave in the orbital wave function is S-wave with small admixture of the D-wave of about 4 %. It has the binding energy of 2.22 MeV and with the size about 4 fm (diameter or relative distance of the proton and neutron). Because the interaction range of OPEP is ∼ 1/m π ∼ 1.4 fm while the deuteron size is sufficiently larger than that, the nucleons in the deuteron spend most of their time without feeling the interaction. This defines loosely bound systems and is the condition for the hadronic molecule.
The main S-wave component of the wave function ψ(r) can be written as those of the free space where µ, B and A are the reduced mass of the nucleon, binding energy and normalization constant. By using this the root mean square distance can be computed as The binding energy and the mass of the nucleon give r 2 1/2 ∼ 4 fm, consistent with the data. Now it is interesting to show that these properties are reproduced by solving the coupled channel Schrödinger equation with only the OPEP included. Explicit form of the coupled channel equations are found in many references, and so we show here only essential results. Employing the axial vector coupling constant for the nucleon g N A ∼ 1.25 and choosing the cutoff parameter at Λ = 837 MeV the binding energy is reproduced. At the same time experimental data for the scattering length and effective range are well reproduced as shown in the third raw of table 4. Note that since g A ∼ 1.25 is fixed, the cutoff Λ is the only parameter here.
Throughout this article, we define that the positive (negative) scattering length stands for the attraction (repulsion) at the threshold.
The cutoff value Λ = 837 MeV is consistent with the intrinsic hadron (nucleon) size. By interpreting the form factor related to the finite structure of the nucleon, we may find the relation The size 0.6 fm corresponds to the core size of the nucleon with the pion cloud removed.
In table 4, results are shown also for those when other meson exchanges are included [54]. By tuning the cutoff parameter Λ around the suitable range as consistent with the nucleon size, low energy properties are reproduced.

OPEP for P ( * )P ( * )
For the interaction of heavy P and P * mesons, we use the Lagrangian and matrix elements of (32) and (33). In deriving the potential, we need to be a bit careful about the normalization of the state; there are 2E particles in unit volume. This requires to divide amplitudes by √ 2E per one external leg as was done for N N ¶ . The OPEP for the P ( * )P ( * ) is given by (87) ¶ In the previous publications by some of the present authors and others the factor √ 2 was missing [54][55][56][57][58][59][60]. It is verified also by the former collaborator (S. Yasui, private communications). In this article this problem has been corrected. Accordingly, it turns out that the OPEP plays an important role for e.g. X(3872), while not so for Z c and Z b as discussed in sections 4 and 5. The baryon systems such asDN will be discussed elsewhere.
Here the axial coupling g A is for P ( * ) (or for the light quark), and ε and S are the spin transition operator between P ↔ P * , and spin one operator for P * , respectively. The polarization vector plays the role of spin transition of P ↔ P * andP ↔P * . The tensor In actual studies for X and Z states, the total isospin of a P ( * )P ( * ) system must be specified. The isospinors of these particles arē and the τ matrices in (84)-(87) are understood to operate these isospin states. When P ( * ) is replaced by P ( * ) , an extra minus sign appears at each vertex reflecting the chargeconjugation or G-parity of the pion as shown in table 3. Finally we note that in (84) and (85), the mass is replaced by an effective one µ by taking into account the energy transfer as discussed in the next subsection.

Effective long-range interaction
In the derivation of P ( * )P ( * ) potential, (86) and (87) we have assumed the static approximation, where the energy transfer p 0 is neglected for the exchanged pion, 1 However, when the masses of the interacting particles changes such as in (84) and (85), the effective mass of the exchanged pion may change from that in the free space due to finite energy transfer. To see how this occurs let us start with the expression where we have included the energy transfer with For heavy particles, we may approximate The ignored higher order terms are of order or less when the molecule size is of order 1 fm or larger. These values can be neglected as compared with the mass differences of MeV. In the integration over q, the pion energy is E π (q) > m π . Therefore, • If m P * − m P < m π which is the case of B * , B mesons, the integrand of (90) is regular but the exchanged pion mass is effectively reduced as where µ is regarded as an effective mass. Therefore the interaction range is extended. Some consequences of this effective long range interactions are discussed in [61].
• If m P * − m P > m π which is marginally the case of D * , D mesons, the integrand of (90) hits the singularity and generate imaginary part. The integral is still performed, and the resulting r-space potential is given by µ → iµ.
The function e iµr /r represents the outgoing wave for the decaying pion with momentum |k| = µ. The plus sign is determined by the boundary condition implemented by +i .

Physical meaning of the imaginary part
The presence of imaginary part implies the instability of a system. For DD * systems, it corresponds to the decay of DD * → DDπ if this process is allowed kinematically. To show this explicitly, let us first consider the matrix element of the complex OPEP (95) by a bound state r|B = ϕ(r), where the momentum wave function is defined bỹ and is normalized as Decomposing the denominator of the interaction by using we find By using the identity 1 where P stands for the principal value integral, the imaginary part of (100) is written as We can now show explicitly that the imaginary part (102) is related to a part of the decay processes of the quasi-bound state ϕ. For illustrative purpose, we consider the three-body decay of DD * → DDπ of isospin symmetric case as shown in figure. 5. There actual small mass differences in the charged and neutral particles are ignored in the right panel.  The three-body decay is computed by the diagrams in the first line (left hand side) of figure. 6, which are for the decay of the quasi-bound state at rest (P = 0) into D(p),D(p ), π(q). Note that there are two possible processes for a given set of the final state momenta p, p , q, whose amplitudes are added coherently. Denoting the interaction vertex ofD * →Dπ as hε · q, where h = g A /2f π the amplitude is written as Here the factor √ 2M B is for the normalization of the initial state; there are 2M B particles in a unit volume. Squaring this and multiplying the three-body phase space the decay rate is computed by Now let us consider the diagrams of the second line of figure. 6. The left diagram is computed by Here we have used the time ordered perturbation theory and taken into account only the terms that contribute to the decay. Similarly we obtain the amplitude for the right diagram by the replacement ϕ(p ) → ϕ(p) in the numerator. Therefore the sum of the diagrams is Picking up the imaginary part, we find The optical theorem says that by writing the S-matrix by S = 1 − iT , Therefore, considering the normalization of this particle (2M B particles in a unit volume), we find that the imaginary part agrees with the decay width. This is nothing but an explicit check of the optical theorem. We see that the off diagonal integral in (107), ϕ * (p) · · · ϕ(p ), agrees with the the potential matrix element (102) modulo a kinematical factor. The difference appears due to different normalization factors in the state decaying into two particles and that of bound states. The diagonal part ϕ * (p) · · · ϕ(p) corresponds to the imaginary part of the self energy ofD * and is not included in the potential matrix element.

The quark model and the hadronic model
Here we define the meson and the exotic states using a quark model. By doing so, we can combine the quark model and the hadron model smoothly into a quark-hadron hybrid model. Such a model enables us to handle the physics with resonances, long range interactions like OPEP, and rather complicated systems, by a hadron model but with the quark degrees of freedom effectively included. Also, by constructing hadrons from the quark degrees of freedom, the charge conjugation of the hadron systems can be defined in a more consistent way as follows.
To obtain observables from the model Lagrangian, one has to choose the initial and/or the final state. A q αqβ meson with a certain spin structure can be defined by using the fermion bilinear as [62] where ψ α is the field operator, Γ stands for the sixteen 4×4 matrices,ñ stands for the spin orientation of the state, and the markñ * corresponds to the complex conjugate ofñ. For a vector meson,ñ * · (ψΓψ) corresponds to 1 , and for a pseudoscalar meson, it corresponds to ( −i √ 2 )(ψiγ 5 ψ). The suffix α or β stands for other quantum numbers, such as color and flavor. φ is a relative motion wave function of the quark and the antiquark. M is the meson mass while m α and m β are the quark and the antiquark masses. The normalization of this state is taken as q α q β |q αqβ = 2M δ α α δ β β (2π) 3 δ 3 (K −K), where K and K are the center of mass momenta of the initial and final qq mesons, which we set to be zero in the following. The above expression can be reduced into where r = r q − rq. Note that in this definition the a s † operator stays always on the left side of the b t † operator, not vice versa, for the qq meson. The charge conjugation, C, changes the creation operators of the quarks to those of the antiquarks: So, the state |q αqβ ;ñΓ is transferred into Here the coordinate r is changed to −r after the charge conjugation because it is defined as r q − rq. The minus sign in the last equation comes from anticommutation of the fermion operators a † and b † . There we also usē For simplicity, let us omit the orbital part of the wave function, φ, and assume = 0. In a nonrelativistic quark model, the higher order term of O(p/m) in the spinors is usually taken care of in operators as the relativistic effects. For further computations here, we follow the convention of [62]. The nonrelativistic spinors are taken as where spin up and down correspond to s = 1 and 2, respectively, for both of quarks and antiquarks. First let us consider the pseudoscalar Qq meson, |Qq; SS z = |Qq; 00 , and its behavior under the charge conjugation. For a pseudoscalar meson we take Γ = iγ 5 and This corresponds to the D meson when Q and q are taken as the charm and the light quarks, respectively. In the last equation, we define |q = a † |0 so that its normalization becomes 1 instead of 2E as in the nonrelativistic quark model. When the charge conjugation is applied to this state, we have This state can be regarded as theD meson, meaning that the charge conjugation of the D meson is theD meson. Or, when both of the Q and q quarks are taken to be the charm quarks, this state corresponds to the η c meson, whose C-parity is +1.
Next we consider the vector meson, |Qq; SS z = |Qq; 11 , and its behavior under the charge conjugation. Now we takeñ to be which corresponds to the D * meson, when Q and q are the charm and nonflavor quarks, respectively. Under the charge conjugation, it becomes This corresponds to −D * meson. Or, when both of the Q and q quarks are taken to be the charm quarks, this state corresponds to the J/ψ meson, whose C-parity is −1. This result and (120) are in accordance with what has been anticipated in (18). Finally, we consider a four-quark system such as D ( * )D( * ) and its charge conjugation. The C-parity can be defined for the neutral systems. The charge conjugation changes D toD, and D * to −D * . Thus the C-parity changes D ( * )D( * ) associated with the orbital relative wave function, ψ L (r), with r = r D − rD, as where L is the orbital angular momentum of the D ( * ) andD ( * ) relative motion, and S is the total spin (= 0, 1, 2). Thus, the C-parity eigenstates of the DD systems are also eigenstates of the parity, Similarly, those with DD * or D * D are For D * D * , the simultaneous eigenstates of the parity and the C-parity relate to the angular momentum L and the total spin S as with In the quark model, the relations concerning the rearrangement between cc-qq and cq-qc are derived in a systematic manner. For this, first we note that there are two color configurations for the qqcc systems: the one where the quark-antiquark pairs qq Table 5. Rearrangement of the qq-cc-type mesons and DD mesons. The definition of [D ( * )D( * ) ] are shown in (130), (135) and (140). [DD * ] ± means [DD * ] + for J P C = 1 +− , and [DD * ] − for 1 ++ . The quark configuration of the η meson is assumed to be One can see in table 5 that the above relation between J P C and the phase in (128)- (140) appears for each C-parity.
Up to now, the D meson corresponds to cq and theD meson to qc. The neutral DD states consist of the isospin 0 and 1 states, ((cūuc) + (cddc))/ In the next section we discuss X(3872), whose J P C = 1 ++ . The rearrangement becomes which, if the state has isospin 0, becomes We will denote above state simply by (DD * − D * D )/ √ 2, or just by DD * in the isospin basis if there is no room for confusion. Or, in the following section on X(3872), when we write D + D * − and D 0D * 0 in the particle basis, they mean 1 respectively. So, by this notation, the isospin eigenstates of I(J P C ) = 0(1 ++ ) are as usual. In many literatures for hadron models, such as [63,64], the charge conjugation of D * is defined asD * , not −D * . This can be realized when theD ( * ) meson is taken to becq, not qc. In such a hadron model, the C-parity ±1 eigenstates are given by (DD * ±D * D )/ √ 2, respectively. This difference in the phase is due only to the definition. An extra factor appears when the observables are calculated, which compensates the above difference.
The observed masses of the X(3872) in the B →D * 0 D 0 K decay mode are (3872.9 +0.6 −0.4
The full width is less than 1.
Since the first observation of the X(3872), it has received much attention because its features are difficult to explain if a simple cc bound state of the quark potential model is assumed [84]. The potential between heavy quark and heavy antiquark is well understood and known that it can be approximately expressed with the Coulomb plus linear terms. This feature is confirmed by the lattice QCD studies [85,86]. The charmonia with J P C = 1 ++ are the χ c1 states. The observed mass of the ground state of the χ c1 (1P ) is 3510.67±0.05 MeV and the quark potential model gives similar mass [87]. The first excited state of the χ c1 is the χ c1 (2P ) and that state cannot be observed so far. The predicted mass of the χ c1 (2P ) in quark potential models is in between 3925 MeV and 3953 MeV [87]. The observed mass of X(3872) is 53 MeV to 81 MeV smaller than these predictions. This is one of the strong ground for the non-cc structure of X(3872).
Typical size of the isospin symmetry breaking is at most a few %. It is interesting to know what is the origin of this strong isospin symmetry breaking. In [97], this problem has been studied by using the chiral unitary model and the effect of the ρ-ω mixing has been discussed in [126]. It was reported that both of the approaches can explain the observed ratio given in (149). In the charmonium-molecule hybrid approach, the difference of the D 0D * 0 and the D + D * − thresholds gives enough amount of the isospin violation to explain the experiments naturally [38,114]. The Friedrichs-model-like scheme can explain the isospin symmetry breaking [127]. Recently, a new isospin one decay mode, X(3872) → π 0 χ c1 has been observed [128]. The production processes have been studied in [69,[129][130][131][132][133][134][135][136][137][138][139][140]. There an unexpectedly large production rates have been observed at large transverse momentum transfer p ⊥ > 10 GeV [141]. These rates are much larger than those of light nuclei such as the deuteron and 3 He, and are about several % of that of ψ(2S). This property is naively explained if X(3872) has some small component structure such as χ c1 . In later subsection we will see that this will be realized in a model of DD * molecular coupled with a cc core.
It is also an important issue that whether the charged partner of X(3872) exists as a measurable peak or not. BABAR has searched such a state in the X(3872) → π − π 0 J/ψ channel and found no signal [2]. The hybrid picture, where the coupling to the cc core is essential to bound neutral X(3872), is consistent with the absence of charged X(3872).

D ( * )D( * ) molecule with OPEP
In this subsection, we demonstrate the analysis of X(3872) as a D ( * )D( * ) molecule with I G (J P C ) = 0 + (1 ++ ). As for the interaction between D ( * ) andD ( * ) mesons, we employ only the OPEP in (84)- (87). In the D ( * )D( * ) coupled channel system, the possible D ( * )D( * ) components with positive charge conjugation are where ( 2S+1 L J ) denotes the total spin S, the orbital angular momentum L and the total angular momentum J [12,57]. We note that the phase convention in (151) is different from the one in the literatures [12,57] as discussed in section 3.7. In this basis, the matrix elements of the OPEP are given by [12,57] V π (r) = where C π = C(r; µ, Λ), T π = T (r; µ, Λ), C π = C(r; m π , Λ), and T π = T (r; m π , Λ), respectively. The functions C(r; m, Λ) and T (r; m, Λ) are defined in (79) and (80). The functions C π and T π with µ emerge because the nonzero energy transfer in the D-D * transition is taken into account in the potential V π PP * -P * P (r) as explained in subsection 3.4. In the charm sector, the mass µ becomes imaginary, i.e. µ 2 = m 2 π −(m D * −m D ) 2 < 0. In this subsection, the hadron masses summarized in table 7 are used, which are the isospin averaged masses. Then, µ 2 = (37.3i) 2 [MeV 2 ] is obtained, and the C π and T π become complex as seen in (95). The explicit expression of the imaginary central and tensor potentials is given in [124]. In this analysis, we consider only the real part of the potential, because the imaginary part is small for small µ. Table 7. Hadron masses used in the numerical calculation. The masses of the pion, pseudoscalar meson P = D, B and vector P * = D * , B * meson are shown in the ud, charm and bottom sectors. The pion mass is given as the averaged mass of π + , π 0 and π − [5]. The P ( * ) mass is the averaged mass of P ( * )0 and P ( * )± [5]. The mass difference ∆M P P * between the masses of P and P * is also shown. The Hamiltonian of the D ( * )D( * ) coupled channel system is where the kinetic term K is given by Here we define for the state of orbital angular momentum , the reduced mass and The D ( * )D( * ) systems are studied by solving the Schrödinger equation with the Hamiltonian H (153). In the potential V π , there are two parameters that are the coupling constant g A and the cutoff parameter Λ. The coupling g A is determined by the D * → Dπ decay as shown in section 2.3. The cutoff Λ is a free parameter, while it can be evaluated by the ratio of the size of hadrons. In [54,55], the cutoff Λ for the heavy meson is determined by the relation Λ/Λ N = r N /r D , with the nucleon cutoff Λ N , and the sizes of the nucleon and D meson r N , and r D , respectively. The nucleon cutoff is determined to reproduce the deuteron properties as discussed in section 3.3, and we use Λ N = 837 MeV. The ratio of the hadron sizes r N /r D = 1.35 is obtained by the quark model in [55]. Thus, Λ = 1.13 GeV is obtained.
To start with, the D ( * )D( * ) system is solved for the standard parameters (g A , Λ) = (0.55, 1.13 GeV). We have found that the OPEP provides an attraction but is not strong enough to generate a bound or resonant state. The resulting scattering length is a = 0.64 fm for the S-wave DD * channel. By changing the parameter set (g A , Λ) by a small amount of value toward more attraction, a bound state is accommodated.
To see better the properties of the interaction, we show parameter regions on the (g A , Λ) plane which allow bound states or not. In figure 7, boundaries of the two regions are plotted for three cases depending on how the system is solved; (i) the full calculations with all coupled-channels of D ( * )D( * ) states included and with energy transfer properly taken into account, (ii) calculations in the full coupled channels but with the energy transfer ignored (static approximation), and (iii) calculations with a truncated coupled channels removing the D * D * states. Those lines indicate the correlation between g A and Λ. If the coupling g A is small, the cutoff Λ should be large to produce the bound state, and vice versa.
The lines for (i) and (ii) are similar, which is a consequence that energy transfer is not very important here. Nevertheless, the dashed line (ii) is slightly on the right side (or upper side) of the solid line (i). When g A = 0.55, Λ = 1.6 GeV on the line (i), while Λ = 1.7 GeV on the line (ii). Hence, introducing the energy transfer produces more attraction due to smaller effective mass or equivalently to longer force range. Even for µ = 0, the result is almost the same as that in the case (i).
The central and tensor forces with various effective pion masses are shown in figure 8. The central force changes drastically when the effective pion mass is changed, although the contribution of the central force is not large. These are because there is an overall factor of effective pion mass m 2 . When m 2 is negative, i.e. m = µ is imaginary, the overall sign of the potential becomes negative. On the other hand, the tensor force does not depend on the effective pion mass strongly as shown in figure 8.
Naively, one would expect that a longer range potential yields more interaction strength, while we do not see it much here. One reason is that the central force has the factor m 2 just as discussed. Another reason is that the tensor force is mostly effective at shorter distances than 1/m, due to the S-D coupling. In momentum space, it is due to the q 2 dependence in the numerator (73) which increases the tensor force for large q 2 .
Turning to figure 7, the line (iii) shows the result without the D * D * channel. This line is far above the lines (i) and (ii), indicating that the attraction is significantly reduced. Since the coupling to D * D * component with the D-wave induces the tensor force as shown in (152), ignoring this component decreases the attraction due to the tensor force significantly. Hence, the full-coupled channel analysis of DD * and D * D * is important when the tensor force of the OPEP is considered.
Finally, the B ( * )B( * ) bound state in the bottom sector is studied. We employ the  same potential as used in the D ( * )D( * ) system (152) because the potential is given as the leading term of the 1/M P expansion and thus the potential form is heavy flavor independent in the heavy quark limit. The cutoff Λ B for the B meson is also evaluated by the hadron size in the similar way to the cutoff Λ D . In [55], the ratio of the hadron size GeV. This value can be the reference point here, while we also vary the cutoff to see the cutoff dependence. The use of different Λ for charm and bottom sectors is to take partly into account 1/(heavy quark mass) corrections due to kinematics, because in the quark model meson size is a function of the reduced mass. In figure 9, the boundary line of the B ( * )B( * ) bound state is shown, where it is compared with the boundary of the D ( * )D( * ) bound state, which is the same as shown in figure 7 (i). The bound region for the B ( * )B( * ) system is larger than that of the D ( * )D( * ) . In the bottom sector, the kinetic term is suppressed by the large B ( * ) meson mass, about 5 GeV, while the D ( * ) meson mass is about 2 GeV. In addition, the small mass difference between B and B * , about 46 MeV, magnifies the mixing rate of the S-D coupled channel due to the tensor force, yielding more attraction. For the parameters (g A , Λ) = (0.55, 1.08 GeV), the bound state is found in the bottom sector, where the binding energy is 6.3 MeV.
Because of the attraction in the bottom sector, the bottom counter partner of X(3872) is also expected to be formed as the B ( * )B( * ) bound state. Verification in experiments is needed.

Admixture of the cc core and the DD * molecule
As discussed in the previous section, the OPEP tensor term induces the D ( * )D( * ) S-Dwave channel mixing, which gives an attraction to the X(3872) system. This attraction is sizable, but seems not large enough to produce a bound state. Another origin of the attraction is discussed in [114], where X(3872) is assumed to be a shallow bound state of the coupled channels of cc, D 0D * 0 and the D + D * − . The cc-DD * coupling occurs in the short range where the light quark pair in the DD * state can annihilate. A nearby cc(1 ++ ) state is χ c1 (2P ), which has not been observed experimentally but was predicted by the quark model [179]. The predicted mass of χ c1 (2P ) is by about 80 MeV above the D 0D * 0 threshold energy according to the quark model. So, the coupling to the cc state pushes the energy of the DD * state downward, toward the threshold. As a result, the coupling provides an attraction for the isospin-0 S-wave DD * components of X(3872). Therefore, in this section, we study X(3872) in a coupled channel model of cc, D ( * )0D( * )0 and D ( * )+ D ( * )− .
To start with, we investigate a simple model of such, a model of coupled channels of DD * and cc where the interaction of DD * takes place only through their coupling to cc channel. We call this cc model, where, in the absence of OPEP, only the S-waves are relevant for DD * channels. It is reported that by assuming a coupling between cc and D 0D * 0 and D + D * − , a shallow bound state appears below the D 0D * 0 threshold; but there is no peak structure found at the D + D * − threshold. The coupling structure is assumed as The coupling strength g cc is taken so as to produce the observed mass of X(3872). The cutoff Λ q is roughly corresponds to the inverse of the size of the region where the qq annihilation occurs, being U ∝ e −Λqr /r. Here we show the results with Λ q = 0.5 GeV ∼ (0.4 fm) −1 [114]. When one uses smaller value for Λ q , e.g. 0.3 GeV, the model gives a sizable enhancement around the mass of χ c1 (2P ), 3950 MeV, in the B → DD * spectrum. Since such a structure is not observed in the B-decay experiments, it can be a constraint to the interaction region from the B-decay experiments that Λ q is more than about 0.5 GeV [114]. The mass of the charged D meson is heavier than the neutral one by 4.822±0.015 MeV and that of D * by 3.41±0.07 MeV [5]. Therefore, the threshold difference between D 0D * 0 and D + D * − is about 8.2 MeV. Since the X(3872) mass is almost on the D 0D * 0 threshold, the major component of the X(3872) is considered to be D 0D * 0 . In such a situation, it is convenient to look into X(3872) in the particle basis rather than in the isospin basis. The wave functions of the S-wave DD * components of the X(3872) obtained by using the cc model are plotted by the long dashed curves in figure 10. In the cc model only the S-waves are relevant. The D 0D * 0 wave function is actually large in size and has a very long tail. Its root mean square distance (rms) is listed in table 8. Note that this number varies rapidly as the binding energy varies because the rms becomes infinite as the binding energy goes zero as seen from (82). The rms of the D + D * − component is much smaller than that of the D 0D * 0 because of the D 0D * 0 -D + D * − threshold difference. As seen from figure 10, the amplitudes of the D 0D * 0 and the D + D * − wave functions are similar in size in the very short range region where the DD * state couples to the cc; the isospin-0 DD * state becomes a dominant component there as shown in figure 11. Probabilities of various components of the bound state are shown in the first line of table 9. As was mentioned in section 4.1, the production rate of X(3872) in the pp collision experiments suggests that the amount of the cc component is expected to be several %. In the present cc model, the admixture is 8.6%. As we will show later, by introducing OPEP between the D andD mesons, this admixture reduces to 5.9 %, which corresponds to the amount just required from the experiments. Table 8. Parameters and the root mean square distance (rms) of the two mesons in X(3872) in the cc model (cc) [114], in the OPEP model (OPEP) [180], and in the cc-OPEP model (cc-OPEP) [180]. The rms 0 (rms ± ) means the rms between D ( * )0 and D ( * )0 (D ( * )+ and D ( * )− ). BE is the binding energy in MeV. (82) stands for the values calculated by using (82).  Table 9. Probabilities of each components of X(3872) in the cc model (cc) [114], in the OPEP model (OPEP) [180], and in the cc-OPEP model (cc-OPEP) [180]. model Now we consider models with OPEP included; the one denoted as OPEP in table 9 is the model with only the D ( * )D( * ) channels included as discussed in section 4.2, and the other one denoted as cc-OPEP is the cc-D ( * )D( * ) coupled channel model with the OPEP and their S-D tensor couplings included [180]. The model space is now taken to be cc, D ( * )0D( * )0 and D ( * )+ D ( * )− found in (151):   where c cc is the amplitude of the cc component, c 0 ( 3 S) is that of the D 0D * 0 ( 3 S 1 ) component, c ± ( 3 S) is that of the D + D * − ( 3 S 1 ) component, and so on.
The OPEP potential among the I(J P C ) = 0(1 ++ ) D ( * )D( * ) states are found in (152). In the particle base calculation, it is convenient to use the expression with the explicitly written isospin factor, where C π , T π , C π , and T π are the same as those defined for (152). The cc-DD * coupling is taken as (158). The parameters are listed in table 8. The OPEP cutoff Λ in the OPEPmodel is taken to be a free parameter to reproduce a bound state with the binding energy, 0.16 MeV. As for the cc-OPEP model, the OPEP cutoff Λ = Λ D is the standard one obtained from the D-meson size as marked in figure 7 in the previous subsection. The cc-DD coupling strength, g cc , in the cc-OPEP model is taken to be a free parameter to fit the binding energy. In table 8, rms of the D 0D * 0 and D + D * − system are listed. The size of X(3872), governed mostly by the binding energy, does not depend much on details of the model. The wave functions of each model are plotted in figure 10 and figure 12. The DD * 3 S 1 wave functions are similar to each other, though they are slightly enhanced at the short distance in the OPEP model. This is due to the tensor force; the S-D coupling causes effectively an attraction in the S-wave channel which contains the square of the D-wave amplitude. In fact, the location of the maximum strength of the D-wave amplitude roughly coincides with where the 3 S 1 wave function is enhanced. In the cc model, the attraction comes from the cc-DD * coupling, while the S-Dwave mixing by the OPEP tensor term provides the attraction in the OPEP model Their effects can be roughly estimated by the amounts of the cc components, |c cc | 2 , and the D-state probabilities, which are listed in table 9. In the cc-OPEP model, where both of these attractions are introduced, |c cc | 2 reduces from 8.6 % to 5.9 %, while the D-state probabilities reduces from 2.0 % to 0.6 %. The former reduces to 2/3, and the latter reduces to 1/3, which are the rough share of the attraction in the cc-DD * coupling model with a reasonable cutoff for the OPEP. The D-state probability and the cc probability depend much on the binding energy, or on slight change of g A , whose value is determined in the heavy quark limit. The size of the cc component can also vary as shown in the next subsection. Therefore, it is difficult to estimate the relative importance of OPEP quantitatively. Qualitatively, however, we can conclude that effects of the cc-DD * coupling and OPEP are comparable in X(3872). One has to consider the coupled channel system of cc, the D 0D * 0 and the D + D * − channels with their mass difference, and the cc-DD * and the OPEP coupling, simultaneously, to understand the feature of X(3872).

The decay spectrum of X(3872)
The observed strong decay modes of X(3872) are J/ψπ + π − , J/ψπ + π − π 0 , and D 0D * 0 . Here we discuss the strong decay of X(3872), especially the following two notable features to understand the X(3872) nature. One is that a large isospin symmetry breaking is found in the final decay fractions: as seen in (149), the decay fractions of X(3872) going into J/ψπ 3 and J/ψπ 2 indicate that amounts of the J/ψω and J/ψρ 0 components in X(3872) are comparable to each other as shown in (150). The other feature we would like to discuss here is that the decay width of the X(3872) is very small for a resonance above the open charm threshold, or for a resonance decaying through the ρ and ω components, which themselves have a large decay width.
In the followings we employ a model which consists of the cc core, D 0D * 0 , D + D * − , J/ψω and J/ψρ 0 . For the discussions of decay properties here, it is sufficient to consider the formation of a loosely bound DD * states, which couple to the cc and to the Jψρ and Jψω with finite decay widths for ρ and ω. We assume effective couplings between cc and DD * , which gives the attraction as we discussed in the previous section, and between DD * and J/ψω(ρ 0 ), which expresses the rearrangements. In this section we do not introduce OPEP; the system is restricted only to the S-waves, and the attraction from the OPEP is effectively taken into account by introducing the central attraction between the D ( * ) andD ( * ) . The widths of the ω and ρ mesons are taken into account as an imaginary part in the J/ψω and J/ψρ propagators. In this way, we consider that the model can simulate essential features of the decay properties of X(3872).
From the quark model point of view, the DD * states of the total charge 0 are the ccuū or ccdd states, which contain also the J/ψω or J/ψρ 0 state with the appropriate color configuration. The observed final J/ψπ 3 and J/ψπ 2 decay modes are considered to come from these components. The rearrangement between DD * and ωJ/ψ or ρ 0 J/ψ occurs at the short distance, where all four quarks exist in the hadron size region. Note that there is no direct channel coupling between the cc channel and the ωJ/ψ or ρ 0 J/ψ channels in the present model setup. They break the OZI rule, and the latter breaks the isospin symmetry.
The model Hamiltonian for the cc, D 0D * 0 , D + D − * , J/ψω and J/ψρ channels is taken as [38]: where m cc is the cc bare mass when the coupling to DD * is switched off. The reduced masses, µ DD0,± and µ J/ψω,ρ , are for the D 0D * 0 , D ± D * ∓ , J/ψω and J/ψρ systems, respectively. The coupling between the cc state and the DD * state is expressed by the transfer potential, U (mc) , which is chosen to be Lorentzian in the momentum space with the strength g cc . The rearrangement between the DD * states and the J/ψω and J/ψρ meson is expressed by a separable potential, v in V (mm) . The base of the matrix expression in (164) and in (167)  observed yet. The strengths g cc and v 0 are free parameters under the condition that the mass of X(3872) can be reproduced. The value of Λ q is the same as the one used in the previous section, Λ q = 0.5 GeV. The parameters are summarized in table 10. The amount of each component in the bound state, X(3872), is also listed in table 10. The bulk feature is similar to the models in the previous section: the dominant component is D 0D * 0 while the D + D * − component is considerably smaller because of the threshold difference. The amount of the cc component is somewhat smaller but still sizable. The J/ψω and J/ψρ components are small comparing to the DD * components. The fact that the J/ψρ and J/ψω components of X(3872) are comparable in size is reproduced in the present model.
As illustrated in figure 13, the mass spectrum of X(3872) from the B meson weak decay is proportional to the sum of the transfer strength from the cc to the two-meson states, f , W (cc → f ; E), which can be expressed as Here G (cc) is the full propagator of the cc state, which can be written by using the self energy Σ cc as: Here we define the free and the full propagators within the two-meson space, G (mm) 0 (E) and G (mm) (E), respectively, with the decay widths as where Γ ω,ρ is the ω or ρ decay width, respectively. The ρ meson width is taken to be energy dependent as discussed in [38]. The widths of D * mesons are neglected. The width Γ is ignored when the bound state energy or the component is calculated above. It, however, is essential to include them when one investigates the decay spectrum.
In order to obtain the decay spectrum of each final two-meson channel separately, we have rewritten the right hand side of (170) as follows. Since the system has only one cc state in the present model, the above G (cc) (E) and Σ cc (E) are single channel functions of the energy E. They become matrices when more than one cc states are introduced, but the following procedure can be extended in a straightforward way. As seen from (171), the imaginary part of G (cc) (E) comes only from the imaginary part of Σ cc (E). Therefore, Im G (cc) = Im (G (cc) * (G (cc) * ) −1 G (cc) ) = Im (G (cc) * Σ cc G (cc) ) = Im (G (cc) * U (mc) † G (mm) U (mc) G (cc) ) (176) where * stands for the complex conjugate. Using the following relation for a real potential and Lippmann Schwinger equation for the propagator, G = (1 + GV )G 0 , we have for ImG (mm) on the right hand side of (176) Thus, (176) can be rewritten as When we apply the plain wave expansion for the G where k f and µ f stand for the three momentum and the reduced mass of the final twomeson state where E = m 1f + m 2f + k 2 f /(2µ f ). Γ f is the the decay width of mesons in the final state f , i.e. 0 if f is DD * , Γ ω or Γ ρ when f is J/ψω or J/ψρ. |f k stands for the plain wave of the channel f with the momentum k. |f k (f k) stands for the distorted wave function of the channel f with the momentum k which is generated from |f k . This can be obtained by the Lippmann Schwinger equation as In the present model, only the DD * channels couple directly to cc. The summation over f in (180) means summation over D 0D * 0 and D + D − * . The final two-meson fraction expressed by f in the above equations can be J/ψω or J/ψρ as well as DD * . For the channels where Γ f is small, the transfer strength becomes The calculated W spectrum is shown for each final state in figure 14. The narrow peak at around the D 0D * 0 threshold as well as the large isospin mixing are successfully reproduced by this picture of the cc, D 0D * 0 , D + D − * , J/ψω and J/ψρ. Let us discuss the mechanism to have the small width of the peak and the isospin mixing in the final fractions, the two notable features of the final J/ψω(ρ) spectrum mentioned before. The small width of X(3872) means that the corresponding pole of G (cc) is close to the real axis. The imaginary part in the denominator of G (cc) comes from the imaginary part of Σ cc as shown (171), which can be expanded by using (172) as where G (DD * ) stands for the full propagator obtained within the DD * space. On the right hand side of (183), the first term is real for the energy below the threshold because there is no bound state without the cc core and because the width of the D * meson is neglected. The second term has an imaginary part which comes from the width of the ω and ρ in the propagator G (J/ψω,ρ) 0 . This term corresponds to the rearrangement, DD * → J/ψω, J/ψρ, which is very much suppressed that enables us to expand Σ cc like (183). One suppression comes from the color factor. Probability to find a color-singlet cc pair in the DD * channel has a factor 1/9 in the color space, which is taken into account by taking a small coupling constant v 0 in the transfer potential v. Another suppression comes from the range of the rearrangement. It occurs with charm quark exchange, which means that the reaction occurs in the hadron size region of the very large object, X(3872). For the final J/ψρ channel, further suppression factor appears from the isospin symmetry breaking.
The suppression factor from isospin symmetry breaking can be evaluated as follows. In the second term of (183), G (DD * ) is the factor which breaks the isospin symmetry. So, the DD * wave function generated from |cc , G (DD * ) u|cc has an isovector component. The transfer potential v, which expresses the rearrangement, is a short range interaction whose momentum dependence is L Λq (q) as in (169). One can evaluate the degrees of the isospin breaking in the relevant region by calculating the overlap between the factor L Λq (q) and the wave function G (DD * ) u|cc . The latter roughly corresponds to the DD * component of the bound state, X(3872), as far as the short range part is concerned. So, we estimate the isospin breaking by taking the ratio of the overlap of the I = 1 wave function and the Lorentzian L Λq (q), over that of I = 0, as which is the suppression factor due to the isospin breaking. Thus, together with the color factor, the originally large width of the light vector mesons reduces gives only a small imaginary part to G (cc) around the D 0D * 0 threshold and hence the sharp peak. In order to have the ratio of the final J/ψω fraction over the J/ψρ fraction, it is necessary to evaluate the G Rough estimate of the amount of the J/ψρ fraction against the J/ψω fraction is 0.024/0.15 ∼ 0.16, which is underestimate, but gives the correct tendency of that calculated from the experimental one, 0.27±0.02, shown in (150) [124]. In figure 15, the J/ψω and J/ψρ peaks are plotted for the cc coupling with four different values of the strength g cc . The corresponding binding energies are shown by the arrows, which, as mentioned before, are calculated without the ρ or ω meson widths. The left two peaks correspond to the ones with stronger coupling than that of the reference [38]. The third one is that of the original value with a binding energy of 0.16 MeV. The most right one corresponds to the one where the coupling g 2 cc is reduced by a factor 0.95; there is no bound state but a virtual state appears. As the coupling g 2 cc is weakened, the peak moves as the bound state moves to the threshold, but stops at the threshold when the pole moves to the second Riemann sheet. Namely, the peak energy of the final J/ψπ n is kept lower than or equal to the D 0D * 0 threshold energy. On the other hand, the peak energy of the final D 0D * 0 fraction is higher than threshold by definition. This means that the X(3872) mass is higher when measured by the final DD * fraction, which is consistent with the observation. Experimentally the X(3872) mass observed by the J/ψ and anything is 3871.69±0.17 MeV with a width of <1.2 MeV [5], while that of the final DD * is 3872.9 +0.6 −0.4

+0.4
−0.5 MeV [74] or 3875.1 +0.7 −0.5 ± 0.5 MeV [75]. The observed threshold mass is 3871.68±0.07 MeV [5]. Thus the mass observed by the J/ψ and anything is consistent with the heavier two peaks in figure 15. This suggests that the X(3872) is either a virtual state or a very shallow bound state. These two states are very hard to distinguish from each other experimentally when one considers the width of the component, such as the ρ meson. In literatures, line shapes were studied by using amplitudes parametrized by effective range method [145] or by Flatte parametrization [182][183][184]. It, however, seems difficult to determine the position of the resonance poles just from the shape of the decay spectrum. To discuss this subject, it will probably be necessary to treat the ρ meson as a resonance of the ππ continuum and perform dynamical analyses. In order to discuss possibilities of the other mechanism for the X(3872) peak as well as the effects of the meson decay width, let us ignore the OZI rule just for now. In strong decay, where the isospin is conserved, J/ψρ does not directly decay from the cc core, but J/ψω may. We estimate the effects of existence of such a process by introducing a direct coupling between J/ψω-cc. Suppose the J/ψω-cc coupling occurs by the same potential, u, there is a bound state with a binding energy of 10 MeV below the D 0D * 0 threshold. But as seen in figure 16 and as expected, almost all the decay fraction is J/ψω and not J/ψρ in that case, which is excluded from the experiments. As far as X(3872) is concerned, no direct coupling seems to occur between the cc core and the J/ψω channel. One interesting point of this trial calculation is that the peak energy approaches to the threshold with less pronounced peak as the width of the ω meson is enlarged by hand. Therefore, exotic hadrons which appears at around a two-meson threshold and which contains a meson with a large decay width should be examined carefully. Let us make one more comment on the direct decay from the cc core. The χ c1 (2P ) peak may not be seen in the DD * decay spectrum, but it may be seen from the radiative decay. There is a selection rule of E1 transition that reduces the fraction χ c1 (2P ) → J/ψγ but not → ψ(2S)γ, which may show clearly the existence of χ c1 (2P ) [171,185], and if so, support the cc admixture of X(3872) discussed in this section.
Let us summarize the features of X(3872) and the B-decay to the final two-meson spectrum obtained in section 4. The two-meson and the cc hadron model which consists of the cc(1 ++ ) core, D ( * )0D( * )0 ( 3 S, 3 D, 5 D), D ( * )± D ( * )∓ ( 3 S, 3 D, 5 D), J/ψω and J/ψρ 0 , the following features are obtained: • X(3872) is a very shallow bound state or a virtual state which is close to the D 0D * 0 threshold, which are very difficult to distinguish each other.
• The state is molecular, mostly D 0D * 0 in the long range region, but has a considerable D ±D * ∓ component at the short distance.
• Two kinds of channel couplings provides attraction for X(3872): one is OPEP tensor, which mixes D ( * )D( * ) S-D-waves, the other is the cc-DD * coupling. These two effects are comparable in size.
• The amount of the cc component is found to be about 6 % in the model which contains both of the OPEP and the cc-couplings, which meets the requirement from the production rate in the pp experiments.
• When considering the whole energy spectrum of the B weak decay, there is one very narrow peak at the D 0D * 0 threshold, but not around the D ±D * ∓ threshold, nor around the cc(1 ++ ) bare mass.
• Among the final products, amounts of DD * 's is the largest, which are produced directly from cc core. There are small amount of J/ψω and J/ψρ final product, which are comparable to each other.
• The spectrum of the J/ψω and the J/ψρ final products makes a very narrow peak at the bound state energy, if a bound state exists, or at the D 0D * 0 threshold, if not.
Z c (4430) was first observed in Belle with the mass of 4433 ± 4 ± 2 MeV [197]. The present world average of the mass is about 50 MeV higher than the original value, namely, 4478 +15 −18 MeV [5] though the name Z c (4430) is continuously used. One of the interesting point is that Z c (4430) is observed in the mode including the first radially excited state of the charmonium ψ(2S), not the ground state charmonium J/ψ. The Dalitz analysis of B → Kπψ(2S) decays was performed in [198] and the full amplitude analysis of B 0 → K + π − ψ(2S) decays was done in [199]. As for the quantum numbers of the Z c (4430), J P = 1 + were favored in [199] and confirmed by LHCb [200]. LHCb also performed Argand diagram analysis and showed its resonance character. The decay patterns of Z c (4430) exhibit interesting features: from Ref. [199], and from Ref. [194]. From the above two results, we obtain the ratio of the branching ratios which indicate that the decay to πψ(2S) is much enhanced than the decay to πJ/ψ. This ordering is against the naive intuition; the decay rate to πψ(2S) should be smaller than the decay to πJ/ψ in terms of the differences of the phase space. Further investigations will be useful for understanding the internal structures of Z c (4430).

Isovector P ( * )P ( * ) molecule with OPEP
In this subsection, we study how the Z c and Z b states are generated as isovector P ( * )P ( * ) molecular state with the OPEP. We focus on the states with J P C = 0 ++ , 1 +− and 1 ++ , where the S-wave P ( * )P ( * ) component is included. Among them, J P C = 1 +− is assigned as the quantum number of Z c (3900), Z c (4200), and Z c (4430). The J P C = 1 ++ state has not been reported, but it is the isospin partner of X(3872).
The components of the isovector P ( * )P ( * ) states for J P C = 0 ++ , 1 ++ , and 1 +− are given by [12,57] The lowest threshold of the J P C = 0 ++ state is PP , while PP * and P * P are the lowest thresholds for J P C = 1 ++ and 1 +− . We note that the phase convention in (189)- (191) is different from the one in the literatures [12,57] as discussed in section 3.7.
As studied in section 4.2, we search the parameter region which gives a bound state by varying the parameters g A and Λ. As a result, the boundary of the isovector D ( * )D( * ) and B ( * )B( * ) bound states in the (g A , Λ) plane is shown in figure 17. The results for J P C = 0 ++ and 1 +− are similar, while we note that the lowest thresholds are different, PP for J P C = 0 ++ , and PP * (P * P ) for J P C = 1 +− . The bound region of the J P C = 1 ++ is slightly larger than the others, and hence the attraction for J P C = 1 ++ is larger than that for J P C = 0 ++ and 1 +− . Comparing the results of the D ( * )D( * ) and B ( * )B( * ) states, the B ( * )B( * ) bound region is wider than the D ( * )D( * ) one. This is because the heavier mass suppresses the kinetic energy, and because the the small BB * mass splitting enhances the attraction from the coupled channel effect. For the reference point (g A , Λ) = (0.55, 1.13 GeV), no bound state is found for the isovector channels of both the charm and bottom sectors. To accommodate bound states, we need larger g A and/or Λ. In fact, our previous choice of the overestimated coupling strength corresponds to the vertical line g A ∼ 0.83 = √ 2 × 0.59 [57] (see the footnote ¶ in page 20) which allowed a shallow bound state for the isovector B ( * )B( * ) channel. When we have only the OPEP, larger g A or Λ is needed to produce a isovector P ( * )P * bound state.
Finally, we compare the results for the isovector and isoscalar channels. In figure 18, the boundary lines of the D ( * )D( * ) bound states for I(J P C ) = I(1 ++ ) (I = 0, 1) are shown, which were obtained in figure 17 for I = 1 and in figure 7 for I = 0. As seen in figure 18, the bound region for I = 0 is obviously larger than that for I = 1. This fact indicates that the attraction in the I = 0 channel is stronger than that in the I = 1 channel. The difference between them comes from the isospin dependence of the OPEP, which is given by the isospin factor τ 1 · τ 2 in (84)-(87); τ 1 · τ 2 = −3 for the isoscalar channel, while τ 1 · τ 2 = +1 for the isovector channel. For the OPEP, the tensor force in the off-diagonal term has the dominant role to produce an attraction rather than the diagonal term. For the off-diagonal term, the sign of the potential is not important, but the magnitude is important because the off-diagonal term contributes as the second order of the perturbation. Thus, the attraction in the isoscalar channel with |τ 1 ·τ 2 | = 3 is larger than that in the isovector one with |τ 1 · τ 2 | = 1 by about a factor 9. In the previous subsection, we have seen that the OPEP contribution is rather small in the isovector channel. In such a situation, the short and middle range interactions may become important. Such interaction includes σ, ρ and ω exchange interactions [301].
In the isoscalar P ( * )P ( * ) channel for X(3872), we have not considered them because of the reason discussed shortly below. In the isovector channel, there is another reason that we expect that the vector meson exchanges are not important; ρ and ω exchange interactions have opposite signs and are mostly canceled. The σ exchange potential may be effective not only for the isovector channel but also for isoscalar channel for X(3872). One of the reasons that we have considered only OPEP in the previous sections is that the effect of the short range interaction including the σ exchange has been effectively taken care of by the suitable choice of the cutoff parameter Λ. To determine the reference value of Λ, we have used the binding energy of the deuteron. The OPEP thus determined for the nucleon-nucleon interaction is extrapolated to that of DD * by assuming hadron structures by constituent quarks and the pion coupling to the light quarks. Strictly, however, we do not know to what extent such extrapolation works. Therefore, in this subsection we consider the role of σ-exchange interaction in some detail. Analysis here also provides an estimate of ambiguities coming from short range interactions. The σ exchange potential is derived by the effective Lagrangian of the heavy and σ mesons [302], From this Lagrangian, the σ exchange potential is obtained as V σ P * P * -P * P * (r) = − where the σ mass m σ = 550 MeV is used. The factor 1/m 2 σ is multiplied because the function C(r; m σ , Λ) includes m 2 σ (see (79)). There is ambiguity in the value of the coupling constant g σ . In [301], g σ = 3.65 is taken, which is determined by the quark model estimation. This value is one-third of the value of the σN N coupling g σN N according to the quark number counting, because the σ meson couples to the scalar charge of hadrons which is additive. Another way to estimate g σ is to use a chiral theory for quarks such as the Nambu-Jona-Lasinio model [25,26,303]. By using the equality of the σqq and πqq couplings, and the Goldberger-Treiman relation for the constituent quark, we have the relation By using the parameter values m q ∼ 300 MeV, f π = 93 MeV, g q A ∼ 0.55, we find g σ ∼ 1.8. Yet, in [302], an even smaller value g σ = g π /2 √ 6 ∼ 0.76 is obtained, where g π = ∆M/f π and ∆M is the mass difference between the 0 + and 0 − heavy-light mesons. ∆M = 349 MeV is used in [302], which is the mass difference between D * + s0 and D + s . In this subsection, we present the results for g σ = 0.76 and 3.65, which are regarded as the lower and upper limits of the attractive contribution due to the sigma meson exchange potential.
Using Λ N = 681 MeV for the πσ potential for the nucleon as shown in table 4, Λ D and Λ B are obtained by 919 MeV and 878 MeV, respectively. In the basis of (189)- (190), the matrix elements of the σ exchange potential for the given J P C are obtained by with the function C σ = C(r; m σ , Λ). In figure 19, the boundaries of the isovector D ( * )D( * ) bound states are shown for the case only with OPEP (V π ) is used and for the case of the πσ exchange potential (V σ ) used in the (g A , Λ) plane. The result with g σ = 0 is corresponds to the one only with the OPEP as shown in figure 17. Since V σ is attractive, the bound region for g σ = 0 is larger than that for g σ = 0. For the small coupling g σ = 0.76, the boundary is close to the one for g σ = 0. Thus the V σ contribution is small, and the OPEP plays the dominant role. For g σ = 3.65, however, the bound region is much larger than that for g σ = 0 and 0.76.
The results for the isovector B ( * )B( * ) state are summarized in figure 20 for J P C = 0 ++ , 1 ++ and 1 +− . As seen in the D ( * )D( * ) state, the bound region of the B ( * )B( * ) state becomes large as the coupling g σ increases. The result for g σ = 0.76 is the similar to the one for g σ = 0. For g σ = 3.65, the bound region is much larger than those for g σ = 0 and 0. 76  In the end, we show g σ -Λ plots to see continuously the change in the role of σ exchange as g σ is varied figure 21. For the charm sector, the left figure indicates an unlikely situation for the molecular states to be generated at the mean value g σ 1.8, where a very large cutoff is needed, Λ ∼ 4 GeV. For bottom sector, the molecular states are not yet generated there, but only slight increase of g σ will do.

Brief summary
In this article we have discussed the hadronic molecule as one of exotic structures of hadrons. It has become possible experimentally to observe various exotic phenomena long after the predictions made about half century ago, which have stimulated the diverse    amount of theoretical works. Ingredients of hadronic molecules are constituent hadrons and their interactions. The constituent hadrons also couple to compact structures. Therefore, we have discussed in detail how the admixture model has been applied to X(3872). From the first principle point of view, such a picture should effectively and conveniently replaces the direct but complicated approach of QCD. In other words, the model should be economized [33,34,304], with its work region, where and how, being under control. As emphasized in introduction, hadronic molecules are expected to appear near threshold regions. Their formation is a consequence of finely tuned hadron dynamics, as their binding or resonant energies are of order MeV which is much smaller than the scale of low energy QCD, Λ QCD ∼ some hundreds MeV. We have discussed that such conditions are better realized for heavy hadrons. Their kinetic energies are suppressed and relatively weak force is sufficient to generate hadronic molecules.
For the interaction, we have emphasized the role of the one pion exchange interaction. The pion dynamics is well established because the pion is the Nambu-Goldstone boson of spontaneous breaking of chiral symmetry, where the pion interaction is dictated by the low energy theorems. In the constituent quark picture the pion interacts with the light u, d quarks by the pseudoscalar Yukawa coupling of σ · q type, whose strength is extracted from the empirically known axial coupling constants of hadrons.
The long range part of hadron interactions is provided by the one pion exchange potential (OPEP). This is what we have discussed in detail in this paper. Because of the spin structure of the Yukawa coupling, the OPEP contributes to the transitions D → D * and D * → D * . It has turned out that they are effective in the formation of a DD * molecule for the X(3872) channel . Therefore, an emphasis has been put on the role of the tensor component of the OPEP which causes mixing of states (channels) with different angular momenta by 2 . Because this transition leads to second order process, the resulting interaction for the relevant low lying channel must be attractive, and more attraction with more channel couplings.
There still remains a question of short range interaction. Ambiguities are in such as sigma (correlated two pion), ρ and ω meson exchanges, and even quark exchange processes. In a chiral symmetry description, a part of such interactions is described by the Weinberg-Tomozawa interaction, application of which to heavy systems is however dubious. Alternatively an interaction has been employed as derived by the vector meson exchanges [305,306], one way of which to formulate is the hidden local symmetry approach [45,46]. For loosely bound or resonant hadrons the less known short range interaction can be effectively parametrized in form factors. The cutoff parameter there is then determined by the known experimental data. In accordance with this strategy, the OPEP with a form factor was shown to work nicely for the low energy properties of the deuteron [53], which is the best established hadronic molecule of the nucleons. Interestingly, the employed cutoff Λ ∼ 900 MeV determined by the deuteron properties is consistent with the size of the nucleon core √ 6/Λ ∼ 0.5 fm. Discussions of the nucleon core have a long history. It is recognized as the repulsive core of the N N interaction, which has been explicitly shown in the study by the quark cluster model [307,308]. It was also discussed in the chiral bag model where the nucleon is expressed as the quark core and pion clouds [44,309].
Having this construction of the interactions, the OPEP provides non-negligible amount of attraction particularly for hidden heavy hadrons such as X(3872) with isospin 0, where the tensor force plays the dominant role, while the central component plays little. Therefore, the inclusion of the channel coupling of SD waves (generally, states with different angular momenta by 2 ) is crucially important. The resulting strength of the attraction, however, turns out not sufficient to generate X(3872) as a molecular state of DD * . Because of this the coupling/mixing with a compact state of cc has been introduced. Such a mixing provides an additional attraction if the compact state has a larger mass as theoretically expected for χ c1 as a charmonium cc. The mixing is also required to explain the large production rates of X(3872) in high energy hadron processes. Quantitative estimates of the production rate, however, are to be done carefully [129-131, 141, 220]. In the present analysis of X(3872), the OPEP and the short-range coupling play roughly in equal manner. In general, however, their relative importance depends on the systems to be studied.
Another possible molecule which we have discussed is Z c (3900). However, the strength of the OPEP for Z c (3900) of isospin one is smaller than for X(3872) of isospin zero by factor three. As a result, the attraction is reduced and the formation of molecular state is less likely. In the remaining part of this article, we have also discussed the above features for the bottom sector.

Pentaquark baryon P c
In this subsection we shortly discuss one important subject which we have not mentioned in this article; hadronic molecules for baryons. The oldest candidate is Λ(1405) which is described as theKN molecule [11]. In the low energy theorem, the interaction is given by the Weinberg-Tomozawa interaction forK. Physically, much strength of attraction is provided betweenK and N by the ω meson exchange. At the quark level, it originates in the interaction between the antiquarkq inK and the quarks q in N . Due to charge conjugation, the sign of the vector type interaction flips from the one between quarks. Such a picture has been shown in the Skyrme model which is one of successful chiral models for baryons [310][311][312][313]. As a result, theKN bound state appears near theKN threshold. By coupling to the lower πΣ channel, the bound state turns into a resonance, which is a typical mechanism of the Feshbach resonance. An experiment for the study of theKN molecule is going at J-PARC and analysis is underway [314][315][316].
Yet an important topic is the hidden charm P c baryons observed by LHCb in the weak decay of Λ b → J/ψpK. In 2019 they reported three narrow peaks near the thresholds [19] with higher statistics than the first analysis in 2015 [317]; one below Σ cD and two below Σ cD * thresholds. This observation has immediately lead to theory discussions of heavy quark multiplets formed by the four combinations of Σ c , Σ * c and D,D * [318][319][320][321]. In the heavy quark limit, the pair of Σ c and Σ * c is considered as spin doublet of J = 1/2 and 3/2 just as the pair ofD andD * is regarded as the one of J = 0 and 1. Their hadronic molecules also form multiplets of heavy quark spin symmetry. In connection with the present discussions, these multiplets provide the direct evidence of the tensor interaction of OPEP [322].
To discuss the unique feature, let us take a look at the two nearby states P c (4440) and P c (4457), which are slightly below the Σ cD * threshold. Because the LHCb observed these states in the J/ψp final state, the isospin of these P c states is I = 1/2. Therefore their masses are, respectively, 23 MeV and 6 MeV below the isospin 1/2 threshold of Σ cD * , 4463 MeV. A relevant question is the origin of the masses, decay widths and quantum numbers of these states. Assuming that the orbital motion of the molecules is dominated by the S-wave, possible total spins are J = 1/2 and 3/2. Now the crucial observation is that the tensor interaction can be effective for both states because the S-D couplings survive for both channels. This is understood by the fact that the sum of S(Σ c ) = 1/2, S(D * ) = 1, L = 0 and the sum of S(Σ c ) = 1/2, S(D * ) = 1, L = 2 can both form the total J = 1/2 and 3/2. This sharply contrasts with the two nucleon system. In the latter, the S-D mixing due to the tensor force is effective only for the spin J = 1 state (corresponding to the deuteron), while not for the spin J = 0 state. Therefore, the two P c 's provide a good opportunity to study the role of the tensor force in the OPEP other than the nucleons in a unique way. Otherwise, they also give an information of spin-spin interactions.
Having said this much, in [322] the role of the tensor force in the OPEP has been discussed in a coupled channel model of Σ ( * ) cD ( * ) , Λ cD ( * ) with OPEP supplemented by the short range interaction which is brought by the coupling of the molecular states with compact five quark states [323]. An interesting observation there is that by adjusting only one model parameter for the short range interaction, they have made predictions of ten states. Three of them correspond to the P c states of the LHCb with good agreement of masses and decay widths. Then their identification of the quantum numbers of P c (4440) and P c (4457) is J = 3/2 and 1/2, respectively. The mechanism of lowering the J = 3/2 state is due to the tensor force. This result is in contrast with other predictions, where the mass splittings are attributed to a spin-spin interaction. In the latter, the spin 1/2 state is naturally expected to appear lower than the 3/2 state because of the spin alignment, antiparallel spins (1/2 ⊕ 1 = 1/2) give attraction while parallel repulsion (1/2 ⊕ 1 = 3/2). The spins of the P c states are not yet known, and hence the determination of them is extremely important to make further understanding of these states.

Resonances or cusps
Here, we briefly mention a question which one would yet like to ask; whether the observed exotic phenomena imply physical resonant states or cusps of virtual states. At this moment, we have not had decisive answer to this question for so far observed signals, while there are many articles discussing the nature of the signals theoretically. Here we just refer only to a few of them in relation with X(3872) [145,[182][183][184]324]. There, amplitude analyses are performed by using parametrizations of Flattè or effective range expansion types. Then an observation was made that by suitably choosing parameters, the line shape X(3872) was shown to emerge as a virtual state cusp at the threshold [145].
To reproduce the very narrow or sharp peak at the D 0D * 0 threshold, however, there must be sufficient amount of attraction between them. If the attraction is larger than the critical value, the peak appears as a resonant state, otherwise a virtual state cusp. The difference between them is subtle because only a small change in the interaction strength may change the nature of the peak. Moreover, in such a situation it is difficult to differentiate them experimentally. But then the important question is; what would be the mechanism to provide that suitable amount of the attraction? In this paper, we have tried to offer an option that a model with the pion exchange interaction does it, supplemented by a coupling to a short distance structure. This is a dynamical approach for the construction of amplitudes that we discuss shortly below.

Hadron interactions and exotics
The last issue that we would like to mention is the dynamical approach for the construction of amplitudes from reliable hadron interactions. For heavy hadrons including charm or bottom quarks, it is formidably difficult to derive interactions from experiments. This is the reason that we have resorted to a model for the study of X(3872) in this paper.
Yet another powerful and promising method is the lattice QCD, which is in principle the first principle method for the strong interaction. In the so called HAL QCD method, hadron interactions are obtained by using the Nambu-Bethe-Salpeter amplitude [325,326]. To obtain hadron interactions, this method is practically more powerful than the widely used Luscher's method [327,328]. An attempt was made for Z c (3900) with coupled channels of DD * , η c ρ, J/ψπ, where they have derived the interactions between these channels and solved the coupled channel problem [259,260]. Unexpected is the finding that there is a rather strong coupling between J/ψπ and DD * channels, which effectively causes attraction for the J/ψπ channel. As a consequence, they have found rather than a resonance but a virtual state pole which contributes to an enhancement near the DD * threshold corresponding to Z c (3900).
For the study of exotic hadrons, an approach based on the coupled channel method with suitable hadron interactions is highly desired. It is a non-trivial program because many channels may couple including those of more than two particles. With complementary approaches of experiments, effective theories and lattice simulations, such an approach can be more elaborated, which enables to elucidate the nature of exotic hadrons.