Can the $\psi (4040)$ explain the peak associated with Y(4008)?

Abstract

We study the well-known resonance $\psi (4040)$, corresponding to a $3^{3}S_{1}$ charm–anticharm vector state $\psi (3S)$, within a QFT approach, in which the decay channels into DD, $D^{*}D$, $D^{*}D^{*}$, $D_{s}D_{s}$ and $D_{s}^{*}D_{s}$ are considered. The spectral function shows sizable deviations from a Breit–Wigner shape (an enhancement, mostly generated by $DD^{*}$ loops, occurs); moreover, besides the $c{\bar{c}}$ pole of $\psi (4040)$, a second dynamically generated broad pole at 4 GeV emerges. Naively, it is tempting to identify this new pole with the unconfirmed state Y(4008). Yet, this state was not seen in the reaction $e^{+}e^{-}\rightarrow \psi (4040)\rightarrow DD^{*}$, but in processes with $\pi ^{+}\pi ^{-}J/\psi$ in the final state. A detailed study shows a related but different mechanism: a broad peak at 4 GeV in the process $e^{+}e^{-}\rightarrow \psi (4040)\rightarrow DD^{*}\rightarrow \pi ^{+}\pi ^{-}J/\psi$ appears when $DD^{*}$ loops are considered. Its existence in this reaction is not necessarily connected to the existence of a dynamically generated pole, but the underlying mechanism – the strong coupling of $c{\bar{c}}$ to $DD^{*}$ loops – can generate both of them. Thus, the controversial state Y(4008) may not be a genuine resonance, but a peak generated by the $\psi (4040)$ and $D^{*}D$ loops with $\pi ^{+}\pi ^{-}J/\psi$ in the final state.

1 Introduction

The understanding of the nature and properties of hadronic states is a substantial challenge for both experimentalists and theorists. Remarkable progress in the field of charmonium spectroscopy was provided in the past decades: while various resonances emerge as conventional charmonium states (ordinary ${\bar{c}}c$ objects), the so-called X,  Y,  and Z states are candidates for exotic hadrons (such as molecules, hybrids, multi-quarks objects or glueballs; see Refs. [1,2,3,4] and refs. therein).

In this work, we shall concentrate on the vector sector in the energy region close to 4 GeV. Here, the well established charmonium vector state $\psi (4040)$ is listed in the Particle Data Group (PDG) [5] (it has $J^{PC}=1^{--}$ where, as usual, P refers to parity and C to charge conjugation). This resonance can be successfully interpreted as a charmonium state with $(n,L,S,J)=(3,0,1,1)$, where n is the principal number, L the angular momentum, S the spin and J the total spin); hence, the nonrelativistic spectroscopic notation reads n $^{2S+1}L_{J}=$ 3  $^{3}S_{1}$ (see e.g. Refs. [6,7,8,9,10,11] and refs. therein).

Very close to 4 GeV, the enigmatic (and not yet confirmed) resonance Y(4008) was also observed as a significant enhancement by the Belle Collaboration when measuring the cross section of $e^{+}e^{-}\rightarrow \pi ^{+}\pi ^{-}J/\psi$ via initial state radiation (ISR) technique [12] and later on confirmed by the same group [13]: its mass was determined as $4008\pm 40_{-28}^{+114}$ MeV and the decay width as $\Gamma =226\pm 44\pm 87$ MeV. Moreover, a broad Y(4008) was also found in the recent analysis of Ref. [14]. However, the statistic at Belle was pretty limited and Y(4008) could not be confirmed by subsequent experiments studying the same production process at BaBar [15] and BESIII [16], making its existence rather controversial. Nevertheless, several possible theoretical assignments on its nature have been suggested, including non-conventional scenarios as $D^{*}{\bar{D}}^{*}$ molecular state [17, 18] (see however also Ref. [19]), tetraquark [20, 21] or even an interference effect with background [22]. Moreover, in Refs. [23, 24] it was proposed to identify Y(4008) as $\psi (3S)$ charmonium state, but this assignment is not favoured, since, as mentioned above, $\psi (4040)$ is well described by a standard $\psi (3S)$ state. The unexplained status of the observed structure corresponding to Y(4008) makes it an interesting subject that deserves clarification, hence we aim to perform a detailed study of the nearby energy region.

To this end, we develop a quantum field theoretical effective model in which a single $c{\bar{c}}$ seed state, to be identified with $\psi (4040),$ couples to DD,  $DD^{*},$ and $D^{*}D^{*}.$ The immediate question is if we can describe both resonances $\psi (4040)$ and Y(4008) at the same time and within a unique setup. The idea that we test is somewhat reminiscent of studies in the light scalar sector, in which the state $a_{0}(980)$ can be seen as a companion pole of the predominantly $q{\bar{q}}$ state $a_{0}(1450)$ [25,26,27,28] as well as the light $\kappa$ state, now named $K_{0}^{*}(700),$ as a companion pole of the $K_{0}^{*}(1430)$ [29]. Quite interestingly, two poles emerge also in the study of the charmonium resonance $\psi (3770)$ [30, 31]. As we shall see, some similarities, but also some important differences, will emerge between those studies and the one that we are going to present.

As a first step of our analysis, we calculate the spectral function of $\psi (4040)$. As expected, it is peaked at about 4.04 GeV, but it cannot be approximated by a standard Breit–Wigner shape; most remarkably, it may develop an additional enhancement below 4 GeV (this is due to the strong coupling of the bare ${\bar{c}}c$ state to the $DD^{*}$ channel). Moreover, two poles appear in the complex plane: one corresponding to the peak in the spectral function (hence to $\psi (4040)$) and an additional companion pole, dynamically generated by meson–meson (mostly $DD^{*}$) quantum fluctuations. A large-$N_{c}$ study shows that $\psi (4040)$ behaves as a conventional quark-antiquark state while the enhancement does not fit into this standard picture.

At a first sight, it appears quite natural to assign the state Y(4008) to this additional dynamically generated pole. Yet, a closer inspection is necessary: the study of the decay chain in which Y(4008) was seen, $e^+e^- \rightarrow \psi (4040)\rightarrow DD^{*}\rightarrow \pi ^{+}\pi ^{-}J/\psi$ (the latter can occur via a light scalar state, most notably $f_{0}(500)$, but not only), which shows that a broad peak at about 3.9 GeV emerges for the cross-section $e^{+}e^{-}\rightarrow \pi ^{+}\pi ^{-}J/\psi$ (also when $e^{+}e^{-}$ comes from a previous ISR process, as observed in experiment). This is due to the fact that the loop contribution of $DD^{*}$ is peaked at about $m_{D} +m_{D^{*}}\simeq 3.9$ GeV. As we shall explain in detail later on, this contribution is multiplied by the modulus squared of the propagator of $\psi (4040)$, centered at 4.04 GeV and about 80 MeV large, hence a sizable overlap is present. As we shall show, the emergent peak at about 3.9 GeV is very far from a Breit–Wigner state, but is rather distorted. Strictly speaking, the very existence of an additional companion pole is not necessary for the emergence of this signal, but both phenomena arise from a strong coupling of the seed state to $DD^{*}$, hence it is rather natural that they both take place at the same time.

The paper is organized as follows: in Sect. 2 we introduce theoretical model, in particular the Lagrangians, the possible decays channels of $\psi (4040)$ resonance with corresponding theoretical expression for decay widths, loop function (hence, the propagator) and spectral function. Moreover, we show in detail the determination of the parameters of our model. In Sect. 3 we present our results. Summary and outlooks are presented in Sect. 4. Additional results for different parameter values are reported in the Appendices.

2 The model

In this section we present the theoretical model used to analyze the energy region close to 4 GeV. In our approach, only a single standard $c{\bar{c}}$ seed state corresponding to $\psi (4040)$ is included.

2.1 Theoretical framework

The resonance $\psi (4040)$ can be described by a relativistic interaction Lagrangian that couples it to its decay products [two pseudoscalar mesons (DD and $D_sD_s$), one vector and one pseudoscalar meson ($DD^{*}$ and $D^*_sD_s$), and two vector mesons ($D^{*}D^{*}$)]:

$$\begin{aligned} {\mathcal {L}}_{\psi (4040)}={\mathcal {L}}_{VPP}+{\mathcal {L}}_{VPV}+{\mathcal {L}} _{VVV} \end{aligned}$$

(1)

where

$$\begin{aligned} {\mathcal {L}}_{VPP}= & {} ig_{\psi DD}\psi _{\mu }\left[ \left( \partial ^{\mu } D^{+}\right) D^{-}+\left( \partial ^{\mu }D^{0}\right) {\bar{D}}^{0} \right. \nonumber \\&\left. +\left( \partial ^{\mu }D_{s}^{+}\right) D_{s}^{-}\right] +h.c. \end{aligned}$$

(2)

$$\begin{aligned} {\mathcal {L}}_{VPV}= & {} ig_{\psi D^{*}D}{\tilde{\psi }}_{\mu \nu }\left[ \partial ^{\mu }D^{*+\nu }D^{-}+\partial ^{\mu }D^{*0\nu }{\bar{D}}^{0}\right. \nonumber \\&\left. +\partial ^{\mu }D_{s}^{*+\nu }D_{s}^{-}\right] +h.c.{,} \end{aligned}$$

(3)

$$\begin{aligned} {\mathcal {L}}_{VVV}= & {} ig_{\psi D^{*}D^{*}}\left[ \psi _{\mu \nu }\left( D^{*+\mu }D^{*-\nu }+D^{*0\mu }{\bar{D}}^{*0\nu } \right. \right. \nonumber \\&\left. \left. +D_{s}^{*+\mu } D_{s}^{*-\nu }\right) \right] +h.c.{} \end{aligned}$$

(4)

The quantities $g_{\psi DD}$, $g_{\psi D^{*}D}$, $g_{\psi D^{*}D^{*} }$ are the coupling constants that are determined by using experimental data from PDG [5], see Sect. 2.2 for details. Moreover $\psi _{\mu \nu }=\partial _{\mu }\psi _{\nu }-\partial _{\nu }\psi _{\mu }$ and ${\tilde{\psi }}_{\mu \nu }=\frac{1}{2}\varepsilon _{\mu \nu \rho \sigma }\psi ^{\rho \sigma }$ are the vector-field tensor and its dual. In particular, the term ${\mathcal {L}}_{VPP}$ describes the decay processes $\psi (4040)\rightarrow D^{+}D^{-}$, $\psi (4040)\rightarrow D^{0}{\bar{D}}^{0}$ and $\psi (4040)\rightarrow D_{s} ^{+}D_{s}^{-}$, the term ${\mathcal {L}}_{VPV}$ the processes $\psi (4040)\rightarrow D^{*0}{\bar{D}}^{0}+h.c.$, $\psi (4040)\rightarrow D^{*+}D^{-}+h.c.$ and $\psi (4040)\rightarrow D_{s}^{*+}D_{s}^{-}+h.c.$, and, finally, the term ${\mathcal {L}}_{VVV}$ the transitions $\psi (4040)\rightarrow D^{*+}D^{*-}$ and $\psi (4040)\rightarrow D^{*0}{\bar{D}}^{*0}$. The masses of the particles are taken from the PDG: $m_{D^+}=m_{D^-}=1869.65 \pm 0.05$ MeV, $m_{D^0}=m_{{\bar{D}}^0}=1864.83 \pm 0.05$ MeV, $m_{D^{*0}}=m_{{\bar{D}}^{*0}}=2006.85 \pm 0.05$ MeV, $m_{D^{*+}}=m_{D^{*-}}=2010.26 \pm 0.05$ MeV, $m_{D^+_s}=m_{D^-_s}=1968.34 \pm 0.07$ MeV and $m_{D_s^{*+}}=m_{D_s^{*-}}=2112.2 \pm 0.4$ MeV. Other decay channels (as for instance $D^*_sD^*_s$) are not considered because they are kinematically forbidden.

As usual, the theoretical expressions for the tree-level decay widths for each type of decay can be obtained from the Feynman rules and read (by keeping the mass of the decaying state as ‘running’ and denoted by m)

$$\begin{aligned}&\Gamma _{\psi \rightarrow D^{+}D^{-}+h.c}(m)\nonumber \\&\quad =\frac{\left[ k(m,m_{D^{+} },m_{D^{-}})\right] ^{3}}{6\pi m^{2}}g_{\psi DD}^{2}F_{\Lambda } (k)~, \end{aligned}$$

(5)

$$\begin{aligned}&\Gamma _{\psi \rightarrow D^{*+}D^{-}+h.c}(m)\nonumber \\&\quad =\frac{2}{3}\frac{\left[ k(m,m_{D^{*+}},m_{D^{-}})\right] ^{3}}{\pi }g_{\psi D^{*}D} ^{2}F_{\Lambda }(k)~, \end{aligned}$$

(6)

$$\begin{aligned}&\Gamma _{\psi \rightarrow D^{*+}D^{*-}}(m)\nonumber \\&\quad =\frac{2}{3}\frac{\left[ k(m,m_{D^{*+}},m_{D^{*-}})\right] ^{3}}{\pi m_{D^{*+}}^{2}}g_{\psi D^{*}D^{*}}^{2}\nonumber \\&\qquad \times \left[ 2 +\frac{\left[ k(m,m_{D^{*+}},m_{D^{*-} })\right] ^{2}}{m_{D^{*+}}^{2}}\right] F_{\Lambda }(k)~. \end{aligned}$$

(7)

The quantity

$$\begin{aligned} k\equiv |\vec {k}|\equiv k(m,m_{A},m_{B})=\frac{\sqrt{\left( m^{2}-m_{A} ^{2}\!-\!m_{B}^{2}\right) ^{2}-4m_{A}^{2}m_{B}^{2}}}{2m}\nonumber \\ \end{aligned}$$

(8)

is the modulus of the three-momentum of one of the outgoing mesons A or B,  with masses $m_{A}$ and $m_{B}$ respectively, in the rest frame of the decaying particle with mass m. The tree-level on-shell decay width for the state $\psi (4040)$ are obtained by setting $m=m_{\psi (4040)}=4.04$ GeV (here and in the following, we use the average mass $4039.6\pm 4.3$ MeV [5] rounded to $4.040\pm 0.004$ GeV).

Another important quantity is the vertex function (or form factor) $F_{\Lambda }(k)$, which assures that each quantity calculated in our model is finite. Note, one could include the vertex function directly in the Lagrangian by making it nonlocal [32,33,34,35,36] (covariance can be also preserved [37]). In our study we employed a Gaussian form factor

$$\begin{aligned} F_{\Lambda }\equiv F_{\Lambda }^{\text {Gauss}}(k)=e^{-2\frac{k^{2}}{\Lambda ^{2} }}\text { ,} \end{aligned}$$

(9)

which emerges in microscopic approaches such as $^{3}P_{0}$ mechanisms (which models the creation of quark-antiquark pairs in the QCD vacuum) used in quark models [38, 39]. However, there are other possibilities of choosing the cutoff function, the basic requirements being a smooth behavior (a step function is not an admissible choice) and a sufficiently fast decrease on the real positive axis. For completeness, another smooth form factor

$$\begin{aligned} F_{\Lambda }\equiv F_{\Lambda }^{\text {Dipolar}}(k)=\left( 1+\frac{k^{4} }{\Lambda ^{4}}\right) ^{-2} \end{aligned}$$

(10)

has been tested here in order to check how the results depend on the choice of this function. As we shall see, there are no substantial changes.

What is rather important is the numerical value of $\Lambda$. We expect a value between $\sim 0.4$ and $\sim 0.8$ GeV. Namely, for the light $\kappa$ meson, $\Lambda \simeq 0.5$ GeV was obtained by a fit to data [29]. In the recent work of Refs. [30, 31], a even smaller value $\Lambda \approx 0.3$ GeV is found (but a value of about 0.4 GeV also delivers results compatible with data). A comparison with the $^{3}P_{0}$ model induces a cutoff of $\Lambda \approx 0.8$ GeV [38, 39] (but that approach was typically not employed to calculate meson–meson loops). In this work, we test how the results vary upon changing $\Lambda$ in the range from 0.4 to 0.8 GeV (and for different form factors), but only up to 0.6 GeV physically acceptable results are obtained.

It should be stressed that our approach is an effective model of QCD, therefore the value of $\Lambda$ does not represent the maximal value for the possible values of the momentum k. When k is larger than $\Lambda$, then that particular decay is suppressed, but this is a physical consequence of the nonlocal interaction between the decaying meson and its decay products (all of them being extended objects). The momentum k can take any value from 0 to $\infty$, even arbitrarily larger than $\Lambda .$ In particular, the normalization of the spectral function (a crucial feature of our approach, see below) involves an integration up to $k\rightarrow \infty$. Of course, even if it is allowed to take k arbitrarily large, this does not imply that the model is physically complete: since we take into account a single resonance, the $\psi (4040)$, our approach can describe some of the features around 4 GeV (and up to about 4.15 GeV). Above that, one should include the resonance $\psi (4160)$ and, even further, the resonance $\psi (4415).$ (For completeness, we have tested the case in which $\psi (4040)$ and $\psi (4160)$ are present at the same time. As we shall comment later on, including $\psi (4160)$ does not substantially change the results for $\psi (4040)$).

Next, the scalar part of the propagator of the vector field $\psi _{\mu }$ is

$$\begin{aligned} \Delta _{\psi }(p^{2}=m^{2})=\frac{1}{m^{2}-M_{\psi }^{2}+\Pi (m^{2} )+i\varepsilon }\text { ,} \end{aligned}$$

(11)

where $M_{\psi }$ is the bare mass of the vector state $\psi (4040)$ (to be identified with the seed ${\bar{c}}c$ mass in absence of loop corrections). The quantity $\Pi (m^{2})={\text {Re}}(\Pi (m^{2}))+i{\text {Im}} (\Pi (m^{2}))$ is the one-particle irreducible self-energy. At the one-loop level, $\Pi (m^{2})$ is the sum of all one-loop contributions:

$$\begin{aligned} \Pi (m^{2})&=\Pi _{D^{+}D^{-}}(m^{2})+\Pi _{D^{0}{\bar{D}}^{0}}(m^{2} )+\Pi _{D_{s}^{+}D_{s}^{-}}(m^{2})\nonumber \\&\quad +\Pi _{D^{*0}{\bar{D}}^{0}+h.c} (m^{2})\nonumber \\&\quad +\Pi _{D^{*+}D^{-}+h.c}(m^{2})+\Pi _{D_{s}^{*+}D_{s}^{-}+h.c} (m^{2})\nonumber \\&\quad +\Pi _{D^{*0}{\bar{D}}^{*0}}(m^{2})+\Pi _{D^{*+}D^{*-}} (m^{2})+\ldots , \end{aligned}$$

(12)

where dots refer to further subleading contributions of other small decay channels.

For future convenience, it is also useful to define the one-loop contributions without the coupling constants at the vertices. For instance, in the $D^{*+}D^{-}$ case, one has:

$$\begin{aligned} \Pi _{D^{*+}D^{-}+h.c.}(m)=g_{\psi D^{*}D}^{2}\Sigma _{D^{*+} D^{-}+h.c.}(m)\text { .} \end{aligned}$$

Similar definitions hold in all other channels.

Moreover, the imaginary part ${\text {Im}}(\Pi (m^{2}))$ reads (optical theorem)

$$\begin{aligned}&{\text {Im}}(\Pi (m^{2})) =m\Bigg (\Gamma _{\psi (4040)\rightarrow DD}(m)+\Gamma _{\psi (4040)\rightarrow D_{s}D_{s}}\nonumber \\&\qquad +\Gamma _{\psi (4040)\rightarrow D^{*}D}(m) +\Gamma _{\psi (4040)\rightarrow D_{s}^{*}D_{s}}(m)\nonumber \\&\qquad +\Gamma _{\psi (4040)\rightarrow D^{*}D^{*}}(m)\Bigg )\text { ,} \end{aligned}$$

(13)

where:

$$\begin{aligned}&\Gamma _{\psi (4040)\rightarrow DD}(m) \nonumber \\&\quad =\Gamma _{\psi (4040)\rightarrow D^{+}D^{-} }(m)+\Gamma _{\psi (4040)\rightarrow D^{0}{\bar{D}}^{0}}(m)\text {,} \end{aligned}$$

(14)

$$\begin{aligned}&\Gamma _{\psi (4040)\rightarrow D^{*}D}(m) \nonumber \\&\quad =\Gamma _{\psi (4040)\rightarrow D^{*+}D^{-}+h.c}(m) +\Gamma _{\psi (4040)\rightarrow D^{*0}{\bar{D}}^{0} +h.c}(m)\text {,} \nonumber \\\end{aligned}$$

(15)

$$\begin{aligned}&\Gamma _{\psi (4040)\rightarrow D^{*}D^{*}} \nonumber \\&\quad =\Gamma _{\psi (4040)\rightarrow D^{*+}D^{*-}}+\Gamma _{\psi (4040)\rightarrow D^{*0}{\bar{D}}^{*0} }\text { .} \end{aligned}$$

(16)

The real part ${\text {Re}}(\Pi (m^{2}))$ is calculated by dispersion relations. For instance, for the decay channel $\psi (4040)\rightarrow D^{+}D^{-}$ one has:

$$\begin{aligned}&{\text {Re}}(\Pi _{D^{+}D^{-}}(m^{2}))\nonumber \\&\quad =-\frac{1}{\pi } PP \int \limits _{2m_{D^{+}}}^{\infty }2m^{\prime }\frac{m^{\prime }\Gamma _{\psi (4040)\rightarrow D^{+}D^{-}}(m^{\prime })}{m^{2}-m^{\prime 2}}\mathrm {dm} ^{\prime }\text { ;} \end{aligned}$$

(17)

(similar expressions hold for all other channels) (Fig. 1).

Fig. 1

Example of one-loop contribution. Here the case of $D^{*0}$ and ${\bar{D}}^{0}$ is shown

The spectral function is connected to the imaginary part of the propagator introduced above as

$$\begin{aligned} d_{\psi }(m)=\frac{2m}{\pi }\left| {\text {Im}}\Delta _{\psi } (p^{2}=m^{2})\right| \text { .} \end{aligned}$$

(18)

The quantity $d_{\psi }(m)\mathrm {dm}$ determines the probability that the state $\psi (4040)$ has a mass between m and $m+\mathrm {dm}$. It must fulfill the normalization condition

$$\begin{aligned} \int \limits _{0}^{\infty }d_{\psi }(m)\mathrm {dm}=1\text { .} \end{aligned}$$

(19)

      The validity of the normalization is a crucial feature of our study, since it guarantees unitarity [40]. It is a consequence of our theoretical approach (for a detailed mathematical proof of its validity, see Ref. [41]). Note, in Eq. (19) the integration is up to $m\rightarrow \infty$ [hence, $k\rightarrow \infty ,$ see Eq. (8)]. In practice, we shall verify numerically that Eq. (19) is fulfilled (we do so by integrating up to 10 GeV, far above the region of interest of about 4 GeV).

In addition, the partial spectral functions read:

$$\begin{aligned} d_{\psi \rightarrow DD}(m)&=\frac{2m}{\pi }\left| \Delta _{\psi } (m^{2})\right| ^{2}m\Gamma _{\psi (4040)\rightarrow DD}(m)\text { ,}\end{aligned}$$

(20)

$$\begin{aligned} d_{\psi \rightarrow D_{s}^{+}D_{s}^{-}}(m)&=\frac{2m}{\pi }\left| \Delta _{\psi }(m^{2})\right| ^{2}m\Gamma _{\psi (4040)\rightarrow D_{s}D_{s} }(m)\text { ,}\end{aligned}$$

(21)

$$\begin{aligned} d_{\psi \rightarrow DD^{*}}(m)&=\frac{2m}{\pi }\left| \Delta _{\psi }(m^{2})\right| ^{2}m\Gamma _{\psi (4040)\rightarrow DD^{*}}(m)\text { ,}\end{aligned}$$

(22)

$$\begin{aligned} d_{\psi \rightarrow D_{s}^{*}D_{s}}(m)&=\frac{2m}{\pi }\left| \Delta _{\psi }(m^{2})\right| ^{2}m\Gamma _{\psi (4040)\rightarrow D_{s} ^{*}D_{s}}(m)\text { ,}\end{aligned}$$

(23)

$$\begin{aligned} d_{\psi \rightarrow D^{*}D^{*}}(m)&=\frac{2m}{\pi }\left| \Delta _{\psi }(m^{2})\right| ^{2}m\Gamma _{\psi (4040)\rightarrow D^{*}D^{*}}(m)\text { .} \end{aligned}$$

(24)

For instance, $d_{\psi \rightarrow DD^{*}}(m)\mathrm {dm}$ is the probability that the resonance $\psi (4040)$ has a mass between m and $m+\mathrm {dm}$ and decays in the channel $DD^{*}$ [42]. Similar interpretations hold for all other channels. These quantities are physically interesting since they emerge when different channels are studied; if, for instance, the process $e^{+}e^{-}\rightarrow DD^{*}$ is considered, the corresponding cross section is proportional to $d_{\psi \rightarrow DD^{*}}(m)$.

2.2 Determination of the parameters

Our model contains five free parameters: the three coupling constants $g_{\psi DD}$, $g_{\psi D^{*}D}$, $g_{\psi D^{*}D^{*}}$ entering Eqs. (2), (3), and (4), the bare mass of the vector state $M_{\psi }$ entering in the propagator (11), and the energy scale (cutoff) $\Lambda$ included in Eqs. (9) or (10).

Table 1 Position of the poles in the complex plane for different parameters used in the model. The coupling constants $g_{\psi DD}$ and $g_{\psi D^{*}D^{*}}$ are dimensionless, $g_{\psi D^{*}D}$ has dimensions GeV$^{-1}$, $\Lambda$ and $M_{\psi }$ are in GeV

We proceed as follows: first, we fix the value of $\Lambda$ in the range 0.4–0.6 GeV. Then, in order to determine the coupling constants three experimental values are needed. We use the following measured ratios of branching fractions reported in PDG [5] (see also Refs. [15, 43,44,45]):

$$\begin{aligned}&\left. \frac{{\mathcal {B}}(\psi (4040)\rightarrow D{\bar{D}})}{{\mathcal {B}} (\psi (4040)\rightarrow D^{*}{\bar{D}})}\right| _{\exp }=0.24\pm 0.05\pm 0.12\text {,} \end{aligned}$$

(25)

$$\begin{aligned}&\left. \frac{{\mathcal {B}}(\psi (4040)\rightarrow D^{*}{\bar{D}}^{*} )}{{\mathcal {B}}(\psi (4040)\rightarrow D^{*}{\bar{D}})}\right| _{\exp }=0.18\pm 0.14\pm 0.03\text {,} \end{aligned}$$

(26)

where the first error is statistical and the second is systematic. Moreover, we employ the experimental value of the total width of the $\psi (4040)$ resonance PDG [5]

$$\begin{aligned} \Gamma _{\psi (4040)}^{\text {tot,exp}}=80\pm 10\;\text {MeV .} \end{aligned}$$

(27)

The vector state $\psi (4040)$ decays into various two-body final states. The decay channels contributing mostly to its total decay width are: DD, $D_{s}D_{s}$, $D^{*}D$, $D_{s}^{*}D_{s}$ and $D^{*}D^{*}$. The corresponding theoretical expression for the total decay width of $\psi (4040)$ state is given by

$$\begin{aligned} \Gamma _{\psi (4040)}^{\text {tot,theory}}= & {} \Gamma _{\psi (4040)\rightarrow DD}^{\text {on}~\text {shell}}+\Gamma _{\psi (4040)\rightarrow D_{s}D_{s}}^{\text {on}~\text {shell}}\nonumber \\&+\,\Gamma _{\psi (4040)\rightarrow D^{*}D}^{\text {on}~\text {shell}}+\Gamma _{\psi (4040)\rightarrow D_{s}^{*}D_{s}}^{\text {on}~\text {shell}}\nonumber \\&+\,\Gamma _{\psi (4040)\rightarrow D^{*}D^{*} }^{\text {on}~\text {shell}}\text { ,} \end{aligned}$$

(28)

where “on-shell” means that the physical PDG mass $m=4.04$ GeV is employed.

Finally, the coupling constants $g_{\psi DD}$, $g_{\psi D^{*}D}$ and $g_{\psi D^{*}D^{*}}$ as well as their errors are obtained upon minimizing the $\chi ^2$ function $F_{E}$ depending on all this three parameters:

$$\begin{aligned}&F_{E}(g_{\psi DD},g_{\psi D^{*}D},g_{\psi D^{*}D^{*}}) \nonumber \\&\quad =\left( \frac{\frac{\Gamma _{\psi \rightarrow DD}^{theory}(g_{\psi DD})}{\Gamma _{\psi \rightarrow D^{*}D}^{theory}(g_{\psi D^{*}D})}-\left. \frac{{\mathcal {B}}(\psi (4040)\rightarrow D{\bar{D}})}{{\mathcal {B}}(\psi (4040)\rightarrow D^{*}{\bar{D}})}\right| _{\exp }}{\delta \left. \frac{{\mathcal {B}}(\psi (4040)\rightarrow D{\bar{D}})}{{\mathcal {B}}(\psi (4040)\rightarrow D^{*}{\bar{D}})}\right| _{\exp }}\right) ^{2} \nonumber \\&\qquad + \left( \frac{\frac{\Gamma _{\psi \rightarrow D^{*}D^{*}}^{theory}(g_{\psi D^{*}D^{*}})}{\Gamma _{\psi \rightarrow D^{*}D}^{theory}(g_{\psi D^{*}D})}-\left. \frac{{\mathcal {B}}(\psi (4040)\rightarrow D^{*}\bar{D}^{*})}{{\mathcal {B}}(\psi (4040)\rightarrow D^{*}{\bar{D}})}\right| _{\exp }}{\delta \left. \frac{{\mathcal {B}}(\psi (4040)\rightarrow D^{*}\bar{D}^{*})}{{\mathcal {B}}(\psi (4040)\rightarrow D^{*}{\bar{D}})}\right| _{\exp }}\right) ^{2} \nonumber \\&\qquad +\left( \frac{\Gamma _{\psi (4040)}^{tot,theory} (g_{\psi DD},g_{\psi D^{*}D},g_{\psi D^{*}D^{*}})-\Gamma _{\psi (4040)}^{tot,exp}}{\delta \Gamma _{\psi (4040)}^{tot,exp}}\right) ^{2} \end{aligned}$$

(29)

The bare mass $M_{\psi }$ was fixed under the requirement that the maximum of the spectral function corresponds to the nominal mass of $\psi (4040)$, hence to 4.04 GeV.

Fig. 2

a Presents the shape of the spectral function of the resonance $\psi (4040)$ of Eq. (18) (blue line) with comparison to the standard Breit–Wigner form (red line). b Presents the partial spectral functions for $D_{s}D_{s}$, $D^{*}D$ and $D^{*}D^{*}$ channels [these channels are described by Eqs. (14), (15) and (16)]. The corresponding parameters are presented in Table 1

The numerical values of $g_{\psi DD}$, $g_{\psi D^{*}D}$, $g_{\psi D^{*}D^{*}}$ and of bare mass $M_{\psi }$ are reported in Table 1 in Sect. 3.1 for given values of the cutoff $\Lambda$. As expected, $g_{\psi D^{*}D}$, $g_{\psi D^{*}D^{*}}$ depend rather mildly on $\Lambda$, but $g_{\psi DD}$ quite strongly. This is so because the DD threshold is the most distant from the on-shell mass of the state $\psi (4040)$ and the corresponding momentum $k_{\psi DD}=\sqrt{\frac{m_{\psi (4040)}^{2}}{4}-m_{D^{+}}^{2}} \simeq 0.76$ GeV is comparable to $\Lambda .$ As a consequence, a variation of $\Lambda$ implies a sizable variation of the corresponding coupling constant. However, the decay width into DD is quite small and affects only slightly the overall picture. For completeness, we report in Appendix A also the partial decay widths for various choices of the cutoff and for different form factors. While the results are basically compatible with each other, future experimental determination of the channel $\psi (4040)\rightarrow D_{s}^{+}D_{s}^{-}$ would be very helpful to constrain our model.

3 Results

In this section we show the results and comment on them. First, in Sect. 3.1 we concentrate on the form of the full spectral function (as well as the partial ones into DD, $D^{*}D$, $D^{*}D^{*}$ channels) of the resonance $\psi (4040)$. Moreover we determine the position of the pole(s) in the complex plane. Then, in Sect. 3.2 we present the discussion of the important process $e^{+}e^{-}\rightarrow J/\psi \pi ^{+}\pi ^{-}$ and the possible generation of a peak for Y(4008).

3.1 Spectral function and pole positions

Since scattering data have still quite large errors, it is not yet possible to determine the value of $\Lambda$ through a fit. Moreover, such a fit would also need to include an unknown background contribution. This is why $\Lambda$ has been varied in a quite large range in Table 1, in which the positions of the poles have also been reported. As already mentioned in the introduction, we always find two poles in the complex plane, one corresponding to the maximum of the spectral function and an additional dynamically generated one. For $\Lambda$ up to $\sim 0.5$ GeV similar results are obtained, but for larger values the second pole (even if always present) appears at higher values. We have also tested values larger than 0.6 GeV, but they do not generate satisfactory results. (This outcome is in agreement with the results of Sect. 3.2 and Appendix B, see later on).

In the following, we choose for the numerical value $\Lambda =0.42,$ since it generates a pole in terms of the variable $\sqrt{s}$ whose imaginary part is 40 MeV, then $\Gamma _{\psi (4040)}^{\text {pole}}=80$ MeV. (Note, through all this work we use $\sqrt{s}$ , hence the definition $\sqrt{s_{\psi (4040)}^{\text {pole}}}=m_{\psi (4040)}^{\text {pole}}-i\Gamma _{\psi (4040)}^{\text {pole}}/2$ holds.) We then use this value for illustration and for the presentation of the plots. (Yet, it should not be considered as a sharp value for the cutoff). The spectral function $d_{\psi }(m)$ defined in Eq. (18) is shown in Fig. 2a together with a standard Breit–Wigner function peaked at 4.04 GeV and with a width of 80 MeV, which serves for comparison.

Only one single peak close to 4.04 GeV corresponding to the standard seed ${\bar{c}}c$ state is present. While the Breit–Wigner function approximates quite well the spectral function close to the peak, sizable deviations close to 3.9 GeV are present. This is due to an enhancement in the energy region below 4 GeV, which is generated by meson–meson loops.

Fig. 3

Study of the changing of the spectral function (a) and pole movement in the complex plane (b) for different values of $\lambda$. The used parameters for the plots are: $g_{\psi DD}=39.84$, $g_{\psi D^{*} D}=3.44$ GeV$^{-1}$, $g_{\psi D^{*}D^{*}}=1.85$, $M_{\psi }=4.007$ GeV and Gaussian form factor $F_{\Lambda }$ with $\Lambda =0.42$ GeV

In Fig. 2b we present the contributions of individual channels (DD, $D^{*}D$ and $D^{*}D^{*}$) to the total spectral function. The $D^{*}D$ channel turns out to be the most important for the deformation on the l.h.s. of the spectral function. In the complex plane we found two poles: one for 4.053–0.039 i GeV, corresponding to $\psi (4040)$ resonance, and one for 3.934–0.030 i GeV. Thus, even if only one single seed state identified with $\psi (4040)$ was included into the calculations, two poles naturally emerge.

At a first sight, it is tempting to identify the additional pole with the controversial state Y(4008). Moreover, a look at Table 1 shows that a second additional pole always exists, and up to values of about 0.5 GeV the dynamically generated pole is not far from 4 GeV. However, care is needed for the following reasons: the pole width of the additional state is too small when compared to the experimental value (about 200 MeV), but it should be stressed that a direct comparison of this pole width and the experiment is misleading, since very different reactions were measured, see later on. If the enhancement in Fig. 1 is real, it should be visible in the cross-section of the channel $e^{+}e^{-}\rightarrow \psi (4040)\rightarrow D^{*}D,$ which is proportional to $d_{\psi \rightarrow DD^{*}}(m)$ (see Fig. 1b); presently, the data have too large errors to resolve such a complicated structure. Quite importantly, the state Y(4008) has been observed in the ISR reaction $e^{+}e^{-}\rightarrow \gamma e^{+}e^{-}\rightarrow \gamma \pi ^{+}\pi ^{-}J/\psi$ and not in the $DD^{*}$ channel, see next section for the discussion of this important point.

Next, we perform a large-$N_{c}$ study of the resonance $\psi (4040)$ (where $N_{c}$ refers to the number of colors in QCD). To this end, we introduce the scaling parameter $\lambda \in (0,1)$, linked to $N_{c}$ as $\lambda =3/N_{c}$, and consider the scaling of the coupling constants as [46]

$$\begin{aligned}&g_{\psi DD}\rightarrow \sqrt{\lambda }g_{\psi DD}\text { , }g_{\psi D^{*} D}\rightarrow \sqrt{\lambda }g_{\psi D^{*}D}\text { , }g_{\psi D^{*} D^{*}}\nonumber \\&\qquad \rightarrow \sqrt{\lambda }g_{\psi D^{*}D^{*}}\text { }. \end{aligned}$$

(30)

Clearly, by setting $\lambda =1$, we reobtain our physical results. In the opposite limiting case, $\lambda =0$, the spectral function reduced to a delta function centered in the seed mass, $\delta (m-M_{\psi })$. In Fig. 3 we test the intermediate values $\lambda =0.8,0.6$ and 0.4 for both the spectral function and the positions of the poles (for the latter, $\lambda =0.5$ is also shown).

The large-$N_{c}$ study shows that for smaller $\lambda$ (hence, larger $N_{c}$), the left enhancement in the spectral function becomes smaller and finally disappears, while the spectral function becomes narrower. For what concerns the pole trajectory, the seed pole corresponding to $\psi (4040)$ resonance moves towards to real energy axis, while the additional companion pole moves away from it. This behavior confirms that the resonance $\psi (4040)$ is a conventional $q{\bar{q}}$ meson while the second pole is dynamically generated.

As a last point, we comment on mixing with other vector state, in particular with the closest quarkonium state $\psi (4160).$ By studying the mix propagator (see Refs. [47, 48] for some formal equation), we tested how the spectral function of $\psi (4040)$ changes when taking into account the processes $\psi (4040)\rightarrow DD^{*}$ $\rightarrow \psi (4160)\rightarrow DD^{*}\rightarrow \psi (4040)$ (this is so because $\psi (4160)$ also couples to $DD^*$; analogous processes with DD and $D^{*}D^{*}$ are possible). The spectral function of $\psi (4040)$ turns out to be only slightly affected in the region of interest, thus the results here presented would hold also in the enlarged scenario in which more ${\bar{c}}c$ states are considered.

3.2 Decay into $J/\psi \pi ^{+}\pi ^{-}$

An important decay channel, in which various Y states have been observed, among which the Y(4008) state is one, is the $e^{+}e^{-}\rightarrow \gamma J/\psi \pi ^{+}\pi ^{-}$ decay, where the photon $\gamma$ comes from ISR. Hence, the reaction may be recasted into two steps: $e^{+}e^{-}\rightarrow \gamma \left( e^{+} e^{-}\right) _{\text {off-shell}}$ and $\left( e^{+}e^{-}\right) _{\text {off-shell}}\rightarrow J/\psi \pi ^{+}\pi .\$

Since the very same fundamental process is involved, for simplicity we consider in the following the process $e^{+}e^{-}\rightarrow J/\psi \pi ^{+} \pi ^{-}$ (we thus ignore that the electron-positron pair is off-shell). In particular, we are interested in the case in which $\psi (4040)$ is an intermediate state of the reaction

$$\begin{aligned} e^{+}e^{-}\rightarrow \psi (4040)\rightarrow J/\psi \pi ^{+}\pi ^{-}\text { .} \end{aligned}$$

(31)

There are basically two ways in which this process can take place. The first one involves the emission of two gluons

$$\begin{aligned} \psi (4040)\equiv & {} c{\bar{c}}\rightarrow c{\bar{c}}+gg\rightarrow J/\psi \nonumber \\&+f_{0}(500)\rightarrow J/\psi +\pi ^{+}\pi ^{-} \end{aligned}$$

(32)

The choice of $f_{0}(500)$ (see [49] for a review) as an intermediate state is motivated by the fact that it is the lightest quantum state with quantum numbers of the vacuum and is in the right kinematic region ($f_{0}(980)$ is at the border, $f_{0}(1370)$ already too heavy). Nevertheless, one can repeat the very same discussion by considering different $f_{0}$ states. This decay mode can be modelled by

$$\begin{aligned} {\mathcal {L}}_{\psi jf_{0}}^{direct}=g_{\psi jf_{0}}^{direct}\psi _{\mu }j^{\mu }f_{0}\text { ,} \end{aligned}$$

(33)

where $j^{\mu }$ is the field corresponding to the $J/\psi$ and $f_{0}$ to $f_{0}(500).$ This term would generate a peak at 4.04 GeV, which has not been seen experimentally (in fact, this would be a “standard” decay of $\psi (4040)$ peaked at its Breit–Wigner mass). It means that $g_{\psi jf_{0}}^{direct}$ should be quite small. We will neglect this channel in the following.

There is, however, a second mechanism

$$\begin{aligned} \psi (4040)\rightarrow DD^{*}\rightarrow J/\psi +f_{0}(500)\rightarrow J/\psi +\pi ^{+}\pi ^{-}, \nonumber \\ \end{aligned}$$

(34)

where the additional vertex is represented by the following four-body interaction

$$\begin{aligned} {\mathcal {L}}_{DD^{*}jf_{0}}=\lambda _{DD^{*}jf_{0}}\left[ \partial ^{\mu }D^{*+\nu }D^{-}+\partial ^{\mu }D^{*0\nu }{\bar{D}}^{0}\right] j_{\mu \nu }f_{0}\text { ,} \nonumber \\ \end{aligned}$$

(35)

where $j_{\mu \nu }=\partial _{\mu }j_{\nu }-\partial _{\nu }j_{\mu }.$

The spectral function in the channel $J/\psi \pi ^{+}\pi ^{-}$ reads

$$\begin{aligned}&d_{\psi (4040)\rightarrow J/\psi \pi ^{+}\pi ^{-}}(m)\nonumber \\&\quad =\frac{2m}{\pi }\left| \Delta _{\psi }(m^{2})\right| ^{2}m\Gamma _{\psi (4040)\rightarrow J/\psi \pi ^{+}\pi ^{-}}(m) \end{aligned}$$

(36)

where

$$\begin{aligned}&\Gamma _{\psi (4040)\rightarrow J/\psi f_{0}(500)}(m)\nonumber \\&\quad =\left| g_{\psi jf_{0}}^{direct}+\lambda _{DD^{*}jf_{0}}g_{\Psi D^{*}D }\left[ \Sigma _{D^{0}D^{*0}}(m^{2}) \right. \right. \nonumber \\&\qquad \left. \left. +\Sigma _{D^{+}D^{*-}}(m^{2})\right] \right| ^{2}\frac{k}{8\pi m^{2}}\left( 3+\frac{k^{2}}{m_{J/\psi }^{2} }\right) e^{-2\frac{k^{2}}{\Lambda ^{2}}}\text { .} \end{aligned}$$

(37)

with $k\equiv k(m,m_{J/\psi },m_{f_{0}(500)}).$ It is then clear that the $DD^{*}$ loop, together with the Lagrangian ${\mathcal {L}}_{DD^{*}jf_{0}}$ of Eq. (35), generates a mass-dependent coupling for the channel $\psi (4040)\rightarrow J/\psi f_{0}(500)$.

Fig. 4

Schematic diagram of the decay $\psi (4040)\rightarrow \pi ^{+}\pi ^{-}J/\psi$ via $D^{+}D^{*-}$ loop

In general, this decay is small, since both mechanisms are suppressed according to the large-$N_{c}$ limit [46, 50]. The exact determination is difficult, since the large-$N_{c}$ properties of the unconventional state $f_{0}(500)$ are not that of a ${\bar{q}}q$ state [49, 51]. If interpreted as a dynamically generated molecular state, it does not exist in the large-$N_{c}$ limit [28, 29, 52], while its behavior as a tetraquark state is disputed. According to [46, 50], tetraquarks are also suppressed for increasing $N_{c}$, but other interpretations are possible (see the recent paper of S. Weinberg [53] and the subsequent discussions in Refs. [54, 55]). Moreover, $f_{0}(500)$ is (most probably) a mixed state and contains a subleading ${\bar{q}}q$ part, which may become more important when $N_{c}$ is large enough. Let us first consider the decay chain of Eq. (32) and let us denote as $g_{f_{0}(500)\text {-}{\bar{q}}q\text { } }$ as the coupling of $f_{0}(500)$ to quark/antiquark. Using the rules listed in Ref. [46] one has the amplitude scaling:

$$\begin{aligned}&A_{\psi (4040)\rightarrow j/\psi \text {-}gg\rightarrow j/\psi \text {-} f_{0}(500)}\nonumber \\&\quad \propto \frac{1}{\sqrt{N_{c}}}\frac{1}{\sqrt{N_{c}}}\left( \frac{1 }{\sqrt{N_{c}}}\right) ^{4}N_{c}^{2}g_{f_{0}(500)\text {-}{\bar{q}}q\text { }} \nonumber \\&\quad =\frac{g_{f_{0}(500)\text {-}{\bar{q}}q\text { }}}{N_{c}}\text { ,} \end{aligned}$$

(38)

where each conventional meson coupled to quarks counts as $1/\sqrt{N_{c}}$ and where the emission of the two gluons from one c-quark and their subsequent conversion into light quarks implies $g_{QCD}^{4}$ with $g_{QCD}\propto N_{c}^{-1/2}.$ Finally, the factor $N_{c}^{2}$ comes from the fact that there are two closed quark lines in this process. For the regular ${\bar{q}}q$ amount into $f_{0}(500),$ one may use $g_{f_{0}(500)\text {-}{\bar{q}} q\text { }}\propto N_{c}^{-1/2},$ hence the amplitude $1/N_{c}^{3/2}$ follows. Within this context, for the final state with two pions, one obtains $A_{\psi (4040)\rightarrow j/\psi \text {-}gg\rightarrow j/\psi \pi ^{+}\pi ^{-}}\propto N_{c}^{-2}.$ For other non-conventional components in $f_{0}(500),$ this decay amplitude should be even more further suppressed. Note, the leading term of the decay width into $j/\psi \pi ^{+}\pi ^{-}$ scales as $N_{c}^{-4},$ thus explaining why such decay is very small.

Fig. 5

Left panel: ${\text {Re}}\Sigma _{D^{0}D^{0*}}$ and ${\text {Re}}\Sigma _{D^{+*}D^{-}}$ as function of m. Both of them have a peak at about 3.9 GeV. Right panel: Plot of the (normalized) cross-section $\sigma _{R}(m)$ of the channel $e^{+}e^{-}\rightarrow \psi (4040)\rightarrow J/\psi \pi ^{+}\pi ^{-}$ defined in Eq. (41). A deformed and broad structure with a peak at about 3.95 GeV is visible

Next, we turn to the diagram of Fig. 4. Here, one has the following scaling of the amplitude:

$$\begin{aligned}&A_{\psi (4040)\rightarrow j/\psi \text {-}DD^{*}\rightarrow j/\psi \text {- }f_{0}(500)}\propto \frac{1}{\sqrt{N_{c}}}\frac{1}{\sqrt{N_{c}}}g_{f_{0}(500) \text {-}{\bar{q}}q\text { }}\nonumber \\&\quad =\frac{g_{f_{0}(500)\text {-}{\bar{q}}q\text { }}}{N_{c} }\text { }, \end{aligned}$$

(39)

which coincides with the previous decay amplitude. (Here, the fact that the meson $\hat{}~3$ coupling goes as $1/\sqrt{N_{c}}$ has been used.) Similar considerations follow: the amplitude is (at least) suppressed as $N_{c}^{-3/2},$ while the amplitude into $j/\psi \pi ^{+}\pi ^{-}$ is suppressed (at least) as $N_{c}^{-2.}$ Thus, from the point of view of the large-$N_{c}$ expansion both mechanisms of Eqs. (32) and (34) are equally suppressed. Yet, the second mechanism is expected to be dominant in our case for various reasons. Namely, the seed state $\psi (4040)$ couples strongly to $D^{*}D$ (this is the dominant decay mode). Moreover, the real part of the $DD^{*}$ loops, depicted in Fig. 5a, has a pronounced peak at $m_{D}+m_{D^{*}}$ at about 3.9 GeV. The transition $D^{*}D\rightarrow J/\psi f_{0}(500)$ via the Lagrangian ${\mathcal {L}}_{DD^{*}jf_{0}}$ is rather natural, since it implies a redistribution of already existing quarks. Moreover, $f_{0}(500)$ couples strongly to ${\bar{u}}u$ and ${\bar{d}}d.$ Hence, in first approximation we shall neglect $g_{\psi jf_{0}}^{direct}$ in the following.

Summarizing, in Eq. (36) the product of two functions is present: $\left| \Delta _{\psi }(m^{2})\right| ^{2},$ peaked at 4.04 GeV, and $\Gamma _{\psi (4040)\rightarrow J/\psi f_{0}(500)}(m),$ peaked at 3.9 GeV. Of course, other channels, such as DD,  would also couple $J/\psi f_{0}(500)$, but (i) the coupling of $\psi (4040)$ to DD is sizably smaller and (ii) the function $\left| \Sigma _{DD}(m^{2})\right| ^{2}$ is peaked at the DD threshold, hence the overlap with $\left| \Delta _{\psi }(m^{2})\right| ^{2}$ is negligible.

For $\sqrt{s}$ in the region of interest, we have

$$\begin{aligned}&\sigma _{e^{+}e^{-}\rightarrow \psi (4040)\rightarrow J/\psi \pi ^{+}\pi ^{-} }\left( m\right) \nonumber \\&\quad =\frac{2\pi }{m}g_{\psi e^{+}e^{-}}^{2}d_{\psi (4040)\rightarrow J/\psi \pi ^{+}\pi ^{-}}(m)\text { .} \end{aligned}$$

(40)

In Fig. 5b we plot

$$\begin{aligned} \sigma _{R}(m)=\frac{\sigma _{e^{+}e^{-}\rightarrow \psi (4040)\rightarrow J/\psi \pi ^{+}\pi ^{-}}\left( m\right) }{\sigma _{e^{+}e^{-}\rightarrow \psi (4040)\rightarrow J/\psi \pi ^{+}\pi ^{-}}\left( m=m_{D^{0}}+m_{D^{*0} }\right) }. \nonumber \\ \end{aligned}$$

(41)

(In this way, the dependence on the unknown coupling $\lambda _{DD^{*} jf_{0}}$ cancels and $\sigma _{R}(m_{D^{0}}+m_{D^{*0}})=1$). The resulting form is quite peculiar and is definitely not a simple Breit–Wigner peak. If the experimental accuracy is not high enough, one may identify this signal as a broad resonance whose peak is centered at about 4 GeV. Thus, we suggest that the ‘state’ Y(4008) is a manifestation of the standard state $\psi (4040),$ which arises due to the decay into $J/\psi \pi ^{+}\pi ^{-}$ through the nearby $DD^{*}$ loop. It is important to stress that this conclusion is independent of the presence of a dynamically generated pole and its precise position, but it is a consequence of the strong coupling to $DD^{*}$ and the fact that the $DD^{*}$ threshold is not far from the peak of $\psi (4040)$. We call the phenomenon leading to Fig. 5b as ‘real-part-loop generated peak’.

This mechanism for the generation of the state Y(4008) does not necessarily correspond to a resonance. Moreover, the width of the dynamical pole is not related to the width of the signal of Fig. 5b. In Appendix B we also present the results for variations of the parameters, and for a cutoff up to 0.5 GeV a very similar outcome is obtained.

4 Summary and discussion

We studied the energy region close to the resonance $\psi (4040)$ in the framework of a QFT effective model. We evaluated its spectral function and found that, besides the expected resonance pole of $\psi (4040)$ (corresponding to peak in the spectral function and to the underlying seed $c{\bar{c}}$ state), an additional companion pole (no peak, but an enhancement in the spectral function) emerges naturally within our approach (see Table 1). Illustrative result in agreement with phenomenology are: $4.053-i0.039$ GeV for the seed state $\psi (4040)$ and $3.934-i0.030$ GeV for the enhancement. A large-$N_{c}$ study confirms that $\psi (4040)$ resonance is predominantly a charm–anticharm state, while the second pole is dynamically generated by meson–meson quantum fluctuations.

The pole itself cannot be directly associated to the $Y(4008).\$Namely, this pole would be mostly visible in a two-peak structure of the cross-section $e^{+}e^{-}\rightarrow \psi (4040)\rightarrow DD^{*}$ (whose data precision is not good enough). Yet, the chain $e^{+}e^{-}\rightarrow \psi (4040)\rightarrow DD^{*}\rightarrow J/\psi +f_{0}(500)\rightarrow J/\psi \pi ^{+}\pi ^{-}$is quite promising: the $DD^{*}$ loop generates a peak at about 3.9 GeV in the cross-section. The strong coupling of $\psi (4040)$ to $DD^{*}$ and the overlap of the $DD^{*}$-loop function with the modulus square of the propagator are responsible for a quite broad peak in the corresponding spectral function, which can be identified with Y(4008). The important point is that the existence of an additional pole corresponding to Y(4008) is possible (and indeed it does exist for the parameters of our model) but is actually not necessary for the process that we describe. The very same mechanism can be also investigated in the future in other channels, as for instance in connection with the states Y(4260) and $\psi (4160)$.

Appendices

A Decay widths on shell

Here we present the results for the on-shell decay widths for both form factors and for different values of the cutoff, respectively. Even if the qualitative picture does not change much, one observes non-negligible variations of the partial decay widths as function of the cutoff $\Lambda .$ In the future, a better determination of the decay $\psi (4040)$ into $D_{s} ^{+}D_{s}^{-}$ (presently only seen) would constitute a useful constraint on our model.

B Decay into $J/\psi \pi ^{+}\pi ^{-}$-variations of the parameter $\Lambda$

As it was discussed in the paper, the decay of $\psi (4040)$ into $J/\psi \pi ^{+}\pi ^{-}$ via $DD^{*}$ loops generates a sizable peak in the cross-section in the energy region close to 4 GeV. This is an important aspect of our theoretical framework, thus we present in Fig. 6 how the results of the cross-section of Eq. (41) depend on the different values of cutoff.

Fig. 6

Study of the cross-section ratio of Eq. (41) upon changing the value of cutoff parameter

Table 2 Partial decay widths on-shell for both types of form factor and different values of $\Lambda$ parameter

When $\Lambda$ varies in the range from 0.38 up to (at most) 0.5 GeV one observes a broad peak at about 3.95 GeV. For $\Lambda =0.5$ one has actually a quite broad structure, but the peak is already at about 4.04 GeV. For larger values of $\Lambda$, there is a single peak close to 4.04 GeV which corresponds to standard $c{\bar{c}}$ seed state $\psi (4040)$. Hence, for the generation of a signal resembling Y(4008) the value of $\Lambda$ should not exceed 0.5 GeV (Table 2).