Decay widths of $^3 P_J$ charmonium to $DD,DD^*,D^*D^*$ and corresponding mass shifts of $^3 P_J$ charmonium

Abstract

In this work, we calculate the amplitudes of the processes $c{{\bar{c}}}({^3P_J}) \rightarrow DD,DD^*, D^*D^* \rightarrow c\bar{c}({^3P_J})$ in the leading order of the nonrelativistic expansion. The imaginary parts of the amplitudes are corresponding to the branch decay widths of the charmonium $c{{\bar{c}}}({^3P_J}) \rightarrow DD,DD^*, D^*D^*$ and the real parts are corresponding to the mass shifts of the charmonium $c{{\bar{c}}}({^3P_J})$ due to these decay channels. After absorbing the polynomial contributions which are pure real, the ratios between the branch decay widths and the corresponding mass shifts are only dependent on the center-of-mass energy. We find the decay widths and the mass shifts of the $^3P_2$ states are exact zero in the leading order. The ratios between the branch decay widths and the mass shifts for the $^3P_0, {^3P_1}$ states are larger than 5 when the center-of-mass energy is above the $DD,DD^*, D^*D^*$ threshold. The dependence of the mass shifts on the center-of-mass energy is nontrivial especially when the center-of-mass energy is below the threshold. The analytic results can be extended to the b quark sector directly.

1 Introduction

The energy spectrum of hadrons is a basic topic in the strong interaction. Up to now, it is still an unsolved problem due to the complex nonperturbative property of QCD. In literatures, many phenomenological models have been developed to study this problem in the quark level, such as the quark model [1,2,3], QCD sum rules [4,5,6,7,8], Bethe–Salpeter equation [9,10,11,12,13,14,15,16,17,18,19,20,21,22], and etc. In these methods, usually the annihilation effects are neglected since they are much smaller than the non-perturbative potential. Physically, if the annihilation effects can be taken as small comparing with the interaction which binds the quarks, then the imaginary part of the annihilation amplitude is corresponding to the branch decay width and the real part is corresponding to the perturbative mass shift. Theoretically such annihilation effects should be considered and estimated carefully when aiming to understand the energy spectrum precisely.

Experimentally, since 2003 many new charmonium-like states are reported by the collaborations of Belle [23,24,25,26,27,28,29], CDF [30, 31], D0 [32], BABAR [33, 34], Cleo-C [35], LHCb [36], BES [37], and CMS [38, 39]. These charmonium-like states cannot be well understood in the traditional quark model and their masses usually lie above the open charm threshold where some new decay modes are opened. In the previous study [40, 41], we studied the mass shifts of $^1S_0$ and $^3P_J$ heavy quarkonia due to the transition $q{\bar{q}}\rightarrow 2g\rightarrow q{\bar{q}}$. Physically, when the masses of the states lie above the threshold of D or $D^*$ pairs, the transitions $c{{\bar{c}}}$ to these mesons’ pairs are opened. It is natural that these opened channels not only result in the visible branch decay widths but also give contributions to the mass shifts of the corresponding charmonium. When the masses of the charmonium lie about the threshold of the meson pairs, one can expect that the nonrelativistic expansion is applicable, which means that one can take the mesons $D,D^*$ like the heavy quark in the nonrelativistic QCD to construct the effective nonrelativistic interactions order by order. In this work, we follow this spirit to calculate the amplitudes of $c{{\bar{c}}} (^3P_J) \rightarrow DD, DD^*, D^*D^* \rightarrow c{{\bar{c}}} (^3P_J)$ in the leading order of non-relativistic expansion. The imaginary parts of the results are corresponding to the branch decay widths which can be used to determine the effective coupling constants. Furthermore, if these annihilation interactions are much smaller than the binding interaction, then the real parts can be used to estimate the corresponding mass shifts.

We organize the paper as follow. In Sect. 2 we describe the basic frame to calculate the amplitudes of $c{{\bar{c}}} (^3P_J) \rightarrow DD, DD^*,D^*D^* \rightarrow {{\bar{c}}} (^3P_J)$ in the leading order of nonrelativistic expansion, in Sect. 3 we give the analytic results for the amplitudes in the leading order of nonrelativistic expansion, in Sect. 4, we present some numerical results to show some properties in detail.

2 Basic formula

When the mass of the charmonium is about $2m_D$ or $2m_{D^*}$ with $m_{D,D^*}$ being the masses of the $D,D^*$ mesons, the three-momenta of the c quarks and the mesons in the decay channels $c{\bar{c}} (^3P_J)\rightarrow DD,DD^*,D^*D^*$ are much smaller than c quarks’ mass $m_c$ or $m_{D,D^*}$. In this case, one can take $m_c\approx m_D \approx m_{D^*}$ as the large scale comparing with $\Lambda _{\text {QCD}}$ and expand the interaction on the small variables $|\overset{\rightharpoonup }{\mathbf {q}}|/m_c$ with $\overset{\rightharpoonup }{\mathbf {q}}$ the three-momenta of the c quarks and the mesons. This nonrelativistic expansion is similar with the spirit of NRQCD where the contact four point interactions are introduced. Different from NRQCD, now there is no hard gluon in the decay channels $c{\bar{c}} (^3P_J)\rightarrow DD,DD^*,D^*D^*$, but only nonrelativistic heavy quarks and heavy mesons. This means that there are only contact interactions between the c quarks and the $D,D^*$ mesons.

To construct the interaction Lagrangian order by order on $|\overset{\rightharpoonup }{\mathbf {q}}|/m_c$, one can use the same spirit with NRQCD [42] where only the three-dimension partial $\overset{\rightharpoonup }{\mathbf {\nabla }}$ on the fields is used and naively one $\overset{\rightharpoonup }{\mathbf {\nabla }}$ contributes one order. A little difference from NRQCD is that now the 0-component of the $D^*$ field $A^{\mu }$ is also a small variable since $\partial _{\mu }A^{\mu }=\partial _{0}A^{0}-\overset{\rightharpoonup }{\mathbf {\nabla }}\cdot \overset{\rightharpoonup }{\mathbf {A}}=0$ for the on shell $D^*$. This means that the inclusion of one $A^{0}$ field also contributes one order like $\overset{\rightharpoonup }{\mathbf {\nabla }}$. Then in the leading order there is no $\overset{\rightharpoonup }{\mathbf {\nabla }}$ or $A^{0}$ in the interaction Lagrangian. For convenience, one can also rewrite the Lagrangian in a covariant form and choose the leading order contribution finally. By this reason the most general interactions with C, P, T invariance in the leading order can be written as follows:

$$\begin{aligned} {{\mathcal {L}}}_1= & {} g_{a} {\overline{\psi }}\psi \phi \phi , \nonumber \\ {{\mathcal {L}}}_2= & {} g_{b}{\overline{\psi }}\gamma ^5 \gamma _{\mu }\psi \phi A^{\mu }+ \textit{h.c.},\nonumber \\ {{\mathcal {L}}}_3= & {} g_{c} {\overline{\psi }}\psi A_{\mu } A^{\mu }, \end{aligned}$$

(1)

where $\psi ,\phi ,A^{\mu }$ are the fields of the c quark, the D meson, and the $D^*$ meson, respectively. Here we do not assume that there is spin asymmetry between the D and $D^*$ mesons since the dynamics of the light quarks inside the D and $D^*$ mesons may break the spin symmetry strongly. This means that the couplings $g_{a,b,c}$ are independent.

By these interactions, the Feynman diagrams for the amplitudes of $c{{\bar{c}}} ({^3P_J})\rightarrow DD, DD^*,D^*D^* \rightarrow c\bar{c}({^3P_J})$ in the leading order are shown in Fig. 1a–c.

Fig. 1

The diagrams for $c{{\bar{c}}}(^3P_J)\rightarrow c{{\bar{c}}}(^3P_J)$ processes where a–d are corresponding to $c{{\bar{c}}} \rightarrow DD \rightarrow c{{\bar{c}}}$, $c{{\bar{c}}} \rightarrow DD^* \rightarrow c{{\bar{c}}}$, $c{{\bar{c}}} \rightarrow D^*D^* \rightarrow c\bar{c}$, and $c{{\bar{c}}}\rightarrow c{{\bar{c}}}$ via contact interactions

Similar with any effective theory, usually the contact interactions are needed to absorb the UV divergences in the loop diagrams. In the practical calculation, we find the following contact interactions are needed to absorb the UV divergences in Fig. 1a–c:

$$\begin{aligned} {{\mathcal {L}}}^{c}_1= & {} g_{10} [{\overline{\psi }}\psi ] [{\overline{\psi }}\psi ]-g_{11}\Big (\partial _\mu \partial ^\mu [{\overline{\psi }}\psi ]\Big ) [{\overline{\psi }}\psi ]\nonumber \\&+g_{12}\Big (\partial _\mu \partial ^\mu \partial _\nu \partial ^\nu [{\overline{\psi }}\psi ]\Big )[{\overline{\psi }}\psi ] , \nonumber \\ {{\mathcal {L}}}^{c}_2= & {} g_{20}[{\overline{\psi }}\gamma ^5 \gamma _{\mu }\psi ][{\overline{\psi }}\gamma ^5 \gamma ^{\mu }\psi ]\nonumber \\&- g_{21}\Big (\partial _\nu \partial ^\nu [{\overline{\psi }}\gamma ^5 \gamma _{\mu }\psi ]\Big )[{\overline{\psi }}\gamma ^5 \gamma ^{\mu }\psi ]. \end{aligned}$$

(2)

The physical reason of such contact interactions is related with the quantum number of $^3P_{0,1}$: ${{\mathcal {L}}}^{c}_1$ is the contact interaction in the scalar channel and ${{\mathcal {L}}}^{c}_2$ is the contact interaction in the axial-vector channel. We want to point out that we just write down such contact interactions here to show the exact cancellations of the UV divergences and the polynomial contributions. In the practical calculation, one can get the same final results even without knowing the form of the contact interactions. The Feynman diagram for the contribution due to these contact interactions is shown in Fig. 1d.

In the center of mass frame, we choose the four external momenta as follows:

$$\begin{aligned}&p_1 \overset{def}{=} \frac{P}{2}+q_i, \quad p_2 \overset{def}{=} \frac{P}{2}-q_i, \nonumber \\&p_3 \overset{def}{=} \frac{P}{2}+q_f, \quad p_4 \overset{def}{=} \frac{P}{2}-q_f. \end{aligned}$$

(3)

For simplicity we define $P \overset{def}{=} (\sqrt{s},0,0,0)$ and use the instantaneous approximation for $q_{i,f}$ which means that we assume $q_i= (0,\mathbf{q} _i)$ and $q_f= (0,\mathbf{q} _f)$, where we use the bold formatting to refer to the three momentum here and in the following.

To project the $c{{\bar{c}}}$ pairs to the $^3P_J$ states we use the project matrices in the on-shell case [43,44,45,46,47] which are defined as follows:

$$\begin{aligned} \sum {\bar{\nu }}(\mathbf{p }_2,s_2)T u(\mathbf{p }_1,s_1)<\frac{1}{2}s_1; \frac{1}{2}s_2|1s_i> \overset{def}{=} \text {Tr}[T\Pi _i(s_i)], \nonumber \\ \sum {\bar{u}}(\mathbf{p }_3,s_3)T\nu (\mathbf{p }_4,s_4)<\frac{1}{2}s_3;\frac{1}{2}s_4|1s_f> \overset{def}{=} \text {Tr}[T\Pi _f(s_f)],\nonumber \\ \end{aligned}$$

(4)

where the Clebsch–Gordan coefficients are the standard ones as in Refs. [44, 45], and the Dirac spinors are normalized as $u^+u=\nu ^+\nu =1$, whose definitions are expressed as

$$\begin{aligned} u(\mathbf{p }_1,s_1) \overset{def}{=} \frac{ {\overline{p}} \!\!\!/_1+m}{\sqrt{E_1(E_1+m)}}\left( \begin{array}{c} \xi ^{s_1} \\ 0 \\ \end{array} \right) , \nonumber \\ \nu (\mathbf{p }_2,s_2) \overset{def}{=} \frac{- {\overline{p}} \!\!\!/_2+m}{\sqrt{E_2(E_2+m)}}\left( \begin{array}{c} 0 \\ \eta ^{s_2} \\ \end{array} \right) , \end{aligned}$$

(5)

with $E_{1,2} =\sqrt{|\mathbf{p }_{1,2}|^2+m_c^2}$, ${\overline{p}}_{1,2}=(E_{1,2},\mathbf{p }_{1,2})$, $\xi ^{1/2}=(1,0)^T$, $\xi ^{-1/2}=(0,1)^T$, $\eta ^{1/2}=(0,1)^T$, and $\eta ^{-1/2}=(-1,0)^T$. Finally the project matrices can be written as

$$\begin{aligned} \Pi _i(s_i)= & {} N_{i} ( {\overline{p}} \!\!\!/_1+m_c)(2E_i+ {\overline{p}} \!\!\!/_1 + {\overline{p}} \!\!\!/_2) \epsilon \!\!\!/_p(s_i)(- {\overline{p}} \!\!\!/_2+m_c),\nonumber \\ \Pi _f(s_f)= & {} N_{f} (- {\overline{p}} \!\!\!/_4+m_c) \epsilon \!\!\!/_p^*(s_f)(2E_f+ {\overline{p}} \!\!\!/_3+ {\overline{p}} \!\!\!/_4)( {\overline{p}} \!\!\!/_3+m_c),\nonumber \\ \end{aligned}$$

(6)

where $E_{i,f}= \sqrt{|\mathbf{q} _{i,f}|^2+m_c^2}$, and

$$\begin{aligned}&\epsilon _p^{\mu }(0) \overset{def}{=} (0,0,0,1), \nonumber \\&\epsilon _p^{\mu }(\pm 1) \overset{def}{=} (0,\mp 1 ,-i,0)/\sqrt{2}, \end{aligned}$$

(7)

and $N_{i,f}$ are the normalized global factors which can be expressed as follows in the nonrelativistic limit

$$\begin{aligned} N_{i,f}=-\frac{1}{8\sqrt{2}E_{i,f}^2(E_{i,f}+m_c)}. \end{aligned}$$

(8)

In principle the form of the project matrix for a bounded $c{\overline{c}}$ pair should be deduced from the Bethe–Salpeter wave function or similar Lorentz covariant matrix element, while in the ultra nonrelativistic limit the above expressions are expected to be correct.

In the leading order of nonrelativistic expansion, the structure of a meson $H(^3P_J)$ can be expressed as follow:

$$\begin{aligned} |H(^3P_J)\rangle \sim \phi (|\mathbf{p }|) \frac{\delta _{ij}}{\sqrt{N_c}}|q^i{\bar{q}}^j(^3P_J)\rangle , \end{aligned}$$

(9)

where $N_c=3$ and $\phi (|\mathbf{p }|)$ is the wave function of $H(^3P_J)$ in the momentum space which is defined as

$$\begin{aligned} \phi (|\mathbf{p }|) Y_{1m}(\Omega _{\mathbf{p }}) \overset{def}{=} \int \frac{d^3\mathbf{r} }{(2\pi )^{3/2}}e^{-i\mathbf{p }\cdot \mathbf{r} }R_1(|\mathbf{r} |)Y_{1m}(\Omega _\mathbf{r }), \end{aligned}$$

(10)

with the normalization condition

$$\begin{aligned} \int d|\mathbf{r} ||\mathbf{r} |^2R^2_1(|\mathbf{r} |)=1. \end{aligned}$$

(11)

Combining the structure of $H(^3P_J)$ and the project matrices, the expression for the amplitudes in the leading order can be expressed as

$$\begin{aligned}&{{\mathcal {M}}}^{(X)}(^3P_J)\nonumber \\&\quad = \int \frac{d |\mathbf{q} _i| d |\mathbf{q} _f|}{(2\pi )^3} |\mathbf{q} _i|^2 |\mathbf{q} _f|^2 \phi (|\mathbf{q} _f|)\phi ^*(|\mathbf{q} _i|){\overline{G}}^{(X)}(^3P_J),\nonumber \\ \end{aligned}$$

(12)

where the index (X) refers to (a, b, c, d) which are corresponding to the contributions from the diagrams (a), (b), (c) and (d) shown in Fig. 1, respectively. ${\overline{G}}^{(X)}(^3P_J)$ are expressed as

$$\begin{aligned} {\overline{G}}^{(X)}(^3P_J)= & {} \sum _{s_i,s_f} \langle JJ_z|1s_f;1m_f\rangle \langle JJ_z|1s_i;1m_i\rangle \nonumber \\&\times \int d\Omega _\mathbf{q _i} d\Omega _\mathbf{q _f} Y_{1m_i}(\Omega _\mathbf{q _i}) \nonumber \\&\times \ Y_{1m_f}^*(\Omega _\mathbf{q _f}) G^{(X)}(s_i,s_f), \end{aligned}$$

(13)

with

$$\begin{aligned}&G^{(a)}(s_i,s_f)\nonumber \\&\quad = -ic_f \mu ^{2\epsilon }\int \frac{d^d k}{(2\pi )^d} \text {Tr}[T_1 \Pi _{i}(s_i)] \text {Tr}[T_1 \Pi _{f}(s_f)] S(k) \nonumber \\&\qquad \times S(p_1+p_2-k), \nonumber \\&G^{(b)}(s_i,s_f)\nonumber \\&\quad = -ic_f \mu ^{2\epsilon }\int \frac{d^d k}{(2\pi )^d} \text {Tr}\big [T_2^{\mu } \Pi _{i}(s_i)\big ] \text {Tr}\big [T_2^{\nu } \Pi _{f}(s_f)\big ]\nonumber \\&\qquad \times D_{\mu \nu }(k) S(p_1+p_2-k), \nonumber \\&G^{(c)}(s_i,s_f)\nonumber \\&\quad = -ic_f \mu ^{2\epsilon }\int \frac{d^d k}{(2\pi )^d} \text {Tr}[T_3^{\mu \rho } \Pi _{i}(s_i)] \text {Tr}[T_3^{\nu \omega } \Pi _{f}(s_f)]\nonumber \\&\qquad \times D_{\mu \nu }(k) D_{\rho \omega } (p_1+p_2-k),\nonumber \\&G^{(d)}(s_i,s_f)= -ic_f \mu ^{2\epsilon }\Big (\text {Tr}[T_4\Pi _{i}(s_i) ] \text {Tr}[\Pi _{f}(s_f)] \nonumber \\&\qquad + \text {Tr}[T_5^{\mu } \Pi _{i}(s_i)] \text {Tr}[\gamma _{5}\gamma _\mu \Pi _{f}(s_f)] \Big ), \end{aligned}$$

(14)

where $d=4-2\epsilon$ is the dimension, $\mu$ is the introduced energy scale, $c_f=\frac{\delta _{ij}}{\sqrt{N_c}}\delta _{ij}\frac{\delta _{i'j'}}{\sqrt{N_c}}\delta _{i'j'}=3$ is the color factor, and

$$\begin{aligned} T_1= & {} ig_a, \nonumber \\ T_2^{\mu }= & {} ig_b \gamma ^5 \gamma ^{\mu },\nonumber \\ T_3^{\mu \rho }= & {} ig_c g^{\mu \rho },\nonumber \\ T_4= & {} i(g_{10}+g_{11}s+g_{12}s^2),\nonumber \\ T_5^{\mu }= & {} i(g_{20} +g_{21}s)\gamma _5 \gamma ^{\mu }, \end{aligned}$$

(15)

and the propagators of the pseudoscalar S and the vector $D_{\mu \nu }$ are defined as

$$\begin{aligned} S(k)= & {} \frac{i}{k^2-m_D^2+i \varepsilon }, \nonumber \\ D_{\mu \nu }(k)= & {} \frac{-i(g^{\mu \nu }- \frac{k^{\mu }k^{\nu }}{m_{D^*}^2})}{k^2-m_{D^*}^2+i \varepsilon }. \end{aligned}$$

(16)

In the practical calculation, the package FeynCalc [48, 49] is used to do the trace in the d dimension. The packages FIESTA [50, 51] and PackageX [52, 53] are independently used to do the loop integration for double check. After the loop integrations, $G^{(X)}(s_i,s_f)$ can be expressed in the following form:

$$\begin{aligned} G^{(X)}(s_i,s_f)= & {} C_1^{(X)} \epsilon _p(s_i)\cdot \epsilon _p^*(s_f) + C_2^{(X)}\epsilon _p(s_i)\cdot {\hat{q}}_i \ \epsilon _p^*(s_f)\nonumber \\&\cdot \, {\hat{q}}_f + C_3^{(X)}\epsilon _p(s_i)\cdot {\hat{q}}_f \ \epsilon _p^*(s_f)\cdot {\hat{q}}_i, \end{aligned}$$

(17)

where $C^{(X)}_i$ can be expressed as

$$\begin{aligned} C^{(X)}_i = \sum _{n=0}^1 C^{(X)}_{in}(|\mathbf{q} _i|,|\mathbf{q} _f|)(\hat{q_i}\cdot {\hat{q}}_f)^n, \end{aligned}$$

(18)

with ${\hat{q}}_{i,f}\overset{def}{=}q_{i,f}/|\mathbf{q} _{i,f}|$, respectively.

To get the coefficients ${\overline{G}}^{(X)}(^3P_J)$, usually the sums of the spins and the integrations of the angles are calculated independently to simplify the expressions [54]. In our calculation, we directly calculate the sums of the spins and the integrations of the angles together after getting the expressions of $C^{(X)}_{in}$. This method is more efficient and has been used in our previous work [40, 41]. The relevant expressions are listed in the Appendix.

3 The energy shift of $^3P_J$ in the leading order

We expand ${\overline{G}}^{(X)}(^3P_J)$ on $|\mathbf{q} _{i}|,|\mathbf{q} _{f}|$ to order 1 as following forms:

$$\begin{aligned} {\overline{G}}^{(a,b,c)}(^3P_J)= & {} 3g^2_{a,b,c}N_{i}N_{f} m_c^4 \Big [ |\mathbf{q} _{i}||\mathbf{q} _{f}|c^{(a,b,c)}_{J} + \text {higher order} \Big ],\nonumber \\ {\overline{G}}^{(d)}(^3P_J)= & {} 3N_{i}N_{f} m_c^4\Big [ |\mathbf{q} _{i}||\mathbf{q} _{f}|c^{(d)}_{J} + \text {higher order} \Big ]. \end{aligned}$$

(19)

Here we want to emphasis that the contributions ${\overline{G}}^{(d)}(^3P_J)$ are used to absorb the UV divergences in ${\overline{G}}^{(a+b+c)}(^3P_J)$ and give no contributions to the decay widths of $^3P_J$ states. The finite parts of the contributions ${\overline{G}}^{(d)}(^3P_J)$ are arbitrary. Actually, they not only absorb the UV divergences but also absorb the polynomial contributions in ${\overline{G}}^{(a+b+c)}(^3P_J)$. This situation is a little different from the results in the $c\bar{c}(^3P_J)\rightarrow 2g \rightarrow c{{\bar{c}}}(^3P_J)$ cases where there are no any contact interactions in the original QCD interaction. The important point is that these absorptions are universal and independent on the processes, and we discuss the details in the following subsection.

3.1 The energy shift of $^3P_0$ state

After the loop integration, the sum of the spins, the integration of the angles, and the Taylor expansion, we get the following results in the $^3P_0$ channel.

$$\begin{aligned} c^{(a)}_0= & {} c^{(a)}_{0,poly}+\frac{256\sqrt{s(s-4m_D^2)}}{\pi s} \ln \nonumber \\&\times \left[ \frac{2m_D^2-s+\sqrt{s(s-4m_D^2)}}{2m_D^2}+i\varepsilon \right] , \nonumber \\ c^{(b)}_0= & {} 0, \nonumber \\ c^{(c)}_0= & {} c^{(c)}_{0,poly}+\frac{64[(s-2m_{D^*}^2)^2+8m_{D^*}^4] \sqrt{s(s-4m_{D^*}^2)}}{\pi s m_{D^*}^4} \ln \nonumber \\&\times \left[ \frac{2m_{D^*}^2-s+\sqrt{s(s-4m_{D^*}^2)}}{2m_{D^*}^2}+i\varepsilon \right] ,\nonumber \\ c^{(d)}_{0}= & {} c^{(d)}_{0,poly}, \end{aligned}$$

(20)

where $c^{(a,c,d)}_{0,poly}$ are some polynomial functions on s which include the UV divergences and are expressed as follows:

$$\begin{aligned} c^{(a)}_{0,poly}= & {} \frac{256}{\pi } \left( 2+\frac{1}{{\overline{\epsilon }}_{\text {UV}}} +\ln \frac{\mu _{\text {UV}}^2}{m_D^2}\right) , \nonumber \\ c^{(c)}_{0,poly}= & {} \frac{64}{\pi m_{D^*}^4} \left[ 4\left( 4+\frac{3}{{\overline{\epsilon }}_{\text {UV}}} +3\ln \frac{\mu _{\text {UV}}^2}{m_{D^*}^2}\right) m_{D^*}^4\right. \nonumber \\&-2\left( 5+\frac{3}{{\overline{\epsilon }}_{\text {UV}}} +3\ln \frac{\mu _{\text {UV}}^2}{m_{D^*}^2}\right) m_{D^*}^2s\nonumber \\&\left. +\left( 2+\frac{1}{{\overline{\epsilon }}_{\text {UV}}} +\ln \frac{\mu _{\text {UV}}^2}{m_{D^*}^2}\right) s^2\right] ,\nonumber \\ c^{(d)}_{0,poly}= & {} \frac{256}{\pi ^3} (g_{10}+g_{11}s+g_{12}s^2), \end{aligned}$$

(21)

with$\frac{1}{{\overline{\epsilon }}_{\text {UV}}}=\frac{1}{\epsilon _{\text {UV}}}-\gamma _E+\log (4\pi )$.

An important property of the two contributions $c^{(a,c)}_{0,poly}$ is that they can be absorbed by the contact interactions $\mathcal{L}_1^{c}$ independently. These contact interactions are independent and give no contributions to the decay widths of the charmonium. This means that their effects can be absorbed by the models which are used to calculate the energy spectrum and do not include the annihilation effects. Here we are only interested in the mass shifts due to the decay modes, then we only focus on the contributions including the imaginary parts due to the loop calculation and neglect the terms $c^{(a,c)}_{0,poly}$. The choices of $g_{10,11,12}$ which can cancel all the polynomial contributions in $c_{0}^{(a,c)}$ can be got directly.

From Eq. (20), one can easily get the imaginary parts as follows:

$$\begin{aligned} \text {Im}[c^{(a)}_0]= & {} \frac{256\sqrt{s(s-4m_D^2)}}{s} \theta (s-4m_D^2), \nonumber \\ \text {Im}[c^{(b)}_0]= & {} 0, \nonumber \\ \text {Im}[c^{(c)}_0]= & {} \frac{64[(s-2m_{D^*}^2)^2+8m_{D^*}^4] \sqrt{s(s-4m_{D^*}^2)}}{s m_{D^*}^4}\nonumber \\&\times \theta (s-4m_{D^*}^2), \nonumber \\ \text {Im}[c^{(d)}_0]= & {} 0. \end{aligned}$$

(22)

Matching the amplitude with the corresponding amplitude in quantum mechanism with a perturbatively potential, one has

$$\begin{aligned} {{\mathcal {M}}}(^3P_J)= -\langle H(^3P_J)|V_{eff}|H(^3P_J)\rangle . \end{aligned}$$

(23)

Finally the decay widths of $^3P_0$ to DD, $DD^*$ and $D^*D^*$ in the leading order are expressed as follows:

$$\begin{aligned} \Gamma (^3P_0\rightarrow DD )= & {} 2 \text {Im}[{{\mathcal {M}}}^{(a)}(^3P_0)]\nonumber \\= & {} {\frac{27g_a^2}{8\pi ^2}}N_iN_fm_c^4 \text {Im}[c_{0}^{(a)}]|R_1^{(1)}(0)|^2, \nonumber \\ \Gamma (^3P_0\rightarrow DD^*)= & {} 2 \text {Im}[{{\mathcal {M}}}^{(b)}(^3P_0)] =0, \nonumber \\ \Gamma (^3P_0\rightarrow D^*D^* )= & {} 2 \text {Im}[{\mathcal {M}}^{(c)}(^3P_0)]\nonumber \\= & {} {\frac{27g_c^2}{8\pi ^2}}N_iN_fm_c^4 \text {Im}[c_{0}^{(c)}]|R_1^{(1)}(0)|^2,\nonumber \\ \end{aligned}$$

(24)

where we have used the relation

$$\begin{aligned} \int \frac{dp}{(2\pi )^{3/2}}\phi (p) p^{2n+3} = (-1)^n \frac{2n+3}{4\pi } R_1^{(2n+1)}(|\mathbf{r} |)\Big |_{|\mathbf{r} |=0}.\nonumber \\ \end{aligned}$$

(25)

The corresponding full mass shift labeled as $\Delta m(^3P_0)$ is expressed as

$$\begin{aligned} \Delta m(^3P_0)= & {} -\text {Re}[{{\mathcal {M}}}^{(a+b+c)}(^3P_0)] \nonumber \\= & {} -\frac{\text {Re}[{\overline{c}}_{0}^{(a)}]}{2\text {Im}[c_{0}^{(a)}]} \Gamma (^3P_0\rightarrow DD)\nonumber \\&-\frac{\text {Re}[{\overline{c}}_{0}^{(c)}]}{2\text {Im}[c_{0}^{(c)}]} \Gamma (^3P_0\rightarrow D^*D^*), \end{aligned}$$

(26)

where ${\overline{c}}_{0}^{(a,c)}=c_{0}^{(a,c)}-c_{0,poly}^{(a,c)}$.

3.2 The energy shift of $^3P_1$ state

In the $^3P_1$ channel, we have the following results

$$\begin{aligned} c^{(a)}_1= & {} 0, \nonumber \\ c^{(b)}_1= & {} c^{(b)}_{1,poly}+\frac{128}{9\pi s^2m_{D*}^2}A(A^2+12sm_{D*}^2)\ln \nonumber \\&\times \left[ \frac{A-s+m_D^2+m_{D^*}^2}{2m_Dm_{D^*}}+i\varepsilon \right] , \nonumber \\ c^{(c)}_1= & {} 0, \nonumber \\ c^{(d)}_1= & {} c^{(d)}_{1,poly}, \end{aligned}$$

(27)

with

$$\begin{aligned} A= & {} \sqrt{\big [ s-(m_D-m_{D*})^2\big ] \big [s-(m_D+m_{D*})^2\big ]}. \end{aligned}$$

(28)

The polynomial terms are expressed as

$$\begin{aligned} c^{(b)}_{1,poly}= & {} \sum \limits _{n=-2}^{1}s^nc_{1;n}^{(b)},\nonumber \\ c^{(d)}_{1,poly}= & {} -\frac{512}{3\pi ^3}(g_{20}+g_{21}s), \end{aligned}$$

(29)

with

$$\begin{aligned} c_{1;-2}^{(b)}= & {} \frac{64}{9\pi m_{D^*}^2}(m_{D^*}^2-m_D^2)^3\ln \frac{m_D^2}{m_{D^*}^2},\nonumber \\ c_{1;-1}^{(b)}= & {} \frac{64}{9\pi m_{D^*}^2}(m_{D^*}^2-m_D^2) \left[ 2(m_{D^*}^2-m_D^2)\right. \nonumber \\&\left. +3(3m_{D^*}^2-m_D^2)\ln \frac{m_D^2}{m_{D^*}^2}\right] ,\nonumber \\ c_{1;0}^{(b)}= & {} -\frac{64}{9\pi m_{D^*}^2} \left[ \left( 2+\frac{6}{{\overline{\epsilon }}_{\text {UV}}} +6\ln \frac{\mu ^2}{m_D^2}\right) m_D^2\right. \nonumber \\&-2\left( 23+\frac{9}{{\overline{\epsilon }}_{\text {UV}}} +9\ln \frac{\mu ^2}{m_{D^*}^2}\right) m_{D^*}^2\nonumber \\&\left. +3(m_D^2-3m_{D^*}^2) \ln \frac{m_D^2}{m_{D^*}^2}\right] ,\nonumber \\ c_{1;1}^{(b)}= & {} \frac{64}{27\pi m_{D^*}^2} \left( 4+\frac{6}{{\overline{\epsilon }}_{\text {UV}}} +6\ln \frac{\mu ^2}{m_D^2}+3\ln \frac{m_D^2}{m_{D^*}^2}\right) .\nonumber \\ \end{aligned}$$

(30)

At first glance, this property is very different from that in the $^3P_0$ channel due to the nonzero values of $c_{1,-2}$ and $c_{1,-1}$. While actually when taking the nonrelativistic approximation $m_D\approx m_{D^*}$, one has $c_{1;-2},c_{1;-2} \approx 0$, this means that these contributions are very small in the nonrelativistic approximation and can be neglected. The numerical calculation also shows such property and we neglect these two terms in the following.

Similarly, the terms $c^{(b,d)}_{1,poly}$ can be neglected when aiming to discuss the contributions from the annihilation effects. The corresponding imaginary parts can be expressed as

$$\begin{aligned} \text {Im}[c_1^{(a)}]= & {} 0,\nonumber \\ \text {Im}[c_1^{(b)}]= & {} \frac{128}{9 s^2m_{D*}^2}A(A^2+12sm_{D*}^2) \theta \big (s-(m_D+m_{D^*})^2\big ),\nonumber \\ \text {Im}[c_1^{(c)}]= & {} 0,\nonumber \\ \text {Im}[c_1^{(d)}]= & {} 0. \end{aligned}$$

(31)

In the leading order, the decay width of $^3P_1$ to DD, $DD^*$ and $D^*D^*$ are expressed as

$$\begin{aligned} \Gamma (^3P_1\rightarrow DD )= & {} 2 \text {Im}[{{\mathcal {M}}}^{(a)}(^3P_1)] = 0, \nonumber \\ \Gamma (^3P_1\rightarrow DD^* )= & {} 2 \text {Im}[{{\mathcal {M}}}^{(b)}(^3P_1)]\nonumber \\= & {} {\frac{27g_b^2}{8\pi ^2}}N_iN_fm_c^4 \text {Im}[c_{1}^{(b)}]|R_1^{(1)}(0)|^2,\nonumber \\ \Gamma (^3P_1\rightarrow D^*D^* )= & {} 2 \text {Im}[\mathcal{M}^{(c)}(^3P_1)] = 0, \end{aligned}$$

(32)

and the corresponding full mass shift labeled as $\Delta m(^3P_1)$ is expressed as

$$\begin{aligned} \Delta m(^3P_1)= & {} -\text {Re}[{{\mathcal {M}}}^{(a+b+c)}(^3P_1)]\nonumber \\= & {} -\frac{\text {Re}[{\overline{c}}_{1}^{(b)}]}{2\text {Im}[c_{1}^{(b)}]}\Gamma (^3P_1\rightarrow DD^*), \end{aligned}$$

(33)

where ${\overline{c}}_{1}^{(b)}=c_{1}^{(b)}-c_{1,poly}^{(b)}$.

3.3 The energy shift of $^3P_2$ state

For $^3P_2$ state, we get

$$\begin{aligned} c^{(a,b,c,d)}_2= & {} 0. \end{aligned}$$

(34)

These results mean that the decay widths $\Gamma (^3P_2\rightarrow DD,DD^*,D^*D^* )$ are exact zero and there are no mass shifts for the $^3P_2$ states in the leading order. This is a strong property which can be tested by the experiments and be used to judge whether a state is a pure $^3P_2$ heavy quarkonium or not.

3.4 Property of the analytic results

Comparing our results with those given by the $^3P_0$ model in Ref. [55] and the Friedrichs model in Ref. [56], one can find that both the two methods give the zero results for $c{\overline{c}}(^3P_0)\rightarrow DD^*$ and $c{\overline{c}}(^3P_1)\rightarrow DD$. The contributions of $c{\overline{c}}(^3P_1)\rightarrow D^*D^*$ and $c{\overline{c}}(^3P_2)\rightarrow DD,DD^*,D^*D^*$ in Ref. [55] are nonzero, and even in the same order with the contribution in $c{\overline{c}}(^3P_1)\rightarrow DD^*$. While in our results, $c{\overline{c}}(^3P_1)\rightarrow D^*D^*$ and $c{\overline{c}}(^3P_2)\rightarrow DD,DD^*, D^*D^*$ are exact zero.

To understand these exact zero results in our method, one can consider the decay widths in the hadronic level at first. The most general form for the amplitudes of $c{\overline{c}}(^3P_J)\rightarrow DD,DD^*,D^*D^*$ with Lorentz invariance can be expressed as follows:

$$\begin{aligned} {{\mathcal {M}}}(^3P_0\rightarrow DD)= & {} \{1\}, \nonumber \\ {{\mathcal {M}}}(^3P_0\rightarrow DD^*)= & {} \{P_{i}^{\mu }\}\epsilon ^*_{\mu }, \nonumber \\ {{\mathcal {M}}}(^3P_0\rightarrow D^*D^*)= & {} \{g^{\mu \nu }, P_{i}^{\mu } P_{j}^{\nu },\epsilon ^{\mu \nu \rho \omega }P_{i\rho }P_{j\omega }\}\epsilon ^*_{1\mu }\epsilon ^*_{2\nu },\nonumber \\ \end{aligned}$$

(35)

and

$$\begin{aligned}&{{\mathcal {M}}}(^3P_1\rightarrow DD)= \{P_{i}^{\mu }\}X_{\mu }, \nonumber \\&{{\mathcal {M}}}(^3P_1\rightarrow DD^*)=\{g^{\mu \nu }, P_{i}^{\mu } P_{j}^{\nu },\epsilon ^{\mu \nu \rho \omega } P_{i\rho }P_{j\omega }\}X_{\mu }\epsilon ^*_{\nu }, \nonumber \\&{{\mathcal {M}}}(^3P_1\rightarrow D^*D^*)\nonumber \\&\quad =\{ P_{i}^{\mu }g^{\nu \rho }, P_{i}^{\nu }g^{\mu \rho }, P_{i}^{\rho }g^{\mu \nu }, P_{i}^{\mu }P_{j}^{\nu }P_k^{\rho }, \nonumber \\&\qquad P_{i}^{\omega }P_{j}^{\sigma }P_k^{\rho }\varepsilon ^{\mu \nu \omega \sigma }, P_{i}^{\omega }P_{j}^{\sigma }P_k^{\nu }\varepsilon ^{\mu \rho \omega \sigma }\} X_{\mu }\epsilon ^*_{1\nu }\epsilon ^*_{2\rho }, \end{aligned}$$

(36)

and

$$\begin{aligned}&{{\mathcal {M}}}(^3P_2\rightarrow DD) = \{P_{i}^{\mu }P_j^{\nu }\}X_{\mu \nu }, \nonumber \\&{{\mathcal {M}}}(^3P_2\rightarrow DD^*)\nonumber \\&\quad =\{P_{i}^{\mu }g^{\nu \rho }, P_{i}^{\nu }g^{\mu \rho }, P_{i}^{\mu }P_{j}^{\nu }P_k^{\rho }, P_{i}^{\omega }P_{j}^{\sigma } P_k^{\nu }\varepsilon ^{\mu \rho \omega \sigma },\nonumber \\&\qquad P_{i}^{\omega }p_{j}^{\sigma }P_k^{\mu }\varepsilon ^{\nu \rho \omega \sigma }\} X_{\mu \nu }\epsilon ^*_{\rho }, \nonumber \\&{{\mathcal {M}}}(^3P_2\rightarrow D^*D^*)\nonumber \\&\quad =\{g^{\mu \rho }g^{\nu \omega },g^{\mu \rho }P_i^{\nu }P_j^{\omega }, g^{\mu \omega }P_i^{\nu }P_j^{\rho },\text {terms with more } P_i\}\nonumber \\&\qquad \times X_{\mu \nu }\epsilon ^*_{1\rho }\epsilon ^*_{2\omega }, \end{aligned}$$

(37)

where $P_i$ are the momenta of $^3P_J$ meson, D meson and $D^*$ meson, $\epsilon _{\mu },X_{\mu }$ and $X_{\mu \nu }$ are the polarization vectors of $D^*$ meson,$^3P_1$ and $^3P_2$ meson, respectively, $\{\ldots \}$ denotes to the linear combination of the terms in the braces. The properties that $X^{\mu }_{\mu }=0$, $X_{\mu \nu }=X_{\nu \mu }$, and $\epsilon _{1\mu }^{*}$ and $\epsilon _{2\nu }^{*}$ are symmetrical have been used.

Due to the P invariance, one can directly get that $\mathcal{M}(^3P_0\rightarrow DD^*)=0$ and ${{\mathcal {M}}}(^3P_1\rightarrow DD)=0$ since the terms in the braces break the P invariance. In principle, other amplitudes are nonzero due to the symmetry.

When taking the nonrelativistic expansion on the three-momenta of the $D,D^*$ mesons, in the leading order all $P_i$ are parallel. This property means that ${{\mathcal {M}}}(^3P_1\rightarrow D^*D^*,~^3P_2\rightarrow DD,DD^*)=0$ in the leading order since $P_i^{\mu }\epsilon ^*_{\mu },P_i^{\mu } X_{\mu }, P_i^{\mu }X_{\mu \nu }$ are zero when all $P_i$ are parallel. Mathematically, all $P_i^{\mu }\epsilon ^*_{\mu },P_i^{\mu } X_{\mu }, P_i^{\mu }X_{\mu \nu }$ are zero or small comparing with the mass $m_D$.

For ${{\mathcal {M}}}(^3P_2\rightarrow D^*D^*)$, when only taking the nonrelativistic expansion on the three-momenta of the $D,D^*$ mesons, there is a nonzero contribution like $X_{\mu \nu }\epsilon _{1}^{*\mu }\epsilon _{2}^{*\nu }$ in the leading order. But when one considers the quark structure of the $^3P_2$ meson and also takes the nonrelativistic expansion on the three-momenta of the quarks, then the interactions at the quark level in the leading order is ${{\mathcal {L}}}_3$ in Eq. (1). This interaction means that the final amplitude is like $\{P_{i}^{\mu }P_j^{\nu }\}X_{\mu \nu }\epsilon _{1\rho }^*\epsilon _{2}^{*\rho }$ after considering both the nonrelativistic expansion on the $D^*$ mesons and the c quarks. This is just zero in the leading order.

In summary, we get the following properties:

• (1) Due to the P invariance, ${{\mathcal {M}}}(^3P_0\rightarrow DD^*$, $^3P_1\rightarrow DD)=0$.

• (2) Due to the nonrelativistic expansion on the three-momenta of the $D,D^*$ mesons, ${{\mathcal {M}}}(^3P_1\rightarrow D^*D^*, ^3P_2\rightarrow DD,DD^*)=0$ in the leading order.

• (3) Due to the nonrelativistic expansion on the three-momenta of the $D,D^*$ mesons and the quarks in $^3P_2$ meson, $\mathcal{M}(^3P_2\rightarrow D^*D^*)=0$ in the leading order.

In the practical calculation, we also take the following next leading order interaction as example

$$\begin{aligned} {{\mathcal {L}}}_{NLO}= & {} g(\partial _{\mu }{\overline{\psi }})(\partial _{\nu }\psi ) A^{\mu } A^{\nu }, \end{aligned}$$

(38)

and we find nonzero result for ${{\mathcal {M}}}(^3P_2\rightarrow D^*D^*)$ in the next leading order.

From the above analysis, one can see that when go to beyond the leading order of nonrelativistic expansion, the amplitudes $\mathcal{M}(^3P_1\rightarrow D^*D^*$, $^3P_2\rightarrow DD,DD^*,D^*D^*)$ will not be zero. The calculation in Ref. [55] is based on the $^3P_0$ model where a light quark pair is dynamically produced in the vacuum and the nonrelativistic wave functions of the mesons are used to estimate the contributions. This means that in Ref. [55] the three-momenta of the quarks and the $D,D^*$ mesons are not expanded order by order and their momenta’ distributions are considered via some wave functions and the quark pair creation mechanism. This is very different from our method which is based on the general model independent interactions under the nonrelativistic expansion. In our calculation, all the dynamics of the light quark and $D,D^*$ meason is absorbed by the coupling constants in the leading order of the nonrelativistic expansion. On another hand, we only consider the contributions due to the annihilation effects and neglect the polynomial contribution since the latter is uncertain.

4 Numerical results and discussion

To show the properties of the above analytic results more clearly, we present some numerical results in this section. Firstly we want to emphasize that the absolute values of $\text {Re}[{\overline{c}}^{(a,b,c)}_{J}]$ and $\text {Im}[c^{(a,b,c)}_{J}]$ can not determine the physical decay widths and the mass shifts directly, since there are some global unknown constant factors. But the ratios of the mass shifts and the decay widths $-\text {Re}[{\overline{c}}^{(a,b,c)}_{J}]/2\text {Im}[c^{(a,b,c)}_{J}]$ are model independent. This means that if the decay widths are measured experimentally, the corresponding corrections to the masses of the heavy quarkonia can be got directly.

In Fig. 2 the dependence of $\text {Im}[c^{(a,b,c)}_{J}]$, $\text {Re}[{\overline{c}}^{(a,b,c)}_{J}]$ and their ratios on $\sqrt{s}$ is presented, respectively, where we take $m_D=1.87~\text {GeV}$ and $m_{D^*}=2.01~\text {GeV}$ as inputs.

Fig. 2

Numerical results for $\text {Im}[c^{(a,b,c)}_{J}]$ vs. $\sqrt{s}$, $\text {Re}[{\overline{c}}^{(a,b,c)}_{J}]$ vs. $\sqrt{s}$ and $-2\text {Im}[c^{(a,b,c)}_{J}]/\text {Re}[{\overline{c}}^{(a,b,c)}_{J}]$ vs. $\sqrt{s}$. The sub figures (a, b, c) are corresponding to $\text {Im}[c^{(a,b,c)}_{J}]$ and $\text {Re}[{\overline{c}}^{(a,b,c)}_{J}]$ vs. $\sqrt{s}$, respectively. The sub figure (d) shows the results for $-2\text {Im}[c^{(X)}_{J}]/\text {Re}[{\overline{c}}^{(X)}_{J}]$ vs. $\sqrt{s}$

The numerical results presented in Fig. 2 show four interesting properties:

• (1) The real parts $\text {Re}[{\overline{c}}^{(a,c)}_{0}]$ and $\text {Re}[{\overline{c}}^{(b)}_{1}]$ which are represented by the solid black curves are always negative. This means that after considering the annihilation effects, the masses of the $^3P_{0,1}$ states move up and the masses of $^3P_2$ states do not move.

• (2) When $\sqrt{s}$ is on the threshold of $DD,DD^*$ or $D^*D^*$ the corresponding mass shifts are exact zero.

• (3) When $\sqrt{s}$ is above the threshold, the mass shifts are much smaller than the corresponding decay widths, the largest mass shift is about 1/5 of the corresponding decay width when $\sqrt{s}\approx 4.5$ GeV which is much larger than the threshold. This property gives a strong constrain on the mass shifts to all the $^3P_{0,1}$ states.

• (4) When $\sqrt{s}$ is below the mass-shell, although the decay widths are exact zero, but the mass-shifts are still nonzero and the dependence of $\text {Re}[{\overline{c}}^{(a,b,c)}_{J}]$ vs. $\sqrt{s}$ shows non-trivial property.

To show the non-trivial dependence of $\text {Re}[{\overline{c}}^{(a,b,c)}_{J}]$ vs. $\sqrt{s}$ more clearly, we present the dependence of $\text {Re}[{\overline{c}}^{(a,b,c)}_{J}(s)]/\text {Re}[{\overline{c}}^{(a,b,c)}_{J}(s_0)]|$ vs. $\sqrt{s}$ with $s_0=3$ GeV in Fig. 3. The curves in Fig. 3 clearly show that when $\sqrt{s}$ increases from 3 to 4.5 GeV the ratios of the mass shifts decrease from 1 to 0 at first and then increase from 0 to 0.5. For the states with the same quantum number, it means that the corresponding mass shifts are non-linear and can not be absorbed by some constants.

Fig. 3

Numerical results for the dependence of $\text {Re}[{\overline{c}}^{(a,b,c)}_{J}(s)]/\text {Re}[{\overline{c}}^{(a,b,c)}_{J}(s_0)]$ vs. $\sqrt{s}$

Experimentally, up to now there are still no definite results for the branch decay widths $\Gamma (^3P_{0,1}\rightarrow DD,DD^*, D^*D^*)$ [57,58,59,60,61,62], this makes it difficult to determine the mass shifts certainly. The experiments reported that the decay widths $\Gamma (X(3915),\chi _{c2}(3930)\rightarrow DD,DD^*, D^*D^*)$ are seen. By our calculation, we expect that the decay widths $\Gamma (^3P_{2}\rightarrow DD,DD^*, D^*D^*)$ are zero in the leading order which suggests that the decay widths $\Gamma (^3P_{2}\rightarrow DD,DD^*, D^*D^*)$ should be much smaller than $\Gamma (^3P_{0}\rightarrow DD, D^*D^*)$ and $\Gamma (^3P_{1}\rightarrow DD^*)$. A relative larger decay width of a resonance to $DD,DD^*, D^*D^*$ suggests that it maybe is not a pure $c{\bar{c}}(^3P_J)$ state. These properties are more reliable in the b quark part and can be tested by the further precise experiments. Furthermore, the similar discussion can be extended to the S wave states and compared with the similar studies in Ref. [63, 64].

In summary, the nonrelativistic asymptotic behaviors of the transitions $c{\bar{c}}(^3P_J) \rightarrow DD,DD^*,D^*D^* \rightarrow c{\bar{c}}(^3P_J)$ with $J=0, 1, 2$ are discussed. We find that the decay widths $\Gamma (^3P_{0}\rightarrow DD^*),\Gamma (^3P_{1}\rightarrow DD,D^*D^*)$ and $\Gamma (^3P_{2}\rightarrow DD,DD^*,D^*D^*)$ are exact zero in the leading order of nonrelativistic expansion. For other channels, the ratios between the branch decay widths and the mass shifts are larger than 5 when the center-of-mass energy is above the threshold. When below the threshold, the mass shifts are dependent on the center-of-mass energy nontrivially and can not be absorbed by a constant.

Appendix: The FIESTA integrations

We define the following functions to refer to the results after summing the spins and integrating the angles:

$$\begin{aligned}&P(J,X,n) \overset{def}{=} \sum _{s_i,s_f} \langle JJ_Z|1s_f;1m_f\rangle \langle JJ_Z|1s_i;1m_i\rangle \nonumber \\&\qquad \times \int d\Omega _\mathbf{q _i} d\Omega _\mathbf{q _f} Y_{1m_i}(\Omega _\mathbf{q _i}) Y_{1m_f}^{*} (\Omega _\mathbf{q _f})\nonumber \\&\qquad \times ({{\hat{q}}}_i\cdot {{\hat{q}}}_f)^n X, \end{aligned}$$

(39)

where X are some functions dependent on ${{\hat{q}}}_{i}, {{\hat{q}}}_{f}, \epsilon _p(s_i)$, and $\epsilon _p^*(s_f)$ with ${{\hat{q}}}_{i,f} \overset{def}{=} q_{i,f}/|\mathbf{q} _{i,f}|$, P(J, X, n) are not dependent on $J_z$. When $J=0,1,2$ and $n=0,1$, we have

$$\begin{aligned}&P(J,\epsilon _p(s_i)\cdot \epsilon _p^*(s_f),1) = \frac{4\pi }{3},\nonumber \\&\quad P(0,\epsilon _p(s_i)\cdot \hat{q_i} \ \epsilon _p^*(s_f)\cdot \hat{q_f},0) = 4\pi , \nonumber \\&P(0,\epsilon _p(s_i)\cdot \hat{q_f} \ \epsilon _p^*(s_f)\cdot \hat{q_i},0) = \frac{4\pi }{3},\nonumber \\&\quad P(1,\epsilon _p(s_i)\cdot \hat{q_f} \ \epsilon _p^*(s_f)\cdot \hat{q_i},0) = -\frac{4\pi }{3}, \nonumber \\&P(2,\epsilon (s_i)\cdot \hat{q_f} \ \epsilon ^*(s_f)\cdot \hat{q_i},0)= \frac{4\pi }{3}, \end{aligned}$$

(40)

and others are zero.