New factorization theory for heavy quarkonium production and decay

TL;DR

This paper introduces a new soft gluon factorization (SGF) approach based on QCD to improve the theoretical description of heavy quarkonium production and decay, addressing limitations of the traditional NRQCD method.

Contribution

The paper proposes a rigorously defined SGF approach that offers better convergence than NRQCD and enhances understanding of quarkonium production and decay mechanisms.

Findings

SGF has better convergence than NRQCD.

SGF simplifies to NRQCD in phenomenological applications.

Provides new insights into quarkonium production mechanisms.

Abstract

The widely used nonrelativistic QCD (NRQCD) factorization theory now encounters some notable difficulties in describing quarkonium production. This may be due to the inadequate treatment of soft hadrons emitted in the hadronization process, which causes bad convergence of velocity expansion in NRQCD. In this paper, starting from QCD we propose a rigorously defined factorization approach, soft gluon factorization (SGF), to better deal with the effects of soft hadrons. After a careful velocity expansion, the SGF can be as simple as the NRQCD factorization in phenomenological studies, but has a much better convergence. The SGF may provide a new insight to understand the mechanisms of quarkonium production and decay.

Equations

$$\begin{array}{rl}(2\pi {)}^{3}2P_H^0\frac{d{\sigma }_H}{d^3{P}_H}=& \sum _n{\mathcal{H}}_n(P_H)⟨{\mathcal{O}}_n^H⟩+\cdots ,\end{array}$$

$$\begin{array}{rl}\frac{d{\sigma }_H}{d{P}_{HT}^2}& \sim \frac{1}{P_T^c}\sim \frac{(1-\lambda {)}^c}{P_{HT}^c},\end{array}$$

$$(2\pi {)}^{3}2P_H^0\frac{d{\sigma }_H}{d^3{P}_H}\approx \sum _n\int \frac{d^{4}P}{(2\pi {)}^4}{\mathcal{H}}_n(P)F_{n\to H}(P,P_H),$$

$$\displaystyle\begin{split}&F_{n\to{H}}(P,P_{H})=\int d^{4}be^{-iP\cdot b}\,\\ &\times\langle 0|[\overline{\Psi}{\cal K}_{{n}}\Psi]^{\dagger}(0)\big{(}a_{H}^{\dagger}a_{H}\big{)}[\overline{\Psi}{\cal K}_{n}\Psi](b)|0\rangle_{\text{S}},\end{split}$$

$$a_H^{†}{a}_H=\sum _X\sum _{J_z^H}|H+X⟩⟨H+X|$$

$${\mathcal{K}}_n(rb)=\frac{\sqrt{M_H}}{M_H+2m}\frac{M_H+/P_H}{2M_H}{\Gamma }_n\frac{M_H-/P_H}{2M_H}{\mathcal{C}}^{[c]},$$

$${\mathcal{C}}^{[1]}$$

$${\mathcal{C}}^{[8]}$$

$${\Phi }^{(A)}(rb)=\mathcal{P}\text{exp}\left\{-i{g}_s{\int }_0^{\infty }d\lambda b_{\ell }\cdot A^{(A)}(rb+\lambda b_{\ell })\right\},$$

$$\begin{array}{rl}& F_{Q\overline{Q}({}^{2S+1}{L}_J^{[c]})\to H}\equiv \sum _{J_z}{F}_{Q\overline{Q}({}^{2S+1}{L}_{J,J_z}^{[c]})\to H}.\end{array}$$

$${\overline{d\sigma }}_n(P)\approx \sum _{n^{\prime }}\int \frac{d^4{P}^{\prime }}{(2\pi {)}^4}{\mathcal{H}}_{n^{\prime }}(P^{\prime })F_{n^{\prime }\to n}(P^{\prime },P),$$

$${\overline{d\sigma }}_n(P)\equiv (2\pi {)}^{3}2P^0\frac{d{\sigma }_n}{d^{3}P}.$$

$$p_Q=P/2+q,p_{\overline{Q}}=P/2-q.$$

$$\text{Nc1c for color singlet},$$

$$\text{Nc2−12ta for color octet},$$

$$\frac{2}{\sqrt{M}(M+2m)}(/p_{\overline{Q}}-m)\frac{M-/P}{2M}{\tilde{\Gamma }}_{S{S}_z}^s\frac{M+/P}{2M}(/p_Q+m),$$

$${\tilde{\Gamma }}_{00}^s$$

$${\tilde{\Gamma }}_{1S_z}^s$$

$$P\cdot q=0,q^2=m^2-P^2/4.$$

$$\sum _{L_z,S_z}⟨L,L_z;S,S_z|J,J_z⟩\int d^2\Omega |\mathbf{q}{|}^{-L}\sqrt{\frac{(2L+1)!!}{4\pi (L!)}}{Y}_L^{\ast L_z}.$$

$$F_{n^{\prime }\to n}^{(0)}(P^{\prime },P)={\delta }_{n^{\prime }n}(2\pi {)}^4{\delta }^4(P^{\prime }-P),$$

$${\overline{d\sigma }}_n$$

$$F_{n^{\prime }\to n}$$

$${\mathcal{H}}_n$$

$${\mathcal{H}}_n^{(0)}(P)=$$

$$\begin{array}{rl}{\mathcal{H}}_n^{(1)}(P)=& {\overline{d\sigma }}_n^{(1)}(P)\\ & -\sum _{n^{\prime }}\int \frac{d^4{P}^{\prime }}{(2\pi {)}^4}{\overline{d\sigma }}_{n^{\prime }}^{(0)}(P^{\prime })F_{n^{\prime }\to n}^{(1)}(P^{\prime },P),\end{array}$$

$$\begin{array}{rl}{\mathcal{H}}_n^{(2)}(P)=& {\overline{d\sigma }}_n^{(2)}(P)\\ & -\sum _{n^{\prime }}\int \frac{d^4{P}^{\prime }}{(2\pi {)}^4}{\overline{d\sigma }}_{n^{\prime }}^{(1)}(P^{\prime })F_{n^{\prime }\to n}^{(1)}(P^{\prime },P)\\ & -\sum _{n^{\prime }}\int \frac{d^4{P}^{\prime }}{(2\pi {)}^4}{\overline{d\sigma }}_{n^{\prime }}^{(0)}(P^{\prime })F_{n^{\prime }\to n}^{(2)}(P^{\prime },P),\end{array}$$

$$(2\pi {)}^{3}2P_H^0\frac{d{\sigma }_H}{d^3{P}_H}\approx \sum _n{\mathcal{H}}_n(P_H)⟨{\tilde{\mathcal{O}}}_n^H⟩,$$

$$⟨{\tilde{\mathcal{O}}}_n^H⟩=\int \frac{d^{4}P}{(2\pi {)}^4}{F}_{n\to H}(P,P_H)$$

$$⟨{\mathcal{O}}_n^H⟩\approx$$

$$(2\pi {)}^{3}2P_H^0\frac{d{\sigma }_H}{d^3{P}_H}\approx \sum _n\int dz{\mathcal{H}}_n(P_H/z)F_{n\to H}(z),$$

$$F_{n\to H}(z)=\int \frac{d^{4}P}{(2\pi {)}^4}\delta (z-\sqrt{P_H^2/P^2})F_{n\to H}(P,P_H).$$

$$(2\pi {)}^{3}2P_H^0\frac{d{\sigma }_H}{d^3{P}_H}\approx \sum _n\int dz{\mathcal{H}}_n(P_H/z)F_{n\to H}^{\hat{e}}(z),$$

Full text

I Introduction

As the simplest bound state of strong interactions, heavy quarkonium is a perfect system to study both perturbative and nonperturbative physics of QCD. Ever since the discovery of the first heavy quarkonium $J/\psi$ in 1974, a lot of efforts have been devoted to interpret the production and decay mechanisms of heavy quarkonium. Among them, the most notable theories include the color-singlet model Ellis:1976fj ; Carlson:1976cd ; Chang:1979nn and the nonrelativistic QCD (NRQCD) factorization theory Bodwin:1994jh 111For quarkonium production only, there are two more well-established theories. One is the color-evaporation model Fritzsch:1977ay ; Halzen:1977rs ; Ma:2016exq , which is a phenomenological model that may not be able to derive from first principles of QCD. The other is QCD collinear factorization Kang:2011mg ; Fleming:2012wy ; Kang:2014tta , which is rigorous but can only describe high momentum quarkonium production..

The NRQCD factorization approach is successful. With the color-octet mechanism, NRQCD can solve the infrared divergence problem encountered in the color-singlet model Bodwin:1992qr , explain the $\psi(2S)$ surplus Braaten:1994vv , and describe the inclusive quarkonium production Boyd:1998km ; Kramer:2001hh ; Klasen:2001cu ; Zhang:2006ay ; Lansberg:2008gk ; Ma:2008gq ; Gong:2009kp ; Butenschoen:2009zy ; Ma:2010vd ; Ma:2010yw ; Butenschoen:2010rq ; Ma:2014mri . In addition, although there is still no convincing all-order proof of NRQCD factorization for quarkonium production, it was found that factorization may hold at least to next-to-next-to-leading order (NNLO) if long-distance matrix elements (LDMEs) are modified to be gauge complete Nayak:2005rt ; Nayak:2005rw ; Nayak:2006fm .

Nevertheless, studies in recent years have shown that NRQCD factorization in describing quarkonium production may encounter some notable difficulties. (1) The polarization puzzle: The leading order calculation in NRQCD implies that $\psi(nS)$ and $\Upsilon(nS)$ produced at hadron colliders are transversely polarized due to ${{}^{3}S_{1}^{[8]}}$ channel dominance Cho:1994ih ; Beneke:1995yb ; Braaten:1999qk . But experimental measurements found these states almost unpolarized Abulencia:2007us ; Abelev:2011md ; Aaij:2013nlm ; Chatrchyan:2013cla ; Chatrchyan:2012woa . Thanks to next-to-leading order calculations Gong:2008hk ; Butenschoen:2012px ; Chao:2012iv ; Gong:2012ug ; Bodwin:2014gia ; Faccioli:2014cqa ; Gong:2013qka ; Han:2014kxa , the observed polarizations of $J/\psi$ and $\Upsilon(nS)$ can be qualitatively explained, but it is still hard to understand the polarization of $\psi(2S)$ Shao:2014yta . (2) The hierarchy problem: The best fit of $J/\psi$ yield data at high transverse momentum in hadronic collisions determines two linearly combined LDMEs, $M_{0}=0.074\text{GeV}^{3}$ and $M_{1}=0.0005\text{GeV}^{3}$ Ma:2010yw (the $J/\psi$ polarization data requires almost the same two combined LDMEs Chao:2012iv ). There is a two-orders difference between the two combined LDMEs. However, velocity scaling rules in NRQCD Bodwin:1994jh expect these LDMEs to be at the same order of magnitude. (3) The universality problem: A necessary condition for NRQCD factorization is that LDMEs are universal, i.e. process independent. Yet phenomenology studies show that $M_{0}$ extracted from hadron colliders Ma:2010yw ; Butenschoen:2010rq ; Gong:2012ug is much larger than the upper bound set by $e^{+}e^{-}$ collisions $M_{0}<0.02\text{GeV}^{3}$ Zhang:2009ym .

In this paper, we show that the velocity expansion in the present NRQCD framework suffers from large high order relativistic corrections due to ignoring the momentum of soft hadrons, which are mainly soft gluons perturbatively, emitted in the hadronization process. Thus, by including only a few low order contributions in relativistic expansion, NRQCD is hard to provide good descriptions for quarkonium production, which may be the reason for the above mentioned difficulties. Starting from QCD, we propose a new factorization approach, called soft gluon factorization (SGF), to better deal with the effects of soft hadrons. The SGF has a much better convergence in the velocity expansion, and it may provide a new insight to fully understand the mechanisms of quarkonium production and decay.

The rest of the paper is organized as follows. In Sec. II, we explain why NRQCD factorization may have bad convergence in velocity expansion, and then define the SGF formula with rigorous operator definition. Relations to other approaches are also discussed. In Sec. III, we do a careful simplification to the SGF so that it will be easy to use in practice, but at the same time it can capture main physics. Following this, we compare the SGF and NRQCD by studying the gluon fragmentation to yield of $J/\psi$ in Sec. IV. As expected, we find that the lowest order NRQCD approximation is not good, which may be the reason why NRQCD encounters many difficulties. The summary and the possibility of solving these difficulties in the SGF framework are given in Sec. V.

II Soft gluon factorization

II.1 Bad convergence of $v^{2}$ expansion in NRQCD

In NRQCD factorization Bodwin:1994jh , the differential cross section of a heavy quarkonium ${H}$ production can be factorized as

[TABLE]

where $P_{H}$ is the momentum of ${H}$, $n$ denotes intermediate $Q{\overline{Q}}$ states, whose quantum numbers are usually expressed in terms of spectroscopic notation ${{}^{{2S+1}}L_{J}^{[c]}}$, with $c=1,8$ denotes color singlet or color octet of the pair, ${\cal H}_{n}(P_{H})$ are perturbative calculable short-distance coefficients which can be expanded order by order in $\alpha_{s}$, $\langle{\cal O}_{n}^{H}\rangle$ are gauge-completed Nayak:2005rt ; Nayak:2005rw ; Nayak:2006fm nonperturbative long-distance matrix elements (LDMEs) which can be classified according to the power counting rules of $v$, and the ellipses denote other relativistic correction terms. It is needed to emphasize that, in NRQCD factorization, one expands the mass of quarkonium $M_{H}$ around twice of heavy quark mass $2m$, which results in $P_{H}^{2}\approx 4m^{2}$ in Eq. (1).

If convergence of velocity expansion in NRQCD factorization is very good, one only needs to consider a very limited number of LDMEs to describe experimental data. Unfortunately, it may not be the case. To see this, for an example let us consider the differential cross section of quarkonium production in hadron colliders with high transverse momentum $P_{{H}T}$. In NRQCD factorization, one calculates the partonic $Q{\overline{Q}}$ production differential cross sections with transverse momentum $P_{T}$, and then expands $P_{T}$ around $P_{{H}T}$. On average, let us set $P_{{H}T}\sim(1-\lambda)P_{T}$, with $\lambda$ being the average transverse momentum fraction carried by soft gluons which can be at the order of $v$ or $v^{2}$. At a specific range of large $P_{{H}T}$, differential cross section behaves as

[TABLE]

with $c$ being usually larger than 4. Expansion of this result with respect to $\lambda$ mimics the relativistic expansion in NRQCD. For a typical choice $\lambda=0.2$ and $c=4$ Mangano:1996kg , we have $(1-\lambda)^{c}=0.4096$, but the expansion of $(1-\lambda)^{c}$ to $O(\lambda^{0})$ gives 1, and expansion to $O(\lambda^{1})$ gives 0.2. The above argument shows a bad convergence of relativistic expansion. We note that the convergence still cannot be improved if one instead is to calculate, e.g., $P_{{H}T}^{4}\frac{d\sigma_{H}}{dP_{{H}T}^{2}}$. Furthermore, as we will see in Sec. IV, there are other relativistic correction terms that make the convergence of relativistic expansion even worse.

Because of the bad convergence of relativistic corrections, theoretic calculations based on only a limited number of LDMEs are sometimes hard to describe experimental data. In literature, e.g. Refs. Beneke:1997qw ; Fleming:2003gt ; Fleming:2006cd ; Leibovich:2007vr , for each problem one can resum a specific subset of LDMEs to improve the theoretic results. The resummation results in many so-called “shape functions”.

A factorization method aiming to improve convergence for an arbitrary problem was proposed in Ref. Beneke:1999gq . However, this method can only be thought of as a model because there is no operator definition for nonperturbative functions. Without an operator definition, the method is not well defined and it is not possible to do rigorous calculation beyond tree level.

II.2 Soft gluon factorization formula

The aim of soft gluon factorization is to resum a subset of relativistic correction terms in NRQCD factorization that are important for a phenomenological purpose. It is convenient to demand the subset to be Lorentz invariant. The formula of SGF for a quarkonium ${H}$ production is

[TABLE]

where ${\cal H}_{n}(P)$ are perturbatively calculable hard parts that, roughly speaking, produce an intermediate state with quantum numbers $n$ and momentum $P$, and $F_{n\to{H}}(P,P_{H})$ are nonperturbative functions, which we call soft gluon distribution functions (SGDs), which describe the hadronization of the intermediate state to physical quarkonium ${H}$. To account for the effect of soft hadrons emission, which are mainly soft gluons perturbatively, the momentum of the observed quarkonium $P_{H}$ is kept different from the momentum of the intermediate state $P$, which is different from the treatment in NRQCD.

For quarkonium decay, one can define a similar formula.

II.3 Intermediate states

Naively, in Eq. (3) we should sum over a complete set of intermediate states that may or may not contain a $Q{\overline{Q}}$ pair. However, if there is no $Q{\overline{Q}}$ pair in the intermediate state, a $Q{\overline{Q}}$ pair must be produced in the hadronization process. In this case, the corresponding $F_{n\to{H}}(P,P_{H})$ has short-distance effects and can be refactorized, which enables us to eliminate this contribution from the factorization formula. Furthermore, if the intermediate state contains other energetic partons in addition to a $Q{\overline{Q}}$ pair, the corresponding $F_{n\to{H}}(P,P_{H})$ can also be refactorized because there is no energetic gluon or light quark in the dominant Fock state of conventional quarkonium. In this consideration, we only need to sum over intermediate states which contain a $Q{\overline{Q}}$ pair and some soft partons. According to physical scales existing in the quarkonium system, energy of these soft partons can be at the order of $mv$ or $mv^{2}$, which are assumed to be much smaller than $m$.

As each conventional quarkonium has a $Q{\overline{Q}}$ state with specific quantum numbers as its leading Fock state, the $Q{\overline{Q}}$ intermediate state with the same quantum numbers has leading contribution. If the intermediate state contains soft partons in addition to a $Q{\overline{Q}}$ pair, which means that soft partons are produced in the hard parts, the contribution is at least suppressed by $v^{2}$ comparing with the contribution from the intermediate state with the same $Q{\overline{Q}}$ but without soft partons. In other words, power counting of the corresponding $F_{n\to{H}}(P,P_{H})$ in the factorization formula is at high order in $v^{2}$. Unlike $v^{2}$ corrections due to kinematic effects of soft gluons emission, contribution of intermediate states with soft partons should have a mild $v^{2}$ corrections. Thus, as the first approximation, we will ignore these contributions and consider only intermediate states constituted by a $Q{\overline{Q}}$ pair in the rest of the paper.

Then $n$ in Eq. (3) denotes intermediate $Q{\overline{Q}}$ states, whose quantum numbers are usually expressed in terms of spectroscopic notation $n=Q{\overline{Q}}({{}^{{2S+1}}L_{J,J_{z}}^{[c]}})$, with $c=1,8$ denoting color singlet or color octet of the pair. We note that the state in the amplitude $n$ can be in general different from that in the complex-conjugate amplitude $\tilde{n}=Q{\overline{Q}}({{}^{{2\tilde{S}+1}}\tilde{L}_{\tilde{J},\tilde{J}_{z}}^{[\tilde{c}]}})$. Color charge, $C$-parity, $P$-parity and angular momentum conservations provide general selection rules Ma:2015yka , which gives that the following relations always hold: $c=\tilde{c}$, $S=\tilde{S}$, and $|L-\tilde{L}|=0,2,4,\cdots$. For the production of a polarization-summed quarkonium, we have further constraints $J=\tilde{J}$, $J_{z}=\tilde{J}_{z}$. Thus, even in the polarization-summed case, there is interference between the ${{}^{{3}}S_{1,J_{z}}^{[c]}}$ state and the ${{}^{{3}}D_{1,J_{z}}^{[c]}}$ state. Selection rules for polarization distributions can be found in Ref. Ma:2015yka . For simplicity, in the rest of the paper we will only discuss polarization-summed quarkonium production and choose the quantum numbers of the intermediate $Q{\overline{Q}}$ pair in the amplitude to be the same as that in the complex-conjugate amplitude, while polarization and interference contributions can be studied similarly.

II.4 Nonperturbative distributions

To make SGDs Lorentz invariant, it is convenient to define them using QCD fields instead of NRQCD fields. We define SGDs for polarization-summed quarkonium ${H}$ production as vacuum expectation values of bilocal operators constructed from QCD fields,

[TABLE]

where

[TABLE]

projects final state to include a polarization-summed ${H}$ with relativistic normalization $\langle{H}(P_{H})|{H}(P_{H}^{\prime})\rangle=(2E_{H})(2\pi)^{3}\delta^{3}({\bm{P}}_{H}-{\bm{P}}_{H}^{\prime})$. The subscript “S” means that, to evaluate the matrix element, one only picks up integration regions where off-shellness of all particles is much smaller than heavy quark mass 222From the point view of method of region Beneke:1997zp , the effect of “S” keeps all regions except the hard region. . ${\cal K}_{n}(rb)$ are projection operators defining the intermediate state $n$,

[TABLE]

where angular momentum operators $\Gamma_{n}$ are defined explicitly in Appendix A, color operators ${\cal C}^{[c]}$ will be defined later, $\frac{M_{H}\pm\not{P}_{H}}{2M_{H}}$ project out large components in velocity expansion, and, as we will see later, the introduction of prefactor $\frac{\sqrt{M_{H}}}{M_{H}+2m}$ makes our SGDs more closely related to NRQCD LDMEs.

If the intermediate $Q{\overline{Q}}$ pair is in color singlet, we define color operator as

[TABLE]

which is nothing but SU(3) color Clebsch-Gordan coefficient $\langle 3i;\bar{3}\bar{i}|00\rangle$ and ${\bm{1}}_{c}$ is the identity matrix in color space. If the intermediate $Q{\overline{Q}}$ pair is in color octet, we define color operator by the multiplication of a color matrix in fundamental representation with a gauge link in adjoint representation,

[TABLE]

where $\sqrt{2}t^{\bar{a}}$ is SU(3) color Clebsch-Gordan coefficient $\langle 3i;\bar{3}\bar{i}|8\bar{a}\rangle$ and the introduction of gauge link $\Phi^{(A)}_{\bar{a}a}(rb)$ is to enable gauge invariance of SGDs.

The next question is what direction of the gauge link should we choose. A natural choice is to define the direction along a light cone, because a gauge link can be obtained by eikonal approximation for soft gluon interaction between the $Q{\overline{Q}}$ pair and jets. Indeed, this is the choice to define gauge-completed NRQCD LDMEs in Refs. Nayak:2005rt ; Nayak:2005rw ; Nayak:2006fm . In the SGF framework, however, this choice will result in uncanceled gauge-link-collinear rapidity divergence because the position $b$ in Eq. (4) is in general not along a light cone. This is the same problem encountered in transverse-momentum-dependent (TMD) factorization where there is a mismatch of rapidity divergence between virtual correction and real correction Collins:2011zzd . To cancel the rapidity divergence, soft factor as while as jet functions may need to be introduced.

If we do not want to introduce additional soft factors 333As we will discuss later, the SGF Eq. (3) is only a conjecture and it is not clear right now whether it is valid to all orders. If Eq. (3) does not hold to all orders, maybe a modification by introducing a gauge invariant soft factor is unavoidable., the gauge links can be defined along the $b$ direction,

[TABLE]

where ${\cal P}$ denotes path ordering, $A^{(A)}$ is gluon field with color matrix in adjoint representation, and $b_{\ell}^{\mu}=b^{\mu}+\varepsilon\ell^{\mu}$. We choose $0<\varepsilon\ll 1$ and a fixed light-like vector $\ell^{\mu}$ so that: when $b^{\mu}$ is finite, $b_{\ell}^{\mu}$ is the same as $b^{\mu}$; while as $b^{\mu}\to 0$, the gauge link is well defined (along the $\ell^{\mu}$ direction).

In the above, we have defined SGDs for $n$ with fixed $J_{z}$. SGDs for polarization-summed intermediate states can be defined accordingly,

[TABLE]

Similar to LDMEs in NRQCD, each SGD has a definite power counting in $v^{2}$, which will be explained later. Thus, for a certain accuracy, we can truncate to use only a few SGDs for each ${H}$ production.

II.5 Perturbative matching

With the definition of SGDs, short-distance hard parts can be obtained order by order in perturbation theory by matching both sides of Eq. (3). As ${\cal H}_{n}(P)$ are independent of nonperturbative physics, they are unchanged if we replace the quarkonium ${H}$ in Eq. (3) by an on-shell $Q{\overline{Q}}$ pair with specific quantum numbers. To make the following equations more convenient to use, we first relabel $n$ and $P$ in Eq. (3) by $n^{\prime}$ and $P^{\prime}$, respectively, and we then replace ${H}$ in both sides of the equation by a state $n=Q{\overline{Q}}({{}^{{2S+1}}{L}_{{J,J_{z}}}^{[{c}]}})$ with momentum $P$, which results in

[TABLE]

where

[TABLE]

Note that, although $\overline{d\sigma_{H}}$ and $F_{n^{\prime}\to{H}}$ are changed to $\overline{d\sigma_{n}}$ and $F_{n^{\prime}\to n}$, respectively, ${\cal H}_{n^{\prime}}(P^{\prime})$ stays unchanged, which is guaranteed by the fact that ${\cal H}_{n^{\prime}}(P^{\prime})$ is perturbatively calculable and thus it is independent of nonperturbative external states. The reason for choosing the on-shell $Q{\overline{Q}}$ pair is that $\overline{d\sigma_{n}}$ and $F_{n^{\prime}\to n}$ are then gauge invariant and can be expanded order by order in $\alpha_{s}$ in any gauge.

The key to determine all ${\cal H}_{n^{\prime}}$ is to project ${H}$ to a complete set of $Q{\overline{Q}}$ states, and then to match coefficients of $\alpha_{s}$ in both sides of Eq. (3). Although any complete set is doable for this purpose, a good choice should satisfy as much as possible that leading order in $\alpha_{s}$ expansions of $F_{n^{\prime}\to n}$ are delta functions.

To define a good complete set, we begin with a $Q{\overline{Q}}$ with momenta

[TABLE]

We project it to color state $c$ by color projectors:

[TABLE]

where $\sqrt{N_{c}^{2}-1}$ in the denominator is to average over color-octet states. We then project it to state with total spin $S$ and $S_{z}$ by replacing spinors of $Q{\overline{Q}}$ by444This expression can be interpreted as covariant spin projectors Kuhn:1979bb ; Guberina:1980dc multiplied by a normalization factor $\sqrt{2/M}$ Fan:2009cj . But for our purpose, this interpretation is unnecessary. The key in our definition is to guarantee the orthogonal relations Eq. (19).

[TABLE]

where $M^{2}=P^{2}$ and

[TABLE]

On-shell conditions $p_{Q}^{2}=p_{{\overline{Q}}}^{2}=m^{2}$ result in

[TABLE]

There are still 2 degrees of freedom of the relative momentum $q$, which can be chosen as spatial angles in the rest frame of the pair. We further do a partial wave expansion in this frame, which removes the $q$ dependence although introducing orbital angular momentum $L$ and $L_{z}$. Finally, we project spin and orbital angular momentum to total momentum $J$ and $J_{z}$. More precisely, in the rest frame of $P$, $d\sigma_{n}$ and $F_{n^{\prime}\to n}$ with quantum number $n=Q{\overline{Q}}({{}^{{2S+1}}{L}_{{J,J_{z}}}^{[{c}]}})$ can be obtained from the corresponding production of $Q{\overline{Q}}$ with fixed $c$, $S$ and $S_{z}$ by the following operation:

[TABLE]

In Appendix B, we show that the above definition can indeed result in orthogonal relations at lowest order

[TABLE]

where the superscript “$(0)$” denotes leading order in $\alpha_{s}$ expansion.

By inserting perturbative expansions

[TABLE]

into Eq. (11) and using the orthogonal relations Eq. (19), coefficients of different powers of $\alpha_{s}$ give the following relations:

[TABLE]

and so on. Based on these relations, ${\cal H}_{n}$ can be obtained by perturbative calculation of $\overline{d\sigma}_{n}$ and $F_{n^{\prime}\to n}$ for on-shell $Q{\overline{Q}}$ pair production. This perturbative calculation relies on the cancellation of IR divergences between $\overline{d\sigma}_{n}$ and $F_{n^{\prime}\to n}$ [like in Eq. (21)] to all orders in perturbation theory, which is a very difficult problem and will be discussed in the next subsection.

It is worth emphasizing that the above perturbative calculation of on-shell $Q{\overline{Q}}$ pair production can only directly determine ${\cal H}_{n}(P)$ with $P^{2}>4m^{2}$. Although on average we may have $P^{2}>4m^{2}$ Ma:2016exq , there is a nonvanishing contribution from the $P^{2}<4m^{2}$ region in Eq. (3). The value of ${\cal H}_{n}(P)$ in the later region can be obtained by analytical continuation of its value in the former region. Analytical continuation here is almost trivial, which means that we use the same functional form of ${\cal H}_{n}(P)$ in all regions of $P^{2}$.

II.6 Justification

The above proposed SGF is only a conjecture, and we cannot provide an all order proof at present. The most dangerous interaction that may ruin the factorization is the elastic scattering between the $Q\bar{Q}$ pair and hard jets by exchanging gluons, the kinematic region of these exchanged gluons are called “Glauber” region in literature. We find that IR divergence from the “Glauber” region at one-loop order is purely imaginary, which cancels with the contribution from its complex conjugate diagram. Therefore, SGF holds at one-loop order. Thanks to the two-loop study of infrared divergences of quarkonium production Nayak:2005rt ; Nayak:2005rw ; Nayak:2006fm , it is hopeful that SGF may also hold to two-loop order, which however still needs further examination.

II.7 Relation to NRQCD factorization

The difference between NRQCD and SGF can be partially understood from the treatment of intermediate momenta. Based on relations in Eq. (17), we can express $q^{\prime}_{0}$ and $\bm{q}^{\prime 2}$ by $P^{\prime}_{0}$ and $\bm{P}^{\prime 2}$, which corresponds to the SGF strategy. Alternatively, we can also use Eq. (17) to express $q^{\prime}_{0}$ and $P^{\prime}_{0}$ by $\bm{q}^{\prime 2}$ and $\bm{P}^{\prime 2}$, which corresponds to the strategy of NRQCD factorization.

To further discuss the relation between SGF and NRQCD, by ignoring velocity corrections we can approximate ${\cal H}_{n}(P)$ by ${\cal H}_{n}(P_{H})$ in Eq. (3). Then the integral over $P$ is applied on SGDs defined in Eq. (4), and thus Eq. (3) becomes

[TABLE]

which has the same form as the NRQCD factorization Eq.(1). In the above formula,

[TABLE]

are closely related to gauge-completed LDMEs $\langle{\mathcal{O}}_{n}^{H}\rangle$ although they are defined by different fields. Considering the differences between normalization in color space, up to lowest order in $v^{2}$ approximation, we have 555Let us first check color-octet states. In NRQCD one usually chooses nonrelativistic normalization for the ${H}$ state, which means that our definition has an extra factor $2E_{H}$ due to normalization of state. For a color operator, in NRQCD one uses $t^{\bar{a}}$ instead of $\sqrt{2}t^{\bar{a}}$, due to which we have an extra factor $2$. Finally, the prefactor in Eq. (6) introduces an extra factor $M_{H}/(M_{H}+2m)^{2}$. The multiplication of these three factors equals $1$ if we choose the ${H}$ rest frame and take $v^{2}\to 0$. Except these differences, all other parts of our definition are identical to the NRQCD definition at the lowest order in $v^{2}$ approximation. For color-singlet states, in NRQCD one uses $\bm{1}_{c}$ as color operator, which results in a difference of factor $2N_{c}$ comparing with the color-octet case.

[TABLE]

Based on NRQCD velocity scaling rules, these relations also tell us $v^{2}$ power counting rules for SGDs.

Beyond the above approximation, the complete SGF formula resums a series of velocity corrections in NRQCD, which are important for phenomenological study. As we know, the TMD factorization is a generalization of collinear factorization with a series of higher twist corrections resummed Collins:2011zzd . In this sense, we can say that SGF is a “TMD” version of NRQCD. We expect that there are a lot of similarities between SGF and TMD factorization.

II.8 Relation to QCD collinear factorization

It is also needed to explain the relation between the SGF method and the QCD collinear factorization method Kang:2011mg ; Fleming:2012wy ; Kang:2014tta . In QCD collinear factorization, one assumes the transverse momentum of the produced quarkonium to be much larger than the heavy quark mass, $P_{HT}\gg m$. Therefore, one can factorize out hard physics at the scale of $P_{HT}$, and leave the physics at the scale of $m$ to be described by input functions, which are called fragmentation functions (FFs). Unlike NRQCD factorization and SGF, QCD collinear factorization has been proved to all orders in perturbation theory at both leading power and next-to-leading power levels. Phenomenologically, the advantage of this method is that large logarithms $\log(P_{HT}/m)$ can be resummed to all orders in perturbation theory using renormalization group equations for single-parton and double-parton FFs Kang:2014tta , and therefore there can be smaller theoretical uncertainties.

Due to the above advantages, when studying high transverse momentum quarkonium production, one should better first use QCD collinear factorization to express cross sections in terms of the convolution of perturbative calculable hard parts Kang:2014pya with FFs. Based on renormalization group equations, these FFs are fully determined by their values at an initial scale. Because of the heavy quark mass, these FFs at the initial scale still have perturbative physics that can be further separated from nonperturbative physics. This can be achieved by using NRQCD or SGF. Using NRQCD factorization, these FFs at an initial scale have been studied extensively Braaten:1993mp ; Braaten:1993rw ; Braaten:1994kd ; Braaten:1995cj ; Ma:1995vi ; Braaten:2000pc ; Bodwin:2003wh ; Hao:2009fa ; Jia:2012qx ; Bodwin:2012xc ; Ma:2013yla ; Ma:2014eja . In Sec. IV, we will see that FFs can also be calculated by using the SGF method, which can improve the convergence of $v^{2}$ expansion comparing with the using of NRQCD factorization.

Note, however, that both NRQCD factorization and SGF can also be used to directly describe low transverse momentum quarkonium production as well as quarkonium decay, where the QCD collinear factorization does not apply.

III Simplification

III.1 Collinear approximation

The factorization formula in Eq. (3) can be further simplified to make it more suitable in practice use. In the rest frame of ${H}$, SGDs in Eq. (3) have support only in the region $P^{\mu}=(M+O(\lambda^{2}/M),O(\lambda),$ $O(\lambda),O(\lambda))$, where $M^{2}=P^{2}$ and $\lambda\sim\text{a few hundreds MeV}<<M$ is the typical energy of emitted gluons. Thus we can expand $O(\lambda)$ and $O(\lambda^{2}/M)$ terms in hard parts. In this expansion, although space components of $P^{\mu}$ are $O(\lambda)$, the convergence should be much better than usual velocity expansion in NRQCD because $\vec{P}$ is integrated out symmetrically around the origin. The leading term in this expansion gives

[TABLE]

where

[TABLE]

In this way, we significantly reduce the complication of input functions from four dimensional to one dimensional. We refer to Eq. (25) as SGF-1d and Eq. (3) as SGF-4d.

Note that a similar formula that relates momentum of $\chi_{cJ}$ to its decaying particle $J/\psi$, $p_{J/\psi}\approx\frac{m_{J/\psi}}{m_{\chi_{cJ}}}p_{\chi_{cJ}}$, was first introduced by us in Ref. Ma:2010vd , and later also used in other publications Gong:2012ug ; Bodwin:2015iua ; Ma:2016exq . The approximation was found to be very good, with only 8% deviation from full results for both yield and polarization measurements Bodwin:2015iua . We thus expect the approximation Eq.(25) to be also very good.

There are also other choices of collinear approximations. For any direction ${\hat{e}}$, we can have an approximation

[TABLE]

with

[TABLE]

In the rest of the paper, we will only take Eq. (25) as an example to show the effect of collinear approximation.

III.2 Statical approximation

As hard gluon emissions are excluded from $F_{n\to{H}}(z)$ and very soft gluon emissions are suppressed by phase space, $F_{n\to{H}}(z)$ should peak around $z=z_{n}\sim 1-O(\lambda/m_{H})$. If the distribution of $F_{n\to{H}}(z)$ is narrow enough, we can approximate it by $\delta(z-z_{n})\,\langle\widetilde{\mathcal{O}}_{n}^{H}\rangle$. Then the SGF becomes

[TABLE]

which we call SGF-0d. While this formula is as simple as NRQCD factorization, we will see later that the lowest order SGF-0d is a much better approximation of real physics than that of the lowest order NRQCD.

III.3 Expansion of $m$

In hard parts of all the above SGFs, there are at least two independent hard scales, $2m$ and $M$. Their difference is, however, much smaller than $2m$. We thus can further simplify hard parts by expanding $m$ around $M/2$, which defines a relativistic correction series in SGFs. We will see that there are no large corrections in this expansion.

IV Example: $J/\psi$ hadroproduction

IV.1 High $p_{T}$ quarkonium production

To compare SGFs with NRQCD factorization, we apply them to the $J/\psi$ hadron production via gluon fragmentation. A charmonium ${H}$ production cross section via gluon fragmentation is given by

[TABLE]

where $d\hat{\sigma}_{g}(p_{T}/x)$ is the well-known hard part that produces a gluon with transverse momentum $p_{T}/x$, and $D_{g\to{H}}(x)$ is the FF of a gluon into an ${H}$ that will be determined in the following.

IV.2 Fragmentation functions calculated in SGF

Using the SGF-4d Eq. (3), we obtain

[TABLE]

where

[TABLE]

which only depends on $x/y$ and $P^{2}$ because of Lorentz symmetry combined with boost invariance along the “+” direction. Using the SGF-1d Eq. (25), we obtain

[TABLE]

In Eqs. (31) and (33), $\hat{D}_{n}(y,P^{2})$ can be calculated perturbatively. Up to $O(\alpha_{s})$, only the $n={{}^{3}S_{1}^{[8]}}$ channel has a nonvanishing contribution machao ,

[TABLE]

with $\Delta=\frac{M}{2m}$. Equation (34) can also be used in NRQCD factorization, where one expresses $\Delta=\sqrt{1+\langle v\rangle^{2}}$ and then expands $\langle v\rangle^{2}$, which gives a normalized series $1-\frac{11}{6}\langle v\rangle^{2}+\frac{191}{72}\langle v\rangle^{4}+\cdots$. The first two terms in this expansion have been calculated before Braaten:1994kd ; Bodwin:2003wh .

IV.3 Model assumptions

For the nonperturbative function $F_{n\to{H}}(P,P_{H})$, we simply employ an existing model in Ref. Beneke:1999gq ,

[TABLE]

where $k=P-P_{H}$ is a timelike momentum with positive energy $k_{0}$. $\Lambda$ is an energy cutoff for emitted soft gluons, for which we will choose $500$MeV 666By assuming the cutoff at the order of $mv^{2}$, this choice corresponds to $v^{2}\approx 0.3$. As we will see, $D^{1d}_{g\to{H}}(x)$ calculated in this model is picked around $x\approx 1-v^{2}/2$.. As we are only interested in the cross section ratio in the following, the value of $a$ is irrelevant. We further set $m=1.55$GeV in studying $J/\psi$ production. We note that, our conclusions are in fact independent of these choices.

IV.4 Numerical results

With these inputs in hand, FFs calculated by using different methods are shown in Fig. 1. In this figure, $D^{4d}_{g\to{H}}(x)$ and $D^{1d}_{g\to{H}}(x)$ are calculated by using the exact value of $m$, while $D^{1d^{(0)}}_{g\to{H}}(x)$ is similar to $D^{1d}_{g\to{H}}(x)$ but calculated by setting $m=M/2$ in Eq. (34). We find that FFs calculated by using SGF-4d and SGF-1d have different shapes, but they have the same accumulated values,

[TABLE]

This is because SGF-4d and SGF-1d are equivalent for the integrated FF. The FF obtained by expanding $m$ to the lowest order, $D^{1d^{(0)}}_{g\to{H}}(x)$, has only a small difference from the complete result $D^{1d}_{g\to{H}}(x)$, which implies that relativistic correction due to the expansion of $m$ is very small.

We also find that $F_{{{}^{3}S_{1}^{[8]}}\to{H}}(z)$ has almost the same shape as that of $D^{1d}_{g\to{H}}(x)$, which is $\delta$-function-like. For cross section, we will find that SGF-0d with $z_{0}=0.86$ can well reproduce SGF-4d.

Using the above FFs, we can calculate the $J/\psi$ cross section based on Eq. (30). Let us assume that the cross section calculated by SGF-4d is “exact”, and examine how good are the SGF-1d, SGF-0d, and NRQCD expansions. We show ratios of $J/\psi$ differential cross section calculated by different methods over that calculated by SGF-4d in Fig. 2. We find that the SGF-1d result is close to the exact value, with largest error less than 6%. This implies that, as expected, the SGF-1d expansion should be very good. By expanding $m$ to the lowest order, SGF-1d*(0)* also provides a good approximation, with deviation smaller than 10%. The SGF-0d can well reproduce SGF-4d if we choose $z_{0}=0.86$, as shown in Fig. 2. So, we expect that the convergence of velocity expansion in SGF is good in general.

On the contrary, the lowest order NRQCD result is larger than the exact value by more than a factor of 4, where NRQCD LDMEs are determined by the approximation Eq. (24). There are two main sources for the large deviation. One comes from the hard part of Eq. (30), which has an approximate scaling behavior $d\hat{\sigma}_{g}(p_{T}/x)\sim(p_{T}/x)^{-4}$ with average value of $x$ being around 0.86. Yet, by ignoring soft gluon emission, NRQCD approximates $x\approx 1$, which enhances the total result by almost a factor of 2. This effect has also been pointed out in Sec.II.1. The other one comes from Eq. (34). By ignoring soft gluon emission and then expanding $M$ around $2m$, NRQCD approximates $\Delta\approx\frac{m_{J/\psi}/0.86}{2m}\approx 1.16$ by 1, which enhances the total result by another factor of 2. Roughly speaking, the lowest order NRQCD approximates $0.86^{9}\approx 1/4\sim(1-v^{2}/2)^{9}$ by 1 in this problem, which is hard to be recovered by traditional relativistic correction in NRQCD Xu:2012am .

For bottomonia production, because $v^{2}$ is smaller, NRQCD approximation can be a little bit better than that for charmonia. According to the above discussion, we expect the lowest order NRQCD approximation overestimates bottomonia production rate by a factor of $(1-v^{2}/2)^{-9}\approx 1.6$.

V Summary and outlook

In this paper, we propose a soft gluon factorization (SGF) approach to describe quarkonium production and decay, which keeps the momentum difference between intermediate $Q{\overline{Q}}$ pair and physical quarkonium. From the point view of NRQCD, the SGF effectively resums a subset of Lorentz-invariant relativistic correction terms in NRQCD factorization. In this sense, the SGF is a generalization of the NRQCD factorization.

By construction, the SGF has a good convergence in velocity expansion, so the lowest order approximation may already capture most physics; whereas for NRQCD factorization, we find that the lowest order result in velocity expansion can deviate from the full result by more than a factor of 4, mainly due to ignoring the momentum taken away by soft gluon emission. With so large deviation, it is not surprising that the NRQCD calculation faces many difficulties.

Hopefully, these difficulties may be resolved or relieved in the SGF framework with well controlled relativistic corrections. Specifically, the universality problem may be due to the fact that, for quarkonium production in different processes, e.g., $e^{+}e^{-}$ collision or $pp$ collision, lowest order NRQCD calculations suffer from large but different relativistic corrections. While in SGF, we do not expect large relativistic corrections. Moreover, considering the large relativistic effect on yield, we may also expect a significant effect on polarization, like color-magnetic dipole transition effect, which may alter the transverse polarization of the ${{}^{3}S_{1}^{[8]}}$ channel. Further examinations should be performed on the above expectations.

Acknowledgements.

We thank G. Bodwin, A.P. Chen, C. Meng, J.W. Qiu, K. Wang, F. Yuan, H. Zhang and Y.J. Zhang for many useful communications and discussions. In particular, we thank J.Z. Li for very helpful discussions in the earlier stage of this work. The work is supported in part by the National Natural Science Foundation of China (Grants No. 11975029, No. 11475005, No. 11875071 and No. 11075002) and the National Key Basic Research Program of China (No. 2015CB856700).

Appendix A Definition of angular momentum operators

Angular momentum operators in Eq. (6) are defined as

[TABLE]

As total spin has only two choices, $S=0$ or $1$, the corresponding spin operators are

[TABLE]

where ${\epsilon}_{S_{z}}^{\mu}$ are polarization vectors. Orbital operators are defined as

[TABLE]

where ${D}_{\mu}$ is the gauge covariant derivative with $\bar{\Psi}\overleftrightarrow{D}_{\mu}\Psi=\bar{\Psi}\left({D}_{\mu}\Psi\right)-\left(D_{\mu}\bar{\Psi}\right)\Psi$ and ${\epsilon}_{L_{z}}^{\mu_{1}\cdots\mu_{L}}$ are $L$-rank polarization tensors. The above polarization tensors ${\epsilon}_{S_{z}}^{\mu}$ and ${\epsilon}_{L_{z}}^{\mu_{1}\cdots\mu_{L}}$ have only spatial components in the rest frame of $P_{{H}\mu}$, which is equivalent to the following relations:

[TABLE]

Furthermore, polarization tensors are chosen to be orthonormal, e.g., ${\epsilon}_{S_{z}}\cdot{\epsilon}_{S_{z}^{\prime}}^{*}=-\delta_{S_{z}S_{z}^{\prime}}$. The combination of the above two properties results in

[TABLE]

Although orbital operators in Eq. (39) are similar to that in NRQCD, their meanings are actually very different. In fact, in addition to these operators, there are an infinite number of relativistic-correction operators in NRQCD, like $(-\frac{i}{2})^{2}\overleftrightarrow{D}^{2}$, which do not show up in SGF. The reason is that, instead of relativistic expansion, the SGF is performing a partial-wave expansion. This can be seen more clearly in momentum space. Then in the rest frame of $P_{H}$, Eq. (39) becomes

[TABLE]

where $q_{\mu}$ is half of the relative momentum between intermediate $Q{\overline{Q}}$, and $\theta$ and $\phi$ are polar angle and azimuthal angle of $\bm{q}$, respectively. Equation (42) can be thought of as the definition of ${\Gamma}^{o}_{LL_{z}}$ in momentum space. As the partial-wave expansion already forms a complete set of operators, we do not need further relativistic-correction operators. Relativistic corrections in SGF are encoded in perturbative calculable hard parts. Our definition of SGDs is to pick up the minimal number of $|{\bm{q}}|$ in each partial wave.

Appendix B Leading order expansion of SGDs

The Feynman diagram which represents the leading order expansion $F_{n\to n^{\prime}}^{(0)}(P,P^{\prime})$ is shown in Fig. 3. It is clear that, if the color of $n$ is different from that of $n^{\prime}$, then the amplitude vanishes. If both $n$ and $n^{\prime}$ are color singlet, we get ${\text{Tr}}_{c}[\frac{1}{\sqrt{N_{c}}}\frac{1}{\sqrt{N_{c}}}]=1$. If both $n$ and $n^{\prime}$ are color octet, for amplitude we have ${\text{Tr}}_{c}[\sqrt{2}t^{a}\sqrt{\frac{2}{N_{c}^{2}-1}}t^{a^{\prime}}]=\frac{\delta_{aa^{\prime}}}{\sqrt{N_{c}^{2}-1}}$, which results in $\frac{\delta_{aa^{\prime}}}{\sqrt{N_{c}^{2}-1}}\frac{\delta_{aa^{\prime}}}{\sqrt{N_{c}^{2}-1}}=1$ after summing over color states. Thus for color factor we always have $\delta_{cc^{\prime}}$.

The angular momentum part of the amplitude gives

[TABLE]

where $(2\pi)^{4}\delta^{4}(p_{Q}-p_{Q}^{\prime})$ presents because the diagram is disconnected and the terms inside of “Tr” project both initial $Q{\overline{Q}}$ and final $Q{\overline{Q}}$ to specific spin states. By setting $P=P^{\prime}$ and using the delta function to integrate out $q$, we get

[TABLE]

where we have used $\int d^{2}\Omega^{\prime}Y_{L}^{L_{z}}(\theta^{\prime},\phi^{\prime})Y_{L^{\prime}}^{*L_{z}^{\prime}}(\theta^{\prime},\phi^{\prime})=\delta_{LL^{\prime}}\delta_{L_{z}L_{z}^{\prime}}$ in the last step. Due to

[TABLE]

the last line of Eq. (44) gives

[TABLE]

where we have used the fact that $P_{H}=P$ at this order and thus $P\cdot\epsilon_{S_{z}}=0$. Therefore,

[TABLE]

Eventually, we have

[TABLE]