1 Introduction

Recently, the LHCb Collaboration reported on evidence for an \(\eta _{c}(1S)\pi ^{-}\) resonance in \(B^{0}\rightarrow K^{+}\eta _{c}(1S)\pi ^{-}\) decays extracted from analysis of pp collisions’ data collected with LHCb detector at center-of-mass energies of \(\sqrt{s}=7,\ 8\) and \(13\ \mathrm {TeV} \) [1]. The mass and width of this new \(Z_{c}^{-}(4100)\) resonance (hereafter \(Z_{c}\)) were found equal to \(m=4096\pm 20_{-22}^{ +18}~ \text {MeV}\) and \(\Gamma =152\pm 58_{-35}^{+60}~\text {MeV}\), respectively. As it was emphasized in Ref. [1], the spin-parity assignments \(J^{P}=0^{+}\) and \(J^{P}=1^{-}\) both are consistent with the data.

From analysis of the decay channel \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\) it becomes evident that \(Z_{c}\) contains four quarks \(cd{\overline{c}} {\overline{u}}\), and it is presumably another member of the family of charged exotic Z-resonances with the same quark content; the well-known axial-vector tetraquarks \(Z_{c}^{\pm }(4430)\) and \(Z_{c}^{\pm }(3900)\) are also built of the quarks \(cd{\overline{c}}{\overline{u}}\) or \(cu{\overline{c}} {\overline{d}}\). The \(Z_{c}^{\pm }(4430)\) were discovered and studied by the Belle Collaboration in B meson decays \(B\rightarrow K\psi ^{\prime }\pi ^{\pm }\) as resonances in the \(\psi ^{\prime }\pi ^{\pm }\) invariant mass distributions [2,3,4]. The decay of \(Z_{c}^{+}(4430)\) to the final state \(J/\psi \pi ^{+}\) also was detected in the Belle experiment [5]. The existence of the \(Z_{c}^{\pm }(4430)\) resonances was confirmed by the LHCb Collaboration as well [6, 7].

Another well-known members of this family are the axial-vector states \( Z_{c}^{\pm }(3900)\), which were detected by the BESIII Collaboration in the process \(e^{+}e^{-}\rightarrow J/\psi \pi ^{+}\pi ^{-}\) as peaks in the \( J/\psi \pi ^{\pm }\) invariant mass distributions [8]. These structures were seen by the Belle and CLEO collaborations as well (see Refs. [9, 10]). The BESIII informed also on observation of the neutral \(Z_{c}^{0}(3900)\) state in the process \( e^{+}e^{-}\rightarrow \pi ^{0}Z_{c}^{0}\rightarrow \pi ^{0}\pi ^{0}J/\psi \) [11].

Various theoretical models and approaches were employed to reveal the internal quark-gluon structure and determine parameters of the charged Z-resonances. Thus, they were considered as hadrocharmonium compounds or tightly bound diquark–antidiquark states, were treated as the four-quark systems built of conventional mesons or interpreted as threshold cusps (see Refs. [12, 13] and references therein).

The diquark model of the exotic four-quark mesons is one of the popular approaches to explain their properties. In accordance with this picture the tetraquark is a bound state of a diquark and an antidiquark. This approach implies the existence of multiplets of the diquark–antidiquarks with the same quark content, but different spin-parities. Because the resonances \( Z_{c}^{\pm }(3900)\) and \(Z_{c}^{\pm }(4430)\) are the axial-vectors, one can interpret them as the ground-state 1S and first radially excited 2S state of the same \([cu][{\overline{c}}{\overline{d}}]\) or \([cd][{\overline{c}} {\overline{u}}]\) multiplets. An idea to consider \(Z_{c}(4430)\) as a radial excitation of the \(Z_{c}(3900)\) state was proposed in Ref. [14], and explored in Refs. [15, 16] in the framework of the QCD sum rule method.

The resonances \(Z_{c}(3900)\), \(Z_{c}(4200)\), \(Z_{c}(4430)\) and \(Z_{c}\) were detected in B meson decays and/or electron-positron annihilations, which suggest that all of them may have the same nature. Therefore, one can consider \(Z_{c}\) as the ground-state scalar or vector tetraquark with \(c {\overline{c}}d{\overline{u}}\) content. The recent theoretical articles devoted to the \(Z_{c}\) resonance are concentrated mainly on exploration of its spin and possible decay channels [17,18,19,20,21]. Thus, sum rule computations carried out in Ref. [17] demonstrated that \(Z_{c}\) is presumably a scalar particle rather than a vector tetraquark. The conclusion about a tetraquark nature of \(Z_{c}\) with quantum numbers \( J^{PC}=0^{++}\) was drawn in Ref. [18] as well. In the hadrocharmonium framework the resonances \(Z_{c}\) and \(Z_{c}^{-}(4200)\) were treated as the scalar \(\eta _{c}\) and vector \(J/\psi \) charmonia embedded in a light-quark excitation with quantum numbers of a pion [19]. In accordance with this picture \(Z_{c}\) and \( Z_{c}^{-}(4200)\) are related by the charm quark spin symmetry which suggests certain relations between their properties and decay channels. The possible decays of a scalar and a vector tetraquark \([cd][{\overline{c}}{\overline{u}}]\) were analyzed also in Ref. [20].

In the present work we treat \(Z_{c}\) as the scalar diquark–antidiquark state \([cd][{\overline{c}}{\overline{u}}]\), since it was observed in the process \( Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\). In fact, for the scalar \(Z_{c}\) this decay is the dominant S-wave channel, whereas for the vector tetraquark \(Z_{c}\) it turns P-wave decay. We are going to calculate the spectroscopic parameters of the tetraquark \(Z_{c}\), i.e., its mass and coupling by means of the two-point QCD sum rule method. The QCD sum rule method is the powerful nonperturbative approach to investigate the conventional hadrons and calculate their parameters [22, 23]. But it can be successfully applied for studying multiquark systems as well. To get reliable predictions for the quantities of concern, in the sum rule computations we take into account the quark, gluon, and mixed vacuum condensates up to dimension ten.

The next problem to be considered in this work is investigating decays of the resonance \(Z_{c}\) and evaluating its total width. It is known, that strong and semileptonic decays of various tetraquark candidates provide valuable information on their internal structure and dynamical features. In the framework of the QCD sum rule approach relevant problems were subject of rather intensive studies [24,25,26,27,28,29,30,31,32,33,34,35,36,37]. The dominant strong decay of the resonance \(Z_{c}\) seems is the channel \( Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\). But its S-wave hidden-charm \(\eta _{c}(2S)\pi ^{-}\), \(J/\psi \rho ^{-}\) and open-charm \(D^{0}D^{-}\) and \( D^{*0}D^{*-}\)decays are also kinematically allowed modes [20].

We calculate the partial width of the dominant S-wave processes and use obtained results to evaluate the total width of the tetraquark \(Z_{c}\). The decays \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\), \(\eta _{c}(2S)\pi ^{-}\), and \(D^{0}D^{-}\) are explored by applying the QCD three-point sum rule method. The quantities extracted from the sum rules are the strong couplings \(g_{Z_{c}\eta _{c1}\pi }\), \(g_{Z_{c}\eta _{c2}\pi }\), and \(g_{Z_{c}DD}\) that correspond to the vertices \(Z_{c}\eta _{c}(1S)\pi ^{-}\), \(Z_{c}\eta _{c}(2S)\pi ^{-}\), and \(Z_{c}D^{0}D^{-}\), respectively. The coupling \( g_{Z_{c}J/\psi \rho }\), which describes the strong vertex \(Z_{c}J/\psi \rho ^{-}\), is found by means of the QCD light-cone sum rule (LCSR) method and technical tools of the soft-meson approximation [38, 39]. For analysis of the tetraquarks this method and approximation was adapted in Ref. [26], and applied to study their numerous strong decay channels. Alongside the mass and coupling of the state \(Z_{c}\) the strong couplings provide an important information to determine the width of the decays under analysis.

This work is organized in the following manner: In Sect. 2 we calculate the mass m and coupling f of the scalar resonance \(Z_{c}\) by employing the two-point sum rule method and including into analysis the quark, gluon, and mixed condensates up to dimension ten. The obtained predictions for these parameters are applied in Sect. 3 to evaluate the partial widths of the decays \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\) and \(\eta _{c}(2S)\pi ^{-}\). The decay \(Z_{c}\rightarrow D^{0}D^{-}\) is considered in Sect. 4, whereas Sect. 5 is devoted to analysis of the decay \(Z_{c}\rightarrow J/\psi \rho ^{-}\). In Sect. 5 we also give our estimate for the total width of the resonance \(Z_{c}\). The Sect. 6 contains the analysis of obtained results and our concluding notes. In the Appendix we write down explicit expressions of the heavy and light quark propagators, as well as the two-point spectral density used in the mass and coupling calculations.

2 Mass and coupling of the scalar tetraquark \(Z_{c}\)

The scalar resonance \(Z_{c}\) can be composed of the scalar diquark \(\epsilon ^{ijk}[c_{j}^{T}C\gamma _{5}d_{k}]\) in the color antitriplet and flavor antisymmetric state and the scalar antidiquark \(\epsilon ^{imn}[{\overline{c}} _{m}\gamma _{5}C{\overline{u}}_{n}^{T}]\) in the color triplet state. These diquarks are most attractive ones, and four-quark mesons composed of them should be lighter and more stable than bound states of diquarks with other quantum numbers [40]. The scalar diquarks were used as building blocks to construct various hidden-charm and -bottom tetraquark states and study their properties (see, for example Refs. [41,42,43]). In the present work for the resonance \(Z_{c}\) we choose namely this favorable structure.

To calculate the mass m and coupling f of the resonance \(Z_{c}\) using the QCD sum rule method, we start from the two-point correlation function

$$\begin{aligned} \Pi (p)=i\int d^{4}xe^{ip\cdot x}\langle 0|{\mathcal {T}}\{J(x)J^{\dagger }(0)\}|0\rangle , \end{aligned}$$
(1)

where J(x) is the interpolating current for the tetraquark \(Z_{c}\). In accordance with our assumption on the structure of \(Z_{c}\) the interpolating current J(x) has the following form

$$\begin{aligned} J(x)=\epsilon {\tilde{\epsilon }}[c_{j}^{T}(x)C\gamma _{5}d_{k}(x)] [{\overline{c}}_{m}(x)\gamma _{5}C{\overline{u}}_{n}^{T}(x)]. \end{aligned}$$
(2)

Here we employ the notations \(\epsilon =\epsilon ^{ijk}\) and \({\tilde{\epsilon }}=\epsilon ^{imn}\), where ijkm and n are color indices, and C is the charge-conjugation operator.

To derive the sum rules for the mass m and coupling f of the ground-state tetraquark \(Z_{c}\) we adopt the “ground-state + continuum” approximation, and calculate the physical or phenomenological side of the sum rule. For these purposes, we insert into the correlation function a full set of relevant states and carry out in Eq. (1) the integration over x, and get

$$\begin{aligned} \Pi ^{\mathrm {Phys}}(p)=\frac{\langle 0|J|Z_{c}(p)\rangle \langle Z_{c}(p)|J^{\dagger }|0\rangle }{m^{2}-p^{2}}+\cdots \end{aligned}$$
(3)

Here we separate the ground-state contribution to \(\Pi ^{\mathrm {Phys}}(p)\) from effects of the higher resonances and continuum states, which are denoted there by the dots. In the calculations we assume that the phenomenological side \(\Pi ^{\mathrm {Phys}}(p)\) can be approximated by a single pole term. In the case of the multiquark systems the physical side, however, receives contribution also from two-meson reducible terms [44, 45]. In other words, the interpolating current J(x) interacts with the two-meson continuum, which generates the finite width \( \Gamma (p^{2})\) of the tetraquark and results in the modification [46]

$$\begin{aligned} \frac{1}{m^{2}-p^{2}}\rightarrow \frac{1}{m^{2}-p^{2}-i\sqrt{p^{2}}\Gamma (p^{2})}. \end{aligned}$$
(4)

The two-meson continuum effects can be properly taken into account by rescaling the coupling f, whereas the mass of the tetraquark m preserves its initial value. But these effects are numerically small, therefore in the phenomenological side of the sum rule we use the zero-width single-pole approximation and check afterwards its self-consistency.

Calculation of \(\Pi ^{\mathrm {Phys}}(p)\) can be finished by introducing the matrix element of the scalar tetraquark

$$\begin{aligned} \langle 0|J|Z_{c}\rangle =fm. \end{aligned}$$
(5)

As a result, we find

$$\begin{aligned} \Pi ^{\mathrm {Phys}}(p)=\frac{f^{2}m^{2}}{m^{2}-p^{2}}+\cdots . \end{aligned}$$
(6)

Because \(\Pi ^{\mathrm {Phys}}(p)\) has trivial Lorentz structure proportional to I, corresponding invariant amplitude \(\Pi ^{\mathrm {Phys}}(p^{2})\) is equal to the function given by Eq. (6).

At the next step one has to find the correlation function \(\Pi (p)\) using the perturbative QCD and express it through the quark propagators, and, as a result, in terms of the vacuum expectation values of various quark, gluon and mixed operators as nonperturbative effects. To this end, we use the interpolating current J(x), contract the relevant heavy and light quark operators in Eq. (1) to generate propagators, and obtain

$$\begin{aligned} \Pi ^{\mathrm {OPE}}(p)= & {} i\int d^{4}xe^{ip\cdot x}\epsilon {\tilde{\epsilon }} \epsilon ^{\prime }{\tilde{\epsilon }}^{\prime }{\mathrm {Tr}}\left[ \gamma _{5} {\widetilde{S}}_{c}^{jj^{\prime }}(x)\gamma _{5}S_{d}^{kk^{\prime }}(x)\right] \nonumber \\&\times \, {\mathrm {Tr}}\left[ \gamma _{5}{\widetilde{S}}_{u}^{n^{\prime }n}(-x)\gamma _{5}S_{c}^{m^{\prime }m}(-x)\right] . \end{aligned}$$
(7)

Here \(S_{c}(x)\) and \(S_{u(d)}(x)\) are the heavy c- and light u(d)-quark propagators, respectively. These propagators contain both the perturbative and nonperturbative components: their explicit expressions are presented in the Appendix. In Eq. (7) we also utilized the shorthand notation

$$\begin{aligned} {\widetilde{S}}_{c(u)}(x)=CS_{c(u)}^{T}(x)C. \end{aligned}$$
(8)

To extract the required sum rules for m and f one must equate \(\Pi ^{ \mathrm {Phys}}(p^{2})\) to the similar amplitude \(\Pi ^{\mathrm {OPE}}(p^{2})\), apply the Borel transformation to both sides of the obtained equality to suppress contributions of the higher resonances and, finally, perform the continuum subtraction in accordance with the assumption on the quark-hadron duality: These manipulations lead to the equality that can be used to get the sum rules. The second equality, which is required for these purposes, can be obtained from the first expression by applying on it by the operator \( d/d(-1/M^{2})\). Then, for the mass of the tetraquark \(Z_{c}\) we get the sum rule

$$\begin{aligned} m^{2}=\frac{\int _{4m_{c}^{2}}^{s_{0}}dss\rho ^{\mathrm {OPE}}(s)e^{-s/M^{2}}}{ \int _{4m_{c}^{2}}^{s_{0}}ds\rho ^{\mathrm {OPE}}(s)e^{-s/M^{2}}}. \end{aligned}$$
(9)

The sum rule for the coupling f reads

$$\begin{aligned} f^{2}=\frac{1}{m^{2}}\int _{4m_{c}^{2}}^{s_{0}}ds\rho ^{\mathrm {OPE} }(s)e^{(m^{2}-s)/M^{2}}. \end{aligned}$$
(10)

where \(M^{2}\) and \(s_{0}\) are the Borel and continuum threshold parameters, respectively. In Eqs. (9) and (10) \(\rho ^{\mathrm { OPE}}(s)\) is the two-point spectral density, which is proportional to the imaginary part of the correlation function \(\Pi ^{\mathrm {OPE}}(p).\) The explicit expression of \(\rho ^{\mathrm {OPE}}(s)\) is presented in the Appendix.

We use the obtained sum rules to compute the mass m and coupling f of the tetraquark \(Z_{c}\). They contain numerous parameters, some of which, such as the vacuum condensates, the mass of the c-quark, are universal quantities and do not depend on the problem under discussion. In computations we utilize the following values for the quark, gluon and mixed condensates:

$$\begin{aligned}&\langle {\bar{q}}q\rangle =-(0.24\pm 0.01)^{3}~\text {GeV}^{3},\nonumber \\&m_{0}^{2}=(0.8\pm 0.1)~\text {GeV}^{2},\quad \langle {\overline{q}}g_{s}\sigma Gq\rangle =m_{0}^{2}\langle {\overline{q}}q\rangle , \nonumber \\&\left\langle \frac{\alpha _{s}G^{2}}{\pi }\right\rangle =(0.012\pm 0.004)~\text {GeV}^{4}, \nonumber \\&\langle g_{s}^{3}G^{3}\rangle =(0.57\pm 0.29)~\text {GeV}^{6}. \end{aligned}$$
(11)

The mass of the c-quark is taken equal to \(m_{c}=1.275_{-0.035}^{+0.025} \text {GeV}\).

The Borel parameter \(M^{2}\) and continuum threshold \(s_{0}\) are the auxiliary parameters and should be chosen in accordance with standard constraints accepted in the sum rule computations. The Borel parameter can be varied within the limits \([M_{\mathrm {min}}^{2},\ M_{\mathrm {max}}^{2}]\) which have to obey the following conditions: At \(M_{\mathrm {max}}^{2}\) the pole contribution (\(\mathrm {PC)}\) defined as the ratio

$$\begin{aligned} {\mathrm {PC}}=\frac{\Pi (M_{\mathrm {max}}^{2},\ s_{0})}{\Pi (M_{\mathrm {max} }^{2},\ \infty )}, \end{aligned}$$
(12)

should be larger than some fixed number. Let us note that \(\Pi (M^{2},\ s_{0})\) in Eq. (12) is the Borel transformed and subtracted invariant amplitude \(\Pi ^{\mathrm {OPE}}(p^{2})\). In the sum rule calculations involving the tetraquarks the minimal value of \({\mathrm {PC}}\) varies between 0.15–0.2. In the present work we choose \({\mathrm {PC}}>0.15\). The minimal value of the Borel parameter \(M_{\mathrm {min}}^{2}\) is fixed from convergence of the sum rules: in other words, at \(M_{\mathrm {min}}^{2}\) contribution of the last term (or a sum of last few terms) to \(\Pi (M^{2},\ s_{0})\) cannot exceed 0.05 part of the whole result

$$\begin{aligned} R(M_{\mathrm {min}}^{2})=\frac{\Pi ^{\mathrm {DimN}}(M_{\mathrm {min}}^{2},\ s_{0})}{\Pi (M_{\mathrm {min}}^{2},\ s_{0})}<0.05. \end{aligned}$$
(13)

The ratio \(R(M_{\mathrm {min}}^{2})\) quantifies the convergence of the OPE and will be used for the numerical analysis. The last restriction on the lower limit \(M_{\mathrm {min}}^{2}\) is the prevalence of the perturbative contribution over the nonperturbative one.

The mass m and coupling f should not depend on the parameters \(M^{2}\) and \(s_{0}\). But in real calculations, these quantities are sensitive to the choice of \(M^{2}\) and \(s_{0}\). Therefore, the parameters \(M^{2}\) and \(s_{0}\) should also be determined in such a way that to minimize the dependence of m and f on them. The analysis allows us to fix the working windows for the parameters \(M^{2}\) and \(s_{0}\)

$$\begin{aligned} M^{2}\in [4,\ 6]~\text {GeV}^{2},\quad s_{0}\in [19,\ 21]~\text {GeV}^{2}, \end{aligned}$$
(14)

which obey all aforementioned restrictions. Thus, at \(M^{2}=6~\text {GeV} ^{2}\) the pole contribution equals to 0.19, and within the region \( M^{2}\in [4,~6]~\text {GeV}^{2}\) it changes from 0.54 till 0.19. To find the lower bound of the Borel parameter from Eq. (13) we use the last three terms in the expansion, i.e. \( \mathrm {DimN}=\mathrm {Dim(8+9+10)}\) [we remind that \(\Pi ^{\mathrm {Dim10}}=0\)]. Then at \(M^{2}=4~\text {GeV}^{2}\) the ratio R becomes equal to \(R(4~ \text {GeV}^{2})=0.02\) which ensures the convergence of the sum rules. At \( M^{2}=4~\text {GeV}^{2}\) the perturbative contribution amounts to \(83\%\) of the full result exceeding considerably the nonperturbative terms.

As it has been noted above, there are residual dependence of m and f on the parameters \(M^{2}\) and \(s_{0}\). In Figs. 1 and 2 we plot the mass and coupling of the tetraquark \(Z_{c}\) as functions of the parameters \(M^{2}\) and \(s_{0}\). It is seen that both the m and f depend on \(M^{2}\) and \(s_{0}\) which generates essential part of the theoretical uncertainties inherent to the sum rule computations. For the mass m these uncertainties are small which has a simple explanation: The sum rule for the mass (9) is equal to the ratio of integrals over the functions \(s\rho ^{\mathrm {OPE}}(s)\) and \(\rho ^{\mathrm {OPE}}(s)\), which reduces effects due to variation of \(M^{2}\) and \(s_{0}\). The coupling f depends on the integral over the spectral density \(\rho ^{\mathrm {OPE} }(s) \), and therefore its variations are sizeable. In the case under analysis, theoretical errors for m and f generated by uncertainties of various parameters including \(M^{2}\) and \(s_{0}\) ones equal to \(\pm \, 3.7\%\) and \(\pm \, 21\%\) of the corresponding central values, respectively.

Our analysis leads for the mass and coupling of the tetraquark \(Z_{c}\) to the following results:

$$\begin{aligned} m= & {} (4080~\pm 150)~\text {MeV}, \nonumber \\ f= & {} (0.58\pm 0.12)\cdot 10^{-2}~\text {GeV}^{4}. \end{aligned}$$
(15)

The mass of the resonance \(Z_{c}\) modeled as the scalar diquark–antidiquark state is in excellent agreement with the data of the LHCb Collaboration.

The scalar tetraquark \(cu{\overline{c}}{\overline{d}}\) with \(C\gamma _{5}\otimes \gamma _{5}C\) structure was studied also in Ref. [47]. The prediction for the mass of this four-quark meson \(m=(3860~\pm 90)~\mathrm {MeV }\) allowed the author to interpret it as the resonance \(X^{*}(3860)\) (actually, its charged partner) observed recently by the Belle Collaboration [48]. This charmoniumlike state was seen in the process \( e^{+}e^{-}\rightarrow J/\psi D{\overline{D}}\), where D refers to either \( D^{0}\) or \(D^{+}\) meson, and was considered there as a \(\chi _{c0}(2P)\) candidate. Comparing our result and that of Ref. [47], one can see the existence of an overlapping region between them, nevertheless the difference between the central values \(200\ ~\text {MeV}\) is sizable. This discrepancy is presumably connected with working regions for \(M^{2}\) and \(s_{0}\), and also with values of the vacuum condensates (fixed or evolved) used in numerical computations.

Fig. 1
figure 1

The mass of the state \(Z_c^{-}(4100)\) as a function of the Borel parameter \(M^2\) at fixed \(s_0\) (left panel), and as a function of the continuum threshold \(s_0\) at fixed \(M^2\) (right panel)

Fig. 2
figure 2

The same as in Fig. 1, but for the coupling f of the resonance \(Z_c^{-}(4100)\)

3 Decays \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\) and \(Z_{c}\rightarrow \eta _{c}(2S)\pi ^{-}\)

The S-wave decays of the resonance \(Z_{c}\) can be divided into two subclasses: The decays to two pseudoscalar and two vector mesons, respectively. The processes \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\) and \( Z_{c}\rightarrow \eta _{c}(2S)\pi ^{-}\) belong to the first subclass of decays. The final stages of these decays contain the ground-state and first radially excited \(\eta _{c}\) mesons, therefore in the QCD sum rule approach they should investigated in a correlated form. An appropriate way to deal with decays \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\) and \(Z_{c}\rightarrow \eta _{c}(2S)\pi ^{-}\) is the QCD three-point sum rule method. Indeed, because we are going to explore the form factors \(g_{Z_{c}\eta _{ci}\pi }(q^{2})\) for the off-shell pion the double Borel transformation will be carried out in the \(Z_{c}\) and \(\eta _{c}\) channels, i.e. over momenta of these particles. This transformation applied to the phenomenological side of the relevant three-point sum rules suppresses contributions of the higher resonances in these two channels eliminating, at the same time, terms associated with the pole-continuum transitions [39, 49]. The elimination of these terms is important for joint analysis of the form factors \(g_{Z_{c}\eta _{ci}\pi }(q^{2})\), because one does not need to apply an additional operator to remove them from the phenomenological side of the sum rules. Nevertheless, there may still exist in the pion channel terms corresponding to excited states of the pion which emerge as contaminations (for the \(NN\pi \) vertex, see discussions in Refs. [50, 51]). To reduce the uncertainties in evaluation of the strong couplings at the vertices and smooth problems with extrapolation of the form factors to the mass-shell, it is possible to fix the pion on the mass-shell and treat one of the remaining heavy states (\(Z_{c}\) or \(\eta _{c}\)) as the off-shell particle. This trick was used numerously to study the conventional heavy-heavy-light mesons’ couplings [52, 53]. Form factors obtained by treating a light or one of heavy mesons off-shell may differ from each other considerably, but after extrapolating to the corresponding mass-shells lead to the same or slightly different strong couplings.

In the framework of the three-point sum rule approach a more detailed representation for the phenomenological side was used in Refs. [31, 36, 37]. This technique generates additional terms in the sum rules and introduces into analysis new free parameters, which should be chosen to obtain stable sum rules with variations of the Borel parameters. In the present work, to calculate \( g_{Z_{c}\eta _{c1}\pi }(q^{2})\) and \(g_{Z_{c}\eta _{c2}\pi }(q^{2})\) we apply the standard three-point sum rule method and choose the pion an off-shell particle. We use this method to study the decay \(Z_{c}\rightarrow D^{0}D^{-}\) as well.

The process \(Z_{c}\rightarrow \) \(J/\psi \rho ^{-}\) belongs to the second subclass of \(Z_{c}\) decays; it is a decay to two vector mesons. We investigate this mode by means of the QCD light-cone sum rule method and soft-meson approximation. The sum rule on the light-cone allows one to find the strong coupling by avoiding extrapolating procedures and express \( g_{Z_{c}J/\psi \rho }\) not only in terms of the vacuum condensates, but also using the \(\rho \)-meson local matrix elements. As for unsuppressed pole-continuum effects that after a single Borel transformation survive in this approach, they can be eliminated by means of well-known prescriptions [49].

To determine the partial widths of the decays \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\) and \(Z_{c}\rightarrow \eta _{c}(2S)\pi ^{-}\) one needs to calculate the strong couplings \(g_{Z_{c}\eta _{c1}\pi }\) and \(g_{Z_{c}\eta _{c2}\pi }\) which can be extracted from the three-point correlation function

$$\begin{aligned} \Pi (p,p^{\prime })= & {} i^{2}\int d^{4}xd^{4}ye^{-ip\cdot x}e^{ip^{\prime }\cdot y} \nonumber \\&\times \, \langle 0|{\mathcal {T}}\{J^{\eta }(y)J^{\pi }(0)J^{\dagger }(x)\}|0\rangle , \end{aligned}$$
(16)

where

$$\begin{aligned} J^{\eta }(y)={\overline{c}}_{a}(y)i\gamma _{5}c_{a}(y),\quad J^{\pi }(0)={\overline{u}}_{b}(0)i\gamma _{5}d_{b}(0) \end{aligned}$$
(17)

are the interpolating currents for the pseudoscalar mesons \(\eta _{c}\) and \( \pi ^{-}\), respectively. The J(x) is the interpolating current for the resonance \(Z_{c}\) and has been introduced above in Eq. (2).

In terms of the physical parameters of the tetraquark and mesons the correlation function \(\Pi (p,p^{\prime })\) takes the form

$$\begin{aligned}&\Pi ^{\mathrm {Phys}}(p,p^{\prime })=\sum _{i=1}^{2}\frac{\langle 0|J^{\eta }|\eta _{ci}\left( p^{\prime }\right) \rangle }{p^{\prime 2}-m_{i}^{2}}\frac{ \langle 0|J^{\pi }|\pi \left( q\right) \rangle }{q^{2}-m_{\pi }^{2}} \nonumber \\&\quad \times \, \frac{\langle \eta _{ci}\left( p^{\prime }\right) \pi (q)|Z_{c}(p)\rangle \langle Z_{c}(p)|J^{\dagger }|0\rangle }{p^{2}-m^{2}} +\cdots , \end{aligned}$$
(18)

where \(m_{\pi }\) is the mass of the pion, and \(m_{i}=m_{1}\), \(m_{2}\) are masses of the mesons \(\eta _{c}(1S)\) and \(\eta _{c}(2S)\), respectively. The four-momenta of the particles are evident from (18). Here by the dots we denote contribution of the higher resonances and continuum states.

To continue we introduce the matrix elements

$$\begin{aligned} \langle 0|J^{\eta }|\eta _{ci}\left( p^{\prime }\right) \rangle =\frac{ f_{i}m_{i}^{2}}{2m_{c}}, \end{aligned}$$
(19)

where \(f_{1}\) and \(f_{2}\) are the decay constants of the mesons \(\eta _{c}(1S)\) and \(\eta _{c}(2S)\), respectively. The relevant matrix element of the pion is well known and has the form

$$\begin{aligned} \langle 0|J^{\pi }|\pi \left( q\right) \rangle =\mu _{\pi }f_{\pi },\quad \mu _{\pi }=-\frac{2\langle {\overline{q}}q\rangle }{f_{\pi }^{2}}, \end{aligned}$$
(20)

where \(f_{\pi }\) is the decay constant of the pion, and \(\langle {\overline{q}} q\rangle \) is the quark condensate. Additionally, the matrix elements of the vertices \(Z_{c}\eta _{c}(1S)\pi ^{-}\) and \(\ Z_{c}\eta _{c}(2S)\pi ^{-}\) are required. To this end, we use

$$\begin{aligned} \langle \eta _{ci}\left( p^{\prime }\right) \pi (q)|Z_{c}(p)\rangle =g_{Z_{c}\eta _{ci}\pi }(p\cdot p^{\prime }). \end{aligned}$$
(21)

Here, the strong coupling \(g_{Z_{c}\eta _{c1}\pi }\) corresponds to the vertex \(Z_{c}\eta _{c}(1S)\pi ^{-}\), whereas \(g_{Z_{c}\eta _{c2}\pi }\) describes \(Z_{c}\eta _{c}(2S)\pi ^{-}\); namely these couplings have to be determined from the sum rules.

Employing Eqs. (19)–(21) for \( \Pi ^{\mathrm {Phys}}(p,p^{\prime })\) we get the simple expression:

$$\begin{aligned} \Pi ^{\mathrm {Phys}}(p,p^{\prime })= & {} \sum _{i=1}^{2}\frac{g_{Z_{c}\eta _{ci}\pi }m_{i}^{2}f_{i}mf}{2m_{c}(p^{\prime 2}-m_{i}^{2})\left( p^{2}-m^{2}\right) } \nonumber \\&\times \frac{\mu _{\pi }f_{\pi }}{q^{2}-m_{\pi }^{2}}(p\cdot p^{\prime })+\cdots . \end{aligned}$$
(22)

The Lorentz structure of the \(\Pi ^{\mathrm {Phys}}(p,p^{\prime })\) is proportional to I therefore the invariant amplitude \(\Pi ^{\mathrm {Phys} }(p^{2},p^{\prime 2},q^{2})\) is given by the sum of two terms from Eq. (22). The double Borel transformation of \(\Pi ^{\mathrm {Phys} }(p^{2},p^{\prime 2},q^{2})\) over the variables \(p^{2}\) and \(p^{\prime 2}\) with the parameters \(M_{1}^{2}\) and \(M_{2}^{2}\) forms one of sides in the sum rule equality.

The QCD side of the sum rule, i.e. the expression of the correlation function in terms of the quark propagators reads

$$\begin{aligned}&\Pi ^{\mathrm {OPE}}(p,p^{\prime })=i^{2}\int d^{4}xd^{4}ye^{-ip\cdot x}e^{ip^{\prime }\cdot y}\epsilon ^{ijk}\epsilon ^{imn} \nonumber \\&\quad \times {\mathrm {Tr}}\left[ \gamma _{5}S_{c}^{aj}(y-x)\gamma _{5}{\widetilde{S}} _{d}^{bk}(-x)\gamma _{5}{\widetilde{S}}_{u}^{nb}(x)\gamma _{5}S_{c}^{ma}(x-y) \right] .\nonumber \\ \end{aligned}$$
(23)

The Borel transformation \({\mathcal {B}}\Pi ^{\mathrm {OPE}}(p^{2},p^{\prime 2},q^{2})\), where \(\Pi ^{\mathrm {OPE}}(p^{2},p^{\prime 2},q^{2})\) is the invariant amplitude that corresponds to the structure \(\sim I\) in \(\Pi ^{ \mathrm {OPE}}(p,p^{\prime })\) constitutes the second component of the sum rule. Equating \({\mathcal {B}}\Pi ^{\mathrm {OPE}}(p^{2},p^{\prime 2},q^{2})\) and the double Borel transformation of \(\Pi ^{\mathrm {Phys}}(p^{2}, p^{\prime 2},q^{2})\) and performing continuum subtraction we find the sum rule for the couplings \(g_{Z_{c}\eta _{c1}\pi }\) and \(g_{Z_{c}\eta _{c2}\pi }\).

The Borel transformed and subtracted amplitude \(\Pi ^{\mathrm {OPE}} (p^{2},p^{\prime 2},q^{2})\) can be expressed in terms of the spectral density \(\rho _{\mathrm {D}}(s,s^{\prime },q^{2})\) which is proportional to the imaginary part of \(\Pi ^{\mathrm {OPE}}(p,p^{\prime })\)

$$\begin{aligned} \Pi ({\mathbf {M}}^{2},\mathbf {\ s}_{0},~q^{2})= & {} \int _{4m_{c}^{2}}^{s_{0}}ds \int _{4m_{c}^{2}}^{s_{0}^{\prime }}ds^{\prime }\rho _{\mathrm {D} }(s,s^{\prime },q^{2}) \nonumber \\&\times e^{-s/M_{1}^{2}}e^{-s^{\prime }/M_{2}^{2}}, \end{aligned}$$
(24)

where \({\mathbf {M}}^{2}=(M_{1}^{2},\ M_{2}^{2})\) and \({\mathbf {s}}_{0}=(s_{0},\ s_{0}^{\prime })\) are the Borel and continuum threshold parameters, respectively.

The obtained sum rule has to be used to determine the couplings \( g_{Z_{c}\eta _{c1}\pi }\) and \(g_{Z_{c}\eta _{c2}\pi }\). A possible way to find them is to get the second sum rule from the first one by applying the operators \(d/d(-1/M_{1}^{2})\) and/or \(d/d(-1/M_{2}^{2})\). But in the present work we choose the alternative approach and use iteratively the master sum rule to extract both \(g_{Z_{c}\eta _{c1}\pi }\) and \(g_{Z_{c}\eta _{c2}\pi }\) . To this end, we fix the continuum threshold parameter \(\sqrt{s_{0}^{\prime }}\) which corresponds to the \(\eta _{c}\) channel just below the mass of the first radially excited state \(\eta _{c}(2S)\). By this manipulation we include \(\eta _{c}(2S)\) into the continuum and obtain the sum rule for the strong coupling of the ground-state meson \(\eta _{c}(1S).\) The physical side of the sum rule (22) at this stage contains only the ground-state term and depends on the coupling \(g_{Z_{c}\eta _{c1}\pi }\). This sum rule can be easily solved to evaluate the unknown parameter \( g_{Z_{c}\eta _{c1}\pi }\)

$$\begin{aligned} g_{Z_{c}\eta _{c1}\pi }({\mathbf {M}}^{2},\mathbf {\ s}_{0}^{(1)},~q^{2})=\frac{ \Pi ({\mathbf {M}}^{2},\mathbf {\ s} _{0}^{(1)},~q^{2})e^{m/M_{1}^{2}}e^{m_{1}^{2}/M_{2}^{2}}}{A_{1}},\nonumber \\ \end{aligned}$$
(25)

where

$$\begin{aligned} A_{1}=\frac{mfm_{1}^{2}f_{1}\mu _{\pi }f_{\pi }}{4m_{c}(q^{2}-m_{\pi }^{2})} ( m^{2}+m_{1}^{2}-q^{2}), \end{aligned}$$

and \({\mathbf {s}}_{0}^{(1)}=(s_{0},\ s_{0}^{\prime }\simeq m_{2}^{2}).\)

At the next step we fix the continuum threshold \(\sqrt{s_{0}^{\prime }}\) at \( m_{2}+(0.5-0.8)\ \text {GeV}\) and use the sum rule that now contains the ground and first radially excited states. The QCD side of this sum rule is given by the expression \(\Pi ({\mathbf {M}}^{2},\mathbf {\ s}_{0}^{(2)},~q^{2})\) with \({\mathbf {s}}_{0}^{(2)}=(s_{0},\ s_{0}^{\prime }\simeq [m_{2}+(0.5-0.8)]^{2})\). By substituting the obtained expression for \( g_{Z_{c}\eta _{c1}\pi }\) into this sum rule it is not difficult to evaluate the second coupling \(g_{Z_{c}\eta _{c2}\pi }\).

The couplings depend on the Borel and continuum threshold parameters and, at the same time, are functions of \(q^{2}\). In what follows we omit their dependence on the parameters, replace \(q^{2}=-Q^{2}\) and denote the obtained couplings as \(g_{Z_{c}\eta _{c1}\pi }(Q^{2})\) and \(g_{Z_{c}\eta _{c2}\pi }(Q^{2})\). For calculation of the decay width we need value of the strong couplings at the pion’s mass-shell, i.e. at \(q^{2}=m_{_{\pi }}^{2}\), which is not accessible for the sum rule calculations. The standard way to avoid this problem is to introduce a fit functions \(F_{1(2)}(Q^{2})\) that for the momenta \(Q^{2}>0\) leads to the same predictions as the sum rules, but can be readily extrapolated to the region of \(Q^{2}<0\). Let us emphasize that values of the fit functions at the mass-shell are the strong couplings \( g_{Z_{c}\eta _{c1}\pi }\) and \(g_{Z_{c}\eta _{c2}\pi }\) to be utilized in calculations.

Expressions for \(g_{Z_{c}\eta _{c1}\pi }(Q^{2})\) and \(g_{Z_{c}\eta _{c2}\pi }(Q^{2})\) depend on various constants, such as the masses and decay constants of the final-state mesons. The values of these parameters are collected in Table 1. Additionally, there are parameters \( {\mathbf {M}}^{2}\) and \({\mathbf {s}}_{0}\) which should also be fixed to carry out numerical analysis. The requirements imposed on these auxiliary parameters have been discussed above and are standard for all sum rule computations. The regions for \(M_{1}^{2}\) and \(s_{0}\) which correspond to the tetraquark \( Z_{c}\) coincide with the working windows for these parameters fixed in the mass calculations \(M_{1}^{2}\in [4,\ 6]\ \text {GeV}^{2},\ s_{0}\in [19,\ 21]\ \text {GeV}^{2}\). The Borel and continuum threshold parameters \(M_{2}^{2},\ s_{0}^{\prime }\) in Eq. (25) are chosen as

$$\begin{aligned} M_{2}^{2}\in [3,\ 4]~\text {GeV}^{2},\quad s_{0}^{\prime }=13~\text {GeV}^{2}, \end{aligned}$$
(26)

whereas in the sum rule for the second coupling \(g_{Z_{c}\eta _{c2}\pi }(Q^{2})\) we employ

$$\begin{aligned} M_{2}^{2}\in [3,\ 4]~\text {GeV}^{2},\quad s_{0}^{\prime }\in [17,\ 19]~\text {GeV}^{2}. \end{aligned}$$
(27)

As it has been emphasized above to evaluate the strong couplings at the mass-shell \(Q^{2}=-m_{\pi }^{2}\) we need to determine the fit functions. To this end, we employ the following functions

$$\begin{aligned} F_{i}(Q^{2})=F_{0}^{i}\mathrm {\exp }\left[ c_{1i}\frac{Q^{2}}{m^{2}} +c_{2i}\left( \frac{Q^{2}}{m^{2}}\right) ^{2}\right] , \end{aligned}$$
(28)

where \(F_{0}^{i}\), \(c_{1i}\) and \(c_{2i}\) are fitting parameters. The performed analysis allows us to find the parameters as \(F_{0}^{1}=0.49\ \text {GeV}^{-1}\), \(c_{11}=27.64\) and \(c_{12}=-34.66\). Another set reads \( F_{0}^{2}=0.39\ \text {GeV}^{-1}\), \(c_{21}=28.13\) and \(c_{22}=-35.24\).

At the mass-shell the strong couplings are equal to

$$\begin{aligned} g_{Z_{c}\eta _{c1}\pi }(-m_{\pi }^{2})= & {} (0.47\pm 0.06)\ \text {GeV}^{-1}, \nonumber \\ g_{Z_{c}\eta _{c2}\pi }(-m_{\pi }^{2})= & {} (0.38\pm 0.05)\ \text {GeV}^{-1}. \end{aligned}$$
(29)

The widths of the decays \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\) and \( Z_{c}\rightarrow \eta _{c}(2S)\pi ^{-}\) can be found by means of the formula

$$\begin{aligned}&\Gamma \left[ Z_{c}\rightarrow \eta _{c}(\mathrm {I}S)\pi ^{-}\right] = \frac{g_{Z_{c}\eta _{ci}\pi }^{2}m_{i}^{2}}{8\pi }\lambda \left( m,m_{i},m_{\pi }\right) \nonumber \\&\quad \times \left[ 1+\frac{\lambda ^{2}\left( m,m_{i},m_{\pi }\right) }{ m_{i}^{2}}\right] ,\quad \mathrm {I}\equiv i=1,2 \end{aligned}$$
(30)

where

$$\begin{aligned} \lambda \left( a,b,c\right) =\frac{1}{2a}\sqrt{a^{4}+b^{4}+c^{4}-2\left( a^{2}b^{2}+a^{2}c^{2}+b^{2}c^{2}\right) }. \end{aligned}$$

For the decay \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\) one has to set \(g_{Z_{c}\eta _{ci}\pi }\rightarrow g_{Z_{c}\eta _{c1}\pi }\) and \(m_{i}\rightarrow m_{1}\), whereas in the case of \(Z_{c}\rightarrow \eta _{c}(2S)\pi ^{-}\) quantities with subscript 2 have to be used.

Using the strong couplings given by Eqs. (29) and (30) it is not difficult to evaluate the partial widths of the decay channels

$$\begin{aligned} \Gamma \left[ Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\right]= & {} (81\pm 17)\ \text {MeV}, \nonumber \\ \Gamma \left[ Z_{c}\rightarrow \eta _{c}(2S)\pi ^{-}\right]= & {} (32\pm 7)~\text {MeV}, \end{aligned}$$
(31)

which are main results of this section.

Table 1 Parameters of the mesons produced in the decays of the resonance \( Z_{c}\)

4 Decay \(Z_{c}\rightarrow D^{0}D^{-}\)

This section is devoted to investigation of the process \(Z_{c}\rightarrow D^{0}D^{-}\), which is the S-wave decay to the two open-charm pseudoscalar mesons. Our starting point is the three-point correlation function

$$\begin{aligned} {\widetilde{\Pi }}(p,p^{\prime })= & {} i^{2}\int d^{4}xd^{4}ye^{-ip\cdot x}e^{ip^{\prime }\cdot y} \nonumber \\&\times \langle 0|{\mathcal {T}}\{J^{D}(y)J^{D^{0}}(0)J^{\dagger }(x)\}|0\rangle , \end{aligned}$$
(32)

where

$$\begin{aligned} \ J^{D}(0)={\overline{c}}_{r}(y)i\gamma _{5}d_{r}(y),\ J^{D^{0}}(y)={\overline{u}}_{s}(0)i\gamma _{5}c_{s}(0), \end{aligned}$$
(33)

are the interpolating currents for the pseudoscalar mesons \(D^{-}\) and \( D^{0} \), respectively.

The correlation function \(\Pi (p,p^{\prime })\) expressed using parameters of the mesons \(D^{0}\) and \(D^{-}\) and tetraquark \(Z_{c}\) has the form

$$\begin{aligned}&{\widetilde{\Pi }}^{\mathrm {Phys}}(p,p^{\prime })=\frac{\langle 0|J^{D}|D^{-}\left( p^{\prime }\right) \rangle }{p^{\prime 2}-m_{D}^{2}} \frac{\langle 0|J^{D^{0}}|D^{0}\left( q\right) \rangle }{q^{2}-m_{D^{0}}^{2}} \nonumber \\&\quad \times \frac{\langle D^{-}\left( p^{\prime }\right) D^{0}(q)|Z_{c}(p)\rangle \langle Z_{c}(p)|J^{\dagger }|0\rangle }{p^{2}-m^{2} }+\cdots , \end{aligned}$$
(34)

where \(m_{D}\) and \(m_{D^{0}}\) are masses of the mesons \(D^{-}\) and \(D^{0}\), respectively. Contribution of the higher resonances and continuum states, as usual, are shown by dots.

We continue by utilizing the matrix elements

$$\begin{aligned} \langle 0|J^{D}|D^{-}\left( p^{\prime }\right) \rangle= & {} \frac{ f_{D}m_{D}^{2}}{m_{c}} \nonumber \\ \langle 0|J^{D^{0}}|D^{0}\left( q\right) \rangle= & {} \frac{ f_{D^{0}}m_{D^{0}}^{2}}{m_{c}} \end{aligned}$$
(35)

and the vertex

$$\begin{aligned} \langle D^{-}\left( p^{\prime }\right) D^{0}(q)|Z_{c}(p)\rangle =g_{Z_{c}DD}(p\cdot p^{\prime }). \end{aligned}$$
(36)

Simple manipulations lead to:

$$\begin{aligned} {\widetilde{\Pi }}^{\mathrm {Phys}}(p,p^{\prime })= & {} \frac{ f_{D^{0}}m_{D^{0}}^{2}f_{D}m_{D}^{2}}{m_{c}^{2}( p^{\prime 2}-m_{D}^{2}) ( q^{2}-m_{D^{0}}^{2}) } \nonumber \\&\times \frac{mf}{p^{2}-m^{2}}(p\cdot p^{\prime })+\cdots . \end{aligned}$$
(37)

Because the Lorentz structure of \(\ {\widetilde{\Pi }}^{\mathrm {Phys} }(p,p^{\prime })\) is trivial and \(\sim I\), the invariant amplitude \( {\widetilde{\Pi }}^{\mathrm {Phys}}(p^{2},p^{\prime 2},q^{2})\) is given by the function from Eq. (37).

The same correlation function written down in terms of the quark propagators is

$$\begin{aligned}&{\widetilde{\Pi }}^{\mathrm {OPE}}(p,p^{\prime })=i^{2}\int d^{4}xd^{4}ye^{-ip\cdot x}e^{ip^{\prime }\cdot y}\epsilon ^{ijk}\epsilon ^{imn} \nonumber \\&\quad \times {\mathrm {Tr}}\left[ \gamma _{5}S_{d}^{rk}(y-x)\gamma _{5}{\widetilde{S}} _{c}^{sj}(-x)\gamma _{5}{\widetilde{S}}_{u}^{ns}(x)\gamma _{5}S_{c}^{mr}(x-y) \right] .\nonumber \\ \end{aligned}$$
(38)

The invariant amplitudes \(\widetilde{\Pi }^{\mathrm {OPE}}(p^{2},p^{\prime 2},q^{2})\) and \(\widetilde{\Pi }^{\mathrm {Phys}}(p^{2},p^{\prime 2},q^{2})\) equated to each other yield the required sum rule. Contributions of higher resonances and continuum states can be suppressed by applying the double Borel transformation, and subtracted in accordance with the quark-hadron duality assumption.

The final sum rule for the strong coupling can be recast to the traditional form

$$\begin{aligned} g_{Z_{c}DD}({\mathbf {M}}^{2},\mathbf {\ s}_{0},~q^{2})=\frac{{\widetilde{\Pi }}( {\mathbf {M}}^{2},\mathbf {\ s}_{0},~q^{2})e^{m/M_{1}^{2}}e^{m_{D}^{2}/M_{2}^{2}} }{B}, \nonumber \\ \end{aligned}$$
(39)

where

$$\begin{aligned} B=\frac{f_{D}m_{D}^{2}f_{D^{0}}m_{D^{0}}^{2}mf}{ 2m_{c}^{2}(q^{2}-m_{D^{0}}^{2})}( m^{2}+m_{D}^{2}-q^{2}). \end{aligned}$$

Here \({\widetilde{\Pi }}({\mathbf {M}}^{2},\mathbf {\ s}_{0},~q^{2})\) is the Borel-transformed and continuum-subtracted amplitude \({\widetilde{\Pi }}^{ \mathrm {OPE}}(p^{2},p^{\prime 2},q^{2})\) given by analogous to Eq. (39) formula:

$$\begin{aligned} {\widetilde{\Pi }}({\mathbf {M}}^{2},{\mathbf {s}}_{0},~q^{2})= & {} \int _{4m_{c}^{2}}^{s_{0}}ds\int _{m_{c}^{2}}^{s_{0}^{\prime }}ds^{\prime } {\widetilde{\rho }}_{\mathrm {D}}(s,s^{\prime },q^{2}) \nonumber \\&\times e^{-s/M_{1}^{2}}e^{-s^{\prime }/M_{2}^{2}}. \end{aligned}$$
(40)

The sum rule for the strong coupling \(g_{Z_{c}DD}\) depends on vacuum condensates, and contains also the masses and decay constants of the mesons \( D^{0}\) and \(D^{-}\), which are shown in Table 1. Constraints imposed on the auxiliary parameters \({\mathbf {M}}^{2}\) and \({\mathbf {s}}_{0}\) are similar to ones discussed above and universal for all sum rules computations.The parameters \(M_{1}^{2}\) and \(s_{0}\) coincide with the working regions for these parameters fixed in the mass calculations (14). The Borel and continuum threshold parameters \(M_{2}^{2},\ s_{0}^{\prime }\) in Eq. (40)

$$\begin{aligned} M_{2}^{2}\in [3,\ 6]\ \text {GeV}^{2},\ s_{0}^{\prime }\in [7,\ 9]\ \text {GeV}^{2}, \end{aligned}$$
(41)

and ones from Eq. (14) lead to stable results for the form factor \(g_{Z_{c}DD}({\mathbf {M}}^{2},\mathbf {\ s}_{0},~q^{2})\) at \(q^{2}<0\). In what follows we denote it \(g_{Z_{c}DD}(Q^{2})\) omitting dependence on \( \mathbf {(M}^{2},\mathbf {\ s}_{0})\) and introducing \(q^{2}=-Q^{2}\).

A sensitivity of the strong coupling \(g_{Z_{c}DD}(Q^{2})\) to the Borel parameters is demonstrated in Fig. 3, which reveals its residual dependence on \(M_{1}^{2}\) and \(M_{2}^{2}\). This dependence of \( g_{Z_{c}DD}(Q^{2})\) as well as its variations generated by the continuum threshold parameters are main sources of ambiguities in the sum rule computations.

Fig. 3
figure 3

The strong coupling \(g_{Z_{c}DD}(Q^{2})\) as a function of the Borel parameters \({\mathbf {M}}^{2}=(M_{1}^{2},\ M_{2}^{2})\) at the fixed \( (s_{0},s_{0}^{\prime })=(20,8)~ \text {GeV}^{2}\) and \(Q^{2}=5~\text {GeV} ^{2}\)

The width of the decay \(Z_{c}\rightarrow D^{0}D^{-}\) depend on the strong coupling at \(D^{0}\) meson’s mass shell. In other words, we need \( g_{Z_{c}DD}(-m_{D^{0}}^{2})\) which cannot be accessed by direct sum rule computations. Therefore,we use the fit function \({\widetilde{F}}(Q^{2})\) that for the momenta \(Q^{2}>0\) coincides with the sum rule results, and can be easily extrapolated to the region of \(Q^{2}<0\). The function (28) with the parameters \({\widetilde{F}}_{0}=0.44\ \text {GeV}^{-1}\), \(\widetilde{ c}_{1}=2.38\) and \({\widetilde{c}}_{2}=-1.61\) meets these requirements. In Fig. 4 we plot \({\widetilde{F}}(Q^{2})\) and the sum rule results for \(g_{Z_{c}DD}(Q^{2})\) demonstrating a very nice agreement between them.

At the mass shell \(Q^{2}=-m_{D^{0}}^{2}\) the strong coupling is

$$\begin{aligned} g_{Z_{c}DD}(-m_{D^{0}}^{2})=(0.25\pm 0.05)\ \text {GeV}^{-1}. \end{aligned}$$
(42)

The width of the decay \(Z_{c}\rightarrow D^{0}D^{-}\) is calculated employing the expression

$$\begin{aligned} \Gamma [Z_{c}\rightarrow D^{0}D^{-}]=\frac{g_{Z_{c}DD}^{2}m_{D}^{2}}{ 8\pi }\lambda \left( 1+\frac{\lambda ^{2}}{m_{D}^{2}}\right) , \end{aligned}$$
(43)

where \(\lambda =\lambda \left( m,m_{D^{0}},m_{D}\right) .\)

The partial width of this decay reads:

$$\begin{aligned} \Gamma [Z_{c}\rightarrow D^{0}D^{-}]=(19\pm 5)~\text {MeV}. \end{aligned}$$
(44)

It will be used below to estimate the total width of the tetraquark \(Z_{c}\).

Fig. 4
figure 4

The sum rule predictions and fit function for the strong coupling \( g_{Z_{c}DD}(Q^{2})\). The star marks the point \(Q^{2}=-m_{D^{0}}^{2}\)

5 Decay \(Z_{c}\rightarrow J/\psi \rho ^{-}\)

The scalar tetraquark \(Z_{c}\) in S-wave can also decay to the final state \( J/\psi \rho ^{-}\). In the QCD light-cone sum rule approach this decay can be explored through the correlation function

$$\begin{aligned} \Pi _{\mu }(p,q)=i\int d^{4}xe^{ipx}\langle \rho (q)|{\mathcal {T}}\{J_{\mu }^{J/\psi }(x)J^{\dagger }(0)\}|0\rangle , \end{aligned}$$
(45)

where

$$\begin{aligned} J_{\mu }^{J/\psi }(x)={\overline{c}}_{i}(x)i\gamma _{\mu }c_{i}(x), \end{aligned}$$
(46)

is the interpolating current for the vector meson \(J/\psi \).

The correlation function \(\Pi _{\mu }^{\mathrm {Phys}}(p,q)\) in terms of the physical parameters of the tetraquark \(Z_{c}\), and of the mesons \(J/\psi \) and \(\rho \) has the following form

$$\begin{aligned} \Pi _{\mu }^{\mathrm {Phys}}(p,q)= & {} \frac{\langle 0|J_{\mu }^{J/\psi }|J/\psi \left( p\right) \rangle }{p^{2}-m_{J/\psi }^{2}}\langle J/\psi \left( p\right) \rho (q)|Z_{c}(p^{\prime })\rangle \nonumber \\&\times \frac{\langle Z_{c}(p^{\prime })|J^{\dagger }|0\rangle }{p^{\prime 2}-m^{2}}+\cdots , \end{aligned}$$
(47)

where \(m_{J/\psi }\) is the mass of the meson \(J/\psi .\) In Eq. (47) by the dots we denote contribution of the higher resonances and continuum states. Here \(p^{\prime }=p+q\) is the momentum of the tetraquark \( Z_{c}\), where p and q are the momenta of the \(J/\psi \) and \(\rho \) mesons, respectively.

Further simplification of \(\Pi _{\mu }^{\mathrm {Phys}}(p,q)\) can be achieved by utilizing explicit expressions of the matrix elements \(\langle 0|J_{\mu }^{J/\psi }|J/\psi \left( p\right) \rangle ,\) \(\langle Z_{c}(p^{\prime })|J^{\dagger }|0\rangle \), and of the vertex \(Z_{c}(p^{\prime })J/\psi \left( p\right) \rho (q)\). The matrix element of the tetraquark \(Z_{c}\) is given by Eq. (5), whereas for the meson \(J/\psi \left( p\right) \) we can use

$$\begin{aligned} \langle 0|J_{\mu }^{J/\psi }|J/\psi \left( p\right) \rangle =m_{J/\psi }f_{J/\psi }\varepsilon _{\mu }, \end{aligned}$$
(48)

where \(f_{J/\psi }\) and \(\varepsilon _{\mu }\) are its decay constant and polarization vector, respectively. We also model the three-state vertex as

$$\begin{aligned} \langle J/\psi \left( p\right) \rho (q)|Z_{c}(p^{\prime })\rangle= & {} g_{Z_{c}J/\psi \rho }\left[ (p\cdot \varepsilon ^{\prime })(q\cdot \varepsilon ^{*})\right. \nonumber \\&-\left. (p\cdot q)(\varepsilon ^{*}\cdot \varepsilon ^{\prime })\right] , \end{aligned}$$
(49)

with \(\varepsilon ^{\prime }\) being the polarization vector of the \(\rho \) -meson. Then \(\Pi _{\mu }^{\mathrm{Phys}}(p,q)\) takes the form:

$$\begin{aligned}&\Pi _{\mu }^{\mathrm {Phys}}(p,q)=g_{Z_{c}J/\psi \rho }\frac{m_{J/\psi }f_{J/\psi }}{p^{2}-m_{J/\psi }^{2}}\frac{mf}{p^{\prime 2}-m^{2}} \nonumber \\&\quad \times \left[ \frac{1}{2}\left( m^{2}-m_{J/\psi }^{2}-q^{2}\right) \varepsilon _{\mu }^{\prime }-p\cdot \varepsilon ^{\prime }q_{\mu }\right] +\cdots . \nonumber \\&\end{aligned}$$
(50)

It contains different Lorentz structures \(\sim \varepsilon _{\mu }^{\prime }\) and \(q_{\mu }\). One of them should be chosen to fix the invariant amplitude and carry out sum rule analysis. We choose the structure \(\sim \varepsilon _{\mu }^{\prime }\) and denote the corresponding invariant amplitude as \(\Pi ^{\mathrm {Phys}}(p^{2},q^{2}).\)

The second component of the sum rule is the correlation function \(\Pi _{\mu }(p,q)\) computed using quark propagators. For \(\Pi _{\mu }^{\mathrm {OPE} }(p,q)\) we obtain

$$\begin{aligned}&\Pi _{\mu }^{\mathrm {OPE}}(p,q)=i^{2}\int d^{4}xe^{ipx}\epsilon {\widetilde{\epsilon }}\left[ \gamma _{5}{\widetilde{S}}_{c}^{aj}(x){}\gamma _{\mu }\right. \nonumber \\&\quad \left. \times {\widetilde{S}}_{c}^{ma}(-x){}\gamma _{5}\phantom {\gamma _{5}{\widetilde{S}}_{c}^{aj}(x){}\gamma _{\mu }}\right] _{\alpha \beta }\langle \rho (q)|{\overline{d}}_{\alpha }^{k}(0)u_{\beta }^{n}(0)|0\rangle , \end{aligned}$$
(51)

where \(\alpha \) and \(\beta \) are the spinor indexes.

The expression for \(\Pi _{\mu }^{\mathrm {OPE}}(p,q)\) can be written down in a more detailed form. For these purposes, we first expand the local operator \({\overline{d}}_{\alpha }^{a}u_{\beta }^{d}\) by means of the formula

$$\begin{aligned} {\overline{d}}_{\alpha }^{a}u_{\beta }^{d}\rightarrow \frac{1}{4}\Gamma _{\beta \alpha }^{j}\left( {\overline{d}}^{a}\Gamma ^{j}u^{d}\right) , \end{aligned}$$
(52)

with \(\Gamma ^{j}\) being the full set of Dirac matrixes

$$\begin{aligned} \Gamma ^{j}=\mathbf {1,\ }\gamma _{5},\ \gamma _{\lambda },\ i\gamma _{5}\gamma _{\lambda },\ \sigma _{\lambda \rho }/\sqrt{2}. \end{aligned}$$

Applying the projector onto a color-singlet state \(\delta ^{ad}/3\) we get

$$\begin{aligned} {\overline{d}}_{\alpha }^{a}u_{\beta }^{d}\rightarrow \frac{1}{12}\Gamma _{\beta \alpha }^{j}\delta ^{ad}\left( {\overline{d}}\Gamma ^{j}u\right) , \end{aligned}$$
(53)

where \({\overline{d}}\Gamma ^{j}u\) are the color-singlet local operators. Substituting the last expression into Eq. (51) we see that the correlation function \(\Pi _{\mu }^{\mathrm {OPE}}(p,q)\) depends on the \(\rho \) -meson’s two-particle local matrix elements. Some of them does not depend on the \(\rho \)-meson momentum,

$$\begin{aligned} \langle 0|{\overline{u}}\gamma _{\mu }d|\rho (q,\lambda )\rangle =\epsilon _{\mu }^{(\lambda )}f_{\rho }m_{\rho }, \end{aligned}$$
(54)

whereas others contain momentum factor

$$\begin{aligned} \langle 0|{\overline{u}}\sigma _{\mu \nu }d|\rho (q,\lambda )\rangle =if_{\rho }^{T}(\epsilon _{\mu }^{(\lambda )}q_{\nu }-\epsilon _{\nu }^{(\lambda )}q_{\mu }). \end{aligned}$$
(55)

There are also three-particle matrix elements that contribute to the correlation function \(\Pi _{\mu }^{\mathrm {OPE}}(p,q)\). They appear due to insertion of gluon field strength tensor G from the c-quark propagators into the local operators \({\overline{d}}\Gamma ^{j}u.\) The \(\rho \)-meson three-particle local matrix element

$$\begin{aligned} \langle 0|{\overline{u}}g{\widetilde{G}}_{\mu \nu }\gamma _{\nu }\gamma _{5}d|\rho (q,\lambda )\rangle =f_{\rho }m_{\rho }^{3}\epsilon _{\mu }^{(\lambda )}\zeta _{4\rho }, \end{aligned}$$
(56)

is a q free quantity. But other matrix elements depend on the \(\rho \) -meson momentum

$$\begin{aligned} \langle 0|{\overline{u}}gG_{\mu \nu }d|\rho (q,\lambda )\rangle= & {} if_{\rho }^{T}m_{\rho }^{3}\zeta _{4}^{T}(\epsilon _{\mu }^{(\lambda )}q_{\nu }-\epsilon _{\nu }^{(\lambda )}q_{\mu }), \nonumber \\ \langle 0|{\overline{u}}g{\widetilde{G}}_{\mu \nu }i\gamma _{5}d|\rho (q,\lambda )\rangle= & {} if_{\rho }^{T}m_{\rho }^{3}\widetilde{\zeta }_{4}^{T}(\epsilon _{\mu }^{(\lambda )}q_{\nu }-\epsilon _{\nu }^{(\lambda )}q_{\mu }). \nonumber \\&\end{aligned}$$
(57)

As a result, the correlation function contains only local matrix elements of the \(\rho \)-meson and depends on the momenta p and q. This is general feature of QCD sum rules on the light-cone with a tetraquark and two conventional mesons. Indeed, because a tetraquark contains four quarks, after contracting two quark fields from its interpolating current with relevant quarks from the interpolating current of a meson one gets a local operator sandwiched between the vacuum and a second meson. The variety of such local operators gives rise to different local matrix elements of the meson rather that to its distribution amplitudes. Then the four-momentum conservation in the tetraquark-meson-meson vertex requires setting \(q=0\) ( for details, see Ref. [26]). In the standard LCSR method the choice \(q=0\) is known as the soft-meson approximation [39]. At vertices composed of conventional mesons in general \( q\ne 0\), and only in the soft-meson approximation one equates q to zero, whereas the tetraquark-meson-meson vertex can be analyzed in the context of the LCSR method only if \(q=0\). An important observation made in Ref. [39] is that the soft-meson approximation and full LCSR treatment of the conventional mesons’ vertices lead to results which numerically are very close to each other. It is worth to note that the full version of the sum rules on the light-cone is applicable to tetraquark-tetraquark-meson vertices [54].

After substituting all aforementioned matrix elements into the expression of the correlation function and performing the summation over color indices we fix the local matrix elements of the \(\rho \) meson that survive the soft limit. It turns out that in the \(q\rightarrow 0\) limit only the matrix elements (54) and (56) contribute to the invariant amplitude \(\Pi ^{\mathrm {OPE}}(p^{2})\) [i.e. to \(\Pi ^{\mathrm {OPE} }(p^{2},\ 0)\)]. These matrix elements depend on the mass and decay constant of the \(\rho \)-meson \(m_{\rho }\), \(f_{\rho }\), and on \(\zeta _{4\rho }\) which normalizes the twist-4 matrix element of the \(\rho \)-meson [55]. The parameter \(\zeta _{4\rho }\) was evaluated in the context of QCD sum rule approach at the renormalization scale \(\mu =1\,\,{\mathrm {GeV }}\) in Ref. [56] and is equal to \(\zeta _{4\rho }=0.07\pm 0.03 \).

The Borel transform of the invariant amplitude \(\Pi ^{\mathrm {OPE}}(p^{2})\) is given by the expression

$$\begin{aligned} \Pi ^{\mathrm {OPE}}(M^{2})=\int _{4m_{c}^{2}}^{\infty }ds{\widetilde{\rho }}^{ \mathrm {OPE}}(s)e^{-s/M^{2}}, \end{aligned}$$
(58)

where \({\widetilde{\rho }}^{\mathrm {OPE}}(s)\) is the corresponding spectral density. In the present work we calculate \(\widetilde{\rho }^{\mathrm {OPE} }(s)\) by taking into account contribution of the condensates up to dimension six. The spectral density has both the perturbative and nonperturbative components

$$\begin{aligned} {\widetilde{\rho }}^{\mathrm {OPE}}(s)=\rho ^{\mathrm {pert.}}(s)+\rho ^{\mathrm { n.-pert.}}(s). \end{aligned}$$
(59)

After some computations for \(\rho ^{\mathrm {pert.}}(s)\) we get

$$\begin{aligned} \rho ^{\mathrm {pert.}}(s)=\frac{f_{\rho }m_{\rho }(s+2m_{c}^{2})\sqrt{ s(s-4m_{c}^{2})}}{24\pi ^{2}s}. \end{aligned}$$
(60)

The nonperturbative part of the spectral density \(\rho ^{\mathrm {n.-pert.} }(s)\) contains terms proportional to the gluon condensates \(\langle \alpha _{s}G^{2}/\pi \rangle \), \(\langle \alpha _{s}G^{2}/\pi \rangle ^{2}\) and \( \langle g_{s}^{3}G^{3}\rangle \): Here we do not provide their explicit expressions. The twist-4 contribution to \(\Pi ^{\mathrm {OPE}}(M^{2})\) reads

$$\begin{aligned} \Pi ^{\mathrm {OPE(tw4)}}(M^{2})=\frac{f_{\rho }m_{\rho }^{3}\zeta _{4\rho }m_{c}^{2}}{8\pi }\int _{0}^{1}d\alpha \frac{e^{-m_{c}^{2}/M^{2}a(1-a)}}{ a(1-a)}. \nonumber \\ \end{aligned}$$
(61)

To derive the expression for the strong coupling \(g_{Z_{c}J/\psi \rho }\) the soft-meson approximation should be applied to the phenomenological side of the sum rule as well. Because in the soft limit \(p^{2}=p^{\prime 2}\), we have to perform the Borel transformation of \(\ \Pi ^{\mathrm {Phys}}(p^{2},0)\) over the variable \(p^{2}\) and carry out calculations with one parameter \( M^{2}\). To this end, we first transform \(\Pi ^{\mathrm {Phys}}(p^{2},0)\) in accordance with the prescription

$$\begin{aligned} \frac{1}{(p^{2}-m_{J/\psi }^{2})\left( p^{\prime 2}-m^{2}\right) } \rightarrow \frac{1}{(p^{2}-{\widetilde{m}}^{2})^{2}}, \end{aligned}$$
(62)

where \({\widetilde{m}}^{2}=( m^{2}+m_{J/\psi }^{2}) /2\), and instead of two terms with different poles get the double pole term. By equating the physical and QCD sides and performing required manipulations we get

$$\begin{aligned}&\left( g_{Z_{c}J/\psi \rho }m_{J/\psi }f_{J/\psi }mf\frac{m^{2}-m_{J/\psi }^{2}}{2}+AM^{2}\right) \nonumber \\&\quad \times \frac{e^{-{\widetilde{m}}^{2}/M^{2}}}{M^{2}}+\cdots =\Pi ^{\mathrm {OPE }}(M^{2}). \end{aligned}$$
(63)

The equality given by Eq. (63) is the master expression which can be used to extract sum rule for the coupling \(g_{Z_{c}J/\psi \rho } \). It contain the term corresponding to the decay of the ground-state tetraquark \(Z_{c}\) and conventional contributions of higher resonances and continuum states suppressed due to the Borel transformation; the latter is denoted in Eq. (63) by the dots. But in the soft limit there are also terms \(\sim A\) in the physical side which remain unsuppressed even after the Borel transformation. They describe transition from the excited states of the tetraquark \(Z_{c}\) to the mesons \(J/\psi \rho ^{-}\). Of course, to obtain the final formula all contributions appearing as the contamination should be removed from the physical side of the sum rule. The situation with the ordinary suppressed terms is clear: they can be subtracted from the correlation function \(\Pi ^{\mathrm {OPE}}(M^{2})\) using assumption on the quark-hadron duality. As a result the correlation function acquires a dependence on the continuum threshold parameter \(s_{0}\), i.e., becomes equal to \(\Pi ^{\mathrm {OPE}}(M^{2},s_{0})\). The treatment of the terms \(\sim A\) requires some additional manipulations; they can be removed by applying the operator [49]

$$\begin{aligned} {\mathcal {P}}(M^{2},{\widetilde{m}}^{2})=\left( 1-M^{2}\frac{d}{dM^{2}}\right) M^{2}e^{{\widetilde{m}}^{2}/M^{2}}, \end{aligned}$$
(64)

to both sides of Eq. (63). Then the sum rule for the strong coupling reads:

$$\begin{aligned} g_{Z_{c}J/\psi \rho }= & {} \frac{2}{m_{J/\psi }f_{J/\psi }mf(m^{2}-m_{J/\psi }^{2})} \nonumber \\&\times {\mathcal {P}}(M^{2},{\widetilde{m}}^{2})\Pi ^{\mathrm {OPE} }(M^{2},s_{0}). \end{aligned}$$
(65)

The width of the decay \(Z_{c}\rightarrow J/\psi \rho ^{-}\) can be calculated using the formula

$$\begin{aligned}&\Gamma \left( Z_{c}\rightarrow J/\psi \rho ^{-}\right) =\frac{ g_{Z_{c}J/\psi \rho }^{2}m_{\rho }^{2}}{8\pi }\lambda \left( m,\ m_{J/\psi },m_{\rho }\right) \nonumber \\&\quad \times \left[ 3+\frac{2\lambda ^{2}\left( m,\ m_{J/\psi },m_{\rho }\right) }{m_{\rho }^{2}}\right] . \end{aligned}$$
(66)

In the sum rule (65) for \(M^{2}\) and \(s_{0}\) we use the working regions given by Eq. (14). For the strong coupling \( g_{Z_{c}J/\psi \rho }\) we get

$$\begin{aligned} g_{Z_{c}J/\psi \rho }=(0.56\pm 0.07)\ \text {GeV}^{-1}. \end{aligned}$$
(67)

Then the width of the decay \(Z_{c}\rightarrow J/\psi \rho ^{-}\) is

$$\begin{aligned} \Gamma \left[ Z_{c}\rightarrow J/\psi \rho ^{-}\right] =(15\pm 3)\ \mathrm { MeV.} \end{aligned}$$
(68)

       In accordance with our investigation, the total width of the resonance \(Z_{c}\) saturated by the dominant decay modes \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\), \(\eta _{c}(2S)\pi ^{-}\), \(D^{0}D^{-}\) and \( Z_{c}\rightarrow \) \(J/\psi \rho ^{-}\) is

$$\begin{aligned} \Gamma =(147\pm 19)\ \mathrm {MeV.} \end{aligned}$$
(69)

This is the second parameter of the resonance \(Z_{c}\) to be compared with the LHCb data; our result for the total with of \(Z_{c}\) is in excellent agreement with existing data \(\Gamma =152_{-68}^{+83}\ \text {MeV}\).

6 Analysis and concluding remarks

We have performed quantitative analysis of the newly observed resonance \( Z_{c}\) by calculating its spectroscopic parameters and total width. In computations we have used different QCD sum rule approaches. Thus, the mass and coupling of \(Z_{c}\) have been evaluated by means of the two-point sum rule method, whereas its decay channels have been analyzed using the three-point and light-cone sum rule techniques.

We have calculated the spectroscopic parameters of the tetraquark \(Z_{c}\) using the zero-width single-pole approximation. But the interpolating current (2) couples also to the two-meson continuum \(\eta _{c}(1S)\pi ^{-}\), \(\eta _{c}(2S)\pi ^{-}\), \(J/\psi \rho ^{-}\), \(D^{0}D^{-}\) and \(D^{*0}D^{*-}\) which can modify the results for m and f obtained in the present work. Effects of the two-meson continuum change the zero-width approximation (4) and lead to the following corrections [46]

$$\begin{aligned} \lambda ^{2}e^{-m^{2}/M^{2}}\rightarrow \lambda ^{2}\int _{{\mathcal {M}} ^{2}}^{s_{0}}dsW(s)e^{-s/M^{2}} \end{aligned}$$
(70)

and

$$\begin{aligned} \lambda ^{2}m^{2}e^{-m^{2}/M^{2}}\rightarrow \lambda ^{2}\int _{{\mathcal {M}} ^{2}}^{s_{0}}dsW(s)se^{-s/M^{2}}, \end{aligned}$$
(71)

where \(\lambda =mf\) and \(\mathcal {M=}m_{D^{*0}}+m_{D^{*-}}\). In Eqs. (70) and (71) we have introduced the weight function

$$\begin{aligned} W(s)=\frac{1}{\pi }\frac{m\Gamma (s)}{\left( s-m^{2}\right) ^{2}+m^{2}\Gamma ^{2}(s)} \end{aligned}$$
(72)

where

$$\begin{aligned} \Gamma (s)=\Gamma \frac{m}{s}\sqrt{\frac{s-{\mathcal {M}}^{2}}{m^{2}-{\mathcal {M}} ^{2}}}. \end{aligned}$$
(73)

Utilizing the central values of the m and \(\Gamma \), as well as \(M^{2}=5\ \text {GeV}^{2}\) and \(s_{0}=20\ \text {GeV}^{2}\), it is not difficult to find that

$$\begin{aligned} \lambda ^{2}\rightarrow 0.903\lambda ^{2}\rightarrow (0.95f)^{2}m^{2}, \end{aligned}$$
(74)

and

$$\begin{aligned} \lambda ^{2}m^{2}\rightarrow 0.91\lambda ^{2}m^{2}\rightarrow (0.955f)^{2}m^{4}. \end{aligned}$$
(75)

As is seen the two-meson effects result in rescaling \(f\rightarrow 0.95f\) which changes it approximately by \(5\%\) relative to its central value. These effects are smaller than theoretical errors of the sum rule computations themselves.

We have saturated the total width of the resonance \(Z_{c}\) by its four dominant decay modes \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\), \(\eta _{c}(2S)\pi ^{-}\), \(D^{0}D^{-}\) and \(J/\psi \rho ^{-}\). To calculate partial widths of these decay channels we used two approaches in the framework of the QCD sum rule method. Thus the decays \(Z_{c}\rightarrow \eta _{c}(1S)\pi ^{-}\), \(Z_{c}\rightarrow \eta _{c}(2S)\pi ^{-}\) and \(Z_{c}\rightarrow D^{0}D^{-}\) have been studied by applying three-point sum rules, whereas the process \(Z_{c}\rightarrow J/\psi \rho ^{-}\) has been investigated using the LCSR method and soft-meson approximation. Predictions obtained for partial widths of these S-wave decay channels have been used to evaluate the total width of the resonance \(Z_{c}\).

Our results for the mass \(m=(4080~\pm 150)~\text {MeV}\) and total width \( \Gamma =(147\pm 19)\ \text {MeV}\) of the resonance \(Z_{c}\) are in a very nice agreement with experimental data of the LHCb Collaboration. This allows us to interpret the new charged resonance as the scalar diquark–antidiquark state with \(cd{\overline{c}}{\overline{u}}\) content and \(C\gamma _{5}\otimes \gamma _{5}C\) structure. It presumably belongs to one of the charged Z-resonance multiplets, axial-vector members of which are the tetraquarks \( Z_{c}^{\pm }(3900)\) and \(Z_{c}^{\pm }(4330)\), respectively. The charged resonances \(Z_{c}^{\pm }(4330)\) and \(Z_{c}^{\pm }(3900)\) were observed in the \(\psi ^{\prime }\pi ^{\pm }\) and \(J/\psi \pi ^{\pm }\) invariant mass distributions, i.e. they dominantly decay to these particles. The neutral resonance \(Z_{c}^{0}(3900)\) was discovered in the process \( e^{+}e^{-}\rightarrow \pi ^{0}\pi ^{0}J/\psi \). Because \(J/\psi \) and \(\psi ^{\prime }\) are vector mesons, and \(\psi ^{\prime }\) is the radial excitation of \(J/\psi \), it is natural to suggest that \(Z_{c}(4330)\) is the excited state of \(Z_{c}(3900)\). This suggestion was originally made in Ref. [14], and confirmed later by sum rule calculations. Then the resonance \(Z_{c}^{-}(4100)\) fixed in the \(\eta _{c}(1S)\pi ^{-}\) channel can be interpreted as a scalar counterpart of these axial-vector tetraquarks. It is also reasonable to assume that the neutral member of this family \( Z_{c}^{0}(4100)\) will be seen in the processes \(e^{+}e^{-}\rightarrow \pi ^{0}\pi ^{0}\eta _{c}(1S)\) with dominantly \(\pi ^{0}\pi ^{0}\) mesons rather than \(D{\overline{D}}\) ones at the final state.