\chi_{cJ}\rightarrow e^{+}e^{-} decays revisited

# χcJ→e+e− decays revisited

N. Kivel  and M. Vanderhaeghen
Helmholtz Institut Mainz, Johannes Gutenberg-Universität, D-55099 Mainz, Germany
Institut für Kernphysik, Johannes Gutenberg-Universität, D-55099 Mainz, Germany
On leave of absence from St. Petersburg Nuclear Physics Institute, 188350, Gatchina, Russia
###### Abstract

We present a calculation of the width for decay. The amplitude of the process is computed within the NRQCD framework. The leading-order contribution is described by two terms associated with the two different integration domains in the electromagnetic loop describing two-photon annihilation of the heavy quark-antiquark pair. The corresponding operators are defined in the framework of NRQCD. The matrix element of one of these operators describes a configuration with an ultrasoft photon and can be associated with the higher Fock state of the heavy meson. In order to compute this contribution we use the heavy hadron chiral perturbation theory. We obtain that this contribution is numerically dominant. The obtained estimates for the decay widths of the and states are eV and eV, respectively.

## 1 Introduction

The leptonic decays of C-even charmonium states into a lepton pair have a very small branching ratio because they can only occur via a two-photon intermediate state . However with the high-luminosity BEPC-II -accelerator operating on the charmonium energy region such measurements of direct production cross section become feasible. The study of the mechanism of the production of -even quarkonium states is especially interesting in view of the production of higher resonances such as the exotic charmonium-like state X(3872).

The decays have already been studied long time ago in Ref.[1]. The authors considered the quarkonium description and a vector dominance model (VDM) in order to describe the decay amplitudes of charmonium states with . It was found that a naive quarkonium description is problematic because of infrared logarithmic divergencies arising in the integrals describing the quark-photon annihilation loop. Such divergencies indicate that the corresponding contribution is also sensitive to long distance physics. The corresponding integrals in Ref.[1] have been regularized by introducing the binding energy MeV. The numerical estimates obtained in this way give very small value of the widths, smaller thÐ°n bounds derived from analyticity and unitarity in the same work. Probably, a more realistic estimate was obtained using generalized VDM which yields a larger numerical value for the widths, consistent with the unitarity constrains:  eV and  eV [1].

Recently the decay rate of the state was again estimated in [2] using the VDM approach, with result  eV. The short distance contributions describing a configuration with highly virtual photons were considered in this framework as unknown contact vertices giving rise to a theoretical uncertainty.

The aim of the present work is to provide a more systematic description of the decay amplitudes for process using the NRQCD factorization framework [4, 5], see also review [6] and references theirin. This technique allows one to perform a systematic description of heavy quarkonium states using the small relative velocity of heavy quarks and the small QCD running coupling at short distances. Within this framework we associate the IR-divergencies found in [1] with a specific quark-photon operator which describes a configuration with an ultrasoft photon. The matrix element of this operator describes the overlapping with the higher Fock state which consists of heavy quark-antiquark and photon. This allows us to perform a systematic separation and description of the different contributions relevant in a leading-order expansion in small velocity .

Our paper is organized as follows. In Sec. 2 we briefly describe our notation and provide definitions of various quantities used in the following. Sec. 3 is devoted to the investigation of the one-loop integral which describes the leading-order contribution. In this section we establish the dominant regions and provide a description of the amplitude within the NRQCD factorization framework. Sec. 4 is devoted to the calculation of the ultrasoft photon matrix element in the heavy hadron chiral effective theory (HHPT). Furthermore, a estimate of the decay rates is given. In Sec. 5 we briefly summarize our results.

## 2 Kinematics and notation

Let us start from the description of the decay kinematics . The initial state momentum can be written as

 P=Mχ ω,  ω2=1, (1)

where denotes the charmonium velocity. In the following, we consider the charmonium rest frame which implies

 ω=(1,→0). (2)

The small relative velocity of heavy quarks in the bound state is denoted as . Neglecting lepton masses, the lepton momenta can be written as

 l1=Mχn2,  l2=Mχ¯n2, (3)

where and denote the light-like vectors which satisfy . Any 4-vector can be expanded as

 (4)

where denotes the components transverse to the light-like vectors : . Similarly, one can also write a decomposition

 Vμ=(V⋅ω)ωμ+Vμ⊤, (5)

where denotes the component which is orthogonal to the velocity : . In the following we assume that in the rest frame

 ω=n2+¯n2. (6)

The momenta of the heavy quark and antiquark with mass which form a quarkonium state can be written as

 p1=12P+Δ,  p2=12P−Δ, (7)

where the relative momentum satisfies

 (Δ⋅ω)=0,  Δ2=−→Δ2. (8)

The heavy quarks which create a bound state are non-relativistic, implying that the relative velocity is quite small: .

The power counting rules of NRQCD has been established in [3, 4]. Following these arguments we assume that the mass is large enough and that the most important scales such as mass , typical three-momentum of the heavy quark and its typical kinetic energy satisfy

 (mv2)2≪(mv)2≪m2. (9)

Integrating out the modes with hard momenta one passes onto the effective theory NRQCD which describes the modes with the soft momenta . If the scale then one can integrate over the soft region together with potential gluons with momenta [7, 8, 5, 9, 10]

 \ p0∼mv2,  →p∼m→v, (10)

After this one obtains a new effective theory which is known as potential NRQCD (pNRQCD). For a more detailed information about these effective theories see Ref.[6] and references therein.

The charm quark mass GeV is not large enough compared to therefore in this case one can only factorize the effects at momentum scales of order . However, in the QED sector one can also consider the possibility to integrate over the soft region too. As we will show further on, such situation is relevant for the factorization of the electromagnetic loop describing the transition .

After factorization of the hard contribution, the nonpertubative QCD dynamics is described by the matrix elements of appropriate operators defined in an effective theory. In the following, we will need the following set of NRQCD operators which describe the matrix elements between the charmonium and the vacuum states:

 Oσ(3S1)=χ†ωγσ⊤ψω, (11)
 O(3P0)=−1√3 χ†ω(−i2)↔Dα⊤γα⊤ψω, (12)
 Oβ(3P1)=12√2 χ†ω↔Dα⊤(−i2)[γα⊤,γβ⊤]γ5ψω, (13)
 Oαβ(3P2)=χ†ω(−i2)↔D(α⊤γβ)⊤ψω, (14)

where the covariant derivative ,  . Furthermore, we use the covariant four-component fields to describe the soft quark and antiquarks within the NRQCD framework. These fields satisfy

 (15)

Using Eq.(2), one can show that the operators in Eqs.(11)-(14) can be reduced to the set of well-known non-relativistic operators constructed from two-component Pauli spinors.

The matrix elements of these operators are well known in the literature and can be written as

 (16) ⟨0|ψ†ωγσ⊤χω∣∣ψ′⟩=ϵσψ′ ⟨O′(3S1)⟩, (17)
 (18)
 (19)
 (20)

The constants and are related to the value of the charmonium wave functions at the origin

 ⟨O(3S1)⟩=√2Nc√2Mψ√14π R10(0), (21) ⟨O′(3S1)⟩=√2Nc√2Mψ′√14π R20(0), (22)
 ⟨O(3P0)⟩=√2Nc√2Mχc0√34πR′21(0), (23)

where is the radial part of the Schrödinger wave function and denotes its derivative. The rhs of Eqs.(18)-20) depends on of the same constant due to the spin symmetry of the leading non-relativistic action [4].  The polarization vectors , and correspond to spin-1 and spin-2 charmonium states, respectively. They are normalized to satisfy

 ∑λϵσX(λ){ϵρX(λ)}∗=−gσρ+PσPρM2X, (24)

with and

 ∑λϵαβχ(λ){ϵα′β′χ(λ)}∗=12Mαα′Mββ′+12Mαβ′Mβα′−13MαβMα′β′, (25)

with .

The decay amplitudes  are  defined as

 ⟨e+e−;out∣∣in; χcJ⟩=i(2π)4δ(l1+l2−P) AJ, (26)

with

 AJ=¯unΓJv¯n TJ, (27)

and where and denotes the spinors of the massles lepton and antilepton, respectively

 (28)

and

 Γ1=ϵσχγ⊥σγ5, Γ2=ϵσρχnργ⊥σ. (29)

The leading-order contribution to these amplitudes arises from the annihilation of heavy quarks into two photons which create the outgoing lepton pair, see Fig. 1. If the one-loop integral is dominated by the hard region where both photons and heavy quark are highly off-shell then one can expect that such process can be described within the NRQCD approach and the amplitude can be factorized into hard and soft parts. In the next section we consider this possibility in more detail.

## 3 Factorization of decay amplitudes in NRQCD

The leading-order in diagrams describing the decay of -even charmonia are shown in Fig. 1. These one-loop diagrams are constructed from the photon, lepton and heavy quark (double lines).

The diagrams in Fig. 1 can be computed in the heavy quark mass limit , performing an expansion in the small parameter . Let us to start from a naive guess that the dominant contribution is only provided by the hard region where the loop momentum , and therefore all propagators are far of off-shell. The leading-order contribution in is provided by projections onto the leading-order operators described in Eqs.(18)-(20). The technical details are well known in the literature, see e.g. [1]. The resulting expressions can be presented as

 (30)
 (31)

where the square brackets for the propagators denote the standard Feynman prescription . The corresponding contribution to the amplitude vanishes and therefore is suppressed by a power of and will not be considered it in this work.

The total structure of the integrands in expressions (30) and (31) can be divided into the lepton and heavy quark parts. The lepton part has in the numerator and includes the photon and lepton propagators in the denominators. The heavy quark part is given by Tr. We introduced the projections onto charmonium states

 (32)
 (33)

 Γαβμ(k)=12m{γμ^DQ(k)+^DQ(k)γμ}+^D′μQ(k), (34)

with

 (35)
 ^D′μQ =i(ieeQ)2[k2−2m(kω)−→Δ2]{γβγμγα+γαγμγβ} (36)

where is the charge of the heavy quark (). The small squared relative momentum which appears in the heavy quark propagator provides an IR-regularization and can be neglected if it is not required. With this regularization the traces and loop integrals are computed in four dimensions with .

The expressions (30) and (31) have been obtained by expanding the heavy quark fields in position space

 (37)

and projecting the soft quark fields and onto leading-order operators (12)-(14). The terms (arising from the multipole expansion of the soft quark field arguments ) lead to the expansion of the integrand with respect to small relative momentum giving the contribution . The terms proportional to give the contribution with . The evaluation of the integrals in Eqs. (30) and (31) gives

 A1=¯unΓ1vn i⟨O(3P0)⟩α2m3e2Q 2√2lnm22→Δ2, (38)
 A2=¯unΓ2vn i⟨O(3P0)⟩ α2m3e2Q2(2ln→Δ2/m2+23(ln2−1+iπ)). (39)

These expressions are in agreement with the results obtained in Ref. [1]. We obtain that both amplitudes depend on the  large logarithm which is sensitive to the soft scale . This shows that the starting assumption about one dominant region is incorrect. There must be at least one more domain where some propagators in the loop integral are soft. One can expect that the additional region is associated with the configuration when one of the photons is soft. In this case the propagator of the heavy quark is also soft and the hard configuration is described by the tree level subdiagram describing the annihilation through one photon.

In order to get an idea about the explicit definition of this region it is useful to investigate the integrals of diagrams in Fig. 1 within the threshold expansion technique worked out in Ref.[5]. According to this analysis the threshold kinematics is described by the following regions

 hard: kμ∼m, (40)
 soft:kμ∼mv, (41)
 potential:k0∼mv2,  →k ∼mv, (42)
 usoft: kμ∼mv2. (43)

The same regions can also be considered for the photon with momentum . These regions can be associated with the fields appearing in the effective Lagrangians, see e.g. Ref.[6].

According to the threshold expansion prescription an integrand is expanded in each domain to a required accuracy and the resulting integral is computed in dimensional regularization. A detailed analysis of the full expressions in Eqs. (30) and (31) is quite similar. To be definite let us consider the integral which enters in Eq.(30)

 J=∫dk e2 ¯unγα(l1−k)γβv¯n[(k−l1)2][k2][(k−P)2]ϵνχ4Tr[P1μνΓαβμ(k)]. (44)

Keeping the denominators of the heavy quark propagators in unexpanded

 (p1−k)2−m2=(12P+Δ−k)2−m2=k2−P0k0+2(→k⋅→Δ)−→Δ2+14P20−m2, (45)
 (p2−k)2−m2=(12P−Δ−k)2−m2=k2−P0k0−2(→k⋅→Δ)−→Δ2+14P20−m2, (46)

where .  In the hard region, the small scalar products with and the term can be neglected resulting in

 [(12P±Δ−k)2−m2]h≃k2−(kP), (47)

which appear in the expressions (35) and (36) (up to small regularization term ). From dimensional counting one immediately finds

 Jh∼¯unΓvnm3, (48)

where denotes the Dirac structure. Computing the hard integral in dimensional regularization one finds the IR poles . These singularities must cancel in the sum with other contribution.

Expanding the integrand (44) in the soft region (41) yields

 (49)

where is given by (34) with

 (50)

Calculating the trace and performing the contractions in the numerator results in

 Js∼1m3¯unΓ1v¯n∫dk 1 [k2][−(kω)]2∼¯unΓ1v¯nm3. (51)

As the integral in (51) is scaleless it therefore vanishes in the dimensional regularization, i.e. .

In the potential region (42), the expansion of the heavy quark propagator reads

 [(12P±Δ−k)2−m2]p≃P20/4−m2−P0k0−(→k±→Δ)2. (52)

The computation of the corresponding integral then yields

 Jp ≃1m¯unΓ1v¯n∫dk  1[−→k2][P20/4−m2−P0k0−(→k+→Δ)2]2+(→Δ→−→Δ)∼¯unΓ1v¯nv−1m3. (53)

However the poles in in the integrand of Eq.(53) lie in the same imaginary half-plane and therefore the integral over vanishes. This observation is also true for the higher order contributions in appearing from this domain. We can therefore conclude that the potential region cannot contribute in this case.

In the ultasoft domain (43), the heavy quark propagators are expanded as

 [(k−12P±Δ)2−m2]us≃P20/4−m2−P0k0−→Δ2. (54)

Performing the expansion of the integrand one gets

 Jus∼1m¯unΓ1v¯n∫dk  1[k2][P20/4−m2−P0k0−→Δ2]2∼¯unΓ1v¯n1m3. (55)

This integral has the same scaling behavior as the hard integral in Eq.(48). One can also see that the integral in Eq.(55) is UV divergent. The similar analysis can also be carried out for the second photon with momentum . Therefore we conclude that the exact integral must be given by sum

 J=Jh+Jus, (56)

where denotes the contributions from the both ultrasoft domains. This conclusion can be checked by explicit calculations. A similar conclusion for the two-photon diagrams in Fig.1 has also been obtained in Ref.[11].

Guided by this consideration we suggest that the additional relevant domain is described by the ultrasoft region. In order to find the description of the appropriate operator in the effective theory one has to integrate out hard and soft photons and leptons. After that the description of QED sector includes only collinear leptons and ultrasoft photons. The integration of the soft photon with the lepton and quark must be described in the framework of the effective theory.

Within the above picture the factorization of the decay amplitudes can be described as a sum of two contributions

 AJ = ¯unΓJv¯n C(J)γγ i⟨O(3P0)⟩ +Cγ⟨e+e−∣∣¯ξn(0)Y†n(0)γσ⊥Y¯n(0)ξ¯n(0) Oσ(3S1) |χcJ⟩. (57)

The first term on rhs of this equation corresponds to the hard domain with the hard photons, denotes the corresponding hard coefficient function.

The second term on rhs of Eq.(57) corresponds to the domain with the ultrasoft photon. The operator is defined in Eq.(11). The outgoing collinear leptons are described by fields and which defined as

 (58)

The photon Wilson lines and describe the interaction of the ultrasoft longitudinal photons with the energetic lepton and antilepton and read

 Y†n(0)=Pexp{ie∫∞0ds n⋅Bus(sn)},   Y¯n(0)=\={P}exp{−ie∫∞0ds ¯n⋅Bus(s¯n)}, (59)

where denotes the ultrasoft photon field. The appearance of these Wilson lines is related with the fact that in a general gauge the tree level diagram with attachments of photon to the collinear field describing the outgoing lepton111We assume electrical charge is measured in proton units (positron is particle and electron is antiparticle) that allows to use the same notation for the covariant derivative and Wilson lines as in QCD. are resummed to the P-ordered exponents

 (60)

The leading-order hard coefficient function is defined by the diagram in Fig.2 and reads

 Cγ=απm2eQ. (61)

The soft and collinear modes in the effective action describing the QED sector are decoupled. This property is well known in the soft-collinear effective theory, see e.g. Refs.[12, 13, 14]. This allows us to contract the lepton fields in the second matrix element in Eq.(57) and rewrite it as

 ⟨e+e−∣∣¯ξn(0)Y†n(0)γσ⊥Y¯n(0)ξ¯n(0) Oσ(3S1) |χcJ⟩=¯unγσ⊥v¯n ⟨0|Oσγ(3S1)|χcJ⟩, (62)

with

 Oσγ(3S1)≡Y†n(0)Y¯n(0)Oσ(3S1). (63)

The presence of the soft scale in Eqs. (38) and (39) can be explained by the contribution with ultrasoft photon. Therefore in order to find the hard coefficient functions we have to perform the matching onto the configuration described by Eq.(57). For that purpose we need to compute the ultrasoft matrix element (62) in the effective theory.

The interaction of ultrasoft photons with quarks are described within the pNRQED. The ultrasoft photons have momentum so that photon field scales as

 Busμ∼mv2. (64)

The scaling of the quark fields reads

 ψω∼(mv)3/2,  →∂iψω∼(mv)ψω,  ∂0ψω∼(mv2)ψω. (65)

Using this counting one finds

 C(J)γγ Oσ(3PJ)∼m−3(mv)4. (66)

At the same time

 CγOσγ(3S1)∼CγOσ(3S1)∼m−2(mv)3. (67)

However the pure quark operator is -odd and therefore it cannot contribute to the matrix element with a -even charmonium state

 (68)

In order to obtain a nontrivial contribution one needs to consider at least one interaction of an ultrasoft photon with the quark in pNRQED. We only need the two-particle sector describing the electromagnetic interactions of quarks (in rest frame )

 Lem0[Bus]=∫d4x ψ†ω(x)γ0(iω⋅∂+i∂⊤⋅i∂⊤2m)ψω(x), (69)
 Lem1[Bus]=∫d4x ψ†ω(x)γ0[ →x⋅∂⊤ eeQ ω⋅Bus(x0)+1meeQBus(x0)⋅i∂⊤] ψω(x), (70)

and analogous contributions with antiquark fields. The arguments of the ultrasoft photon field are expanded because the space components of the quark fields varies at , the measure scales as . With these rules one finds that and . The leading-order term (69) provides the soft quark propagator

 Δω(k)=i(ωk)−→k2/2m+iε. (71)

A nontrivial contribution to the matrix element   can be obtained from -product

 T{Oσγ(3S1),Lem1[Bus]}∼mv4, (72)

which is of the same order as the hard contribution in Eq.(66). Calculation of this -product gives diagrams shown in Fig.3. The dashed lines can be associated with the collinear leptons or equivalently with the ultrasoft Wilson lines (59).

These diagrams induce a mixing of the operators and due to electromagnetic interaction in the framework of pNRQED.

In order to perform the matching onto operators according to formula (57) one has also to compute the contribution of the diagrams in Fig.3. The simplest way to proceed is to follow the same technique as we used above for diagrams in Fig.1.

Let us consider as initial state. In this calculation we set and only keep the relative momentum . Then the sum of all four diagrams gives

 ⟨e+e−∣∣CγT{Oσγ(3S1),Lem1[Bus]}|χc2⟩=¯unΓ2v¯n i⟨O(3P0)⟩ 8Cγe2eQ2m Jus. (73)

Computing these diagrams we project the soft quarks fields on the operator and substitute the corresponding matrix element which gives the factor , the coefficient arises from the sum of the four diagrams shown in Fig.3, the ultrasoft loop integral reads

 Jus=(−i)∫dk 1[k2]1[−(ωk)−→Δ2/2m]2. (74)

This integral coincides with the ultrasoft integral of Eq.(55) obtained within the threshold expansion approach up to term which vanishes because we set . The integral in Eq.(74) is UV-divergent and we use dimension regularization in order to compute it. The result reads

 (75)

where is the factorization scale. The pole is the UV-pole which describes UV-mixing of the operators and , schematically

 [Oγ(3S1)]R=Oγ(3S1)+ZJ O(3PJ), (76)

where on the lhs of Eq.(76) denotes the renormalized operator. Furthermore, is the corresponding renormalization constant. Assuming -subtraction scheme one finds

 [Jus]R=14π2ln→Δ2mμF. (77)

Hence we obtain

 ⟨e+e−∣∣CγT{Oσγ(3S1),Lem1[Bus]}|χc2⟩R=¯unΓ2v¯n i⟨O(3P0)⟩ α2m3e2Q 4ln→Δ2mμF. (78)

The soft matrix element for the can be computed in the same way, resulting in

 ⟨e+e−∣∣CγT{Oσγ(3S1),Lem1[Bus]}|χc1⟩R=¯unΓ1v¯n i⟨O(3P0)⟩ α2m3e2Q 2√2lnmμF→Δ2. (79)

The hard coefficients are given by

 (80)

where the expressions for are given by Eqs.(38) and (39). The important check of the factorization formula (57) is the cancellation of the ultrasoft scale in the expressions for obtained from Eq.(80). Substituting the computed expressions in Eq.(80) we obtain

 C(1)γγ=α2m3e2Q √2lnm24μ2F, (81)
 C(2)γγ= α2m3e2Q2{lnμ2Fm2+23(ln2−1+iπ)}. (82)

These expressions are the main result of this section. We observe that the soft scale cancel in Eqs.(81) and (82) as it must be. Hence the factorization formula described by Eq.(57) describes properly the ultrasoft region of the one-loop diagram.

The coefficient function has an imaginary part which originates from the two-photon cut. Such mechanism can not work for state, therefore is real.

The hard coefficient functions depend on the factorization scale . Therefore

 AJ=¯unΓJv¯n C(J)γγ(μF) i⟨O(3P0)⟩+Cγ ¯unγσv¯n⟨0|Oσγ(3S1)|χcJ⟩(μF), (83)

and the independence of the amplitude on yields the evolution equation

 (84)

The solution of this equation depends on the initial condition defined at some scale . Performing numerical estimates one has to fix a value of this scale. By derivation this scale separates the hard region (two hard photons) from the ultrasoft region (hard and ultrasoft photons). Therefore it is natural to associate this scale with the virtuality of the ultrasoft photon and to set to be of order MeV. Then the matrix element of the operator on the rhs of Eq.(83) describes only the ultrasoft nonperturbative contribution which can be only estimated within some low-energy effective theory or model. Similar to the well known color octet mechanism, the operator can also be associated with the electromagnetic mechanism. The corresponding matrix element can be interpreted as an overlap with the higher Fock state which includes a dynamical photon while the matrix elements of the operators describe the coupling to the dominant quark-antiquark state. Therefore the full description of the leptonic decay requires a knowledge on the subleading structure of the quarkonium state.

In the large mass limit one can consider a specific situation known as the Coulomb limit when the binding energy is larger then the typical hadronic scale . In this case the strong coupling is quite small and ultrasoft contribution can be estimated within the pNRQCD. Then one has to compute the diagram as in Fig.4 resumming the interactions with Coulomb gluons.

Such calculation has been carried out for the radiation function in Ref.[15]. Perhaps, such calculation might also be interesting here in order to get an idea about the relative value of this matrix element in the Coulomb limit. In present paper we will obtain an estimate of the ultrasoft matrix element using the so-called heavy hadron chiral perturbation theory (HHPT) framework in the next section.

## 4 Phenomenology

### 4.1 Calculation of the ultrasoft matrix element in the heavy hadron chiral perturbation theory

In order to provide a numerical estimate of the decay rate we need to estimate the ultrasoft matrix elements

 ⟨0|Oσγ(3S1)|χcJ⟩ (85)

which describes an overlap with the higher Fock component of the charmonium state in which a dynamical photon is present. One can expect that the soft photon has already quite large wavelength and therefore it interacts with heavy charmonium as with a point-like source. Then it is natural to expect that the relevant dynamical degrees of freedom in this case are associated with mesonic fields and the corresponding low energy dynamics is described by the most generic effective action compatible with the symmetries of NRQCD. Such an approach is known as heavy hadron chiral perturbation theory in Refs.[16, 17] for the heavy-light mesons and then generalized on quarkonia in Refs.[18, 19, 20]. This framework can also be used for the calculation of the matrix element in Eq.(85).

For our purpose we need only the electromagnetic sector of the HHPT described by the effective action which includes the kinetic terms for and states and the vertices describing the electromagnetic vertices and 222We are grateful to Maxim Polyakov for discussion of the contribution with the virtual state . . As before we assume the rest frame for the initial state . The kinetic Lagrangian reads

 Lkin(x)=122Mχ ψ(ω)μ(x){i(ω∂)−ΔM}ψ(ω)μ(x)+122Mχ ψ′(ω)μ(x){i(ω∂)−Δ′M}ψ′(ω)μ(x), (86)

with the residual masses and . The fields and describes the residual motion of the heavy and particles and satisfy .

 LemSP=12eeQfγ Tr[γ0~{}J†Sγ0JμP]Fμνων+LemSP+12eeQf′γ Tr[γ0~{}J′†Sγ0JμP]Fμνων+h.c. (87)

with

 (88) (89)

and

 (90)

The currents and describe particles from - and -wave multiplets, respectively. In Eq.(87) we introduced the dimensionless couplings and . The fields describe charmonium states . Computing the trace in Eq.(87) one finds

 LemSP=−eeQfγ χμα2ψ(ω)αFμνων−eeQfγ √2iεμαβρψ(ω)αχ1βωρFμνων (91) −eeQ f′γ χμα2ψ′(ω)αFμνων−eeQf′γ √2iεμαβρψ′(ω)αχ1βωρFμνων+ ... (92)

where we show only the relevant terms.

Our calculations involve operators and which have also to be matched onto physical quarkonium fields. The spin symmetry in the heavy quark limit yields

 (93)
 [Oμ(2s+1PJ)]αβ=⟨O(3P0)⟩[Jμ]αβ, (94)

where are spinor indices. Taking the matrix element and computing the traces one can see that Eqs.(93) and (94) reproduce correctly the matrix elements (16)-(20). Using these results one finds

 Oσγ(3S1)≃{⟨O(3S1)⟩ψ(ω)σ(0)+⟨O′(3S1)⟩ψ′