Dynamical Approach to Decays of XYZ States

: We review the existing results on the exotic XYZ states and their decays obtained within the conﬁned covariant quark model. This dynamical approach is based on a non-local Lagrangian of hadrons with quarks, has built-in quark conﬁnement, and is suited well for the description of different multiquark states, including the four quark ones. We focus our analysis on the various decay modes of ﬁve exotic states, X ( 3872 ) , Z c ( 3900 ) , Y ( 4260 ) , Z b ( 10610 ) , and Z (cid:48) b ( 10650 ) , aiming to clarify their internal quark structures. By considering mostly branching fractions and decay widths using the molecular-type or the tetraquark-type interpolating currents, conclusions about the nature of these particles are drawn: the molecular structure is favored for Z c ( 3900 ) , Z b ( 10610 ) , and Z (cid:48) b ( 10650 ) and the tetraquark for X ( 3872 ) and Y ( 4260 )


Introduction
The concept of multiquark states composed of more then three quarks hypothesized decades ago [1] was for the first time confirmed in 2003 where multiquark state candidates were measured by the BES [2], BaBar [3], and Belle [4] experiments. The latter observation, seen in the π + π − J/ψ invariant mass spectrum, was the first observation of a charmonium-like state X(3872), which did not fit expectations of existing quark models for any conventional hadronic particle. The reason was mainly its measured mass 3872 MeV, not predicted by models, and also the difficulty in interpreting it as an excited charmonium ψ : its eventual decay into ρJ/ψ is strongly suppressed because of isospin violation. In the following years, other heavy quarkonium-like states X, Y, Z were discovered, where Y usually denotes electrically neutral exotic (i.e., non-cc) charmonia having quantum numbers J PC = 1 −− , Z is used for charged states, and X labels any non-Y and non-Z cases. With the aim to report on the results and achievements of the confined covariant quark model, we narrow our review of experimental outcomes to a relevant subset of the whole exotic meson family.
The study of the Y(4260) decay channel J/ψπ + π − by BESII [17] and Belle [25] in 2013 led to the discovery of the charged Z c (3900) resonance in the invariant mass distribution of J/ψπ ± . The Z ± c particle was in the same year observed also by the CLEO-c detector [26]. In addition, the latter experiment provided the first evidence of the neutral member of the Z c isotriplet, the Z 0 c state, discovered in the π 0 J/ψ channel. A state Z c (3885) was seen in the DD * spectrum of the e + e − → π ± (DD * ) ∓ reaction at BESIII in 2014 [27]. Assuming it can be identified with the Z c (3900) particle, the measurement provided arguments in favor of J P = 1 + quantum numbers. The same experiment reaffirmed in 2015 the existence of the neutral Z c state [28], in 2017 confirmed with high significance the J P = 1 + assignment [29], and in 2019 provided the evidence for the ρ ± η c decay channel [30]. The D0 collaboration published the observation of the Z c (3900) state in pp collision data in 2018 [31] and studied its mass and width in [32] (2019). The current Z c parameters are [13] m Z c (3900) = 3887.2 ± 2.3 MeV, Γ Z c (3900) = 28.2 ± 2.6 MeV and I G (J PC ) = 1 + (1 +− ). Z c (3900) as a charmonium-like state with an electric charge is a prominent candidate for an exotic multiquark state and is largely discussed in the existing literature.
Two narrow bottomonium-like four quark state candidates were detected in the Belle detector [33] in 2012. They were labeled Z b (10610) and Z b (10650) and were observed as peaks in the mass spectra of π ± Υ(ns), (n = 1, 2, 3) and π ± h b (ms), (m = 1, 2). The same experiment published two other papers dedicated to these exotics. In [34], the evidence was given for the quantum number assignment I G (J P ) = 1 + (1 + ) for both of the states. In [35], they were observed in different decay channels Z b (10610) → [BB * + cc] ± and Z b (10650) → [B * B * + cc] ± , where one can notice the proximity of the two states to the corresponding B ( * )B * thresholds. These decays dominated the studied final states, which besides two bottom mesons, included also a pion and for which the Born cross-section was given. The decay into BB was found to be suppressed with respect to the two previous final states, and an upper limit was given. The masses and widths are m Z ± b (10610) = 10,607.2 ±2.0 MeV, m Z b (10650) = 10,652.2 ±1.5 MeV, Γ Z ± b (10610) = 18.4 ± 2.4 MeV, and Γ Z b (10650) = 11.5 ± 2.2 MeV. Growing evidence suggests that the mentioned and also other, unmentioned exotic heavy quarkonium-like states observed since 2003 cannot be described as simple hadrons in the usual quark model. The effort to understand their nature combined with the non-applicability of the perturbative approach in the low energy domain of quantum chromodynamics (QCD) resulted in a large number of more or less model dependent strategies. In existing reviews [36][37][38][39][40][41][42][43][44][45][46][47][48][49], different ideas are analyzed. The proximity of the X, Y, Z masses to meson pair thresholds naturally leads to a popular concept of the hadronic molecule, more closely reviewed in different contexts. In [50], the authors studied the implications of the heavy quark flavor symmetry on molecular states. The authors of [51] argued in favor of a molecular picture using an isospin-exchange model, and a nice review of the molecular approach was given in [52]. A frequent treatment of four quark states is represented also by QCD sum rules [53][54][55] and different quark models. A dynamical approach based on a relativistic quark model with a diquark-antidiquark assumption was proposed in [56,57], where tetraquark masses were computed. A non-relativistic screened potential model, presented in [58], was used to compute the masses, electromagnetic decays, and E1 transitions of charmonium states. Treatment of tetraquarks as compact dynamical diquark-antidiquark systems in [59] had the ambition to explain why some of the exotic states preferred to decay into excited charmonia. Several hypotheses (molecular description, tetraquark description, hadro-charmonium picture) for different exotic states were investigated in [60] using tools based on the heavy quark spin symmetry: besides drawing conclusions for some XYZparticles, also possible discovery channels were given. The hybrid and tetraquark interpretation for several exotic states were discussed in paper [61] using the Born-Oppenheimer approximation.
A very complete review of exotic states with some emphasis on the chromomagnetic interaction was provided in a recent publication [62]. The ideas of coupled channels ( [63]) and heavy quark limit ( [64]) are also often seen in the context of the exotic quarkonia. Finally, one has to mention the possibility of peaks in invariant mass distributions being explained by the kinematic effect. This was investigated in detail in a recent text [65]. The arguments for X, Y, and Z states not being purely kinematic effects were given in [66].
In the present paper, we want to review the description of the exotic heavy quarkonia-like states by the confined covariant quark model (CCQM). The model [67][68][69] was proposed and developed as a practical and reliable tool for the theoretical description of exclusive reactions involving the mesons, baryons, and other multiquark states. It was based on a non-local interaction Lagrangian, which introduces a coupling between a hadron and its constituent quarks. The Lagrangian guarantees a full frame independence and the computations relay on standard quantum field theory techniques where matrix elements are given by the set of quark-loop Feynman diagrams according to the 1/N c expansion. Earlier, a confinement was not implemented in the model, and thus, it was not suited for heavy particles (with baryon mass exceeding those of the constituent quarks summed). This was changed in [69], where a smart cutoff was introduced for integration over the space of Schwinger parameters. Since then, arbitrary heavy hadrons could be treated by the CCQM. The CCQM represents a framework where the hadron and the quarks coexist, which raises questions about the proper description of bound states and the double counting. They are solved using a so-called compositeness condition. It guarantees, by setting the hadron renormalization constant Z H to zero, that the dressed state and the bare one have a vanishing overlap. In order to describe radiative decays, one also needs to introduce gauge fields properly in a non-local theory such as CCQM. This was done by the formalism developed in [70] where the path integral of gauge field appeared in the quark-field transformation exponential. One should also mention that the model had no gluons: their dynamics was effectively taken into account by the quark-hadron vertex functions, which depended on one hadron size related parameter. The model has a limited number of free parameters; besides the hadron related ones, it has six "global" parameters: five constituent quark masses and one universal cutoff. The model was applied to with success light and heavy mesons and baryons (e.g., [71][72][73][74][75][76]) and also to exotic four quark states [77][78][79][80][81][82][83]. The latter will be reviewed in the rest of this article.
All sketched features of the CCQM (interaction Lagrangian, confinement, compositeness condition, implementation of electromagnetic interaction) are addressed in more detail in Section 2. Section 3 is dedicated to the X(3872) state and its decays to J/ψ + ρ andD + D * . Its radiative decays are analyzed in Section 4. In Section 5, molecular and tetraquark hypotheses for the nature of Z c (3900) are put in place and the results compared with experimental data. The exotic to exotic reaction Y(4260) → Z c (3900) ± + π ∓ and the decay of Y(4260) to open charm are presented in Section 6. Decays of the bottomonium-like states Z b (10610) and Z b (10650) to several different final states are studied within the molecular picture in Section 7. We close the text by a summary and conclusion given in Section 8.

Interaction Lagrangian
The dynamical description of hadrons in the CCQM follows from the interaction Lagrangian: where the hadronic field is coupled to a non-local quark current. The latter takes different forms for different hadrons: for the mesons, for the baryons, and for the tetraquarks. Here, C stands for the charge conjugation matrix C = γ 0 γ 2 with C = C † = C −1 = −C T and Γ is an appropriate Dirac matrix (or string of Dirac matrices) to describe the spin quantum numbers of the hadron. One has CΓ T C −1 = Γ for the (pseudo)scalar and axial-vector fields and CΓ T C −1 = −Γ for vectors and tensors. The color indices are denoted by superscripts a i , and F H (x; x 1 , . . . , x n ) represents a non-local vertex function, which characterizes the quark distribution inside the hadron. We assume it takes the form: The first factor reflects the natural expectation that the barycenter of the quark system corresponds to the position of the hadron, and the second term has a general form dependent on the relative quark coordinates. Obviously, the vertex function is invariant under translations: for any four-vector a. In principle, any form of the function Φ H is allowed as long as it has an appropriate fall-off behavior in the Euclidean momentum space to guarantee the ultraviolet convergence of the Feynman diagrams. Various alternatives of the vertex function for non-local quark currents were analyzed in [84], and it was found that the dependence of the observables on different choices was small. Because of convenience of performing calculations, the exponential form for the Fourier transform of the function Φ H was adopted: where K 2 is the combination of the loop and external momenta. The minus sign indicates that we are working in the Minkowski space, and the wicked-rotated argument K 2 → −K 2 E makes explicit the appropriate fall-off behavior in the Euclidean region. Λ H is an adjustable parameter of the CCQM, which can be related to the hadron size. Additional free parameters are the constituent quark masses and a universal infrared cutoff (discussed in more detail later). Their values, summarized in Table 1, were determined by adjusting the model predictions to experimental data.

Compositeness Condition
In the Lagrangian of the CCQM, quarks and hadrons are treated equally. However in nature, hadrons are made of quarks. Therefore, questions about an appropriate description of the bound states and double counting arise. The issue is resolved by imposing the so-called compositeness condition [85,86], which requires the renormalization constant of the hadron field to vanish. Since the renormalization constant Z 1/2 H can be interpreted as the matrix element between the physical state and the corresponding bare state, Z 1/2 H = 0 implies that the physical state has no overlap with the bare state and is therefore described as a bound state. For a spin-one particle, the compositeness condition reads: where Π H is the derivative of the scalar part of the vector-meson mass operator: The condition Z 1/2 H = 0 also effectively removes the constituent degrees of freedom from the space of physical states and so eliminates the double counting. A general tetraquark self-energy diagram to be used for the compositeness condition is show in Figure 1.
One should also notice that the application of the compositeness condition lowers the number of model parameters because its fulfillment is reached by tuning the coupling constant value. Equation (4) thus fixes the coupling and increases the predictive power of the CCQM over the wide range of hadronic data. The determination of g H for all participating hadrons by means of the compositeness condition is the first step in the application of the CCQM. It should be remarked that the compositeness condition can be interpreted also in terms of the normalization of the electric form factor at q 2 = 0, as shown in [69].

Infrared Confinement
If the mass of a hadron reaches the limit defined by the sum of the masses of constituent quarks, then in a model without a confinement, the hadron becomes unstable and decays into its constituents. In order to correct this unphysical behavior and enlarge the applicability of the model also to the (increasing) experimental data on heavy hadrons, the confinement of quarks was introduced in [69]. Its implementation assumes the Schwinger representation of quark propagators: with the subsequent cutoff in the upper integration limit applied, in a clever way, to the whole structure of a Feynman diagram. The latter, containing l loops, m vertices, and n propagators, can be schematically written as: where Φ symbolizes vertex functions, p denotes external momenta, v linear combinations of external momenta, k represents loop momenta, and κ linear combinations of the latter. The expression in curly brackets is the argument of the vertex function in the momentum representation. The second line makes explicit the integration over the space of the Schwinger parameters with the whole structure of the first line catch in the F symbol. The next step is to go from the integration over the Schwinger parameters to the integration over the n − 1 simplex of the dimensionless Feynman parameters combined with a one-dimensional integral over a dimension variable t by using the simple insertion of unity in the above expression: One has: Note that a dimension variable t is analogous of the Fock-Schwinger proper time. By performing an analytical continuation of the kinematical variables to the Minkowski space, one can encounter the branch points, which, in particular, correspond to the quark unitary thresholds. The appearance of the imaginary parts in Equation (7) is witnessed on the quark production in the physical spectrum, i.e., on the absence of the quark confinement. One possibility to resolve this problem is to cut the upper limit of the integration over the proper time t, i.e., ∞ → 1/λ 2 with λ being the "infrared" cutoff parameter. This allows one to remove all singularities of the diagram related to the local quark propagators. The integral becomes smooth and always convergent. For clarity, one can demonstrate the approach on a scalar one-loop two-point function. One starts with the loop integral in the Euclidean space: By using the above transformations, one gets: where p 2 = −p 2 E . The expression has a branch point at p 2 = 4m 2 as follows from the vanishing of the first term inside the exponential at α = 1/2. However, this singularity is removed by the cutoff.
We take the value of the cutoff parameter λ presented in Table 1 as universal for all processes we describe.

Electromagnetic Interactions
Inclusion of the electromagnetic (EM) interaction into the non-local CCQM in a gauge invariant way requires a dedicated approach. Our main interest will be in the radiative decays of neutral particles (i.e., the X(3872) tetraquark; see [79]), and so, we will focus on the EM interactions of quarks. The free part of the quark Lagrangian is gauged using the standard minimal coupling prescription: with e u = 2/3, e d = −1/3 in units of the proton charge. This defines the first part of the interaction Lagrangian of quarks with photons: A second term comes from gauging the non-local quark-hadron interaction Lagrangian (1). First, one multiplies each quark field by an exponential expression, which depends on the gauge field: with I being defined as the integral: over the path that connects the hadron and quark positions. It can be readily seen that the local gauge transformations: leave the Lagrangian unchanged for any arbitrary function f . Indeed, the gauge field induced modification of the path integral I( is canceled by the contribution coming from the quark transformations. The exact form of the gauged non-local Lagrangian L EM (2) depends on the quark current structure (i.e., hadron quantum numbers), and we will write down an explicit form of it in the dedicated section. To use the gauged Lagrangian in perturbative calculations, one expands the gauge exponentials into the powers of the coupling constant (and thus, powers of A µ ) up to a desired order. The expansion contains only the derivatives of the path integral I, and using the approach proposed in [70,87], one can define them in a path independent manner: Here, P denotes a path derived from P by extending P from its endpoint by dx. This definition gives: where the independence of the derivative on the path P becomes explicit.

Selected Computational Aspects
To proceed with the calculations, it is convenient to use the following representation for the correlation function: It may be easily obtained by using the Jacobi coordinates. The Gaussian function form of the correlation function Φ X in Equation (3) can be joined with the exponents coming from the Schwinger representation of quark propagators given by Equation (5) into a single exponential function. Its argument takes a Gaussian form in loop momenta: where a is a 3 × 3 matrix, r = (r 1 , r 2 , r 3 ) is a vector constructed from external momenta, and the constant R behaves as a quadratic form of external momenta. As a result, one observes that, with respect to loop momenta, the general expression (6) is a product of a polynomial P (originating from evaluation of traces) with an exponential function. The tensorial loop momenta integration is then performed using the differential identity: which allows us the move the k independent differential operator in front of the integral. The action of the latter operator on the result of the integration is further simplified by applying a second operator identity: which permits commuting the differential operator with the exponential. The next steps are automated using a FORM program [88]. It repeatedly performs the differentiation using the chain rule, thus effectively commuting the differential operator to the right (and eventually making it vanish by acting on a constant). At last, one is left with an integral over the space of the Schwinger parameters (see Section 2.3). The latter is computed numerically with the help of a FORTRAN code. Most of the time, one is interested in the q 2 dependent hadronic form factors: for the purposes of this text, the CCQM should be seen as a smart and effective tool that provides these form factors from the assumed quark currents as inputs.

Strong Decays of X(3872)
The controversy raised by the discovery of the X(3872) state can be best seen in the large number of publications it provoked (with many different interpretations). The proximity of the D * 0 D 0 threshold: naturally suggests the idea of a loosely bound charm meson molecule. This idea was studied in several texts: implications of the molecular hypothesis for interference and binding effects were discussed in [89]; the authors of [90] found support for the molecular interpretation within a non-relativistic quark model; a published text [91] analyzed the molecular assumption in an effective field theory approach; and further works [92][93][94] based their analyses on an effective field theory with pion exchange, Monte Carlo simulations, and heavy quark spin symmetry. A rather strong support for the molecular picture was given in [95] (line shapes study) and [96] (potential model). The lattice study [97] found an explanation for X(3872) in both the molecular and tetraquark scenario. An important group of analyses focused on charmonium [98][99][100] or mixed charmonium [101][102][103][104] explanations. Further arguments in favor of a charmonium structure followed from the Flatté analysis performed in [105], and both molecular and charmonium hypotheses were discussed in [106]. Several works [107][108][109][110] disfavored the molecular description. The authors of [107] based their conclusion on a non-relativistic quark model with the pion exchange, and the analysis presented in [108] favored the charmonium picture instead, while the conclusions in [109] were based on the pion and sigma exchange model. More rare were approaches based on the glueball picture [111] and chromomagnetic interaction [112]. The authors of [113] put in question the existence of a bound state at all. A hybrid hypothesis was considered in [114] and [115] (here, together with the molecular and charmonium one). Lattice computations in relation to X(3872) were used in [116,117], QCD sum rules in [118,119], and the coupled channel approach in [120][121][122]. One should also mention the studies based on quark models [56,123,124] and other strategies [125][126][127].
The description of the X(3872) state by the CCQM was presented in [78]. There, one assumed a tetraquark structure, and within this assumption, decays X → J/ψ + 2π(3π) and X →D 0 + D 0 + π 0 , proceeding through the off-shell ρ/ω and D * states respectively, were computed. In addition, possible implications of the X(3872) dominance in the s-channel dissociation of J/ψ were discussed.
When describing the X(3872) state, one follows the suggestions of [123,128], where a symmetric spin distribution of this J PC = 1 ++ state was proposed: A non-local generalization of this diquark-antidiquark current is written as: with simplified weights resulting from only two quark flavors being present: The strong isospin violation observed by comparing the ρ and ω vector meson mediated decays: experimentally established by Belle [129] suggested a mixed nature of the physical states X l , X h : where θ is the mixing angle. The state X u breaks the isospin symmetry maximally: The mixing angle is to be adjusted to fit the branching fraction ratio (24). The first step in our calculation is to determine the coupling constant g X by using the so-called compositeness condition discussed before. The derivative of the tetraquark mass operator needed for this can be written as: q γ 5 tr S [3] c γ µ S [13] q γ ν + w q tr S [12] c γ 5 S c γ µ S [13] q γ ν + w q tr S [12] c γ 5 S where the short notations for the quark propagators and loop momenta are: The evaluation of this expression is related to the determination of the size parameter Λ X value and allows us to study the Λ X dependence of the results.
Because the X(3872) mass lies close to the studied thresholds: the off-mass-shell character of the ρ, ω, and D * vector mesons has to be taken into account when evaluating the transition amplitudes X → J/ψ + ρ(ω) and X → D * 0D0 . The Feynman diagrams to be considered within the CCQM are depicted in Figure 2. In what follows, we use the notation for the light vector mesons v 0 = ρ, ω. The amplitude of the decay X u →D + D * is written as: where the argument of the X-vertex function is equal to: The amplitude of the decay X u → J/ψ + v 0 is written as: where the argument of the X-vertex function is equal to: In the latter expression, the number of Lorentz structures is reduced to six when X and J/ψ are on the mass-shell because, in that case, one has µ (q µ 1 + q µ 2 ) = 0 and ν q ν 1 = 0. Obvious relations: allow expressing all amplitudes of physical states transitions in terms of the X u ones: The differential decay rate in the narrow-width approximation is written as [130]: where p * (q 2 ) = λ 1/2 (m 2 X , m 2 J/ψ , q 2 )/2m X is the momentum of the v 0 in the X rest frame. The allowed kinematic range is given by: where n = 2 for the ρ meson and n = 3 for the ω meson. The masses, decay widths, and branching fractions appearing in (28) were taken from PDG [13]. In addition to the model parameter values presented in Table 1, further model parameters are needed, namely the size parameters of the appearing mesons. Their values have been settled earlier and are presented in Table 2. Table 2. Size parameters for selected mesons in GeV.
Two adjustable parameters remain, the size parameter Λ X and the mixing angle θ. It was found out that the dependence of the branching fraction: on the size parameter Λ x is in the CCQM small and close to 1/4. Using this observation and the central value of the experimental ratio in Equation (24), one can deduce the mixing angle from: The latter equation yields θ ≈ 18.4 • for X l and θ ≈ −18.4 • for X h . When not considering the ratio, the sensitivity of the decay widths on the size parameter is more important. One may expect the size parameter value to be close to those of the charmonia Λ J/ψ and Λ η c , i.e., to be in the range 3 GeV < Λ X < 4 GeV. This range was scanned, and the behavior of the decay width is depicted in Figure 3. One can conclude that the predicted values in the interval 2.5 ≤ q 2 ≤ 3.5 GeV lie in the range 0.05 MeV < Γ X(3872) < 0.23 MeV, which is in agreement with the upper limit of 1.2 MeV.
The differential rate of the decay X(3872) →D 0 D 0 π 0 in the narrow-width approximation is written as: The matrix element M µν was defined above by Equation (26). One has to note that the allowed kinematic range: is very narrow. Taking the masses, widths, and branching fractions of appearing D * mesons from [13,89,[131][132][133][134], we can calculate the decay width: and study its dependence on the size parameter Λ X . This is shown in Figure 3. By using the experimental data from PDG [13] for the ratio: one finds: The latter is to be compared to the CCQM prediction: where the uncertainty of the result reflects the uncertainty on Λ X . One can see that the two numbers agree within errors.

Implications of X(3872) in the Charm Dissociation Process by Light Mesons
It is interesting to check the significance of X(3872) in the reaction of the charm dissociation process J/ψ + ρ (ω) → X(3872) →DD * , which plays an important role in heavy ion physics. This state will contribute to the s channel of the process. The X-addition to the full cross-section is written as: where p = p 1 + p 2 = q 1 + q 2 . The ∓ sign in the first equation is negative for the ρ meson and positive for ω. A Breit-Wigner propagator is used with Γ X = 1 MeV, and the size parameter value is fixed to Λ X = 3.5 GeV. With this setting, the dependence of the cross-section on the energy E = √ s is shown in Figure 4.  One can compare the predicted behavior to available results for the charged D-mesons: At E = 4.0 GeV, a theoretical evaluation [135] predicts σ(J/ψ + π → D +D * ) = 0.9 mb, and the work in [136] predicted σ(J/ψ + ρ → D +D * ) = 2.9 mb at E = 3.9 GeV. In the case of X(3872), the cross-section reaches the maximum of approximately 0.32 mb at E = 3.88 GeV, and one can conclude that the expected contribution of X(3872) in the charm dissociation is non-negligible.

Radiative Decays of X(3872)
The first experimental evidence for the radiative decay of the X(3872) particle was given in [129] by the Belle experiment. From the measured branching fraction product: the partial width ratio was deduced: This finding was supported by the BaBar observation [137]: which had a limited significance of 3.4 σ. The same experiment reaffirmed the observation in 2009 [138] with smaller errors: from which one can deduce [36]: BaBar also presented a result related to ψ(2s): B(B ± → XK ± ) · B(X → γ + ψ(2S)) = (9.5 ± 2.7 (stat) ± 0.6 (syst)) × 10 −6 .
In 2011, the Belle collaboration published measurements with J/ψ and ψ(2s) in the final state [139]: The first result was in good agreement with the previous one from the same experiment (36); however, the second number brought some tension when compared to BaBar and a later LHCb measurement [140]: The theoretical study of radiative X(3872) decays includes several different approaches. Such decays were analyzed in [98] in the charmonium picture. The authors studied excited 1D and 2P states and their decays in relation with the electric dipole radiation and provided implications for quantum number assignments. The molecular hypothesis was considered in [106]. There, the authors argued that the validity of the molecular picture could be determined from the study of several X(3872) decay channels (including some with the photon emission). The work in [141] was dedicated to radiative decays with two D mesons in the final state. It was claimed that the discrimination between the molecular and charmonium picture could be obtained via analysis of the photon spectrum. Several decay modes, which also included J/ψ + γ, were examined in [142] within a phenomenological Lagrangian approach. The predicted value of the radiative decay width depended on the model parameters and varied from 125 KeV to 250 KeV. In [143], X(3872) was described as a mixture of charmonium and exotic molecular states and treated using QCD sum rules. The predicted radiative decay width ratio Γ X (J/ψγ)/Γ X (J/ψπ + π − ) = 0.19 ± 0.13 was in agreement with experimental measurements. The excited charmonium hypothesis and study of E1 decay widths within the relativistic Salpeter method was presented in [144]. A description based on a charmonium-like picture with high spin 2 −+ using a light front quark model was proposed in [145]. Later works [146][147][148][149] were mostly interested in the puzzling Γ X (ψ(2s)γ)/Γ X (J/ψγ) ratio (43) and analyzed it with different approaches (quark potential model, single-channel approximation, coupled-channel approach, charmonium-molecule hybrid model, and an effective theory framework).
Here, we focus on the J/ψ decay channel, which was studied using the CCQM in [79]. The non-local quark current for the X(3872) hadron was given in the previous section; see Equation (22). The J/ψ quark current is written as: The related size parameter was established in earlier works and has the value of Λ J/ψ = 1.738 GeV.
The knowledge of the quark currents enables us to give more details concerning the interaction with photons, addressed before in Section 2.4. The second part of the electromagnetic interaction Lagrangian stands: where J µ 4q and J µ 2q correspond to the parts of usual currents (22), (44) not containing the vertex function. In order to make use of the definition (15), it is convenient to switch to the Fourier transforms of the vertex functions and quark fields: so that the differential operator can be placed in front of the path integrals: where the long derivatives are defined as D with ω i being combinations of the integration four-vectors p i and mass parameters w q and w c . Next, the identity involving the operator function action on the path integral [150] is applied: Its validity extends to all functions F analytic at zero. The result for X(3872) reads: , l ij = w ij (w ij r + 2 ω j ), (i = 1, . . . , 4; j = 1, . . . , 3), For J/ψ, one obtains: The amplitude evaluation requires evaluation of four Feynman diagrams displayed in Figure 5. The corresponding expression stands: where T µρν (q 1 , q 2 ) can be expanded in terms of appropriate Lorentz structures. Using the on-mass shell condition, gauge invariance, and Schouten identities [151], one can show that only two independent structures remain: The functions W A/B are to be extracted from the expression following from the CCQM computation: where the separate contributions are written down: One evaluates the traces and the loop momenta integrals, and the expression is re-arranged in two terms following the mentioned Lorentz structure. The behavior of coefficient functions W A/B is predicted using a numerical integration over the Schwinger parameters: The decay width is expressed as: where H i denote the helicity amplitudes: with | q 2 | = m 2 X − m 2 J/ψ /(2m X ). The dependence of the predicted decay width on the size parameter Λ X is shown in Figure 6.   If we follow the approach from the previous section and take Λ X = 3.0 ± 0.5 GeV, then the model predicts: which is to be compared with the experimental results Equation (37) and Equation (40). One may conclude that the bound-tetraquark description of the X(3872) state by the CCQM is in an agreement with the experimental observations.

Nature of Z c (3900)
As stated in the Introduction, the detected Z c (3900) decays include both the π ± J/ψ and D * D final states (assuming Z c (3885) and Z c (3900) are the same particle). The ratio of the decay widths of these two channels was measured by BESIII [27]: and represents a quantitative observation to be explained by the theorists. There are many different theoretical approaches that are trying to understand the nature of this state. The tetraquark interpretation was intensively discussed within QCD sum rules [152][153][154] and also in the color flux-tube model [155]. The molecular scenario seems to be more abundant in the literature and is discussed or preferred in several theoretical frameworks. A light front theory description was presented in [156]; an effective field theory description was proposed in [157]; and QCD sum rules were used in [158,159]. The molecular interpretation was also supported by the quark model developed in [160]. The authors of [161] made a proposal for BESIII and forthcoming Belle II measurements by using also the molecular scenario. Further molecular picture oriented works can be found in [162] (constituent quark model, coupled channels) and in [163] (quark interchange model). It is interesting to note that most of the lattice QCD based studies obtained different results from previous ones: some did not see (within the approach they used) a bound state at all [164][165][166][167], invoked a threshold cusp explanation [168,169], or indicated that the understanding of Z c within the lattice QCD was only approaching [170]. For completeness, one can mention the charmonium hybrid interpretation studied in [171], the hadro charmonium picture presented in [172] with the tetraquark and molecular interpretation and the color magnetic interaction [173]. Further ideas can be found in [174][175][176][177][178][179][180][181][182][183][184].
The description of Z c (3900) in the framework of the CCQM was presented in [80]. Two options were tested: the molecular interpretation and the tetraquark hypothesis. For each option, the strong decays into J/ψπ + , η c ρ + ,D 0 D * + , andD * 0 D + were computed and compared to available experimental data. First, we investigate the tetraquark hypothesis. In this scenario, the non-local Z c current is written as: The tetraquark mass operator looks like: where the momenta are defined by: The matrix elements of the decays Z + c → J/ψ + π + and Z + c → η c + ρ + are written down: where the argument of the Z c -vertex function is given by: The notations used are as follows: The amplitudes of the Z + c →D 0 + D * + and Z + c →D * 0 + D + decays are: with the argument of the Z c -vertex function being: Now, the notation used is The decay width for the 1 + (p) → 1 − (q v ) + 0 − (q s ) transition is given by: where H denotes the helicity amplitudes and q v is the three-momentum of the final state vector particle q The helicity amplitudes can be related to the invariant amplitudes A 1 and A 2 , which parametrize the matrix element in terms of the Lorentz structures: by means of the relations: From the comparison of Equation (64) with Equations (58)-(61), one can express A 1,2 as a function of A xy , B xy . The results are importantly influenced by the fact that the amplitudes AD D * and A D * D (Formulas (60) and (61)) vanish exactly within the CCQM description AD D * = A D * D = 0, and the contributions from the non-zero B amplitudes are strongly suppressed by the |q v | 5 factor. Before arriving at the numerical predictions, the size parameters need to be specified, and a strategy with respect to the choice of Λ Z c value has to be settled. The numerical values of the size parameters were in [80] (i.e., the herein presented Z c analysis) re-adjusted with respect to those in [78] and are shown in Table 3. Table 3. The size parameters for selected mesons in GeV used in the Z c (3900) analysis.
As concerns the Λ Z c parameter, first, it is taken as Λ Z c = 2.24 ± 0.10 GeV to make the predicted value of the decay width Γ(Z + c → J/ψ + π + ) close to the one from [152,176]. One obtains: These outputs contradict the experimental number (see Equation (55)), which indicates a larger coupling to DD * than to the J/ψπ mode. If trying to adjust the Λ Z c parameter to a more realistic value, the results do not become any better. Assuming Λ Z c = 3.3 ± 1.1 GeV, one gets: These predictions suggest that the tetraquark picture is not appropriate for the Z c (3900) state.
The molecular description of Z c (3900) appears as a natural alternative. In such a scenario, the non-local interpolation quark current is written as [53]: By using similar steps as in the tetraquark analysis, one writes down the Fourier transformed Z c mass operator in the form: in order to pin down the Λ Zc dependence of the coupling g Z c . Next, the transition amplitudes are constructed: where the argument of the function Φ Z c is given by: The meaning of all other letters and symbols is the same as was in the previous paragraph dedicated to the tetraquark description. The decay widths are also evaluated in a fully analogous way. However, the parameter Λ Z c needs to be adjusted independently. Tuning its value in such a way so as to provide the best description of the BESIII measurement [27], one gets Λ Z c = 3.3 ± 1.1 GeV with the following values for the decay widths: One can see that the obtained results at this time are in agreement with the experimental observations by showing an enhancement of the DD * sector and are in agreement with the observed branching fraction ratio in Equation (55) within the errors. One can conclude that the CCQM supports the molecular picture of the Z c (3900) state.

The Nature of Y(4260)
The distinctive characteristics of the Y(4260) are its mass, which does not fit any charmonium in the same mass region, the suppression of open charm decays with respect to the J/ψπ + π − final state, and the appearance of the exotic charmonium Z c (3900) among its decay products. This interesting mix of properties is addressed in quite a few theoretical works, and like in other cases, the molecular, tetraquark, and several other explanations are invoked.
A support for the molecular picture was provided by the QCD lattice computations in [185], by QCD sum rules in [186], by a meson exchange model in [187], and also by the authors of [188], which favored it over the hadro-charmonium interpretation. Further arguments for Y(4260) being a molecule were based on the line shape study in [189], and the authors of [190] proposed an unconventional state with a large, but not completely dominant molecular component. An interesting paper [191] came up with a baryonic molecule concept, and the molecular hypothesis was also analyzed in [192][193][194][195].
On the contrary, the molecular scenario is strongly disfavored in [196] because of reasons related to the heavy quark spin symmetry and the molecular scenario was rejected in [197] in favor of a charmonium hybrid one. Here, the crux of the argument lies in an important separation between Y(4260) mass and its decay threshold. Further arguments to support the charmonium or hybrid-charmonium picture were given in the publications [198][199][200].
One should also mention different quark models [201][202][203][204] with some of them favoring the tetraquark description of Y(4260). The tetraquark hypothesis was also analyzed in the QCD sum rules study [205], and the coupled channels approach combined with the three-particle Faddeev equations was used to describe Y(4260) in [206].
The analysis of Y(4260) is within the CCQM [207] done in a similar way to the Z c case: its decay modes are analyzed in both the molecular and tetraquark scenario. With quantitative measurements related to Y(4260) not being very numerous, one can analyze the partial decay widths to J/ψπ + π − and open charm final states and see whether the latter ones are suppressed. The Feynman diagrams describing the studied transitions are drawn in Figure 7. The considered open charm final states include DD, DD * , D * D , and D * D * . As follows from the previous section, Z c (3900) is described as a molecular state (67). The molecular-type non-local interpolating current for Y(4260) is written as: with: The matrix element corresponding to the open charm production is given by: where: and the momenta are defined as: The decay into Z c + π involves a three-loop diagram, and the corresponding matrix element is: where: and the momenta are defined as: In the tetraquark scenario, the non-local Y(4260) current takes the form: The matrix element of the decay into DD is expressed as: with the momenta: The matrix element of the decay into Z c π is given by: with the momenta: Here, the summation over Γ is defined by: The considered decays comprise different combinations of pseudoscalar, vector, and axial-vector particles in the final state. The relevant expressions for the matrix elements and decay widths are written down: The value of Λ Z c is set to 3.3 GeV, and guided by our experience, we assume that Λ Y(4260) = 3.3 ± 0.1 GeV. The numerical evaluation leads to the results presented in Table 4.
In both scenarios, the open charm decays are suppressed with respect to the J/ψπ decay channel. The discrimination between them is provided by the total decay width Γ[Y(4260)] = 55 ± 19 MeV, which is in contradiction with the molecular description. Thus, one can conclude that the CCQM approach favors the tetraquark structure of Y(4260).

Bottomonium-Like States Z b (10610) and Z b (10650)
Exotic quarkonia states appear also in the bottomonium sector: Z b (10610) and Z b (10650) are two examples. Even though the exotic bottomonia masses tend to be significantly higher than the charmonia ones, the underlying dynamics is similar, and one finds the molecular, tetraquark, and other hypotheses in theoretical approaches that describe them.
Z b (10610) and Z b (10650) were seen as molecules in the boson exchange model of [208], and the molecular picture was also favored in [209], where the spin structure of these two particles was analyzed. Further support of the molecular scenario came from the quark model based on a phenomenological Lagrangian used by the authors of [210] and also from other analyses preformed in [211] (QCD multipole expansion), [212] (effective field theory), [213] (pion exchange model), [214] (QCD sum rules, only Z b (10650) included), [215] (heavy quark spin symmetry and coupled channels analysis),and [216] (coupled channels approach with pion exchange model). A different set of works supports, with various intensity, the tetraquark structure of the two bottomonia states. In [217], the conclusion followed from an effective diquark-antidiquark Hamiltonian combined with meson-loop induced effects. The authors of [218] based their analysis on the QCD sum rules and interpreted Z b and Z b as axial-vector tetraquarks. The two works [219,220] also drew their conclusions from the QCD sum rules and allowed the tetraquark and molecular scenario. The former work suggested that Z b and Z b could have both the diquark-antidiquark and molecular components (following from a mixed interpolating current). The latter one excluded neither the tetraquark nor molecular the interpretation of Z b (10610), and the idea of a mixed current appeared also. The mentioned analyses could be supplemented by numerous other works  where further ideas and approaches were exploited.
The theoretical analysis of the Z b (10610) and Z b (10650) states by the CCQM was performed in [81]. The work assumed a molecular-type interpolating current, which is favored by most theoretical approaches when interpreting the experimental results. It is a natural choice reflecting the proximity of the particle masses to the corresponding thresholds: The quantum numbers of the two states I G (J PC ) = 1 + (1 +− ) lead to the choice of (local) interpolating currents: which guarantees that, when considering the transitions into B ( * )B( * ) , the Z b state can decay only to the [B * B + c.c.] pair, while the Z b state can decay only to aB * B * pair. Decays into the BB channels are not allowed. Further decay channels include a bottomonium particle accompanied with a charged light meson. Taking into account the G parity, which is conserved in strong interactions and kinematic considerations, only three possible bottomonium-meson decay channels are available: All mentioned Z + b transition can be arranged into three groups with respect to the spin kinematics: The classification of the bottomonia particles based on their quantum numbers is shown in Table 5. Table 5. The bottomonium states 2S+1 L J . We use the notation The expressions for matrix elements and decay widths depend on the spin structure and are for the three cases as follows.
• For 1 + → 1 − + 0 − transitions, the matrix element can be parameterized with two Lorentz structures: The invariant amplitudes A and B can be combined into the helicity amplitudes: which are practical to express the decay width. For the derivation of the latter, it is useful to work in the rest frame of the initial particle, where |q 1 | = λ 1/2 (M 2 , M 2 1 , M 2 2 )/2M is the three-momentum and E 1 = (M 2 + M 2 1 − M 2 2 )/2M is the energy of the final state vector. Furthermore, the on-mass-shell character of the initial and final state particles is taken into account by p 2 = M 2 , q 2 1 = M 2 1 , q 2 2 = M 2 2 , and p µ ε µ = 0. One arrives at: • The matrix element for the 1 + → 1 + + 0 − transitions is expressed through one covariant term only: The decay with the formula can be written as: where one can note the p-wave suppression factor |q 1 | 3 .
With all the above theoretical expressions, one can proceed to the numerical evaluation of the decay widths. The first step is the adjustment of the size parameters Λ Z and Λ Z . They are tuned so as to respect the observables measured by the Belle collaboration [35]: leading to: With the decays into B pairs dominating all other decay channels, we approximate the total decay width as the sum of all herein evaluated channels. The CCQM gives: which is in fair agreement with (108). The predicted partial decay widths of Z b (10610) and Z b (10650) particles are summarized in Table 6. The Z b and Z b decays are dominated [13] by Γ Z b (B + B * 0 +B * + B 0 ) = (85.6 +2.1 −2.9 )% and Γ Z b (B * + B * 0 ) = (74 +4 −6 )%, respectively, meaning that the bottomonia modes should not exceed 15 and 25 percent. This is observed for the h b (1P)π + final state; the other bottomonia channels are suppressed, but not so much as seen in the data: The model also allows us to make predictions: One can conclude that the CCQM provides, within a molecular picture, a fair description of Z b (10610)/Z b (10650) states and related decay observables and catches the tendencies seen in experimental data. Some deviations are observed when the fraction of bottomonium in final states is considered.

Summary and Conclusions
The confined covariant quark model is an approach based on a non-local interaction Lagrangian of quarks and hadrons. It has many appealing features: a full Lorentz invariance, confinement, large applicability range (from mesons to exotic hadrons), inclusion of the electromagnetic interaction, and a limited number of free parameters. As a practical tool, it allows overcoming the difficulties related to the non-applicability of the perturbative approach for bound states in QCD. In this text, we used it to describe four quark exotic states X(3872), Z c (3900), Y(4260), Z b (10610), and Z b (10650). We demonstrated that the CCQM had enough predictive power to make the distinctions between various hypothesis, with respect to the exotic quarkonia mostly related to their structure (molecular versus tetraquark one). At the same time, the model provides a good description of experimental data without large deviations and predictions for future measurements. Concerning the structure of the studied particles, the molecular picture is favored for Z c (3900), Z b (10610), and Z b (10650) and the tetraquark one for X(3872) and Y(4260). These conclusions follow from the measured decay characteristics of the considered exotic states and the related model description: with the expected increase in the number and quality of experimental data, one may hope the quarkonia-structure puzzle will be solved in the years to come.