Novel implementation of the multipole expansion to quarkonium hadronic transitions

TL;DR

This paper introduces a novel implementation of the multipole expansion for calculating hadronic transitions in heavy quarkonium, avoiding the twist expansion and using a different counting scheme based on hybrid interpolating fields.

Contribution

It proposes an alternative counting method for hadronic transitions that enables the use of chiral theorems without relying on the twist expansion, providing new insights into quarkonium decay processes.

Findings

Same predictions for two-pion decay rates as the twist expansion approach.

Effective theory implementation using potential NRQCD.

Phenomenological impact of the new counting scheme explored.

Abstract

We compute hadronic transitions between heavy quarkonium states with two, or one, pion/eta particles in the final state. We use the multipole expansion but not the twist expansion. The latter cannot be justified for the energy release of hadronic transitions between heavy quarkonium states with different principal quantum numbers. Instead, we use a counting based on the dimension of the interpolating field of the hybrid. This alternative counting allows us to still use chiral low-energy theorems to compute the pion production by local gluonic operators. We explore the phenomenological impact of this counting. Remarkably enough, for the two-pion transitions, we obtain the same predictions for the normalized differential decay rate as those obtained assuming the twist expansion. We implement this computational scheme using the hadronic representation of the effective theory potential NRQCD. We assume that the inverse Bohr radius of the heavy quarkonium is much larger than $\Lambda_{\rm QCD}$ but do not impose any constraint on the relative size of $\Lambda_{\rm QCD}$ and the typical kinetic energy of the bound state.

Figures

Tables

Table 1. Table 1: Transformation properties of the heavy quarkonium and gluonic fields, and the projection vectors, under discrete symmetries. The octet field O O {\rm O} has the same transformation properties as S S {\rm S} . The Ψ k i subscript superscript Ψ 𝑖 𝑘 \Psi^{i}_{k} transform as O O {\rm O} combined with the k P ​ C superscript 𝑘 𝑃 𝐶 k^{PC} of the gluelump. The transformations of the projected fields can be obtained by further adding those of the projection vectors; however, these are actually not relevant for the construction of the Lagrangian since the projection vectors always appear in pairs, one explicit in the operator and another implicit in Ψ k ​ λ = 𝒓 ^ λ ⋅ 𝚿 k subscript Ψ 𝑘 𝜆 ⋅ subscript ^ 𝒓 𝜆 subscript 𝚿 𝑘 \Psi_{k\lambda}=\hat{\bm{r}}_{\lambda}\cdot{\bf\Psi}_{k} . For this reason, we give the transformation properties of the unprojected fields Ψ k i subscript superscript Ψ 𝑖 𝑘 \Psi^{i}_{k} . Note that the difference in the transformation of 𝒓 ^ λ subscript ^ 𝒓 𝜆 \hat{\bm{r}}_{\lambda} with respect to Ref.  Berwein  et al. ( 2015 ) are due to the different definition of 𝒓 ^ + subscript ^ 𝒓 \hat{\bm{r}}_{+} .

	P	T	C
$S(t,𝒓,𝑹)$	$-S(t,-𝒓,-𝑹)$	${\sigma }_{2}S(-t,𝒓,𝑹){\sigma }_2$	${\sigma }_2{S}^{\top }(t,-𝒓,𝑹){\sigma }_2$
${\Psi }_{1^{+-}}^i(t,𝒓,𝑹)$	$-{\Psi }_{1^{+-}}^i(t,-𝒓,-𝑹)$	$-{\sigma }_2{\Psi }_{1^{+-}}^i(-t,𝒓,𝑹){\sigma }_2$	$-{\sigma }_2{\left({\Psi }_{1^{+-}}^i\right)}^{\top }(t,-𝒓,𝑹){\sigma }_2$
${\Psi }_{1^{--}}^i(t,𝒓,𝑹)$	${\Psi }_{1^{--}}^i(t,-𝒓,-𝑹)$	${\sigma }_2{\Psi }_{1^{--}}^i(-t,𝒓,𝑹){\sigma }_2$	$-{\sigma }_2{\left({\Psi }_{1^{--}}^i\right)}^{\top }(t,-𝒓,𝑹){\sigma }_2$
$𝑬(t,𝑹)$	$-𝑬(t,-𝑹)$	$𝑬(-t,𝑹)$	$-{𝑬}^{\top }(t,𝑹)$
$𝑩(t,𝑹)$	$𝑩(t,-𝑹)$	$-𝑩(-t,𝑹)$	$-{𝑩}^{\top }(t,𝑹)$
${\widehat{𝒓}}_{\lambda }$	$-{\widehat{𝒓}}_{-\lambda }$	${\widehat{𝒓}}_{\lambda }^{\ast }$	$-{\widehat{𝒓}}_{-\lambda }$

Table 2. Table 2: Transformation properties of the basic chiral building blocks under discrete symmetries.

	$u$	$u_{\mu }$	${\chi }_{\pm }$	$D_{\mu }$
P	$u^{†}$	$-u^{\mu }$	$\pm {\chi }_{\pm }$	$D^{\mu }$
C	$u^{\top }$	$u_{\mu }^{\top }$	${\chi }_{\pm }^{\top }$	$D_{\mu }^{\top }$
T	$u$	$u^{\mu }$	${\chi }_{\pm }$	$D^{\mu }$
h.c	$u^{†}$	$u_{\mu }$	$\pm {\chi }_{\pm }$	$D_{\mu }$

Table 3. Table 3: Quarkonium spectrum obtained from the potential in Eq. ( 102 ) fitted to the lattice data Ref.  Bazavov  et al. ( 2018 ) . The origin of energies is chosen in order to reproduce the spin average of the ground state. All dimension-full entries are in GeV.

nl	$J^{PC}$	$c\overline{c}$			$b\overline{b}$
		$m_H$	$⟨1/r⟩$	$E_{kin}$	$m_H$	$⟨1/r⟩$	$E_{kin}$
1S	$\{0^{-+},1^{--}\}$	3.068	0.66	0.31	9.442	1.15	0.28
2S	$\{0^{-+},1^{--}\}$	3.686	0.43	0.48	9.916	0.67	0.36
1P	$\{1^{+-},{(0,1,2)}^{++}\}$	3.464	0.37	0.39	9.779	0.59	0.30
1D	$\{2^{-+},{(1,2,3)}^{--}\}$	3.775	0.27	0.48	10.009	0.42	0.35

Table 4. Table 4: Quarkonium spectrum obtained from the potential in Eq. ( 15 ) fitted to the lattice data Ref.  Bazavov  et al. ( 2018 ) . The origin of energies is chosen in order to reproduce the spin average of the ground state. All dimension-full entries are in GeV.

nl	$J^{PC}$	$c\overline{c}$			$b\overline{b}$
		$m_H$	$⟨1/r⟩$	$E_{kin}$	$m_H$	$⟨1/r⟩$	$E_{kin}$
1S	$\{0^{-+},1^{--}\}$	3.068	0.73	0.38	9.442	1.36	0.38
2S	$\{0^{-+},1^{--}\}$	3.847	0.51	0.67	10.015	0.76	0.45
1P	$\{1^{+-},{(0,1,2)}^{++}\}$	3.553	0.42	0.50	9.873	0.64	0.35
1D	$\{2^{-+},{(1,2,3)}^{--}\}$	3.949	0.32	0.65	10.145	0.46	0.41

Table 5. Table 5: Heavy hybrid spectrum for k P ​ C = 1 − − superscript 𝑘 𝑃 𝐶 superscript 1 absent k^{PC}=1^{--} obtained from the fit of Eq. ( 102 ) to the lattice data of Ref.  Juge  et al. ( 2003 ) . All dimension-full entries are in GeV. P Π = ∫ d 3 ​ r ​ | ψ + ⁣ / − ( m ) | 2 subscript 𝑃 Π superscript 𝑑 3 𝑟 superscript superscript subscript 𝜓 absent 𝑚 2 P_{\Pi}=\int d^{3}r|\psi_{+/-}^{(m)}|^{2} is the integral over the square of the wave function associated with the Π g subscript Π 𝑔 \Pi_{g} potential. It can be interpreted as the probability to find the hybrid in a Π g subscript Π 𝑔 \Pi_{g} configuration, thus it gives a measure of the mixing effects.

multiplet	l	$J^{PC}$	${\Lambda }_{\eta }^{\sigma }$	$c\overline{c}$				$b\overline{b}$
				$m_H$	$⟨1/r⟩$	$E_{kin}$	$P_{\Pi }$	$m_H$	$⟨1/r⟩$	$E_{kin}$	$P_{\Pi }$
$H_1\left(1𝒫\right)$	1	$\{1^{+-},{(0,1,2)}^{++}\}$	${\Sigma }_g^{+\prime }$-${\Pi }_g$	4.51	0.42	0.15	0.57	11.14	0.55	0.09	0.44
$H_1^{\prime }\left(2𝒫\right)$				4.86	0.27	0.33	0.20	11.33	0.40	0.19	0.27
$H_2\left(1𝒫\right)$		$\{1^{-+},{(0,1,2)}^{--}\}$	${\Pi }_g$	4.74	0.27	0.22	1.00	11.30	0.38	0.13	1.00
$H_2^{\prime }\left(2𝒫\right)$				5.01	0.24	0.37	1.00	11.51	0.34	0.26	1.00
$H_3\left(1𝒮\right)$	0	$\{0^{-+},1^{--}\}$	${\Sigma }_g^{+\prime }$	4.63	0.29	0.25	0.00	11.19	0.40	0.15	0.00
$H_3^{\prime }\left(2𝒮\right)$				5.01	0.25	0.40	0.00	11.41	0.35	0.24	0.00
$H_4\left(1𝒟\right)$	2	$\{2^{-+},{(1,2,3)}^{--}\}$	${\Sigma }_g^{+\prime }$-${\Pi }_g$	4.70	0.28	0.24	0.50	11.24	0.38	0.14	0.39
$H_4^{\prime }\left(2𝒟\right)$				5.06	0.23	0.42	0.24	11.44	0.32	0.25	0.26
$H_5\left(1𝒟\right)$		$\{2^{+-},{(1,2,3)}^{++}\}$	${\Pi }_g$	4.92	0.22	0.32	1.00	11.40	0.30	0.18	1.00
$H_5^{\prime }\left(2𝒟\right)$				5.27	0.19	0.47	1.00	11.61	0.28	0.28	1.00

Table 6. Table 6: Heavy hybrid spectrum for k P ​ C = 1 − − superscript 𝑘 𝑃 𝐶 superscript 1 absent k^{PC}=1^{--} obtained from the fit of Eq. ( 15 ) to the lattice data of Ref.  Juge  et al. ( 2003 ) . All dimension-full entries are in GeV. P Π = ∫ d 3 ​ r ​ | ψ + ⁣ / − ( m ) | 2 subscript 𝑃 Π superscript 𝑑 3 𝑟 superscript superscript subscript 𝜓 absent 𝑚 2 P_{\Pi}=\int d^{3}r|\psi_{+/-}^{(m)}|^{2} is the integral over the square of the wave function associated with the Π g subscript Π 𝑔 \Pi_{g} potential. It can be interpreted as the probability to find the hybrid in a Π g subscript Π 𝑔 \Pi_{g} configuration, thus it gives a measure of the mixing effects.

multiplet	l	$J^{PC}$	${\Lambda }_{\eta }^{\sigma }$	$c\overline{c}$				$b\overline{b}$
				$m_H$	$⟨1/r⟩$	$E_{kin}$	$P_{\Pi }$	$m_H$	$⟨1/r⟩$	$E_{kin}$	$P_{\Pi }$
$H_1\left(1𝒫\right)$	1	$\{1^{+-},{(0,1,2)}^{++}\}$	${\Sigma }_g^{+\prime }$-${\Pi }_g$	4.54	0.44	0.15	0.59	11.18	0.57	0.08	0.58
$H_1^{\prime }\left(2𝒫\right)$				4.87	0.28	0.32	0.20	11.36	0.37	0.18	0.19
$H_2\left(1𝒫\right)$		$\{1^{-+},{(0,1,2)}^{--}\}$	${\Pi }_g$	4.79	0.31	0.28	1.00	11.32	0.41	0.15	1.00
$H_2^{\prime }\left(2𝒫\right)$				5.24	0.28	0.51	1.00	11.57	0.37	0.28	1.00
$H_3\left(1𝒮\right)$	0	$\{0^{-+},1^{--}\}$	${\Sigma }_g^{+\prime }$	4.68	0.28	0.23	0.00	11.26	0.37	0.12	0.00
$H_3^{\prime }\left(2𝒮\right)$				5.05	0.25	0.41	0.00	11.46	0.33	0.23	0.00
$H_4\left(1𝒟\right)$	2	$\{2^{-+},{(1,2,3)}^{--}\}$	${\Sigma }_g^{+\prime }$-${\Pi }_g$	4.74	0.29	0.26	0.46	11.29	0.39	0.14	0.46
$H_4^{\prime }\left(2𝒟\right)$				5.06	0.23	0.42	0.25	11.47	0.30	0.23	0.25
$H_5\left(1𝒟\right)$		$\{2^{+-},{(1,2,3)}^{++}\}$	${\Pi }_g$	5.02	0.25	0.39	1.00	11.44	0.33	0.22	1.00
$H_5^{\prime }\left(2𝒟\right)$				5.47	0.23	0.62	1.00	11.69	0.31	0.34	1.00

Equations

$$L_{\mathrm{pNRQCD}}=$$

$$-\frac{1}{4}{G}_{\mu \nu }^a{G}^{\mu \nu a}+\sum _{i=1}^{n_f}{\bar{q}}_i(i/D-m_i)q_i$$

$$\displaystyle+gV_{A}(r){\rm Tr}\left[{\rm O}^{\dagger}{\bf r}\cdot{\bf E}\,{\rm S}+{\rm S}^{\dagger}{\bf r}\cdot{\bf E}\,{\rm O}\right]+g\frac{V_{B}(r)}{2}{\rm Tr}\left[{\rm O}^{\dagger}\left\{{\bf r}\cdot{\bf E},{\rm O}\right\}\right]\Biggr{\}}\,.$$

$$h_s=$$

$$h_o=$$

$$H$$

$$+\int d^3\mathbf{R}\left(\frac{1}{2}\left({\mathbf{E}}^a\cdot {\mathbf{E}}^a+{\mathbf{B}}^a\cdot {\mathbf{B}}^a\right)-\sum _{i=1}^{n_i}{\bar{q}}_i[i\mathbf{D}\cdot \mathbf{\gamma }-m_i]q_i\right).$$

$$H{O}^{a†}(\mathbf{R},\mathbf{r})G_k^{ia}(\mathbf{R})|0⟩=(V_o^{(0)}+{\Lambda }_k)O^{a†}(\mathbf{R},\mathbf{r})G_k^{ia}(\mathbf{R})|0⟩,$$

$$⟨0|G_k^{ia†}({\mathbf{R}}^{\prime })O^a({\mathbf{R}}^{\prime },{\mathbf{r}}^{\prime })O^{b†}(\mathbf{R},\mathbf{r})G_{k^{\prime }}^{jb}(\mathbf{R})|0⟩={\delta }^{ij}{\delta }_{k{k}^{\prime }}\delta ({\mathbf{R}}^{\prime }-\mathbf{R})\delta ({\mathbf{r}}^{\prime }-\mathbf{r}).$$

$$|\mathbf{R},\mathbf{r};k,\lambda ⟩=P_{k\lambda }^i{O}^{a†}\left(\mathbf{R},\mathbf{r}\right)G_k^{ia}(\mathbf{R})|0⟩,$$

$$\int d^{3}r{d}^{3}R\sum _{k\lambda }|\mathbf{R},\mathbf{r};k,\lambda ⟩{\Psi }_{k\lambda }(t,\mathbf{r},\mathbf{R}),$$

$$P_{10}^i$$

$$P_{1\pm 1}^i$$

$$U=e^{i\Phi/F}\,,\quad\Phi=\left(\begin{array}[]{ccc}\pi^{0}+\frac{\eta}{\sqrt{3}}&\sqrt{2}\pi^{+}&\sqrt{2}K^{+}\\ \sqrt{2}\pi^{-}&-\pi^{0}+\frac{\eta}{\sqrt{3}}&\sqrt{2}K^{0}\\ \sqrt{2}K^{-}&\sqrt{2}\bar{K}^{0}&-\frac{2}{\sqrt{3}}\eta\\ \end{array}\right)\,,$$

$$u_{\mu }=$$

$${\chi }_{\pm }=$$

$$\displaystyle L_{pNRQCD}^{\rm had}=\int d^{3}Rd^{3}r\,{\rm Tr}\Bigl{[}S^{\dagger}\left(i\partial_{t}-V_{s}(r)+\frac{{\bm{\nabla}}^{2}_{r}}{m_{Q}}\right){\rm S}+\frac{F^{2}}{4}\left(\langle u_{\mu}u^{\mu}\rangle+\langle\chi_{+}\rangle\right)$$

$$\displaystyle+\sum_{k^{PC}=1^{+-},\,1^{--}}\sum_{\lambda\lambda^{\prime}}\Psi^{\dagger}_{k\lambda}\biggl{\{}(i\partial_{t}-(V_{o}^{(0)}+\Lambda_{k}))\delta_{\lambda\lambda^{\prime}}+\hat{\bm{r}}^{i\dagger}_{\lambda}\frac{{\bm{\nabla}}^{2}_{r}}{m_{Q}}\hat{\bm{r}}^{i}_{\lambda^{\prime}}\biggr{\}}\Psi_{k\lambda^{\prime}}$$

$$+\left(\mathbf{r}\cdot {\hat{\mathbf{r}}}_{\lambda }{S}^{†}{\Psi }_{1^{--}\lambda }+\text{h.c}\right)t^{(r{1}^{--})}+\mathbf{r}\cdot \delta \mathcal{L}({\Psi }_{k\lambda },{\Psi }_{k^{\prime }{\lambda }^{\prime }})$$

$$\displaystyle+\left(\bm{r}\cdot\hat{\bm{r}}_{\lambda}S^{\dagger}\Psi_{1^{--}\lambda}+\text{h.c}\right)\left(t^{(r1^{--})}_{d0}F^{2}\langle u_{0}u_{0}\rangle+t^{(r1^{--})}_{di}F^{2}\langle u_{i}u^{i}\rangle+t^{(r1^{--})}_{m}F^{2}\langle\chi_{+}\rangle\right)\Bigr{]}\,.$$

$$V_{{\Sigma }_g^{+}}(r)$$

$$V_{k\lambda {\lambda }^{\prime }}(r)$$

$$=\left({\Lambda }_k+V_o^{(0)}+b_{k\lambda }{r}^2+\cdots \right){\delta }_{\lambda {\lambda }^{\prime }}+\mathcal{O}(1/m_Q).$$

$$⟨0|{\Psi }_{\lambda }(t,{\mathbf{r}}^{\prime },{\mathbf{R}}^{\prime })S^{†}(0,\mathbf{r},\mathbf{R})|0{⟩}_{amp.}=i\mathbf{r}\cdot {\hat{\mathbf{r}}}_{\lambda }^{†}{t}^{(r{1}^{--})},$$

$$Z_E⟨0|{\hat{\mathbf{r}}}_{\lambda }^{†}\cdot {\mathbf{E}}^{a†}(t,{\mathbf{R}}^{\prime })O^a\left(t,{\mathbf{r}}^{\prime },{\mathbf{R}}^{\prime }\right)S^{†}(0,\mathbf{r},\mathbf{R})|0{⟩}_{amp.}=i\sqrt{\frac{T_F}{N_c}}{\hat{\mathbf{r}}}_{\lambda }^{†i}{\mathbf{r}}^j⟨0|G_{1^{--}}^{ai†}(\mathbf{R})g{\mathbf{E}}^{aj}(\mathbf{R})|0⟩,$$

$$G_{1^{--}}^a=Z_E^{-1/2}g{\mathbf{E}}^a+Z_{D\times B}^{-1/2}{\left(\mathbf{D}\times g\mathbf{B}\right)}^a+\cdots .$$

$$t^{(r{1}^{--})}=\sqrt{\frac{T_F{Z}_E}{N_c}}.$$

$$\int d^4{x}_{+}{d}^4{x}_{-}{e}^{i{p}_{+}\cdot x_{+}}{e}^{i{p}_{-}\cdot x_{-}}⟨0|{\pi }^{+}(x_{+}){\pi }^{-}(x_{-})S(t,\mathbf{r},\mathbf{R}){\Psi }_{1^{--}\lambda }^{†}(0,\mathbf{r},\mathbf{R})|0{⟩}_{amp.}$$

$$=i4\mathbf{r}\cdot {\hat{\mathbf{r}}}_{\lambda }\left(-t_{d0}^{(r{1}^{--})}{p}_{+}^0{p}_{-}^0+t_{di}^{(r{1}^{--})}{\mathbf{p}}_{+}\cdot {\mathbf{p}}_{-}-t_m^{(r{1}^{--})}{m}_{\pi }^2\right),$$

$$⟨{\pi }^{+}(p_{+}){\pi }^{-}(p_{-})|S(t,\mathbf{r},\mathbf{R}){\hat{\mathbf{r}}}_{\lambda }\cdot G_{1^{--}}^a(0,\mathbf{R})O^{a†}(0,\mathbf{r},\mathbf{R})|0{⟩}_{amp.}$$

$$=ig\sqrt{\frac{T_F}{N_c}}{\hat{\mathbf{r}}}_{\lambda }\cdot \mathbf{r}⟨{\pi }^{+}(p_{+}){\pi }^{-}(p_{-})|\mathbf{E}(\mathbf{R})\cdot G_{1^{--}}(\mathbf{R})|0⟩$$

$$=\frac{i}{3}\sqrt{\frac{T_F}{N_c{Z}_E}}{\hat{\mathbf{r}}}_{\lambda }\cdot \mathbf{r}\frac{8{\pi }^2}{{\beta }_0}\left(\left(2-\frac{9}{2}\kappa \right)p_{+}^0{p}_{-}^0-\left(2+\frac{3}{2}\kappa \right){\mathbf{p}}_{+}\cdot {\mathbf{p}}_{-}+3m_{\pi }^2\right),$$

$$t_{d0}^{(r{1}^{--})}=-\frac{2{\pi }^2}{3{\beta }_0}\sqrt{\frac{T_F}{N_c{Z}_E}}\left(2-\frac{9}{2}\kappa \right),$$

$$t_{di}^{(r{1}^{--})}=-\frac{2{\pi }^2}{3{\beta }_0}\sqrt{\frac{T_F}{N_c{Z}_E}}\left(2+\frac{3}{2}\kappa \right),$$

$$t_m^{(r{1}^{--})}=-\frac{2{\pi }^2}{3{\beta }_0}\sqrt{\frac{T_F}{N_c{Z}_E}}3.$$

$${\mathcal{A}}_{\chi }$$

$$i\mathcal{A}$$

$${\beta }_{r,Q}^{(n^{\prime }n)}=\sum _m⟨S_{n^{\prime }}|{\hat{\mathbf{r}}}_{\lambda }^{†}\cdot \mathbf{r}|{\Psi }_m⟩\left(\frac{1}{m_n-m_m}+\frac{1}{m_{n^{\prime }}-m_m}\right)⟨{\Psi }_m|{\hat{\mathbf{r}}}_{\lambda }^{†}\cdot \mathbf{r}|S_n⟩.$$

Full text

1. I Introduction

2. II Hadronization of the pNRQCD Lagrangian

3. III $\bar{Q}Q(2S)\rightarrow\bar{Q}Q(1S)\pi\pi$ hadronic transitions

4. III.1 Line-shape analysis

5. III.2 Total decay widths

6. IV $\bar{Q}Q(2S)\rightarrow\bar{Q}Q(1P)\pi$ hadronic transitions

7. V $\bar{Q}Q(2S)\rightarrow\bar{Q}Q(1S)\pi^{0}(\eta)$ hadronic transitions

8. VI Ratios

9. VII Conclusions

10. A Differential decay widths

11. B Singlet and hybrid bound states

12. C Matrix Elements

13. D Hadronization of local gluonic operators

14. D.1 Two-pion production through the anomaly of the energy-momentum tensor

15. D.2 Single pion creation from the axial anomaly

16. D.3 Single $P$-wave pion anomaly

I Introduction

Hadronic transitions between heavy quarkonium states have been studied since the middle seventies. Particular attention has been devoted to the transitions with one or two pions in the final state. The relatively small energy of the outgoing pions makes the analysis of these transitions ideal for the implementation of (the nonlinear realization of) chiral symmetry. This was first implemented by Cahn and Brown in Ref. Brown and Cahn (1975). Using current algebra, they parametrized the amplitude of the two-pion transition in the strict chiral limit. The observation of the small variation of the width distribution with respect to the helicity angle of one of the pions led them to an approximated simple one-parameter description of the decay width spectrum. The implementation of chiral symmetry using chiral Lagrangians and the incorporation of the leading light-quark mass corrections was made in Ref. Mannel and Urech (1997). This gave a good description of the width spectrum, and also of the pion helicity angle distribution Yan et al. (1999); Bai et al. (2000).

Alternatively, in Ref. Gottfried (1978), Gottfried realized that these transitions could be though as a two-step process: first a short-distance gluon emission by the heavy quarks, and then the hadronization of the gluons into the light-quark hadrons at a relatively long distance. This picture follows from the use of the multipole expansion of the gluon-heavy quark interaction. This allowed Gottfried to give selection rules and decay width rate estimates beyond the approaches based solely on chiral symmetry. The work of Ref. Yan (1980) connected these nonlocal gluonic matrix elements with chiral symmetry by parametrizing them according to chiral symmetry. This was used to obtain constrains for higher angular momentum channels.

At leading order (LO) in the multipole expansion, the intermediate heavy quarkonium state is in a color-octet configuration. In Ref. Voloshin (1979), Voloshin introduced an additional expansion consisting of an operator product expansion of the nonlocal heavy quarkonium color-octet two-point function. As a result, the transition amplitude can be written as a series of matrix elements of local operators, which we will refer to as a twist expansion. As before, these local matrix elements can be hadronized and parametrized using chiral symmetry but the implementation of the axial Gross et al. (1979); Novikov et al. (1980) and energy-momentum tensor anomalies Voloshin and Zakharov (1980); Novikov and Shifman (1981); Chivukula et al. (1989) constraints their general structure. Since then it has become quite customary to describe the two-pion transitions using the multipole, twist, and the chiral expansion. See, for instance, Ref. Chen and Savage (1998), where the local gluonic matrix elements were obtained up to next-to-leading order (NLO) in the chiral expansions. Unfortunately, this computational scheme is not a model independent derivation of QCD. The reason is that there is no kinematic regime where the twist expansion can be justified, as shown in Ref. Luty and Sundrum (1993), because the transfer energy between heavy quarkonium states in these transitions, $E$, is of order $m_{Q}v^{2}$, whereas the twist expansion requires that $E\ll m_{Q}v^{2}$, as well as $\Lambda_{\rm QCD}\ll mv^{2}$.

We want to retake this discussion within the context of effective field theories (EFTs), which allows us to make a systematic analysis of the scales involved in the problem. We will do the analysis using the weak-coupling version of potential NRQCD (pNRQCD) Pineda and Soto (1998); Brambilla et al. (2000). This allows us to obtain the multipole expansion in a controlled way, since we assume that $m_{Q}v\gg\Lambda_{\rm QCD}$. The resulting EFT has the multipole expansion built in in the Lagrangian. The LO Lagrangian will also have spin symmetry, which follows from the heavy quark mass expansion. To simplify the problem further, we will also organize the computation within a $1/N_{c}$ expansion. If we hadronize this EFT, we obtain a Lagrangian in terms of the singlet (standard heavy quarkonium), hybrids and pion fields. $B/D$ mesons and possible tetraquarks are, a priori, subleading in the large $N_{c}$, as we will discuss in more detail later in the paper. Therefore, we will not consider these degrees of freedom in this paper and their possible incorporation will be relegated to future work. As their effect could be important for hadronic transitions of states close or above to open flavor thresholds, we focus on states below threshold in this paper.

We then write the most general hadronic representation of the weakly coupled pNRQCD Lagrangian made of the singlet, hybrids and pions fields in a combined expansion in the chiral counting and the multipole, $1/m_{Q}$ and $1/N_{c}$ expansions. It is formally possible to obtain the coefficients of the Lagrangian by matching to suitable Green functions in the pNRQCD theory in terms of quarks and gluons. Nevertheless, these coefficients endure a complicated relation with the elementary fields of the theory, and to determine them would require quite costly lattice simulations. Nevertheless, no all hope is lost. In Ref. Brambilla et al. (2000), it was observed a correlation between the dimensionality of the interpolating operator of the hybrid/gluelump and its position in the spectrum. Though not equivalent, we hypothesize in this paper that there is also a correlation between the dimensionality of the interpolating operator and the strength of the interpolation with the hybrid, such that higher dimension operators are subleading (in this respect, it would also be interesting to study whether this approximation can be related with the analysis made in Fitzpatrick et al. (2013)). We will see that such hypothesis plus the chiral low-energy theorems (generated by the axial and energy-momentum tensor anomalies) lead to the same predictions for the normalized differential decay rates of the two-pion transitions as using the twist expansion. Moreover, we do not need to impose extra conditions on $\Lambda_{\rm QCD}$, except the one we already imposed for the multipole expansion $m_{Q}v\gg\Lambda_{\rm QCD}$. In other words, our computational framework would still be valid even if $\Lambda_{\rm QCD}\gg m_{Q}v^{2}$. We explore the implications of this computational framework for a series of observables. Finally, we want to mention that the pure quarkonium hybrid sector of this theory has already been developed in Refs. Berwein et al. (2015); Brambilla et al. (2018); Oncala and Soto (2017); Brambilla et al. (2019) for small energy fluctuations (much smaller than the energy transitions we consider in this paper).

We organize the paper as follows: in Sec. II, we discuss the hadronization of the pNRQCD Lagrangian. In Sec. III, we apply it to describe the $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi\pi$ transitions and compare the results for the decay width spectrum with experiment and the purely chiral description. Extending our Lagrangian beyond LO we use our approach to study one-pion transitions, $Q\bar{Q}(2S)\to Q\bar{Q}(1P)\pi$ in Sec. IV, and $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi$ in Sec. V. Certain uncertainties of our approach cancel out for specific ratios of the decay widths. We compute and study them in Sec VI. Finally, we give our conclusions in Sec. VII.

II Hadronization of the pNRQCD Lagrangian

The pNRQCD Lagrangian at LO in $1/m_{Q}$ (except for the kinetic term) and at NLO in the multipole expansion reads

[TABLE]

${\rm S}$ and ${\rm O}$ are the quark singlet and octet fields, respectively, normalized with respect to color as $S=S\bm{1}_{c}/\sqrt{N_{c}}$ and ${\rm O}=O^{a}T^{a}/\sqrt{T_{F}}$. They should be understood as functions of $t$, the relative coordinates $\bm{r}$, and the center of mass coordinates $\bm{R}$ of the heavy quarks. The trace should be understood as a double trace in color and spin. The singlet, octet and hybrid fields in the Lagrangians that appear in this paper are organized in $SU(2)$ spin multiplets. For instance, $S=\frac{1}{\sqrt{2}}({\bf S}\cdot{{\bm{\sigma}}}+S_{\eta}{\cal I})$. All the fields of the light degrees of freedom in Eq. (1) are evaluated at $\bm{R}$ and $t$; in particular, $G^{\mu\nu\,a}\equiv G^{\mu\nu\,a}(\bm{R},\,t)$, $q_{i}\equiv q_{i}(\bm{R},\,t)$, and $iD_{0}O\equiv i\partial_{0}O-g\left[A_{0}(\bm{R},\,t),O\right]$. $h_{s}$ and $h_{o}$ are the singlet and octet Hamiltonian densities. They read as

[TABLE]

where $V_{s}(r)$ and $V_{o}(r)$ are computed in perturbation theory. Note that we have spin symmetry. Unless stated otherwise we will work in the isospin limit $m_{i}=\hat{m}\equiv\frac{m_{u}+m_{d}}{2}$ with $i=u,d$. At leading log Pineda and Soto (2000), and next to leading log Brambilla et al. (2009) accuracy $V_{A}=1$.

We next aim to construct the hadronic version of the above Lagrangian. We first need to characterize the hadronic degrees of freedom relevant to our case. We first work in the static limit. In the short heavy-quark-antiquark distance limit the gluonic excitations can be characterized by the so-called gluelump operators. They organize themselves in irreducible representations of the $O(3)\otimes$C group. The LO Hamiltonian density in the $1/m_{Q}$ and multipole expansions corresponding to the Lagrangian in Eq. (1) is given by

[TABLE]

For later convenience, we will restrict the discussion to $L=1$ gluelump states. We define the gluelump operators, $G_{k}^{ia}$, as the color-octet gluonic operators that generate the eigenstates of $H$ in the presence of a local heavy-quark-antiquark octet source:

[TABLE]

where $a$ is the color index, $k$ labels the quantum $J^{PC}$ numbers of the gluelump, and $i$ labels its vector components. At this stage, we do not have to make explicit the spin content of $O^{a}$. We normalize the gluelump operators as

[TABLE]

Going beyond the LO in the multipole expansion the system is no longer spherically symmetric, instead it is cylindrically symmetric around the heavy-quark-antiquark axis.111The symmetry group is $D_{\infty h}$, with P replaced by CP. Representations of the cylindrical symmetry group can be constructed by projecting the gluelump operators on various directions with respect to the heavy-quark-antiquark axis. Therefore, we work with states with good transformation properties under the cylindrical symmetry group

[TABLE]

where summation over index $i$ is implied. $P^{i}_{k\lambda}$ is a projector that acts onto the gluelump angular momentum and projects it into an eigenstate of $\mathbf{K}\cdot\hat{\mathbf{r}}$ (where $\mathbf{K}$ is the angular momentum operator for the gluelump) with eigenvalue $\lambda$.

It is useful to project the pNRQCD Lagrangian onto the Fock subspace spanned by the $|{\bf R},{\bf r};k,\,\lambda\rangle$ states

[TABLE]

where $\Psi_{k\lambda}(t,\,\bm{r},\,\bm{R})$ will represent the hybrid field in the hadronic version of the EFT. As we have already mentioned, the case of most interest to us is that of the spin $1$ gluelumps. The projectors for this case are

[TABLE]

and the discrete symmetry transformations for the ${\rm S}$, ${\rm O}$, and $\Psi^{i}_{\kappa}$ fields are given in Table 1.

Besides the singlet and hybrid fields, we will incorporate pions to our Lagrangian. As a basic building block for the Goldstone bosons, we use the unitary matrix $U(t,\,\bm{R})$, which [for $SU(3)$] may be taken as

[TABLE]

although final results for observable quantities do not depend on this specific choice. $F=92.419$ MeV is the pion decay constant. Under chiral symmetry, $U$ transforms as $U\stackrel{{\scriptstyle g}}{{\to}}RUL$, where $R\in SU_{R}(N)$ and $L\in SU_{L}(N)$. Related useful matrices are $u$, defined from $u^{2}=U$, and

[TABLE]

where $\chi=2B{\rm diag}(\hat{m},\hat{m},m_{s})$, with $B$ being related to the vacuum quark condensate. In the isospin limit, the pion mass is $m_{\pi}^{2}=2B\hat{m}$. The transformation properties can be found in Table 2.

We now construct the hadronic Lagrangian. It is fixed by the degrees of freedom, the symmetries, and the parameter expansions we have: $1/N_{c}$, $r$, $E$, $1/m_{Q}$, $m_{i}$. We emphasize that, at this level, we do not integrate out extra degrees of freedom when going to the hadronic representation of Eq. (1). Therefore, it is not a different EFT but the very same pNRQCD, including the same degrees of freedom and scales. Instead what we do is to write the most general Lagrangian consistent with the symmetries made out of the heavy quarkonium, hybrids, and pions. We write this Lagrangian at LO in the chiral counting and the $1/m_{Q}$ expansion (except for the kinetic term) and at NLO in the multipole expansion. For the interaction between heavy quarkonium, hybrids, and pions, we also incorporate the large $N_{c}$ expansion, and consider only the leading terms in it. Strictly speaking, one should also include glueballs, which may interact with the hybrids at LO in the multipole expansion, but at NLO in the $1/N_{c}$ expansion. Such effects would be subleading in our computation of the decays. Therefore, we neglect them. The hadronic version of the pNRQCD Lagrangian projected onto the subspace of Eq. (8) reads

[TABLE]

Note that for each hybrid channel one should also include the excitations. The fields $S$ and $\Psi_{k\lambda}$ should be understood as depending on $t$, $\bm{r}$, and $\bm{R}$. The pion fields depend on $t$ and $\bm{R}$. Note that Tr[] now only stands for the trace over spin indices. $\langle A\rangle$ stands for the trace of $A$ in the isospin index.

Since at NLO in the multipole expansion the singlet can only mix with $k=1$ gluelumps, we only include hybrid states that can be generated by such gluelumps in the Lagrangian. According to the operator analysis of Table II in Ref. Brambilla et al. (2000), they correspond to $\Sigma_{g}^{+^{\prime}}$, $\Pi_{g}$, $\Sigma_{u}^{-}$, $\Pi_{u}$, and associated excitations with bigger gluelump masses. The mixing between hybrids at NLO in the multipole expansion is more complicated. It is encoded in ${\bf r}\cdot\delta{\bf{\cal L}}(\Psi_{k\lambda},\Psi_{k^{\prime}\lambda^{\prime}})$, where $\delta{\bf{\cal L}}(\Psi_{k\lambda},\Psi_{k^{\prime}\lambda^{\prime}})$ is bilinear in the hybrid fields and transforms as a $1^{--}$ vector. Note that in general $k\not=1$ states may contribute to this term. Fortunately, we will not need the details of this interaction for the analysis of this paper. Note also that at NLO in the multipole expansion we do not have hybrid bound states but plane waves. We will need to iterate the ${\cal O}(r)$ vertices to obtain bound states in the hybrid sector.

The last line in Eq. (14) encodes the interaction of the hybrids with the singlet and pions at NLO in the multipole expansion and at the leading nonvanishing order in the large $N_{c}$ and chiral expansions. The interaction with two pions scales with $N_{c}$ as $1/N_{c}^{2}$. One may consider the incorporation of subleading operators in the chiral counting but still leading in the $1/N_{c}$ expansion. Since for the process we consider in this paper, the typical energy of the pions will be of ${\cal O}(mv^{2})$ or smaller (which in general we will consider them to be smaller than $\Lambda_{\rm QCD}$); this will not be necessary. This is also the reason we do not include interactions of the pions with the hybrids that are not suppressed by the multipole expansion. The contribution of these operators to heavy quarkonium hadronic transitions would be suppressed by a power of $mv^{2}/\Lambda_{\rm QCD}$ compared with the contributions considered in this paper.

In the above discussion we have not included tetraquarks (like the $Z_{b}$ or $Z_{c}$) , nor $Q\bar{q}$-$\bar{Q}q$ states, in the physical spectrum. One may wonder whether we should do so. Here we will guide our discussion by the $1/N_{c}$ counting. The inclusion of tetraquarks is delicate. They are not stable in the large $N_{c}$. If the decay width grows like $N_{c}$, then $\Gamma\gg E$ and diagrams with intermediate tetraquarks become effectively local and $1/N_{c}$ suppressed relative to the accuracy for the physical processes we consider in this paper. The same happens for the $Q\bar{q}$-$\bar{Q}q$ loops. The suppression of tetraquark effects can only be bypassed if, for some reason, $\Gamma$ is much smaller than what is expected by the $N_{c}$ counting. This could happen if the channels that would contribute at leading order in $N_{c}$ are closed because they are below threshold. This is the situation discussed in Ref. Weinberg (2013). Then, indeed, the scaling in $N_{c}$ of $\Gamma$ is, at most, of order $1/N_{c}$. However, we are not in this situation. For instance, $Z_{b}$ states are above threshold, albeit very close to it. Therefore, there is little phase space free, effectively producing that $\Gamma$ is very small. In this scenario the mixing with $Q\bar{q}$-$\bar{Q}q$ loops is expected to be very large and it does not make much sense to consider one without the other. In this respect, it is interesting to note that in Refs. Chen et al. (2016, 2017) it has been advocated that tetraquarks and $Q\bar{q}$-$\bar{Q}q$ loops may play an important role for some observables like the $Q\bar{Q}(3S/4S)\rightarrow Q\bar{Q}(1S)$ decays. Nevertheless, in this paper, we want to specifically study the effect associated to the inclusion of hybrids within a context where the multipole expansion can be applied. Therefore, the possible incorporation of tetraquark and $Q\bar{q}$-$\bar{Q}q$ states will be postponed to future work.

In this paper, we will generally consider that the quarkonium binding energies fulfill $m_{Q}v^{2}\ll\Lambda_{\rm QCD}$. This allows us to write the static singlet potential up to ${\cal O}(r^{2})$ as

[TABLE]

The coefficient $b_{\Sigma^{+}_{g}}$ can actually be determined in terms of the coefficients $t^{(r1^{--})}$ from the Lagrangian in Eq. (14), since ${\cal O}(r^{2})$ terms in the Lagrangian do not contribute to the static energies of the singlet.

In the case of hybrid bound states, $m_{Q}v^{2}\ll\Lambda_{\rm QCD}$ is also the natural hierarchy between the scales since the bound states are small energy fluctuation around the minima of the hybrid static energies. Therefore, we can also generically write

[TABLE]

Nevertheless, in this case, we cannot guaranty that ${\cal O}(r^{2})$ terms in the Lagrangian will not contribute to the hybrid potential. Therefore, we cannot determine $b_{k\lambda}$ from the ${\cal O}(r)$ coefficients of the hadronic Lagrangian alone. We will determine this coefficient by fitting the potential to the lattice data for the static energies.

The LO singlet and hybrid spectrum and wave functions will be obtained from solving the Schrödinger equation with different variants of the singlet and hybrid potential based on Eqs. (15) and (16). The details of these solutions can be found in Appendix B.

Let us now discuss the vertices of the hadronic Lagrangian. To obtain the coupling in the hadronic Lagrangian we match equivalent correlation functions computed in both representations of the theory. The equalities are obtained in the static limit. The singlet-hybrid mixing term in the hadronic theory is given by

[TABLE]

which is matched to the following correlator of pNRQCD in term of quarks and gluons:

[TABLE]

where amp. signals that only amputated contributions are considered [overall $\delta({\bf r}^{\prime}-{\bf r})$ are also factored out]. To evaluate the matrix element, we consider the interpolating field for $G_{1^{--}}^{a}$ to be given by a sum of all possible local gluonic operators with the same quantum numbers,

[TABLE]

Now, we then hypothesize that there is a correlation between the dimensionality of the interpolating operator and the strength of the interpolation with the hybrid, such that higher dimension operators are subleading, so the series can be truncated at LO. Though not rigorous, one may consider counting $Z_{O}^{-1/2}\sim\Lambda^{1-{\rm dim}(O)}_{\rm QCD}$, where dim(O) is the dimension of the operator O, and the size of the operator to be given by $\sim\left(mv^{2}\right)^{{\rm dim}(O)}$, so the series can be truncated at LO up to corrections of $\mathcal{O}(mv^{2}/\Lambda_{\rm QCD})$. Using this truncation and Eq. (6) in Eq. (18), we obtain

[TABLE]

The operators with an even number of pions in Eq. (14) are matched in a similar way. In the hadronic pNRQCD,

[TABLE]

and the corresponding correlator in the partonic pNRQCD reads

[TABLE]

where in the last step we use Eq. (19) truncated at the first term and hadronize the resulting gluonic matrix element using Eq. (158), which uses the anomaly relation of the energy-momentum tensor of QCD Voloshin (2008); Novikov and Shifman (1981). The derivation of this equation is reviewed in Appendix D.1.

Comparing Eqs. (21) and (22), we obtain

[TABLE]

III $\bar{Q}Q(2S)\rightarrow\bar{Q}Q(1S)\pi\pi$ hadronic transitions

The two-pion quarkonium transitions are well described using chiral perturbation theory and spin symmetry. The most general amplitude for the two-pion transitions at ${\cal O}(p^{2})$ in the chiral counting reads

[TABLE]

where $p_{+}$, $p_{-}$ are the momentum of the $\pi^{+}$ and $\pi^{-}$, respectively. The coefficients $a_{i}$ can be thought as linear combinations of Wilson coefficients of an effective chiral Lagrangian made only by heavy quarkonium and pions. Such effective Lagrangian can be found in Eqs. (4) and (6) of Ref. Mannel and Urech (1997). Imposing heavy quark spin symmetry sets the $g_{2}$ low-energy constant of Ref. Mannel and Urech (1997) to zero.

Alternatively, we can use the computational scheme developed in this paper. Using Eq. (14), we can compute the leading nonvanishing contribution to the two-pion transitions between quarkonium states. The diagrams are drawn in Fig. 1. One vertex mixes singlet and hybrid fields whereas the other produces a two-pion emission vertex. The amplitude for this process reads

[TABLE]

with

[TABLE]

The index $m$ sums over all states solution of the Schrödinger equation and also for the gluelump excitations. For the second term in the brakets of Eq. (27) (corresponding to the right diagram in Fig. 1), we have used that $m_{n}=m_{n^{\prime}}+p^{0}_{+}+p^{0}_{-}$. It is remarkable that the normalization factors $Z_{E}$ cancels out, which allows us to completely evaluate Eq. (27) except for the parameter $\kappa$ that appears in the couplings in Eqs. (23)-(25) and has its origin in the hadronization of the gluonic operator. In the next section we will fix $\kappa$ by fitting the normalized differential decay width spectrum.

For the case at hand, the quarkonia states are $n^{\prime}=1^{3}S_{1}$ and $n=2^{3}S_{1}$. The intermediate hybrid has $1^{--}$ quantum numbers. Expressions for the hybrid-quarkonia matrix elements that appear in Eq. (28) are given in Appendix C. The wave functions of the bound states that appear in those matrix elements can be found in Appendix B. The matrix elements in Eq. (28) are diagonal in the spin state of the heavy quarks. Thus, only spin triplet hybrid states are allowed as intermediate states. Furthermore, the matrix elements in Eq. (28) only receive contribution from the $\lambda=0$ component of the hybrid state, which leads to the selection rule $l=\ell$ for the angular momentum eigenvalues. Therefore, the only possible hybrid intermediate states correspond to $m^{3}\mathcal{S}_{1}$ states.

Equation (27) yields the following prediction for the chiral parameters $a_{i}$:

[TABLE]

We now confront the above predictions with experiment. Since the experimental data on $\frac{d\Gamma}{dm_{\pi\pi}}$ (where $m_{\pi\pi}=(p_{+}+p_{-})^{2}$ is the dipion invariant mass) are normalized to an unknown constant, it is convenient to fit the theoretical expressions to the normalized differential decay width

[TABLE]

As we will see there is also a strong theoretical motivation to consider this ratio. This object is only sensitive to ratios of the theory parameters. The overall normalization can be obtained from the total decay width. Therefore, in the following subsections, we fit the normalized differential decay width and total decay width using the formulas discussed above. Detailed formulas can be found in Appendix A. Here we will analyze the transitions $\psi(2S)\to J/\Psi\pi^{+}\pi^{-}$ and $\Upsilon(2S)\to\Upsilon(1S)\pi^{+}\pi^{-}$. The experimental data for the former are taken from Ref. Aaboud et al. (2017), and for the latter from Ref. Alexander et al. (1998).222In the first version of this paper the more recent data from Ref. Guido et al. (2017) was used. However this data appears to have some problems. This reflects in an strange behavior of the angular distributions. We acknowledge conversations on this issue with the authors. The total decay widths are taken from Ref. Tanabashi et al. (2018).

III.1 Line-shape analysis

The chiral fit to $\frac{1}{\Gamma}\frac{d\Gamma}{dm_{\pi\pi}}$ produces the following ratios of parameters:

[TABLE]

where the superindex $c$ and $b$ label the results for the $\psi(2S)\to J/\Psi\pi^{+}\pi^{-}$ and $\Upsilon(2S)\to\Upsilon(1S)\pi^{+}\pi^{-}$ transitions respectively. These fits are performed using the relativistic kinematics [see Eqs. (91) and (92)]. If instead we use nonrelativistic kinematics (see Eqs. (93) and (94)), fits of very similar quality are obtained. The results are shown in Fig. 2 together with the experimental data. Note the small $\chi_{\text{d.o.f}}^{2}$ obtained in the fits to both transitions. The uncertainties quoted in Eqs. (33) and (34) correspond to the range of parameter values with $\chi^{2}_{\text{d.o.f}}\leq 1$, which is shown in Fig. 3. In this figure, it can be appreciated that there is a strong correlation between $a^{Q}_{1}/a^{Q}_{2}$ and $a^{Q}_{3}/a^{Q}_{2}$. Therefore, it would be wrong to consider the error of $a^{Q}_{1}/a^{Q}_{2}$ and $a^{Q}_{3}/a^{Q}_{2}$ independently. Instead, the allowed region of parameter space is given by the dashed region in Fig. 3. This corresponds to a particular linear combination of $a^{Q}_{1}/a^{Q}_{2}$ and $a^{Q}_{3}/a^{Q}_{2}$. We conclude then that, with the present data on the width spectrum, we cannot simultaneously fit both parameter ratios with high accuracy .

If we compare these fits with the ones in Ref. Mannel and Urech (1997), we observe some differences. The authors in Ref. Mannel and Urech (1997) report that the best fit yields $a_{3}=0$ ($g_{3}$ in their notation). This does not coincide with our best fit, but is consistent with the uncertainty. The ratio $a_{1}/a_{2}$ corresponds to $1+g_{1}/(2g)$ in the notation of Ref. Mannel and Urech (1997) and their best fits correspond to $a^{c}_{1}/a^{c}_{2}\sim 1.2$ and $a^{b}_{1}/a^{b}_{2}\sim 1.1$. This is compatible with our fits considering the uncertainty. One should keep in mind that the experimental data source in Ref. Mannel and Urech (1997) is Ref. Albrecht et al. (1987), whereas we use more recent data Alexander et al. (1998)333Very similar results to those in Ref. Mannel and Urech (1997) were obtained in Ref. Guo et al. (2005) using more recent experimental data Bai et al. (2000); Alexander et al. (1998) including pion final state interactions through unitarization of $\chi$PT..

We now use the amplitude in Eq. (27) to obtain the normalized differential decay width and fit it to the experimental data. Equation (27) is equivalent to Eq. (26) with the parameters $a_{i}$ taken as in Eqs. (29)-(31). The normalized decay width is independent of $\beta^{(12)}_{r,Q}$, and its functional form depends only on the parameter $\kappa$. Thus, the hadronic pNRQCD, which incorporates the multipole expansion, and a dimensional counting for the overlap of the hybrids with gluonic operators, is more predictive than the EFT relying on chiral and spin symmetry only, since the number of free parameters is reduced from two to one. Note that our approach yields the same normalized decay width line shape, as the one obtained using the twist expansion.

We fit the line shapes of the charmonium and bottomoniun data independently. Since the value of $\kappa$ should be independent of the heavy quarkonium dynamics, we also perform a simultaneous fit to both data sets. Using the nonrelativistic kinematics, we obtain the following values for $\kappa$:

[TABLE]

The range of values correspond to the range with $\chi^{2}-\chi_{\rm min}^{2}\leq 1$. If we use the relativistic kinematics, the fits yield slightly different values,

[TABLE]

The differences with the nonrelativistic fit are of the order of the difference between the charmonium and bottomonium fit. This is reasonable, as they both measure relativistic effects. We will take the difference between Eqs. (37) and (40) as a measure of the size of subleading effects. We then combine it in quadrature with the statistical error quoted in Eq. (40) and give

[TABLE]

as our default value. In any case, it is remarkable that all fits yield similar values for $\kappa$, which we take as a confirmation of the independence of $\kappa$ of the heavy quarkonium dynamics. For illustration, we show the plots of the fit using relativistic kinematics for charmonium and bottomonium data in Fig. 4. Actually, we can also observe in Figs. 2 and 4 the similarity of the experimental data for the line shapes of the charmonium and bottomonium spectra. This can be taken as a reflection of the independence of this observable on the heavy quarkonium dynamics, which is a prediction of the effective theory.

Previous fits of the decay width spectrum of $\psi(2S)\to J/\psi\pi\pi$ have been carried out in Ref. Bai et al. (2000) using the transition amplitude from Ref. Novikov and Shifman (1981), and in Ref. Chen (2019) in which the pion final state interactions were taken into consideration through an Omnès function. The reported values are $\kappa_{c}=0.186(3)$ and $\kappa_{c}=0.135(5)$, respectively. We have checked that the main source of this discrepancy is that the transition amplitudes used in those references do not include the complete ${\cal O}(p^{2})$ pion mass contribution: in Ref. Chen (2019), the coefficient $a_{3}$ is set to zero, and in Ref. Bai et al. (2000), terms proportional to the pion mass have been set to zero.

III.2 Total decay widths

The expressions for the total decay width in terms of the chiral coefficients $a_{i}$ can be found in Eqs. (97) and (98). Inserting the values from Eq. (33) or Eq. (34), the remaining free parameter, $a_{2}$, can be adjusted to reproduce the experimental value of the total decay width. In fact, since the total width is a quadratic function in $a_{2}$, two solutions are possible,

[TABLE]

In our EFT, the total decay width is proportional to $(\beta^{(12)}_{r,Q})^{2}$ (note that it is at this level where there is a difference with predictions using the twist expansion). This object is dependent on the wave functions of heavy quarkonium and hybrids, as well as on the energy difference among them. Typically, this quantity will suffer from rather large uncertainties, as we are forced to make strong approximations to compute this object. We neglect the effect of higher gluelump excitations with the same quantum numbers. We expect that higher gluelumps will give smaller contributions, as they are suppressed by larger energy differences. Still, this is an approximation. In some circumstances, it can be compulsory to sum all of them to recover some high energy logarithms. Nevertheless, compared with other uncertainties this effect will be small. Therefore, in this paper, we neglect the error associated to neglecting higher gluelump channels. We account for the error associated to the energy splitting between singlet and hybrid states using the error of $\Lambda_{1}$, the lowest gluelump mass (for further details, see Appendix B). For the lowest hybrid, we compute $\beta^{(12)}_{r}$ with the $m^{3}\mathcal{S}_{1}$ hybrid intermediate states $m=1,...,4$ (we observe that the effect of introducing three or four hybrid states is comparatively small compared with other uncertainties). For the static potential of the singlet and hybrid, we take two possible parametrization that we explain in further detail in Appendix D.1. One aims for a good description at all distances of the lattice-evaluated static singlet and hybrid potentials, but constrained to have the right behavior at short distance. The other parametrization keeps the shape of the potentials predicted by the multipole expansion in the whole fitted range. Throughout the paper, we will take the former for our central values and the difference with the latter as an estimate of the error associated to the approximate knowledge of the singlet and hybrid wave functions and binding energies. We obtain the following values for $\beta^{(12)}_{r,Q}$:

[TABLE]

The uncertainties are labeled according to the source: $\Lambda_{1}$ (the lowest lying gluelump mass), and s.p. (the different parametrization for the singlet and hybrid static potentials). At this point, it is worth mentioning that $\beta^{(11)}_{r,Q}$, relevant for heavy quarkonium-nucleon interaction and the chiral-associated heavy quarkonium energy shift, has also been analyzed using hybrid potentials in Ref. Lakhina and Swanson (2004).

Introducing these values of $\beta^{(12)}_{r,Q}$, and $\kappa=0.247(20)$, in Eqs. (29)-(31), the expressions for the total width in Eqs. (97) and (98) read (the experimental values are taken from Tanabashi et al. (2018))

[TABLE]

Our numbers differ from the experimental ones by about a factor 2. One should keep in mind however that our estimates suffer from large uncertainties. We find a significant dependence on variations of the wave function of the hybrid and singlet states. The error generated by the uncertainty on the energy difference between singlet and hybrid states is somewhat smaller.These error estimates are of the right magnitude, though not large enough, to completely account for the difference with experiment. One should keep in mind, however, that besides those errors already estimated, one error that has not been incorporated in this analysis is due to the uncertainties associated to the hadronization of the local operator: on the one hand, we have ${\cal O}(\alpha_{s})$ corrections to Eqs. (157) and (158) due to the beta function. These ${\cal O}(\alpha_{s})$ corrections are generated at a low-energy scale, which makes their evaluation not feasible. On the other hand, these effects factor out and could be reabsorbed in $(\beta^{(12)}_{r,Q})^{2}$ if we let this object to be a free parameter, not fixed by theory. Alternatively, these effects are independent of the heavy quarkonium dynamics and would cancel in the ratio $\Gamma_{\psi(2S)\to J/\psi\pi^{+}\pi^{-}}/\Gamma_{\Upsilon(2S)\to\Upsilon(1S)\pi^{+}\pi^{-}}$. We will discuss this ratio later in Sec. VI. Other corrections to Eqs. (157) and (158) are due to ${\cal O}(p^{4})$ chiral corrections. We expect those not to be very important due to the limited phase space available. Other source of error comes from neglecting the anomalous dimension of the light-quark mass. These two sources of error would affect the determination of the line shapes. As we have obtained a pretty good fit for them, we will neglect these sources of error in the following.

The other error we have not incorporated in this analysis, nor in those we will perform later, is the error due to the working hypothesis we use in this paper (saturation of the interpolating field by those with smaller dimensionality). The reason is that we want to see whether such hypothesis is feasible, and if so what is the expected error, by comparing our predictions with experiment.

IV $\bar{Q}Q(2S)\rightarrow\bar{Q}Q(1P)\pi$ hadronic transitions

The $\bar{Q}Q(2S)\to\bar{Q}Q(1P)\pi$ hadronic transition444Let us note that this type of transition was first considered in Ref. Kuang and Yan (1981) for the $\Upsilon(2S)\to h_{b}(1P)\pi^{0}$ decay. However, this particular case turned out not to be kinematically allowed. is zero with the LO Lagrangian in Eq. (1) considered this far. The first nonzero contribution is generated by spin-isospin breaking effects. The leading spin-dependent operators originate from the following pNRQCD operators:

[TABLE]

We remind that the singlet, octet, and hybrid fields in the Lagrangians that appear in this paper are organized in $SU(2)$ spin multiplets. An alternative representation of Eq. (48) in terms of the spins of the quark and antiquark can be found in Eq. (105) of Brambilla et al. (2003). This last representation is the one we will customarily use for the computation of the matrix elements.

The renormalization group improved expression for the matching coefficient $c_{F}$ is known with next-to-leading logarithmic accuracy Amoros et al. (1997); Czarnecki and Grozin (1997). In order to include the leading isospin violation effects, we should no longer consider the light-quark masses degenerate in Eq. (1): we take $m_{u}=2.118$ MeV and $m_{d}=4.690$ MeV instead Tanabashi et al. (2018).

Adding the operator in Eq. (48) to the leading pNRQCD Lagrangian in Eq. (1) and considering isospin violation effects produces the following new terms in the hadronic Lagrangian at LO555There is another possible pseudoscalar operator of the same order as $\langle\chi_{-}\rangle$: $\langle D_{\mu}u^{\mu}\rangle$. However they are both related through the leading order equations of motion $D_{\mu}u^{\mu}=i\left(\chi_{-}-\langle\chi_{-}\rangle/2\right)/2$.:

[TABLE]

Terms that break spin and isospin symmetry simultaneously are not considered, since they produce subleading contributions to the transition we are considering.

Let us match the parton and hadronic description of the spin-dependent mixing operator. In the hadronic EFT, the correlator reads

[TABLE]

and in the partonic version of pNRQCD, we have

[TABLE]

where amp. signals that only amputated contributions are considered. Now we consider $G_{1^{+-}}^{a}$ to be given by a sum of all possible gluonic operators with the same quantum numbers

[TABLE]

Similarly, to the previous section we hypothesize that there is a correlation between the dimensionality of the interpolating operator and the strength of the interpolation with the hybrid, such that higher dimension operators are subleading, so the series can be truncated at LO. Using this truncation and Eq. (6) in Eq. (51), and matching to Eq. (50), we arrive at

[TABLE]

Let us now match the operators with an odd number of pions in the Lagrangian in Eq. (49). In the hadronic theory, we have

[TABLE]

and in pNRQCD,

[TABLE]

where in the last step we have made use of the results of Appendix D.2 for the hadronization of the gluonic operator through the axial anomaly. Comparing Eqs. (54) and (55), we arrive at

[TABLE]

We are interested in investigating the decay $\psi(2S)\to h_{c}(1P)\pi^{0}$, which is given in the hadronic pNRQCD EFT by the diagrams in Fig. 5. The matrix element reads

[TABLE]

with

[TABLE]

In the first term, the $\psi(2S)$ mixes into a hybrid state with $J^{PC}=1^{--}$ but in a spin singlet, then the hybrid decays into a $h_{c}$ conserving the spin but changing the angular momentum. The intermediate hybrid states fulfilling these conditions are $m^{1}\mathcal{P}_{1}$. This term turns out to be small because the second vertex only has contributions from the $\lambda=0$ piece of the hybrid state, which is small for ${}^{1}\mathcal{P}_{1}$ states. In the second term, the order is switched; in this case the allowed hybrid intermediate states are $m^{3}\mathcal{S}_{1}$. Note that all the hybrid states appearing as intermediate states in this transition are associated to the $1^{+-}$ gluelump. As in the two-pion decays, the normalization factors $Z_{B}$ cancels out. The matrix elements are given in Appendix C.

We compute the amplitude for up to $m=4$ for the ${}^{3}\mathcal{S}_{1}$ and up to $m=2$ for the ${}^{1}\mathcal{P}_{1}$ hybrid intermediate states. We use the renormalization group improved expression of $c_{F}(1\leavevmode\nobreak\ {\rm GeV})$=1.12155 up to next-to-leading logarithmic order, where we have used the values $\alpha_{s}(1\leavevmode\nobreak\ {\rm GeV})$ = 0.4798 and $\alpha_{s}(1.496\leavevmode\nobreak\ {\rm GeV})$ = 0.3522. The decay width can be obtained using the amplitude of Eq. (57) in Eq. (99),

[TABLE]

The error analysis has been performed similarly to the previous section. In this case, we do not have error associated to $\kappa$, but still have a dependence on $\Lambda_{1}$ (the lowest lying gluelump mass), and s.p. (the different parametrization for the singlet and hybrid static potentials). We also estimate the error due to the value of the light-quark (l.q.) mass ratio, which we take as

[TABLE]

with the value of $m_{d}/m_{u}$ from the PDG.

This observable is interesting. Unlike previous decays, it does not suffer from the uncertainties associated to the hadronization of the local operator: The axial anomaly does not get ${\cal O}(\alpha)$ corrections, nor there are ${\cal O}(p^{4})$ chiral corrections. We also observe a very weak dependence on variations of the wave function of the hybrid and singlet states. The major error is generated by the uncertainty on the energy difference between singlet and hybrid states. Once this error is taken into account, the result is roughly compatible with experiment.

In principle, for $\bar{Q}Q(2S)\rightarrow\bar{Q}Q(1P)\pi$ transitions, the energy release is small so one can use the twist expansion with no fear if one assumes that $\Lambda_{\rm QCD}\ll mv^{2}$. We then would have an alternative determination with which one can compare. In this respect, let us note that the only previous theoretical estimate was $\Gamma_{\psi(2S)\to h_{c}(1P)\pi^{0}}\sim 15$ eV, from Ref. Voloshin (2008), using the twist expansion ($\Lambda_{\rm QCD}\ll mv^{2}$ and $E\ll mv^{2}$).

V $\bar{Q}Q(2S)\rightarrow\bar{Q}Q(1S)\pi^{0}(\eta)$ hadronic transitions

We now turn our attention to $\bar{Q}Q(2S)\rightarrow\bar{Q}Q(1S)\pi^{0}(\eta)$ hadronic transitions. Since these transitions break spin symmetry, they are zero at LO in the EFT. The leading contribution to these decays is generated by

[TABLE]

which has to be added to Eq. (1). In the hadronic EFT, these operators correspond to

[TABLE]

To determine the Wilson coefficients, we compute the transition amplitude both in the hadronic and partonic versions of the effective theory. Nevertheless, there is one extra subtlety to be taken into account. The $\pi^{0}$ and $\eta$ fields in Eq. (11) mix, and to obtain the physical states the mass matrix needs to be diagonalized. The physical states correspond to

[TABLE]

with $\epsilon$ the mixing angle

[TABLE]

To compute the decay amplitudes, the physical states must be considered. In the transitions with $\eta$ emission, the contribution due to the mixing is subleading, but in the case of $\pi^{0}$ emission both contributions are of the same order. In practice, this amounts to an extra factor $3/2$ in the $\pi^{0}$ decay amplitude. In the hadronic EFT, we have

[TABLE]

and in the partonic pNRQCD,

[TABLE]

where in the last step we have used Eq. (19) truncated to the first term and the anomaly relation in Appendix D.3 [note that this relation does not suffer from ${\cal O}(\alpha_{s})$ corrections]. Matching Eqs. (66) and (67), we obtain

[TABLE]

We are interested in investigating the decays $\psi(2S)\to J/\psi\pi^{0}$, $\Upsilon(2S)\to\Upsilon(1S)\pi^{0}$, as well as the decays $\psi(2S)\to J/\psi\eta$, $\Upsilon(2S)\to\Upsilon(1S)\eta$. The diagrams involved in the decay are drawn in Fig. 5. These contain the $t^{(r1^{--})}$ vertex from Eq. (14) mixing the quarkonia into a hybrid associated to the $1^{--}$ gluelump, and the $P$-wave pion emission vertices from Eq. (62). The former vertex conserves the spin state and is only nonvanishing for $l=\ell$. Therefore, the intermediate hybrid states for the transitions we are considering must be ${}^{3}\mathcal{S}_{1}$. The amplitude for a pion emission then reads

[TABLE]

Note that the factor involving the sum over intermediate hybrid states: $\beta^{(12)}_{r,Q}$, given in Eq. (28), is the same to the one appearing in the two-pion transitions. $\epsilon_{mS}$ stands for the polarization of a $m^{3}S_{1}$ quarkonium state. The decays amplitude to one $\eta$ is

[TABLE]

We have computed the decay widths by summing over final polarizations and averaging over the incoming ones the square of the amplitudes in Eqs. (71) and (72), and inserting the result in Eq. (99). We use relativistic kinematics for the outgoing quarkonium state in the phase space calculation and the renormalization group improved expression of $c_{F}$ with next-to-leading logarithmic accuracy. The expression for charmonium can be found in the previous section. For bottomonium, we have $c_{F}(1\leavevmode\nobreak\ {\rm GeV})$=0.87897 using $\alpha_{s}(4.885\leavevmode\nobreak\ {\rm GeV})$ = 0.2148. $\beta^{(21)}_{r,Q}$ is computed considering up to four hybrid intermediate states. The results obtained for the decays are the following:

[TABLE]

The subindices label the source of the uncertainty. The error analysis is equal to the one performed in the previous section. Again we do not have error associated to $\kappa$, but still have a dependence on $\Lambda_{1}$ (the lowest lying gluelump mass), and s.p. (the different parametrization for the singlet and hybrid static potentials). We also estimate the error due to the value of the light-quark mass ratio.

As in the previous section, these decays do not suffer from the uncertainties associated to the hadronization of the local operator: the axial anomaly does not get ${\cal O}(\alpha_{s})$ corrections, nor there are ${\cal O}(p^{4})$ chiral corrections. On the other hand, unlike in the previous section, we find a significant dependence on variations of the wave function of the hybrid and singlet states. Indeed, the dependence is the same as the one we had in Sec. III.2. The error generated by the uncertainty on the energy difference between singlet and hybrid states is of the same order. Nevertheless, these error estimates are not large enough to account for the difference with experiment, particularly for the charmonium case. We later retake this issue when we consider ratios of decay rates.

VI Ratios

So far we have seen that the hadronic $S$- to $P$-wave heavy quarkonium transitions are roughly compatible with theory within errors. For decays to $S$-wave heavy quarkonium, we have seen that the overall magnitude of our predictions is smaller than experiment, specially for the $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi^{0}(\eta)$ decays. Nevertheless, the normalization of the decays are the most uncertain object in our predictions. Therefore, we expect several uncertainties to cancel for ratios. In fact, we have already seen in Sec. III.1 that the normalized differential decay rates are well described by theory. Going further in this direction, we may try to explore the magnitude of the corrections associated to the different approximations we have made in this paper by studying different ratios. This we do in the following.

The ratios of the $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi^{0}(\eta)$ transitions do not depend on $\beta^{(12)}_{r,Q}$, and are completely determined at LO by chiral symmetry. Indeed the same result is obtained using a pure chiral approach or by using the twist expansion (though the use of the latter suffers from the same drawback as its use in $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi\pi$ transitions Luty and Sundrum (1993)), as studied previously in Refs. Ioffe and Shifman (1981); Donoghue and Wyler (1992); Casalbuoni et al. (1993). What changes is the overall coefficient, which cancels in the ratio. The only uncertainties affecting the ratios are the ones associated to the light-quark mass values, and possible chiral corrections to the amplitude affecting the ratio $F_{\eta}/F_{\pi}$ and the mixing angle $\epsilon$. Generically, we expect these to be of $\mathcal{O}((m_{\pi},m_{K})/\Lambda_{\chi})\sim 14-50\%$. We can check this by comparing the theoretical and experimental ratios (for these and for the experimental ratios below we add the error quadratically)

[TABLE]

Since the theoretical values are compatible experiment we conclude that chiral corrections are not needed at this level of precision.

It is also interesting to compare the theoretical and experimental ratios of the two-pion transitions of Sec. III over the one-pion or eta transition computed in Sec. V, since this set of ratios is also independent of $\beta_{r,Q}^{(12)}$, or, in other words, of the heavy quarkonium dynamics (which is also the case with the twist expansion). The hadronization of the two-pion production using the energy-momentum tensor anomaly has corrections of $\mathcal{O}(\alpha_{s})$ at a low-energy scale, as well as contributions of the anomalous dimension of the light-quark mass operator, both of which could be potentially large.666In principle there are also ${\cal O}(p^{4})$ chiral corrections. Nevertheless, we expect those to be comparatively small due to the limited phase space available. On the other hand the determination of the one single pion or eta matrix element through the axial anomaly only has chiral corrections, as discussed in the previous paragraph, and, in the case of one pion in the final states, the uncertainty on $m_{u}/m_{d}$. Hence, we explore whether comparing the ratios of the two- and one-pion (or eta) transitions to their experimental values allows us asses the size of the corrections to Eqs. (157) and (158).

[TABLE]

The differences with experiment are large, of the order of 75%-60% for the charmonium and bottomonium, respectively. In principle, we expect the most important uncertainties of these ratios to be due to the effect of higher order operators in the interpolating function of the hybrids and of the neglected ${\cal O}(\alpha_{s})$ corrections in the hadronization of the local operators (this latter source of error only applies to the $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi\pi$ transitions). In this respect, it is interesting to note that we expect the ${\cal O}(\alpha_{s})$ corrections in the hadronization of the local operators to be largely independent of the bound state dynamics. Along this line of thought, it is rewarding that we can roughly get agreement with experiment (within errors) both for charmonium and botttomoniun using the same correcting factor $\sim 1/3.5$. This factor could be understood as evaluating the $\beta(\alpha_{s})$ function at a low scale: $-\frac{2\pi\beta(\alpha_{s})}{\beta_{0}\alpha_{s}^{2}}\sim 3.5$. Still, it is premature to draw any conclusion out of this. Note also that these estimates are significantly affected by the uncertainties of $\kappa$, and, for the $\pi^{0}$ case, by the uncertainty of $m_{u}/m_{d}$. On the other hand, the strong dependence on the light-quark masses makes these observables interesting for possible determination of the light-quark masses.

For the ratios considered so far, Eqs. (77)-(82), the prediction of the twist expansion is the same to the one we have found here. This is not so for the following ratio:

[TABLE]

In principle, we expect that for this ratio most of the neglected ${\cal O}(\alpha_{s})$ corrections in the hadronization of the local operators vanish. This observable can be considered a rough measure of $\beta_{r,b}^{(21)}/\beta_{r,c}^{(21)}$. The agreement with experiment is remarkable: below 20%, and could be accounted for by the quoted errors.

We can also consider the following ratios:

[TABLE]

Again, the twist expansion would yield a different prediction for these ratios, and, again, these ratios can be considered a measure of $\beta_{r,b}^{(21)}/\beta_{r,c}^{(21)}$. The agreement with experiment is quite reasonable, with difference of the order of 30%. In order to assess whether this error and the error of Eq. (83) are really related to $\beta_{r,b}^{(21)}/\beta_{r,c}^{(21)}$ or to something else, one may consider the double ratios

[TABLE]

These double ratios are independent of $\beta_{r,b}^{(21)}/\beta_{r,c}^{(21)}$. We also expect several other uncertainties to cancel: the dependence on the light-quark masses cancels; in principle, we would also expect most of the effect due to unaccounted for ${\cal O}(\alpha_{s})$ corrections in the hadronization of the local operators vanish, since the phase space is similar. An indirect reflection of this is that the dependence on $\kappa$ nearly vanishes. Nevertheless, this should be studied with more detail to make this statement more quantitative. For instance, for the $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\eta$ transition, there is little phase space free. As a result, the phase space computation yields very different values for the bottomonium and charmonium transitions. Overall, we consider these double ratios as the cleanest objects to perform dedicated studies of the computational scheme discussed in this paper (and of the twist expansion, which yields the same prediction for these double ratios). For now, we just want to highlight that we find remarkably good agreement with experiment. This may indicate that the overlap of higher dimensional operators with the hybrids is small. In this respect, the possible incorporation of $B_{c}$ transitions to these analyses could be of much help.

VII Conclusions

We have studied one- and two-pion transitions between quarkonium states below threshold. We have used the weakly coupled version of pNRQCD, which assumes that $mv\gg\Lambda_{\rm QCD}$. On the other hand, our analysis does not need to constrain the relative size of $\Lambda_{\rm QCD}$ with respect to $mv^{2}$. For definiteness, we have generically considered the situation $mv^{2}\ll\Lambda_{\rm QCD}$. This EFT of QCD for heavy quark-antiquark systems incorporates the $1/m_{Q}$ and multipole expansions systematically. We then write the hadronic representation for this EFT in terms of the singlet, representing the heavy quarkonium, hybrids and pion fields. In weakly coupled pNRQCD, the hybrid states consist of a color-octet heavy quark-antiquark field and a gluonic excitation operator, called gluelump, characterized by a $J^{PC}$. In order to explicitly compute the matching between both versions of the theory, we assume that the gluelump operators overlap predominantly with the lowest dimension gluonic operator with the same $J^{PC}$. The matching is completed by making use of low-energy theorems generated by the axial and energy-momentum tensor anomalies to obtain the local matrix elements for pion production by gluonic operators.

In this framework, we have computed the $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi\pi$ transitions and compared with the description obtained solely from chiral symmetry, as well as with experimental data. We do so for the normalized decay width spectra and for the total widths of charmonium and bottomonium. In the case of the normalized decay width spectra, both our approach and the purely chiral description fit the data well. However, the chiral description depends on two parameters that cannot be strongly constrained from the experimental data, whereas our approach only depends on one parameter $\kappa$. Our best fit to the combined charmonium and bottomonium data yields the value quoted in Eq. (41) for $\kappa$.

Our computational scheme produces the same theoretical expression as the twist expansion Voloshin (1979) for the normalized decay width spectra computed in Sec. III.1. Note, however, that the twist expansion requires that $\Lambda_{\rm QCD}\ll mv^{2}$, and that the energy release $E\ll mv^{2}$, the latter condition is never fulfilled for two-pion transitions. Lower values for $\kappa$ than the ones in Eq. (41) are found in Bai et al. (2000); Chen (2019). The main source for the discrepancy is due to ${\cal O}(m_{\pi}^{2})$ terms not included in these analysis.

Our prediction for the total two-pion transition decay width depends on $\kappa$ [which we take from Eq. (41)], and on $\beta_{r,Q}^{(12)}$. The latter coefficient suffers from large uncertainties. It involves a double sum: one over the gluelump states with the same quantum numbers and, for each gluelump, other sum over the states solution of its associated Schrödinger equation. In this paper, we have made a first estimate considering only the first gluelump and, for it, summing the first few states solution of the Schrödinger equation. Moreover, the states solution of the Schrödinger equation have a significant overlap with long distances. This makes the result quite dependent on the shape of the potentials. We have estimated this dependence considering different fit functions for the singlet and hybrid potentials. Overall, our estimates for the total decay widths $\Gamma_{\psi(2S)\to J/\psi\pi^{+}\pi^{-}}$, and $\Gamma_{\Upsilon(2S)\to\Upsilon(1S)\pi^{+}\pi^{-}}$ can be found in Eqs. (46) and (47).

Enlarging our Lagrangian to include spin-dependent and isospin breaking operators we can use the same procedure to compute the width for $Q\bar{Q}(2S)\to Q\bar{Q}(1P)\pi^{0}$ transitions. Our prediction for the total width $\Gamma_{\psi(2S)\to h_{c}(1P)\pi^{0}}$ can be found in Eq. (59). This observable is interesting. Unlike previous decays, it does not suffer from the uncertainties associated to the hadronization of the local operator: the axial anomaly does not get ${\cal O}(\alpha)$ corrections, nor there are ${\cal O}(p^{4})$ chiral corrections. We also observe a very weak dependence on variations of the wave function of the hybrid and singlet states: this object depends on $\beta_{\sigma}^{(12)}$. The major error is generated by the uncertainty on the energy difference between singlet and hybrid states. Once this error is taken into account the result is roughly compatible with experiment.

Finally, we considered $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi^{0}(\eta)$ transitions. Our results for the total decay widths $\Gamma_{\psi(2S)\to J/\psi\pi^{0}}$, $\Gamma_{\psi(2S)\to J/\psi\eta}$, $\Gamma_{\Upsilon(2S)\to\Upsilon(1S)\pi^{0}}$, and $\Gamma_{\Upsilon(2S)\to\Upsilon(1S)\eta}$ can be found in Eqs. (73)-(76).

Overall, we find that the $S$-wave hadronic transitions to P-wave heavy quarkonium are roughly compatible with theory within errors. For decays to $S$-wave heavy quarkonium, the magnitude of our predictions is smaller than experiment, especially for the $Q\bar{Q}(2S)\to Q\bar{Q}(1S)\pi^{0}(\eta)$ decays.

We have also computed the ratios of the above transition rates in Sec. VI. We expect a more solid prediction for them. Depending on the ratio, different qualitative information on the theory can be obtained. The ratios of Eqs. (77) and (78) can be considered to be a test of chiral symmetry, as they are independent of the bound state dynamics, and there is no error associated to the hadronization of local gluonic operators. Good agreement with experiment is found. This may indicate a good behavior of chiral symmetry for these decays. Alternative ratios one may consider are Eqs. (79)-(82). These ratios are also independent of the bound state dynamics but on the other hand suffer from unquantified errors due to the hadronization of the local gluonic operators. Large differences with experiment were found, which however could be accommodated by a single constant for bottomonium and charmonium. This agrees with expectations that the main corrections are independent of the bound state dynamics. These ratios and the ratios Eqs. (77) and (78) are independent of the bound state dynamics, i.e., they are independent of $\beta_{r,Q}^{(12)}$ and yield the same result as the twist expansion. On the other hand, the overall coefficients which, in principle, are calculable in both cases, are different for the total transition rates.

One may consider ratios that yield different predictions than the twist expansion. These are Eqs. (83)-(85). We expect these ratios to be quite independent of the error associated to the hadronization of the local gluonic operator. On the other hand they depend on the bound state dynamics through the ratio: $\beta_{r,b}^{(12)}/\beta_{r,c}^{(12)}$. For Eq. (83), good agreement is found (below 20%), and is quite reasonable for the other two ratios: of the order of 30%. This may indicate that our evaluation of $\beta_{r,b}^{(12)}/\beta_{r,c}^{(12)}$ is quite reasonable. Finally, we have considered double ratios in Eq. (86). In principle, these are the cleanest objects to compute [leaving aside Eqs. (77) and (78)]. They are independent of the bound state dynamics, and we also expect them to be rather independent of the uncertainties associated to the hadronization of the local gluonic operator. The agreement with experiment is quite good for them.

We believe there is room for improvement, and leave for future work more dedicated analysis, in particular of the ratios and, specially, of the double ratio discussed in Eq. (86). This last object may allow a quantitative analysis of the validity of our computational scheme (and alternatively of the twist expansion). It would also be interesting to apply our computational scheme to the two- (or one) pion hadronic transition of the $B_{c}$. On a later stage, the possible incorporation of tetraquarks and B-meson loops for heavy quarkonium states near, or above threshold, should also be studied in greater detail. Actually, if enough precision is reached, discrepancies with experiment may hint at the need of incorporating those states.

Acknowledgments

We thank E. Guido and R. Mussa for sharing the experimental data of Ref. Guido et al. (2017) with us. We thank P. Petreczky for sharing the lattice data of Ref. Bazavov et al. (2018) with us. We thank Feng-Kun Guo for providing us the data of Alexander et al. (1998). This work was supported in part by the Spanish Grants No. FPA2017-86989-P and No. SEV-2016-0588 from the Ministerio de Ciencia, Innovación y Universidades, and the Grant No. 2017SGR1069 from the Generalitat de Catalunya. J.T.C acknowledges the financial support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska–Curie Grant No. 665919.

Appendix A Differential decay widths

The most general amplitude for the two-pion transitions at ${\cal O}(p^{2})$ in the chiral counting and using spin symmetry is given in Eq. (26). It is convenient to write it as

[TABLE]

where $m_{\pi\pi}^{2}=(p_{+}+p_{-})^{2}$ is the dipion invariant mass, and $p_{+}$, $p_{-}$ refers to the momentum of the $\pi^{+}$ and $\pi^{-}$, respectively. In the reference frame of the decaying quarkonia, we find

[TABLE]

with

[TABLE]

where $m_{H_{i}}$ and $m_{H_{f}}$ are the masses of initial and final quarkonium respectively. In the nonrelativistic approximation of the final quarkonium momentum, the above expressions reduce to

[TABLE]

Using Eq. (89) in Eq. (88), we can write

[TABLE]

The differential decay width is

[TABLE]

and after integrating $\theta$ it reads

[TABLE]

To obtain the total decay width we integrate numerically

[TABLE]

In the one-pion transition, the momenta are fixed by momentum conservation, in particular the final pion momentum is $|\bm{p}_{\pi}|=\rho(m_{\pi})$, and the decay width is given by

[TABLE]

with $\rho(m_{\pi})$ is given in Eq. (92) replacing $m_{\pi\pi}$ by $m_{\pi}$.

Appendix B Singlet and hybrid bound states

The heavy quarkonium wave function reads

[TABLE]

where $\phi^{(n)}(r)$ is the solution of

[TABLE]

It is not our aim in this paper to make a detailed analysis and optimization of the solutions of the bound states but rather to make a qualitative analysis. For the numerical evaluations of the hybrid and quarkonium spectra, we use the values $m_{c}(1{\rm GeV})=1.496$ GeV and $m_{b}(1{\rm GeV})=4.885$ GeV Peset et al. (2018) for the heavy quark masses in the RS*′* scheme Pineda (2001). To explore the sensitivity of the matrix elements to the shape of the potential, we have considered two versions of the potentials. The first is the potential in Eq. (15) with the nonperturbative constant $b_{\Sigma_{g}^{+}}$ obtained by fitting to the lattice data from Ref. Bazavov et al. (2018) up to $r=1$ fm together with a free energy offset, $b^{\rm offset}_{\Sigma_{g}^{+}}$. The second version of the potential is a fit of the same lattice data to a function constrained to reproduce Eq. (15) in the short distance and a to have a linear behavior in the long distance, while overall adjusting well to the data. The function chosen is

[TABLE]

For the perturbative static singlet potential we take the $\mathcal{O}(\alpha^{3}_{s})$ result777The ${\cal O}(\alpha_{s})$ term was computed in Fischler (1977), the ${\cal O}(\alpha_{s}^{2})$ in Schroder (1999); Peter (1997), the ${\cal O}(\alpha^{3}_{s})$ logarithmic term in Brambilla et al. (1999), the light-flavor finite piece in Smirnov et al. (2008), and the pure gluonic finite piece in Anzai et al. (2010); Smirnov et al. (2010). evaluated at $\nu=5.6$ GeV (the shortest available scale). The parameters obtained from the fits are the following

[TABLE]

In Fig. 6, we plot the lattice data together with the fitted Eq. (15) (with the offset) and Eq. (102) potentials as dashed and solid lines, respectively. The spectra obtained are displayed in Tables 3 and 4 for the potential in Eqs. (15) and in (102), respectively. The origin of energies is chosen in order to reproduce the experimental spin average of the ground state.

The hybrid wave functions are obtained following the procedure described in Ref. Berwein et al. (2015), by solving the coupled Schrödinger equations involving the potentials $V^{(0)}_{\Sigma^{-}_{u}}(r)-V^{(0)}_{\Pi_{u}}(r)$ or $V^{(0)}_{\Sigma^{+\prime}_{g}}(r)-V^{(0)}_{\Pi_{g}}(r)$ generated by the $1^{+-}$ and $1^{--}$ gluelump at short distances, respectively [see Eq. (16)]. There are two types of solution of $\Psi_{1\lambda}$ corresponding to states with opposite parity $\left(\Psi^{mjm_{j}\ell s}_{+}\right)_{\lambda}$ and $\left(\Psi^{mjm_{j}\ell s}_{-}\right)_{\lambda}$,

[TABLE]

The parity and charge conjugation of these states corresponding to the $\Psi_{\pm}$ solutions are

[TABLE]

with $P_{G}$ and $C_{G}$ the parity and charge conjugation of the gluelump associated to the states.

The angular wave functions of the hybrid states are the eigenfunctions of

[TABLE]

which for $\lambda=0$ corresponds to the usual differential equation for spherical harmonics and, therefore, $v_{l\,m_{l}}^{0}=Y_{l\,m_{l}}$. Note that we have chosen to represent the angular momentum quantum numbers of the hybrids and quarkonia by $\ell$ and $l$, respectively to highlight the difference in angular wave functions. When using spectroscopic notation to specify the quarkonium states we use the usual $S,P,D,\dots$ for $l=0,1,2,\dots,$ and for hybrid states, we use $\mathcal{S},\mathcal{P},\mathcal{D},\dots$ for $\ell=0,1,2,\cdots$.

The $\lambda^{2s+1}\ell_{j}$ spin-angular hybrid wave functions can be shown to be

[TABLE]

where $\hat{e}$ are the polarization vectors

[TABLE]

and the tensors $h^{ij}_{2m_{J}}$ are traceless, completely symmetric, and normalized as

[TABLE]

The ${}^{2s+1}l_{j}$ spin-angular wave functions for quarkonium correspond to the $\lambda=0$, ${}^{2s+1}\ell_{j}$ hybrids wave functions.

The radial wave function in Eq. (111) is obtained by (numerically) solving the following coupled radial Schrödinger equations:

[TABLE]

For the special case $\ell=0$ the equations decouple and in fact $\psi_{+}^{(m)}$ becomes irrelevant since $v^{\pm 1}_{00}$ does not exist. The radial wave function in Eq. (118) is obtained from the uncoupled radial Schrödinger equation

[TABLE]

The hybrid static potentials, $V_{\Pi_{u}}$, $V_{\Sigma^{-}_{u}}$, $V_{\Pi_{g}}$, and $V_{\Sigma^{+\prime}_{g}}$, match to the NRQCD heavy quark-antiquark static energies, which have been computed on the lattice Juge et al. (2003); Bali et al. (2000); Capitani et al. (2019) and are given up to NLO in the multipole expansion by the potential in Eq. (16). As in the static potential for the singlet, we have considered two versions of the potentials in order to assess the sensitivity of the matrix elements to the shape of the potential. The first is the potential in Eq. (16) with the nonperturbative constant $b_{k\lambda}$ obtained by fitting to the lattice data from Ref. Juge et al. (2003) up to $r=1$ fm together with a free energy offset. The second version of the potential is a fit of the same lattice data up to $r\leq 2$ fm by a function constrained to reproduce Eq. (16) in the short distance and a to have a linear behavior in the long distance while overall adjusting well to the lattice data. The function is

[TABLE]

The perturbative static octet potential is taken at $\mathcal{O}(\alpha^{3}_{s})$ from Ref. Kniehl et al. (2005), and evaluated at $\nu=2.6$ GeV (the shortest available scale). The results for the hybrid spectrum for $V_{\Pi_{u}}$-$V_{\Sigma^{-}_{u}}$ associated to the $1^{+-}$ gluelump are discussed at length in Ref. Berwein et al. (2015).888In Ref. Berwein et al. (2015), the choice of Eq. (132) was different but the resulting potentials are completely equivalent in practice. The parameters obtained from the fits are the following:

[TABLE]

The fits of the potentials $V_{\Pi_{g}}$-$V_{\Sigma^{+\prime}_{g}}$, associated to the $1^{--}$ gluelump, are shown in Fig. 7 in solid and dashed lines and the corresponding spectra in Tables 5 and  6 for the potentials in Eqs. (132) and (16), respectively. The origin of energies is chosen, as in Ref. Berwein et al. (2015), such that the short distance limit of $V_{\Pi_{u}}$ and $V_{\Sigma^{-}_{u}}$ tend to $V_{o}^{0}+\Lambda_{1^{+-}}$ with the value of $\Lambda_{1^{+-}}=0.87(15)$ GeV computed in Ref. Bali and Pineda (2004).

Appendix C Matrix Elements

We now compute the matrix elements that appear in the hadronic transitions. We suppress the labels $k^{PC}=1^{+-}$ and $k^{PC}=1^{--}$ on the hybrid fields since it does not affect the calculation. For simplicity, we denote the quarkonia and hybrid wave functions just as $S_{n}$ and $\Psi_{m}$, respectively:

[TABLE]

where we have used that $\hat{\bm{r}}^{\dagger}_{\lambda}\cdot\bm{r}=r\delta_{\lambda 0}$. For $\lambda=0$ $v^{0}_{lm}=Y_{l\,m_{l}}$ so we can use the spherical harmonics orthogonality relations

[TABLE]

Let us work out some specific matrix elements necessary for the $\psi(2S)\to h_{c}\pi^{0}$ transition,

[TABLE]

[TABLE]

where we have used that hybrid states in $\ell=0$ must have $\lambda=0$ and also that

[TABLE]

For $P$-wave pion (or eta) emission, we need the following matrix element:

[TABLE]

and due to the symmetry of the angular integral with respect to the indices [see Eq. (141)], we can also derive

[TABLE]

Appendix D Hadronization of local gluonic operators

D.1 Two-pion production through the anomaly of the energy-momentum tensor

The traceless part of the energy momentum tensor for QCD is

[TABLE]

The trace of the energy momentum tensor reads

[TABLE]

where the $\gamma_{i}$ are the anomalous dimension of the $\bar{q}_{i}q_{i}$ operator and $\beta$ is the QCD $\beta$-function

[TABLE]

The first term in Eq. (145) is generated by the conformal anomaly, and the second by the explicit breaking of the scale invariance due to the quark masses.

The matrix element $\langle\pi^{+}(p_{+})\pi^{-}(p_{-})|\theta_{\mu\nu}|0\rangle$ can be determined up to quadratic order in $p_{+}$, $p_{-}$ and $m_{\pi}$ from Voloshin (2008); Novikov and Shifman (1981)

• •

Conservation on the mass shell: $(p_{+}+p_{-})^{\mu}\langle\pi^{+}(p_{+})\pi^{-}(p_{-})|\theta_{\mu\nu}|0\rangle_{p^{2}_{+}=p^{2}_{-}=m^{2}_{\pi}}=0$.

• •

Normalization: $\langle\pi^{+}(p)\pi^{-}(-p)|\theta_{\mu\nu}|0\rangle_{p^{2}=m^{2}_{\pi}}=2p_{\mu}p_{\nu}$.

• •

Adler zero condition: $\langle\pi^{+}(p)\pi^{-}(0)|\theta_{\mu\nu}|0\rangle_{p^{2}=m^{2}_{\pi}}=0$.

and reads as

[TABLE]

The square of chromoelectric and chromomagnetic fields can be written as

[TABLE]

with $v_{\mu}=(1,\,\bm{0})$ and

[TABLE]

We can rewrite the chromoelectric and chromomagnetic fields in terms of $\theta^{\mu}_{\mu}$,

[TABLE]

The problem of computing the matrix elements $\langle\pi^{+}(p_{+})\pi^{-}(p_{-})|\bm{E}^{2}|0\rangle$ and $\langle\pi^{+}(p_{+})\pi^{-}(p_{-})|\bm{B}^{2}|0\rangle$ is transformed into the computation of $\langle\pi^{+}(p_{+})\pi^{-}(p_{-})|\theta^{g}_{\mu\nu}|0\rangle$, $\langle\pi^{+}(p_{+})\pi^{-}(p_{-})|\theta^{\mu}_{\mu}|0\rangle$, and $\langle\pi^{+}(p_{+})\pi^{-}(p_{-})|\sum_{i}m_{i}\bar{q}_{i}q_{i}|0\rangle$. These can be evaluated in $\chi$PT except for a possible normalization, and they correspond to the traceless part of the energy-momentum tensor, the trace of the energy-momentum tensor, and the scalar current, respectively. Let us start by the trace of the energy-momentum tensor

[TABLE]

If we compare with Eq. (148), we determine that $V_{1}(\mu)=1$. The matrix element of the scalar current reads

[TABLE]

where one can think the parameter $B$ ($\chi=2B\hat{m}\mathbbm{1}$) as taking the role of the normalization parameter. Finally, for $\theta^{g}_{\mu\nu}$, we have no way to determine the normalization

[TABLE]

Adding everything together we arrive at

[TABLE]

where we have defined $\kappa=\alpha_{s}\beta_{0}V_{2}(\mu)/(6\pi)$ as in Ref. Novikov and Shifman (1981), and we arrive at the results of Refs. Voloshin (2008); Brambilla et al. (2016). We emphasize that in this relation we have neglected higher order $\alpha_{s}$ corrections to the beta function. These $\alpha_{s}$ corrections are computed at a low-energy scale, which makes their computation unfeasible. Still, they could be parametrized as an overall coefficient (redefining $\kappa$) of the equation. In Eq. (158), we have also neglected the anomalous dimension of the $\bar{q}q$ operator. Otherwise, the coefficient proportional to $m^{2}_{\pi}$ becomes arbitrary. Whereas this is customary done, this should be further investigated.

D.2 Single pion creation from the axial anomaly

The axial current matrix elements are defined as

[TABLE]

with $\lambda$ the Gell-Mann matrices and $q=(u,\,d,\,s)$. The divergence of the axial current reads

[TABLE]

with

[TABLE]

and the same holds replacing $u$ by $d$ or $s$. Then, we have

[TABLE]

Since the diagonal flavor current does not interpolate the $\pi^{0}$ we have

[TABLE]

and using Eq. (162) in Eq. (163), we obtain

[TABLE]

which is the result from Ref. Gross et al. (1979). Finally, we note that

[TABLE]

A similar analysis for the $\eta$ produces

[TABLE]

We remark that, compared with the prediction of the previous section, these results do not have neither ${\cal O}(\alpha_{s})$ nor ${\cal O}(p^{4})$ corrections.

D.3 Single $P$-wave pion anomaly

We first consider the matrix elements of following gluonic operators

[TABLE]

where we have used the Schouten identity

[TABLE]

Due to the Jacobi identity

[TABLE]

we have

[TABLE]

and similarly $X^{\prime}=2Y^{\prime}$.

The sum of the two matrix elements is proportional to the pion momentum:

[TABLE]

The sum of the matrix elements in Eqs. (167)-(168) is

[TABLE]

where we have used Eq. (169) to reexpress Eq. (168). In order for Eq. (173) to be compatible with Eq. (172), the relation $Y=Y^{\prime}$ must hold. Therefore, we arrive at

[TABLE]

Contracting both sides of the last equality with $\epsilon^{\mu\nu\lambda\sigma}/2$,

[TABLE]

the matrix element is precisely the one given in Eq. (164).

Now we can write the matrix elements we are interested in,

[TABLE]

We note that

[TABLE]

differs from Refs. Voloshin (2003, 2008) by a factor $1/2$.