Novel implementation of the multipole expansion to quarkonium hadronic transitions

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.


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 non-linear realization of) chiral symmetry. This was first implemented by Cahn & Brown in Ref. [1]. 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 lead 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. [2]. This gave a good description of the width spectrum, and also of the pion helicity angle distribution [3,4].
Alternatively, in Ref. [5], 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. [6] connected these non-local gluonic matrix elements with chiral symmetry by parameterising 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. [7], Voloshin introduced an additional expansion consisting of an operator product expansion of the non-local 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 [8,9] and energy-momentum tensor anomalies [10][11][12] 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. [13], 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. [14], 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 m Q v 2 , as well as Λ QCD mv 2 .
We want to retake this discussion within the context of effective field theories (EFT's), 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) [15,16].
This allows us to obtain the multipole expansion in a controlled way, since we assume that 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. [16], 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 be also be interesting to study whether this approximation can related with the analysis made in [17]). We will see that such hypothesis plus the chiral low-energy theorems (generated by the axial and energy-momentum tensor anomalies) leads 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 Λ QCD , except the one we already imposed for the multipole expansion: In other words, our computational framework would still be valid even if 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. [18][19][20][21] 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 QQ(2S) → QQ(1S)ππ 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, QQ(2S) → QQ(1P )π in Sec. IV, and QQ(2S) → QQ(1S)π 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 S and O are the quark singlet and octet fields respectively normalized with respect to color as S = S1 c / √ N c and O = O a T a / √ T F . They should be understood as functions of t, the relative coordinates r, and the center of mass coordinates R of the heavy quarks. The trace should be understood as a double trace in colour 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 = 1 √ 2 (S · σ + S η I). All the fields of the light degrees of freedom in Eq. (1) are evaluated at R and t; in particular, G µν a ≡ G µν a (R, t), q i ≡ q i (R, t) and 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 =m ≡ mu+m d 2 with i = u, d. At leading log [22] and next to leading log [23] 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.
For later convenience, we will restrict the discussion to L = 1 gluelump states. We define the gluelump operators, G ia k , as the color-octet gluonic operators that generate the eigenstates of H in the presence of a local heavy-quark-antiquark octet source: where a is the color index, k labels the quantum J P C 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 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 1 .
Representations of the cylindrical symmetry group can be constructed by projecting the 1 The symmetry group is D ∞h , with P replaced by CP.
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: where summation over index i is implied. P i kλ is a projector that acts onto the gluelump angular momentum and projects it into an eigenstate of K ·r (where K is the angular momentum operator for the gluelump) with eigenvalue λ.
It is useful to project the pNRQCD Lagrangian onto the Fock subspace spanned by the |R, r; k, λ states where Ψ kλ (t, r, 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 and the discrete symmetry transformations for the S, O and Ψ i κ fields are given in Table I. 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, R), which (for SU (3)) may be taken as where R ∈ SU R (N ) and L ∈ SU L (N ). Related useful matrices are u, defined from u 2 = U , and Table I. Transformation properties of the heavy quarkonium and gluonic fields, and the projection vectors, under discrete symmetries. The octet field O has the same transformation properties as S. The Ψ i k transform as O combined with the k P C 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λ =r λ · Ψ k . For this reason we give the transformation properties of the unprojected fields Ψ i k . Note that the difference in the transformation ofr λ with respect to Ref. [18] are due to the different definition ofr + . where χ = 2Bdiag(m,m, m s ), with B being related to the vacuum quark condensate. In the isospin limit, the pion mass is m 2 π = 2Bm. The transformation properties can be found in table II.
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 Note that for each hybrid channel one should also include the excitations. The fields S and Ψ kλ should be understood as depending on t, r and R. The pion fields depend on t and R.
Note that Tr[] now only stands for the trace over spin indices. A 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. [16], they correspond to Σ + g , Π g , Σ − u , Π 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 r · δL(Ψ kλ , Ψ k λ ), where δL(Ψ kλ , Ψ k λ ) is bilinear in the hybrid fields and transforms as a 1 −− vector. Note that in general k = 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 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 2 c . 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 O(mv 2 ) or smaller (which in general we will consider them to be smaller than Λ QCD ) this will not be necessary.
In the above discussion we have not included tetraquarks (like the Z b or Z c ) , nor Qq-Qq 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. threshold. This is the situation discussed in Ref. [24]. Then, indeed, the scaling in N c of Γ 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 Γ is very small. In this scenario the mixing with Qq-Qq 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. [25,26] it has been advocated that tetraquarks and Qq-Qq loops may play an important role for some observables like the QQ(3S/4S) → QQ(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 Qq-Qq states will be postponed to future work.
In this paper we will generally consider that the quarkonium binding energies fulfill m Q v 2 Λ QCD . This allows us to write the static singlet potential up to O(r 2 ) as The coefficient b Σ + g can actually be determined in terms of the coefficients t (r1 −− ) from the Lagrangian in Eq. (14), since 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 Λ 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 Nevertheless, in this case, we cannot guaranty that O(r 2 ) terms in the Lagrangian will not contribute to the hybrid potential. Therefore, we cannot determine b kλ from the 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 Schroedinger equation with different variants of the singlet and hybrid potential based on Eq. (15) and Eq. (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 singlethybrid mixing term in the hadronic theory is given by which is matched to the following correlator of pNRQCD in term of quarks and gluons: where amp. signals that only amputated contributions are considered (overall δ(r − r) are also factored out). To evaluate the matrix element, we consider the interpolating field for G a 1 −− to be given by a sum of all possible local gluonic operators with the same quantum numbers: Now, we then hypothesize that there is a correlation between the dimensionality of the The operators with an even number of pions in Eq. (14) are matched in a similar way.
In the hadronic pNRQCD and the corresponding correlator in the partonic pNRQCD reads where in the last step we use Eq. (19) truncated at the first term and hadronize the resulting gluonic matrix element using Eq. (D15) which uses the anomaly relation of the energymomentum tensor of QCD [11,27]. The derivation of this equation is reviewed in Appendix D 1.

III.QQ(2S) →QQ(1S)ππ 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 O(p 2 ) in the chiral counting reads where p + , p − are the momentum of the π + and π − 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. [2]. Imposing heavy quark spin symmetry sets the g 2 low-energy constant of Ref.
with β (n n) 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 + 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 κ that appears in the couplings in Eq. (23)-Eq. (25) and has its origin in the hadronization of the gluonic operator. In the next section we will fix this parameter by fitting the normalized differential decay width spectrum. Eq. (27) yields the following prediction for the chiral parameters a i : We now confront the above predictions with experiment. Since the experimental data on dΓ dmππ (where m ππ = (p + + p − ) 2 is the dipion invariant mass) is normalized to an unknown constant, it is convenient to fit the theoretical expressions to the normalized differential decay width: 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 ψ(2S) → J/Ψπ + π − and Υ(2S) → Υ(1S)π + π − . The experimental data for the former is taken from Ref. [28], and for the latter from Ref. [29]. The total decay widths are taken from Ref. [30].

A. Line-shape analysis
The chiral fit to 1 Γ dΓ dmππ produces the following ratios of parameters: where the superindex c and b label the results for the ψ(2S) → J/Ψπ + π − and Υ(2S) → Υ(1S)π + π − transitions respectively. These fits are performed using the relativistic kinemat-  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. [2] we observe some differences. The authors in Ref. [2] 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 for the charmonium transition, though not for the bottomonium transition. The ratio a 1 /a 2 corresponds to 1 + g 1 /(2g) in the notation of Ref. [2] and their best fits correspond to a c 1 /a c 2 ∼ 1.2 and a b 1 /a b 2 ∼ 1.1. This is consistent with our uncertainty in the first case but slightly outside 1σ for the second.
One should keep in mind that the experimental data source in Ref. [2] is Ref. [31], whereas we use more recent data [29] We now use the amplitude in Eq. (27) to obtain the normalized differential decay width and fit it to the experimental data. Eq. (27) is equivalent to Eq. (26) with the parameters a i taken as in Eqs. (29)- (31). The normalized decay width is independent of β (12) r,Q and its 2 Very similar results to those in Ref. [2] were obtained in Ref. [32] using more recent experimental data [4,33] including pion final state interactions through an unitarization of χPT.
functional form depends only on the parameter κ. 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 κ 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 κ: The range of values correspond to the range with χ 2 − χ 2 min ≤ 1. If we use the relativistic 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 Eq. (37) and Eq. (40) as our default value. In any case, it is remarkable that all fits yield similar values for κ, which we take as a confirmation of the independence of κ 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. Intriguingly, older data for these transitions, see Ref. [31], showed even closer line shapes between charmonium and bottomonium. This reflects on sizable differences between older and more recent data the origin of which should be better understood.
Previous fits of the decay width spectrum of ψ(2S) → J/ψππ has been carried out in Ref. [4] using the transition amplitude from Ref. [11], and in Ref. [34] in which the pion final state interactions were taken into consideration through an Omnès function. The reported values are κ c = 0.186 (3) and κ 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 O(p 2 ) pion mass contribution: In Ref. [34] the coefficient a 3 is set to zero, and in Ref. [4] terms proportional to the pion mass have been set to zero. x is defined as

B. Total decay widths
The expressions for the total decay width in terms of the chiral coefficients a i can be found in Eqs. (A10) and (A11). 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 In our EFT the total decay width is proportional to (β (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 Λ 1 , the lowest gluelump mass (for further details see Appendix B). For the lowest hybrid we compute β Introducing the value of (β (12) r,Q ) 2 , obtained as described above, and κ = 0.301, in Eq. (29)-Eq. (31), the expressions for the total width in Eq. (A10) and Eq. (A11) read (the experimental values are taken from [30]) The 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 (β (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 Γ ψ(2S)→J/ψπ + π − /Γ Υ(2S)→Υ(1S)π + π − .
We will discuss this ratio later in Sec. VI. Other corrections to Eqs. (D14) and (D15) are due to 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.QQ(2S) →QQ(1P )π HADRONIC TRANSITIONS TheQQ(2S) →QQ(1P )π hadronic transition 3 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 oper- 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. (46) in terms of the spins of the quark and antiquark can be found in Eq. (105) of [36]. 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 [37,38]. 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 [30].
Adding the operator in Eq. (46) to the leading pNRQCD Lagrangian in Eq. (1) and considering isospin violation effects produces the following new terms in the hadronic Lagrangian at LO 4 : Terms that break spin and isospin symmetry simultaneously are not considered, since they produce subleading contributions to the transition we are considering. 3 Let us note that this type of transition was first considered in Ref. [35] for the Υ(2S) → h b (1P )π 0 decay.
However, this particular case turned out not to be kinematically allowed. 4 There is another possible pseudoscalar operator of the same order as χ − : D µ u µ . However they are both related through the leading order equations of motion D µ u µ = i (χ − − χ − /2) /2.
Let us match the parton and hadronic description of the spin-dependent mixing operator.
In the hadronic EFT the correlator reads and in the partonic version of pNRQCD we have where amp. signals that only amputated contributions are considered. Now we consider G a 1 +− to be given by a sum of all possible gluonic operators with the same quantum numbers: 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. (49), and matching to Eq. (48), we arrive at t (S1 Let us now match the operators with an odd number of pions in the Lagrangian in Eq. (47). In the hadronic theory we have and in pNRQCD: 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 Eq. (52) and Eq. (53) we arrive at We are interested in investigating the decay ψ(2S) → h c (1P )π 0 , which is given in the hadronic pNRQCD EFT by the diagrams in Fig. 5. The matrix element reads with In the first term the ψ(2S) mixes into a hybrid state with J P C = 1 −− but in a spin singlet, then the hybrid decays into a h c conserving the spin but changing the angular momentum.
The error analysis has been performed similarly to the previous section. In this case we do not have error associated to κ, but still have a dependence on Λ 1 (the lowest lying gluelump mass), and s.p. (the different parameterization for the singlet and hybrid static potentials).
We also estimate the error due to the value of the light-quark mass ratio (l.q.), which we take as In principle, forQQ(2S) →QQ(1P )π transitions, the energy release is small so one can use the twist expansion with no fear if one assumes that Λ QCD 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 Γ ψ(2S)→hc(1P )π 0 ∼ 15 eV, from Ref. [27], using the twist expansion (Λ QCD mv 2 and E mv 2 ).

V.QQ(2S) →QQ(1S)π 0 (η) HADRONIC TRANSITIONS
We now turn our attention toQQ(2S) →QQ(1S)π 0 (η) 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 which has to be added to Eq. (1). In the hadronic EFT these operators correspond to To determine the Wilson coefficients we compute the transition amplitude both in the hadronic and partonic version of the effective theory. Nevertheless, there is one extra subtlety to be taken into account. The π 0 and η fields in Eq. (11) mix, and to obtain the physical states the mass matrix needs to be diagonalized. The physical states correspond to π 0 phys. = π 0 + η , with the mixing angle .
To compute the decay amplitudes the physical states must be considered. In the transitions with η emission the contribution due to the mixing is subleading, but in the case of π 0 emission both contributions are of the same order. In practice this amounts to an extra factor 3/2 in the π 0 decay amplitude. In the hadronic EFT we have and in the partonic pNRQCD, 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 O(α s ) corrections).
Matching Eq. (64) and Eq. (65) we obtain Therefore, the intermediate hybrid states for the transitions we are considering must be 3 S 1 .
The amplitude for a pion emission then reads Note that the factor involving the sum over intermediate hybrid states: β (12) r,Q , given in Eq. (28), is the same to the one appearing in the two-pion transitions. mS stands for the polarization of a m 3 S 1 quarkonium state. The decays amplitude to one η is  state. Indeed the dependence is the same as the one we had in Sec. III B. 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 QQ(2S) → QQ(1S)π 0 (η) 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 A 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 QQ(2S) → QQ(1S)π 0 (η) transitions do not depend on β (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 QQ(2S) → QQ(1S)ππ transitions [14]), as studied previously in Refs. [39][40][41]. 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 η /F π and the mixing angle . Generically, we expect these to be of O((m π , m K )/Λ χ ) ∼ 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) 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 β The differences with experiment are large, of the order of 70%-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 O(α s ) corrections in the hadronization of the local operators (this latter source of error only applies to the QQ(2S) → QQ(1S)ππ transitions). In this respect, it is interesting to note that we expect the O(α 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 get agreement with experiment (within errors) both for charmonium and botttomoniun using the same correcting factor: ∼ 1/3. This factor could be understood as evaluating the β(α s ) function at a low scale: − 2πβ(αs) Still, it is premature to draw any conclusion out of this. Note also that these estimates are significantly affected by the uncertainties of κ, and, for the π 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. (75)-(80), the prediction of the twist expansion is the same to the one we have found here. This is not so for the following ratio: R exp bc,ππ =5.59(0.50) × 10 −2 .
In principle, we expect that for this ratio most of the neglected O(α s ) corrections in the hadronization of the local operators vanish. This observable can be considered a rough measure of β (21) r,b /β (21) r,c . The agreement with experiment is remarkable: below 20%, and could be accounted for by the quoted errors.
We can also consider the following ratios: Again, the twist expansion would yield a different prediction for these ratios, and, again, these ratios can be considered a measure of β (21) r,b /β (21) r,c . The agreement with experiment is quite reasonable, with difference of the order of 30%. In order to asses whether this error and the error of Eq. (81) is really related to β (21) r,b /β (21) r,c or to something else, one may consider the double ratios: These double ratios are independent of β (21) r,b /β (21) r,c . 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 O(α 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 κ nearly vanishes. Nevertheless, this should be studied with more detail to make this statement more quantitative. For instance, for the QQ(2S) → QQ(1S)η 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 Λ QCD .
On the other hand our analysis does not need to constrain the relative size of Λ QCD with respect to mv 2 . For definiteness, we have generically considered the situation mv 2 Λ 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 weaklycoupled pNRQCD the hybrid states consist of a color-octet heavy quark-antiquark field and a gluonic excitation operator, called gluelump, characterized by a J P C . 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 P C . 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 QQ(2S) → QQ(1S)ππ 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: κ. Our best fit to the combined charmonium and bottomonium data yields the value quoted in Eq. (41) for κ.
Our computational scheme produces the same theoretical expression as the twist expansion [7] for the normalized decay width spectra computed in Sec. III A. Note, however, that the twist expansion requires that Λ QCD mv 2 , and that the energy release E mv 2 , the latter condition is never fulfilled for two-pion transitions. Lower values for κ than the ones in Eq. (41) are found in [4,34]. The main source for the discrepancy is due to O(m 2 π ) terms not included in these analysis.
Our prediction for the total two-pion transition decay width depends on κ (which we take from Eq. (41)), and on β (12) r,Q . 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 Γ ψ(2S)→J/ψπ + π − , and Γ Υ(2S)→Υ(1S)π + π − can be found in Eqs. (44) and (45).
Enlarging our Lagrangian to include spin-dependent and isospin breaking operators we can use the same procedure to compute the width for QQ(2S) → QQ(1P )π 0 transitions. Our prediction for the total width Γ ψ(2S)→hc(1P )π 0 can be found in Eq. (57). 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 O(α) corrections, nor there are O(p 4 ) chiral corrections. We also observe a very weak dependence on variations of the wave function of the hybrid and singlet state: This object depend on β (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.
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, specially for the QQ(2S) → QQ(1S)π 0 (η) 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 Eqs. (75) and (76) 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. (77)-(80). 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 accomodated 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. (75) and (76) are independent of the bound state dynamics, i.e. they are independent of β (12) r,Q 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. (81)-(83). 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: β (12) r,b /β (12) r,c . For Eq. (81), good agreement is found (below 20%), and quite reasonable for the other two ratios: of the order of 30%. This may indicate that our evaluation of β (12) r,b /β (12) r,c is quite reasonable. Finally, we have considered double ratios in Eq. (84). In principle, these are the cleanest objects to compute (leaving aside Eqs. (75) and (76)). 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 analyses, in particular of the ratios and, specially, of the double ratio discussed in Eq. (84).
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.
where m 2 ππ = (p + + p − ) 2 is the dipion invariant mass, and p + , p − refers to the momentum of the π + and π − respectively. In the reference frame of the decaying quarkonia we find 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 Using Eq. (A2) in Eq. (A1) we can write The differential decay width is and after integrating θ it reads To obtain the total decay width we integrate numerically In the one pion transition the momenta are fixed by momentum conservation, in particular the final pion momentum is |p π | = ρ(m π ), and the decay width is given by with ρ(m π ) is given in Eq. (A5) replacing m ππ by m π .

Appendix B: Singlet and hybrid bound states
The heavy quarkonium wave function reads where φ (n) (r) is the solution of 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 (1GeV) = 1.496 GeV and m b (1GeV) = 4.885 GeV [43] for the heavy quark masses in the RS scheme [44].
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 Σ + g obtained by fitting to the lattice data from Ref. [42] up to r = 1 fm together with a free energy offset, b offset Σ + 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 For the perturbative static singlet potential we take the O(α 3 s ) result 6 evaluated at ν = 5.6 GeV (the shortest available scale). The parameters obtained from the fits are the following In Fig. 6 we plot the lattice data together with the fitted Eq. (15) (with the offset) and Eq. (B3) potentials as dashed and solid lines respectively. The spectra obtained are displayed in tables III and IV for the potential in Eq. (15) and in Eq. (B3) respectively. The origin of energies is chosen in order to reproduce the experimental spin average of the ground state. 6 The O(α s ) term was computed in [45], the O(α 2 s ) in [46,47], the O(α 3 s ) logarithmic term in [48], the light-flavor finite piece in [49], and the pure gluonic finite piece in [50,51].  Table III. Quarkonium spectrum obtained from the potential in Eq. (B3) fitted to the lattice data Ref. [42]. The origin of energies is chosen in order to reproduce the spin average of the ground state. All dimension-full entries are in GeV.
The hybrid wave functions are obtained following the procedure described in Ref. [18], by solving the coupled Schrödinger equations involving the potentials V  Table IV. Quarkonium spectrum obtained from the potential in Eq. (15) fitted to the lattice data Ref. [42]. The origin of energies is chosen in order to reproduce the spin average of the ground state. All dimension-full entries are in GeV.
The parity and charge conjugation of these states corresponding to the Ψ ± solutions are 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 which for λ = 0 corresponds to the usual differential equation for spherical harmonics and therefore v 0 l m l = Y l m l . Note that we have chosen to represent the angular momentum quantum numbers of the hybrids and quarkonia by 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, . . . for l = 0, 1, 2, . . . and for hybrid states we use S, P, D, . . . for = 0, 1, 2, · · · .
The λ 2s+1 j spin-angular hybrids wave functions can be shown to be whereê are the polarization vectorsê and the tensors h ij 2m J are traceless, completely symmetric and normalized as The 2s+1 l j spin-angular wave functions for quarkonium correspond to the λ = 0, 2s+1 j hybrids wave functions.
The radial wave function in Eq. (B6) are obtained by (numerically) solving the following coupled radial Schrödinger equations, For the special case = 0 the equations decouple and in fact ψ The hybrid static potentials, V Πu , V Σ − u , V Πg and V Σ + g , match to the NRQCD heavy quarkantiquark static energies, which have been computed on the lattice [52][53][54] 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λ obtained by fitting to the lattice data from Ref. [52] 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 ≤ 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 The perturbative static octet potential is taken at O(α 3 s ) from Ref. [55], and evaluated at ν = 2.6 GeV (the shortest available scale). The results for the hybrid spectrum for V Πu -V Σ − u associated to the 1 +− gluelump were discussed at length in Ref. [18] 7 . The parameters obtained from the fits are the following Πg ≡ c The fits of the potentials V Πg -V Σ + g , associated to the 1 −− gluelump, are shown in Fig. 7 in solid and dashed lines and the corresponding spectra in tables V and VI for the potentials in Eq. (B21) and Eq. (16) respectively. The origin of energies is chosen, as in Ref. [18], such that the short distance limit of V Πu and V Σ − u tend to V 0 o + Λ 1 +− with the value of Λ 1 +− = 0.87 (15) GeV computed in Ref. [56].
We remark that, compared with the prediction of the previous section, these results do not have neither O(α s ) nor O(p 4 ) corrections.