**1. 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 *π*+*π*<sup>−</sup> *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 *JPC* = 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 first observation of the *X*(3872) mentioned in the previous paragraph was later confirmed in the *pp*¯ collisions by the CDF [5] and D0 [6] experiments in 2004, by the LHCb experiment [7] in 2011, and also by the BESIIIcollaboration [8] in 2014. Further experimental investigations [9–12] increased the mass measurement precision, established the quantum numbers, and put limits on several decay related observables. As of now [13], *X*(3872) is a particle with the mass *mX*(3872) = 3871.69 ± 0.17 MeV, width Γ*X*(3872) < 1.2 MeV, and quantum numbers *IG*(*JPC*) = 0+(1++), mostly decaying to *<sup>D</sup>*<sup>∗</sup> <sup>0</sup> (<sup>→</sup> *<sup>D</sup>*<sup>0</sup> *π*0)*D*0.

Charmonium-like state *Y*(4260) was for the first time observed by the BaBar experiment [14] in 2005 in the *J*/*ψπ*+*π*<sup>−</sup> mass distribution. Its existence was further confirmed by the CLEO [15] (2006), Belle [16] (2007), and BESII [17] (2013) collaborations. Later investigations by BaBar [18] and BESIII [19] provided further updates on the mass and width parameters. With mass above the *DD*¯ threshold, the *Y*(4260) was also searched for in the open charm decay channels, however with negative results [20–24]. The *Y*(4260) is [13] an *IG*(*JPC*) = 0−(1−−) state with the mass and width *mY*(4260) = 4230 ± 8 MeV, Γ*Y*(4260) = 55 ± 19 MeV.

The study of the *Y*(4260) decay channel *J*/*ψπ*+*π*<sup>−</sup> by BESII [17] and Belle [25] in 2013 led to the discovery of the charged *Zc*(3900) resonance in the invariant mass distribution of *J*/*ψπ*±. The *Z*± *<sup>c</sup>* 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 *Zc* isotriplet, the *Z*<sup>0</sup> *c* state, discovered in the *π*<sup>0</sup> *J*/*ψ* channel. A state *Zc*(3885) was seen in the *DD*¯ <sup>∗</sup> spectrum of the *<sup>e</sup>*+*e*<sup>−</sup> <sup>→</sup> *<sup>π</sup>*±(*DD*¯ <sup>∗</sup>)<sup>∓</sup> reaction at BESIII in 2014 [27]. Assuming it can be identified with the *Zc*(3900) particle, the measurement provided arguments in favor of *J<sup>P</sup>* = 1<sup>+</sup> quantum numbers. The same experiment reaffirmed in 2015 the existence of the neutral *Zc* state [28], in 2017 confirmed with high significance the *J<sup>P</sup>* = 1<sup>+</sup> assignment [29], and in 2019 provided the evidence for the *ρ*±*η<sup>c</sup>* decay channel [30]. The D0 collaboration published the observation of the *Zc*(3900) state in *pp*¯ collision data in 2018 [31] and studied its mass and width in [32] (2019). The current *Zc* parameters are [13] *mZc* (3900) <sup>=</sup> 3887.2 <sup>±</sup> 2.3 MeV, <sup>Γ</sup>*Zc* (3900) <sup>=</sup> 28.2 <sup>±</sup> 2.6 MeV and *<sup>I</sup>G*(*JPC*) = <sup>1</sup>+(1+−). *Zc*(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 *Zb*(10610) and *Z*- *<sup>b</sup>*(10650) and were observed as peaks in the mass spectra of *π*±Υ(*ns*), (*n* = 1, 2, 3) and *π*±*hb*(*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 *IG*(*JP*) = 1+(1+) for both of the states. In [35], they were observed in different decay channels *Zb*(10610) <sup>→</sup> [*BB*¯ <sup>∗</sup> <sup>+</sup> *cc*] ± and *Z*- *<sup>b</sup>*(10650) <sup>→</sup> [*B*∗*B*¯ <sup>∗</sup> <sup>+</sup> *cc*] ±, where one can notice the proximity of the two states to the corresponding *B*(∗)*B*¯ <sup>∗</sup> 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 *mZ*<sup>±</sup> *<sup>b</sup>* (10610) = 10,607.2 <sup>±</sup>2.0 MeV, *mZ*- *<sup>b</sup>*(10650) = 10,652.2 <sup>±</sup>1.5 MeV, <sup>Γ</sup>*Z*<sup>±</sup> *<sup>b</sup>* (10610) <sup>=</sup> 18.4 <sup>±</sup> 2.4 MeV, and <sup>Γ</sup>*Z*- *<sup>b</sup>*(10650) <sup>=</sup> 11.5 <sup>±</sup> 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–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–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–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/*Nc* 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 *ZH* 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–76]) and also to exotic four quark states [77–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*/*ψ* + *ρ* and *D*¯ + *D*∗. Its radiative decays are analyzed in Section 4. In Section 5, molecular and tetraquark hypotheses for the nature of *Zc*(3900) are put in place and the results compared with experimental data. The exotic to exotic reaction *Y*(4260) → *Zc*(3900)<sup>±</sup> + *π*<sup>∓</sup> and the decay of *Y*(4260) to open charm are presented in Section 6. Decays of the bottomonium-like states *Zb*(10610) and *Z*- *<sup>b</sup>*(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.

#### **2. Confined Covariant Quark Model**

#### *2.1. Interaction Lagrangian*

The dynamical description of hadrons in the CCQM follows from the interaction Lagrangian:

$$\mathcal{L}\_{\text{int}} = \lg\_H \cdot H(\mathbf{x}) \cdot \mathbf{J}\_H(\mathbf{x}), \tag{1}$$

where the hadronic field is coupled to a non-local quark current. The latter takes different forms for different hadrons: 

$$J\_M(\mathbf{x}) = \int d\mathbf{x}\_1 \int d\mathbf{x}\_2 \, F\_M(\mathbf{x}; \mathbf{x}\_1, \mathbf{x}\_2) \cdot \eta\_1^a(\mathbf{x}\_1) \, \Gamma\_M \, q\_2^a(\mathbf{x}\_2)$$

*Symmetry* **2020**, *12*, 884

for the mesons,

$$J\_B(\mathbf{x}) = \int d\mathbf{x}\_1 \int d\mathbf{x}\_2 \int d\mathbf{x}\_3 \, F\_B(\mathbf{x}; \mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) \cdot \Gamma\_1 \, q\_1^{a\_1}(\mathbf{x}\_1) \left( q\_2^{a\_2}(\mathbf{x}\_2) \mathbb{C} \, \Gamma\_2 \, q\_3^{a\_3}(\mathbf{x}\_3) \right) \cdot \varepsilon^{a\_1 a\_2 a\_3}$$

for the baryons, and

*J μ <sup>T</sup>*(*x*) = - *dx*<sup>1</sup> ... - *dx*<sup>4</sup> *FT*(*x*; *x*1,..., *x*4) · *q a*1 <sup>1</sup> (*x*1) *<sup>C</sup>*Γ<sup>1</sup> *<sup>q</sup>a*<sup>2</sup> <sup>2</sup> (*x*2) · *q*¯ *a*3 <sup>3</sup> (*x*3) Γ2*C q*¯ *a*4 <sup>4</sup> (*x*4) · *ε a*1*a*2*c ε a*3*a*4*c*

for the tetraquarks. Here, *C* stands for the charge conjugation matrix *C* = *γ*0*γ*<sup>2</sup> with *C* = *C*† = *<sup>C</sup>*−<sup>1</sup> <sup>=</sup> <sup>−</sup>*C<sup>T</sup>* and <sup>Γ</sup> is an appropriate Dirac matrix (or string of Dirac matrices) to describe the spin quantum numbers of the hadron. One has *C*Γ*TC*−<sup>1</sup> = Γ for the (pseudo)scalar and axial-vector fields and *<sup>C</sup>*Γ*TC*−<sup>1</sup> <sup>=</sup> <sup>−</sup><sup>Γ</sup> for vectors and tensors. The color indices are denoted by superscripts *ai*, and *FH*(*x*; *x*1,..., *xn*) represents a non-local vertex function, which characterizes the quark distribution inside the hadron. We assume it takes the form:

$$F\_H(\mathbf{x}; \mathbf{x}\_1, \dots, \mathbf{x}\_n) = \delta^{(4)}\left(\mathbf{x} - \sum\_{i=1}^n w\_i \mathbf{x}\_i\right) \Phi\_H\left(\sum\_{i$$

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:

$$F\_H(\mathbf{x} + a; \mathbf{x}\_1 + a\_\prime \dots \mathbf{x}\_n + a) = F\_H(\mathbf{x}; \mathbf{x}\_1, \dots, \mathbf{x}\_n).$$

for any four-vector *a*. In principle, any form of the function Φ*<sup>H</sup>* 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 Φ*<sup>H</sup>* was adopted:

$$\check{\Phi}\_H(-\mathcal{K}^2) = \exp\left(\frac{\mathcal{K}^2}{\Lambda\_H^2}\right) \tag{3}$$

where *K*<sup>2</sup> 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 *<sup>K</sup>*<sup>2</sup> → −*K*<sup>2</sup> *<sup>E</sup>* makes explicit the appropriate fall-off behavior in the Euclidean region. Λ*<sup>H</sup>* 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.

**Table 1.** Constituent quark masses and universal cutoff *λ* in GeV.


#### *2.2. 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 *<sup>H</sup>* can be interpreted as the matrix element between the physical state and the corresponding bare state, *Z*1/2 *<sup>H</sup>* = 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:

$$Z\_H = 1 - \mathcal{g}\_H^2 \Pi\_H'(m\_H^2) = 0,\tag{4}$$

where Π- *<sup>H</sup>* is the derivative of the scalar part of the vector-meson mass operator:

$$\Pi\_{H}^{\mu\nu}(p) = \mathcal{G}^{\mu\nu}\Pi\_{H}(p^{2}) + p^{\mu}p^{\nu}\Pi\_{H}^{(1)}(p^{2}),$$

$$\Pi\_{H}(p^{2}) = -\frac{1}{3}\left(\mathcal{g}\_{\mu\nu} - \frac{p\_{\mu}p\_{\nu}}{p^{2}}\right)\Pi\_{H}^{\mu\nu}(p).$$

The condition *Z*1/2 *<sup>H</sup>* = 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 *gH* 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*<sup>2</sup> = 0, as shown in [69].

**Figure 1.** General confined covariant quark model (CCQM) tetraquark self-energy diagram.

#### *2.3. 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:

$$S(k) = \frac{(m+k)}{m^2 - k^2} = (m+k)\int\_0^\infty d\beta \exp\{-\beta(m^2 - k^2)\}\tag{5}$$

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:

$$\begin{split} \Pi(p\_1, \ldots, p\_m) &= \int (d^4 k)^l \prod\_{i=1}^m \Phi\_{i+n} \left\{-\sum\_j (\kappa\_{i+n}^{(j)} + v\_{i+n}^{(j)})^2 \right\} \cdot \prod\_{k=1}^n \mathcal{S}\_k(\kappa\_k + v\_k) \\ &= \int\_0^\infty d^n \beta \, F(\beta\_1, \ldots, \beta\_n), \end{split} \tag{6}$$

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:

$$1 = \int\_0^\infty dt \delta\left(t - \sum\_{i=1}^n \beta\_i\right).$$

One has:

$$\Pi = \int\_0^\infty dt \, t^{n-1} \int d^n a \, \delta \left( 1 - \sum\_{i=1}^n a\_i \right) F(ta\_1, \dots, ta\_n). \tag{7}$$

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 *<sup>t</sup>*, i.e., <sup>∞</sup> <sup>→</sup> 1/*λ*<sup>2</sup> with *<sup>λ</sup>* 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:

$$\Pi(p^2) = \int \frac{d^4k\_E}{\pi^2} \frac{\exp\left(-sk\_E^2\right)}{[m^2 + (k\_E + p\_E/2)^2][m^2 + (k\_E - p\_E/2)^2]}.\tag{8}$$

By using the above transformations, one gets:

$$\Pi(p^2) = \int\_0^\infty dt \, \frac{t}{(s+t)^2} \int\_0^1 da \varepsilon^{-t[m^2 - a(1-a)p^2] + \frac{pt}{s+t}(a-1/2)^2 p^2} \,\tag{9}$$

where *<sup>p</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>*p*<sup>2</sup> *<sup>E</sup>*. The expression has a branch point at *<sup>p</sup>*<sup>2</sup> = <sup>4</sup>*m*<sup>2</sup> 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.

#### *2.4. 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:

$$
\partial^{\mu}q \to (\partial^{\mu} - i e\_{q} A^{\mu}) q, \qquad \partial^{\mu}\bar{q} \to (\partial^{\mu} + i e\_{q} A^{\mu}) \bar{q}, \tag{10}
$$

with *eu* = 2/3, *ed* = −1/3 in units of the proton charge. This defines the first part of the interaction Lagrangian of quarks with photons:

$$\mathcal{L}\_{\text{int}}^{EM(1)} = \sum\_{q} \varepsilon\_{q} A\_{\mu}(\mathbf{x}) J\_{q}^{\mu}(\mathbf{x}), \quad \text{with} \quad J\_{q}^{\mu}(\mathbf{x}) = \overline{q}(\mathbf{x}) \gamma^{\mu} q(\mathbf{x}). \tag{11}$$

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:

$$q(\mathbf{x}\_i) \to \exp\{-ie\_q I(\mathbf{x}\_i, \mathbf{x}, P)\} q(\mathbf{x}\_i), \qquad \bar{q}(\mathbf{x}\_i) \to \exp\{ie\_q I(\mathbf{x}\_i, \mathbf{x}, P)\} \bar{q}(\mathbf{x}\_i), \tag{12}$$

with *I* being defined as the integral:

$$I(\mathbf{x}\_{i\prime}, \mathbf{x}\_{\prime}P) = \int\_{\mathbf{x}}^{\mathbf{x}\_{i}} A^{\mu}(z)dz\_{\mu}$$

over the path that connects the hadron and quark positions. It can be readily seen that the local gauge transformations:

$$q(\mathbf{x}\_i) \to \quad \exp\{ie\_q f(\mathbf{x}\_i)\} q(\mathbf{x}\_i), \quad \not q(\mathbf{x}\_i) \to \exp\{-ie\_q f(\mathbf{x}\_i)\} \not q(\mathbf{x}\_i),$$

$$A^\mu(z) \to A^\mu(z) + \partial^\mu f(z), \tag{13}$$

leave the Lagrangian unchanged for any arbitrary function *f* . Indeed, the gauge field induced modification of the path integral *I*(*xi*, *x*, *P*) → *I*(*xi*, *x*, *P*) + *f*(*xi*) − *f*(*x*) is canceled by the contribution coming from the quark transformations. The exact form of the gauged non-local Lagrangian <sup>L</sup>*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:

$$\lim\_{dx^{\mu}\to 0} dx^{\mu} \frac{\partial}{\partial \mathbf{x}^{\mu}} I(\mathbf{x}, \mathbf{y}, P) = \lim\_{dx^{\mu}\to 0} \left[ I(\mathbf{x} + dx, \mathbf{y}, P') - I(\mathbf{x}, \mathbf{y}, P) \right]. \tag{14}$$

Here, *P*denotes a path derived from *P* by extending *P* from its endpoint by *dx*. This definition gives:

$$\frac{\partial}{\partial \mathbf{x}^{\mu}} I(\mathbf{x}, \mathbf{y}, \mathbf{P}) = A\_{\mu}(\mathbf{x}), \tag{15}$$

where the independence of the derivative on the path *P* becomes explicit.

#### *2.5. Selected Computational Aspects*

To proceed with the calculations, it is convenient to use the following representation for the correlation function:

$$\Phi\_X\left(\sum\_{1\le i$$

It may be easily obtained by using the Jacobi coordinates. The Gaussian function form of the correlation function Φ *<sup>X</sup>* 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:

$$\exp(kak + 2kr + R),\tag{17}$$

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:

$$P\left(k\_i^{\mu}\right)\varepsilon^{2kr} = P\left(\frac{1}{2}\frac{\partial}{\partial r\_{i\mu}}\right)\varepsilon^{2kr},\tag{18}$$

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:

$$P\left(\frac{1}{2}\frac{\partial}{\partial r\_i}\right)e^{-ra^{-1}r} = e^{-ra^{-1}r}P\left(\frac{1}{2}\frac{\partial}{\partial r\_i} - [a^{-1}r]\_i\right),\tag{19}$$

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*<sup>2</sup> 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.

#### **3. Strong Decays of X**(**3872**)
