# Novel pentaquark picture for singly heavy baryons from chiral symmetry

Daiki Suenaga<EMAIL_ADDRESS>Research Center for Nuclear Physics,
Osaka University, Ibaraki, 567-0048, Japan  Atsushi Hosaka
<EMAIL_ADDRESS>Research Center for Nuclear Physics, Osaka
University, Ibaraki, 567-0048, Japan  Advanced Science Research Center, Japan
Atomic Energy Agency (JAEA), Tokai 319-1195, Japan

###### Abstract

We propose a new type of structure for singly heavy baryons of $Qqq\bar{q}q$
in addition to the conventional one of $Qqq$. Based on chiral symmetry of the
light quarks, we show that the $Qqq\bar{q}q$ baryon offers a novel picture for
heavy quark spin-singlet and flavor-antisymmetric baryons. By making use of
the effective Lagrangian approach, we find $\Lambda_{c}(2765)$ and
$\Xi_{c}(2967)$ are mostly $Qqq\bar{q}q$ while $\Lambda_{c}(2286)$ and
$\Xi_{c}(2470)$ are mostly $Qqq$. The masses of negative-parity baryons are
predicted. We also derive a sum rule and the extended Goldberger-Treiman
relation that the masses of the baryons satisfy. Furthermore, a mass
degeneracy of parity partners of the baryons with the restoration of chiral
symmetry is discussed. These are the unique features that the conventional
picture of radial excitation in the quark model does not accommodate. Our
findings provide useful information not only for future experiments but also
for future lattice simulations on diquarks.

##  I Introduction

Investigation of heavy baryons has been attracting great attention with recent
development of experimental observation at, e.g., KEK, LHC, and SLAC. Singly
heavy baryons contain only one heavy quark, the mass of which is larger than
the typical energy scale of Quantum Chromodynamics (QCD). There the heavy
quark behaves not only as a color source but also as a spectator for the
remaining light quarks governed by nonperturbative QCD Manohar and Wise
(2000). Therefore, singly heavy baryons provide useful testing ground toward
elucidation of colorful light-quark dynamics inside hadrons.

In this regard, understanding of hadrons from chiral symmetry is important,
since it is one of the fundamental symmetries of QCD. In fact, the dynamics of
light pseudo-scalar mesons and nucleons are systematically formulated by the
spontaneous breakdown of chiral symmetry. In addition, the expected mass
degeneracy of chiral partners of light hadrons at extreme conditions provides
a key signal of the restoration of chiral symmetry for experiments and lattice
simulations Hatsuda and Kunihiro (1994); Harada and Yamawaki (2003). In this
paper we discuss singly heavy baryons from the chiral symmetry point of view
which cannot be easily implemented in the quark model Copley  _et al._ (1979);
Yoshida  _et al._ (2015). In particular, we investigate unique features that
are derived from the linear representations of chiral symmetry.

For baryons of light flavors, the mirror representation has been proposed for
the negative parity nucleon $N^{*}(1535)$, referred to as the mirror nucleon.
This picture contrasts with the quark model description of orbital excitation
for $N^{*}(1535)$. In the chiral representation approach, not only
spectroscopies but also the properties at the extreme conditions such as
change in masses and a degeneracy reflecting the chiral partner structurecan
be studied Detar and Kunihiro (1989); Jido  _et al._ (2001); Gallas  _et al._
(2010); Yamazaki and Harada (2019). The latter has been strongly suggested by
the lattice simulations Detar and Kogut (1987); Aarts  _et al._ (2017). The
quark content of the mirror nucleon is understood as an exotic pentaquark
state ($qqq\bar{q}q$).

Figure 1: A sketch of conventional picture (left) and novel picture proposed
in this paper (right) for describing singly heavy baryons.

Turning to the heavy hadrons, the linear representations of chiral symmetry
have been studied in various contexts Nowak  _et al._ (1993); Bardeen and Hill
(1994); Dmitrasinovic (2012); Ma and Harada (2018); Kawakami and Harada
(2019); Dmitrašinović and Chen (2020). Recently $\Xi_{c}(2967)$ was studied
experimentally and its spin and parity were determined to be
$J^{P}=\frac{1}{2}^{+}$ Moon  _et al._ (2021). The properties of
$\Lambda_{c}(2765)$ are similar to $\Xi_{c}(2967)$, so that the spin and
parity of $\Lambda_{c}(2765)$ are likely to be $J^{P}=\frac{1}{2}^{+}$ Arifi
_et al._ (2020a). We can safely assume that both the ground state
$\Lambda_{c}(2286)$ [$\Xi_{c}(2470)$] and the excited state
$\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] are heavy quark spin-singlet and flavor-
antisymmetric, since no spin partner baryons have been observed Zyla  _et al._
(2020), i.e., the total angular momentum of the light quarks in these heavy
baryons is zero: $j_{l}=0$. Recently observed $\Lambda_{b}(6072)$ baryon also
seems to share the same properties Aaij  _et al._ (2020). Within the
conventional quark model, $\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] can be
understood as a radial excitation. In this paper, as an analogue of the mirror
nucleon for $N^{*}(1535)$ we propose a new description where these heavy
baryons are described by the mirror representation of chiral symmetry. As
anticipated from the previous paragraph, the mirror representation for heavy
baryons is realized by exotic multi-quark contents. Due to the presence of one
heavy quark, it is achieved by the mirror diquark, which is actually a
tetraquark $d^{\prime}\sim qqq\bar{q}$ as displayed in the right panel in Fig.
1. The significance of the mirror diquark is implied when we implement the
chiral representation for diquarks. Introduction of the mirror diquark enables
us to understand $\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] as well as
$\Lambda_{c}(2286)$ [$\Xi_{c}(2470)$] in a unified way based on chiral
symmetry.

The novel state shown in the right panel of Fig. 1 is distinct from the
hadronic molecule state Chen  _et al._ (2016); Guo  _et al._ (2018), but
should be regarded as a compact pentaquark one because $d^{\prime}$ is a
colorful object binding the remaining heavy quark strongly.

In the quark model Copley  _et al._ (1979); Yoshida  _et al._ (2015),
$\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] is considered as analogous states of the
Roper resonance $N(1440)$ Roper (1964) which have been studied as a radial
excitation Burkert and Roberts (2019).111In the quark model approach, the
properties of this Roper-like baryons $\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$]
such as its mass and decay width can be reproduced by employing the
$\,{}^{3}P_{0}$ model Chen  _et al._ (2017), or by including relativistic
corrections Arifi  _et al._ (2021). In contrast, in our approach by
introducing the mirror diquark with chiral symmetry we will find various
unique features that cannot be seen in the quark model. Namely
$\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] and $\Lambda_{c}(2286)$
[$\Xi_{c}(2470)$] are described as superpositions of the three-quark and
pentaquark states. Moreover general relations that the masses of the baryons
satisfy such as the sum rule and the extended Goldberger-Treiman relation are
derived. Furthermore, our approach gives a natural explanation for the strong
suppression of direct decay process of $\Lambda_{c}(2765)$ by the two pion
emission Arifi  _et al._ (2020a). In addition to those findings, our present
approach can provide useful information on hadron properties such as a mass
degeneracy of the parity partners Suenaga  _et al._ (2015); Sasaki (2014);
Suenaga  _et al._ (2017); Buchheim  _et al._ (2018); Ishikawa  _et al._
(2019a); Montaña  _et al._ (2020) at extreme conditions, e.g., finite
temperature and/or density where chiral symmetry is (partially) restored.

This paper is organized as follows. In Sec. II we present an effective
Lagrangian for the heavy baryons formed by the conventional and mirror
diquarks. In Sec. III we determine the parameters of our model and show the
resultant mass spectrum of the baryons. In Sec. IV we discuss the mass
degeneracy of the parity partners with the restoration of chiral symmetry. In
Sec. V we conclude our present work.

##  II Model

The interpolating fields of the conventional and mirror diquarks contributing
to $\Lambda_{c}(2286)$ [$\Xi_{c}(2470)$] and $\Lambda_{c}(2765)$
[$\Xi_{c}(2967)$] are

$\displaystyle(d_{R})^{a}_{i}$ $\displaystyle\sim$
$\displaystyle\epsilon_{ijk}\epsilon^{abc}(q_{R}^{T})^{b}_{j}C(q_{R})^{c}_{k}\
,$ $\displaystyle(d_{L})^{a}_{i}$ $\displaystyle\sim$
$\displaystyle\epsilon_{ijk}\epsilon^{abc}(q_{L}^{T})^{b}_{j}C(q_{L})^{c}_{k}\
,$ $\displaystyle(d^{\prime}_{R})_{i}^{a}$ $\displaystyle\sim$
$\displaystyle\epsilon_{jkl}\epsilon^{abc}(q_{R}^{T})^{b}_{k}C(q_{R})^{c}_{l}[(\bar{q}_{L})^{d}_{i}(q_{R})^{d}_{j}]\
,$ $\displaystyle(d^{\prime}_{L})_{i}^{a}$ $\displaystyle\sim$
$\displaystyle\epsilon_{jkl}\epsilon^{abc}(q_{L}^{T})^{b}_{k}C(q_{L})^{c}_{l}[(\bar{q}_{R})^{d}_{i}(q_{L})^{d}_{j}]\
,$ (1)

where the right- (left-)handed light quark $q_{R(L)}$ is defined by
$q_{R(L)}=\frac{1\pm\gamma_{5}}{2}q$ with $q=(u,d,s)^{T}$. The subscript
“$i,j,\cdots$” and the superscript “$a,b,\cdots$” indicate the chiral and
color indices, respectively. We have introduced scalar diquarks whose total
angular momentum is zero: $j_{l}=0$, here because $\Lambda_{c}(2286)$
[$\Xi_{c}(2470)$] and $\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] are heavy quark
spin-singlet. Equation (1) shows that $d_{R}$, $d_{L}$, $d_{R}^{\prime}$, and
$d_{L}^{\prime}$ belong to the chiral representation of

$\displaystyle d_{R}\sim({\bm{1}},\bar{\bm{3}})\ ,\ \
d_{L}\sim(\bar{\bm{3}},{\bm{1}})\ ,$ $\displaystyle
d^{\prime}_{R}\sim(\bar{\bm{3}},{\bm{1}})\ ,\ \
d^{\prime}_{L}\sim({\bm{1}},\bar{\bm{3}})\ .$ (2)

Therefore, singly heavy baryons formed by a heavy quark $Q$ and a diquark in
Eq. (1) transform as

$\displaystyle B_{R}\to B_{R}g^{\dagger}_{R}\ ,\ \ B_{L}\to
B_{L}g_{L}^{\dagger}\ ,$ $\displaystyle B^{\prime}_{R}\to
B^{\prime}_{R}g^{\dagger}_{L}\ ,\ \ B^{\prime}_{L}\to
B^{\prime}_{L}g_{R}^{\dagger}\ ,$ (3)

with $g_{R(L)}\in SU(3)_{R(L)}$, where $B_{R}$, $B_{L}$, $B^{\prime}_{R}$, and
$B_{L}^{\prime}$ are given by

$\displaystyle B_{R,i}\sim Q^{a}(d_{R})_{i}^{a}\ ,\ \ B_{L,i}\sim
Q^{a}(d_{L})_{i}^{a}\ ,$ $\displaystyle B^{\prime}_{R,i}\sim
Q^{a}(d^{\prime}_{R})_{i}^{a}\ ,\ \ B^{\prime}_{L,i}\sim
Q^{a}(d^{\prime}_{L})_{i}^{a}\ ,$ (4)

respectively. The interpolating fields in Eq. (4) imply that the heavy quark
only plays a role of a spectator for light-quark degrees of freedom.

Heavy baryon fields in Eq. (4) allow us to construct an $SU(3)_{L}\times
SU(3)_{R}$ symmetric effective Lagrangian within the heavy baryon effective
theory (HBET) as

$\displaystyle{\cal L}_{\rm eff}$ $\displaystyle=$
$\displaystyle\sum_{\chi=L,R}\Big{(}\bar{B}_{\chi}iv\cdot\partial{B}_{\chi}-\mu_{1}\bar{B}_{\chi}B_{\chi}$
(5) $\displaystyle+$ $\displaystyle\bar{B}^{\prime}_{\chi}iv\cdot\partial
B^{\prime}_{\chi}-\mu_{2}\bar{B}^{\prime}_{\chi}B^{\prime}_{\chi}\Big{)}$
$\displaystyle-$
$\displaystyle\mu_{3}(\bar{B}_{R}B^{\prime}_{L}+\bar{B}^{\prime}_{L}B_{R}+\bar{B}_{L}B^{\prime}_{R}+\bar{B}^{\prime}_{R}{B}_{L})$
$\displaystyle-$ $\displaystyle
g_{1}(\bar{B}_{L}\Sigma^{*}B_{R}+\bar{B}_{R}\Sigma^{T}B_{L})$ $\displaystyle-$
$\displaystyle
g_{2}(\bar{B}^{\prime}_{R}\Sigma^{*}B^{\prime}_{L}+\bar{B}^{\prime}_{L}\Sigma^{T}B_{R})$
$\displaystyle-$ $\displaystyle
g_{3}(\bar{B}^{\prime}_{R}\Sigma^{*}B_{R}+\bar{B}_{L}\Sigma^{*}B^{\prime}_{L}+{\rm
h.c.})\ .$

Here $\Sigma=S+iP$ consists of scalar ($S$) and pseudo-scalar ($P$) meson
nonets and transforms as $\Sigma\to g_{L}\Sigma g_{R}^{\dagger}$. Also,
$v^{\mu}$ is a velocity of the heavy baryons. In Lagrangian (5), $\mu_{1}$,
$\mu_{2}$, and $\mu_{3}$ terms are responsible for the part of heavy baryon
mass which are of order $\Lambda_{\rm QCD}$, since a common mass parameter
$m_{B}$ for all the heavy baryons is subtracted for defining the HBET. The
remaining $g_{1}$, $g_{2}$, and $g_{3}$ terms are for the interactions between
light quarks and light mesons. These terms also contribute to the masses of
heavy baryons when chiral symmetry is spontaneously broken. The terms
proportional to $\mu_{3}$, $g_{1}$ and $g_{2}$ in Eq. (5) break the $U(1)_{A}$
axial symmetry allowed by the quantum anomaly. We note that the mass term of
$\bar{B}^{(\prime)}_{R}B_{R}^{(\prime)}+\bar{B}^{(\prime)}_{L}B_{L}^{(\prime)}$
is allowed in Eq. (5) because the property of spinor for the baryons is
determined by the heavy quark which is irrelevant to the chiral representation
of the diquarks. The physical baryons as parity eigenstates of $\pm$ are given
by

$\displaystyle B_{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(B_{R}\mp B_{L})\ ,$ $\displaystyle
B_{\pm}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(B_{L}^{\prime}\mp B_{R}^{\prime})\ ,$ (6)

with $B^{(\prime)}_{R}\leftrightarrow-B^{(\prime)}_{L}$ under parity
transformation Harada  _et al._ (2020).

Under the spontaneous breakdown of chiral symmetry, $\Sigma$ acquirers its
vacuum expectation values (VEVs) as $\langle\Sigma\rangle={\rm
diag}(\sigma_{q},\sigma_{q},\sigma_{s})$. When we ignore the $u$ and $d$ quark
masses, $\sigma_{q}$ becomes identical to the pion decay constant. We will use
$\sigma_{q}=93$ MeV as one of inputs for model parameters, whereas
$\sigma_{s}$ is not determined by the decay constants due to corrections from
the non-negligible $s$ quark mass $m_{s}\sim 100$ MeV. By substituting the
parity eigenstates $B_{\pm}$ and $B_{\pm}^{\prime}$ into the Lagrangian (5),
and by reading off the mass terms together with the VEVs
$\langle\Sigma\rangle$, mass eigenvalues of the heavy baryons are obtained as

$\displaystyle M(B_{+,i}^{H/L})$ $\displaystyle=$ $\displaystyle
m_{B}+\frac{1}{2}\Big{[}m_{+,i}+m_{+,i}^{\prime}$
$\displaystyle\pm\sqrt{(m_{+,i}-m_{+,i}^{\prime})^{2}+4\tilde{m}^{2}_{+,i}}\
\Big{]}\ ,$ $\displaystyle M(B_{-,i}^{H/L})$ $\displaystyle=$ $\displaystyle
m_{B}+\frac{1}{2}\Big{[}m_{-,i}+m_{-,i}^{\prime}$ (7)
$\displaystyle\pm\sqrt{(m_{-,i}-m_{-,i}^{\prime})^{2}+4\tilde{m}^{2}_{-,i}}\
\Big{]}\ ,$

with

$\displaystyle m_{\pm,i}$ $\displaystyle=$ $\displaystyle\mu_{1}\mp
g_{1}\sigma_{i}\ ,$ $\displaystyle m_{\pm,i}^{\prime}$ $\displaystyle=$
$\displaystyle\mu_{2}\mp g_{2}\sigma_{i}\ ,$ $\displaystyle\tilde{m}_{\pm,i}$
$\displaystyle=$ $\displaystyle\mu_{3}\mp g_{3}\sigma_{i}\ .$ (8)

Here $i=1,2,3$ stands for the flavor index where $i=1,2$ corresponds to
$\Xi_{c}$’s while $i=3$ to $\Lambda_{c}$’s. In Eq. (7) we have introduced a
common heavy mass parameter $m_{B}$ to define the HBET Manohar and Wise
(2000). The choice of $m_{B}$ is arbitrary, thus we take the average mass of
all the heavy baryons which will be explicitly given latter. For the VEVs, we
have defined $\sigma_{i=1}=\sigma_{i=2}\equiv\sigma_{q}$ and
$\sigma_{i=3}\equiv\sigma_{s}$. The superscript $H$ $(L)$ represents the
higher (lower) mass eigenstate, i.e. $H$ ($L$) corresponds to “$+$” (“$-$”)
sign in front of the square root in right-hand-side. In obtaining Eq. (7), we
have defined mass eigenstates by

$\displaystyle\left(\begin{array}[]{c}B_{\pm,i}^{L}\\\ B_{\pm,i}^{H}\\\
\end{array}\right)=\left(\begin{array}[]{cc}\cos\theta_{B_{\pm,i}}&\sin\theta_{B_{\pm,i}}\\\
-\sin\theta_{B_{\pm,i}}&\cos\theta_{B_{\pm,i}}\\\
\end{array}\right)\left(\begin{array}[]{c}B_{\pm,i}\\\ B_{\pm,i}^{\prime}\\\
\end{array}\right)\ ,$ (15)

with the mixing angle satisfying
$\tan\theta_{B_{\pm,i}}=(2\tilde{m}_{\pm,i})/(m_{\pm,i}-m^{\prime}_{\pm,i})$.
Hence the physical states are superpositions of three-quark ($B_{\pm,i}$) and
pentaquark ($B_{\pm,i}^{\prime}$) states as depicted in Fig. 1. This is the
novel picture of singly heavy baryons proposed in this paper.

##  III Mass spectrum

To start with, we explain how the remaining parameters $\mu_{1}$, $\mu_{2}$,
$\mu_{3}$, $g_{1}$, $g_{2}$, $g_{3}$, and $\sigma_{s}$ are fixed. First, we
employ masses of experimentally observed states of $J^{P}=\frac{1}{2}^{+}$;
$\Lambda_{c}(2286)$, $\Lambda_{c}(2765)$, $\Xi_{c}(2470)$, and $\Xi_{c}(2967)$
are used as inputs for Zyla  _et al._ (2020)

$\displaystyle M(B_{+,i=3}^{L})=2286\,{\rm MeV}\ ,\
M(B_{+,i=3}^{H})=2765\,{\rm MeV}\ ,$ $\displaystyle
M(B_{+,i=1,2}^{L})=2470\,{\rm MeV}\ ,\ M(B_{+,i=1,2}^{H})=2967\,{\rm MeV}\ .$

Next, we use the masses of the conventional diquarks which are considered to
be decoupled from the mirror ones. These masses were estimated by the lattice
QCD Bi  _et al._ (2016) as

$\displaystyle M(d_{+,i=1,2})=906\,{\rm MeV}\ ,\ M(d_{-,i=1,2})=1142\,{\rm
MeV}\ ,$ $\displaystyle M(d_{+,i=3})=725\,{\rm MeV}\ ,\
M(d_{-,i=3})=1265\,{\rm MeV}\ .$ (17)

It should be noted that the last one in Eq. (17) is estimated by a chiral
effective theory Harada  _et al._ (2020) together with the simulation Bi  _et
al._ (2016). Here “$\pm$” and “$i$” stand for the parity and flavor indices of
the diquark, respectively. These diquark masses are then related to the masses
of baryons of three-quark states $B_{\pm,i}$ with the pentaquark states
switched off: $B^{\prime}_{\pm,i}=0$ in the Lagrangian (5), leading to
$M(B_{\pm,i})=m_{B}+m_{\pm,i}$. The mass difference of the heavy baryons may
equal that of diquarks, i.e.

$\displaystyle M(B_{-,i})-M(B_{+,i})=2g_{1}\sigma_{i}=M(d_{-,i})-M(d_{+,i})\
.$ (18)

which can be used to fix $g_{1}$ and $\sigma_{s}$.

In addition to the above inputs, we employ a quark model prediction as a
useful reference for the demonstration. Namely, we take the mass of the
lightest heavy quark spin-singlet $\Lambda_{c}$ baryon carrying
$J^{P}=\frac{1}{2}^{-}$ predicted in Ref. Yoshida  _et al._ (2015) as another
input: $M(B_{-,i=3}^{L})=2890$ MeV. We note that the observed
$J^{P}=\frac{1}{2}^{-}$ baryons of $\Lambda_{c}(2595)$ and $\Xi_{c}(2790)$ are
heavy quark spin-doublet which are not treated in this paper.

Figure 2: Mass spectrum of $\Lambda_{c}$’s and $\Xi_{c}$’s with $J^{P}=\frac{1}{2}^{+},\frac{1}{2}^{-}$ obtained in our model. The asterisk ($*$) stands for the inputs. The ratios indicated under the bars correspond to the components of three-quark and pentaquark states in each baryon: $Qqq:Qqq\bar{q}q$, in which the upper and lower ratios correspond to the parameter sets (I) and (II), respectively. |  $\mu_{1}$ [MeV] |  $\mu_{2}$ [MeV] |  $\mu_{3}$ [MeV] | $g_{1}$ | $g_{2}$ | $g_{3}$ |  $\sigma_{q}$ [MeV] |  $\sigma_{s}$ [MeV]
---|---|---|---|---|---|---|---|---
Set (I) | $-247$ | $247$ | $\mp 91.0$ | $1.27$ | 1.94 | $\pm 0.34$ | $93^{*}$ | 212
Set (II) | $94.1$ | $-94.1$ | $\pm 246$ | $1.27$ | 1.94 | $\pm 0.34$ | $93^{*}$ | 212
Table 1: Two parameter sets (I) and (II). The mass parameter $m_{B}$ is fixed to be $m_{B}=2868$ MeV for both the sets. The asterisk ($\ast$) stands for the inputs.  | $\theta_{B_{+,i=1,2}}$ | $\theta_{B_{-,i=1,2}}$ | $\theta_{B_{+,i=3}}$ | $\theta_{B_{-,i=3}}$
---|---|---|---|---
Set (I) | $\pm 14.8^{\circ}$ | $\pm 6.01^{\circ}$ | $\pm 21.4^{\circ}$ | $\pm 1.67^{\circ}$
Set (II) | $\pm 29.9^{\circ}$ | $\pm 38.6^{\circ}$ | $\pm 23.2^{\circ}$ | $\pm 43.0^{\circ}$
Table 2: The fixed mixing angles for parameter sets (I) and (II), with the
sign corresponding to the ones in the text. Note that $\theta_{B_{\pm,i=3}}$
and $\theta_{B_{\pm,i=1,2}}$ stand for the mixing angles of $\Lambda_{c}$
baryons and those of $\Xi_{c}$ baryons, respectively.

Now we can fix all the parameters. Since the mass eigenvalues in Eq. (7)
include square roots, the coupled equations yield four solutions for the
parameter sets. Physically these four sets are classified into two sets (I)
and (II) as displayed in Table 1. Both the parameter sets (I) and (II) predict
the remaining negative-parity baryon masses as $M(B_{-,i=1,2}^{L})=2732$ MeV,
$M(B_{-,i=1,2}^{H})=3302$ MeV, and $M(B_{-,i=3}^{H})=3529$ MeV. We show the
resultant mass spectrum of the heavy baryons in Fig. 2 where the asterisk
($*$) stands for the inputs. Also, the mixing angles defined in Eq. (15) are
determined as shown in Table 2, where $\theta_{B_{\pm,i=3}}$ and
$\theta_{B_{\pm,i=1,2}}$ stand for the mixing angles of $\Lambda_{c}$’s and
those of $\Xi_{c}$’s, respectively. Table 2 indicates that, for
$J^{P}=\frac{1}{2}^{+}$ baryons the lower state $\Lambda_{c}(2286)$
[$\Xi_{c}(2470)$] is dominated by $Qqq$ while the higher one
$\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] by $Qqq\bar{q}q$ for both the parameter
sets. On the other hand, for $J^{P}=\frac{1}{2}^{-}$ baryons the ratio is
largely dependent on the parameters. The ratio of the $Qqq$ and $Qqq\bar{q}q$
components for each baryon is also shown under the bars in Fig. 2, in which
the upper and lower ratios correspond to the parameter sets (I) and (II),
respectively. The ratio is estimated by the square of each coefficient in Eq.
(15), i.e. $\cos^{2}\theta_{B_{\pm,i}}$ or $\sin^{2}\theta_{B_{\pm,i}}$.

---
Figure 3: The left panel shows the $\sigma_{s}$ dependence of
$M(B_{+,i=3}^{L})$, $M(B_{-,i=3}^{L})$, $M(B_{+,i=3}^{H})$, and
$M(B_{-,i=3}^{H})$, while the right one shows the $\sigma_{q}$ dependence of
$M(B_{+,i=1,2}^{L})$, $M(B_{-,i=1,2}^{L})$, $M(B_{+,i=1,2}^{H})$, and
$M(B_{-,i=1,2}^{H})$. The results are identical for the parameter sets (I) and
(II) in Table 1.

The $\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] baryon have been found to be mainly
$Qqq\bar{q}q$ with $J^{P}=\frac{1}{2}^{+}$ where the $\bar{q}q$ constituent in
the mirror diquark $qq\bar{q}q$ requires a $P$-wave excitation. Although
adding such a $P$-wave $\bar{q}q$ pair costs an energy of order 1 GeV in the
quark model, the above analysis shows that the mirror diquark based on chiral
symmetry costs only about 0.5 GeV to form $\Lambda_{c}(2765)$
[$\Xi_{c}(2967)$]. This finding is similar to the small mass of the light
scalar meson $\bar{q}q$ in the chiral model.

One of the most important consequences of our model is the sum rule of the
masses:

$\displaystyle\sum_{p=\pm,n=H,L}M(B_{p,i=1,2}^{n})=\sum_{p=\pm,n=H,L}M(B_{p,i=3}^{n})\
,$ (19)

which can be derived from Eq. (7). Namely the sum of mass of the four
$\Lambda_{c}$’s coincides with that of the four $\Xi_{c}$’s. The mass formula
(7) also yields the extended Goldberger-Treiman relation

$\displaystyle\sum_{n=H,L}M(B_{-,i}^{n})-\sum_{n=H,L}M(B_{+,i}^{n})=2(g_{1}+g_{2})\sigma_{i}$
(20)

for each $i$, which gives a constraint on the mass difference between parity
partners and the coupling constant of one pion (kaon) emission. In deriving
Eq. (20) the higher and lower masses $M^{H}_{p,i}$ and $M^{L}_{p,i}$ have been
summed up to cancel out the square roots in Eq. (7).

At the end of this section we give comments on decay properties of the excited
baryons in our present model. In this paper, we have investigated mostly the
mass spectrum of the baryons as the first step towards understanding of the
importance of the newly introduced pentaquark from chiral symmetry. In
addition to the masses, our Lagrangian (5) can derive various couplings among
the light mesons and heavy baryons. Namely, decay properties of the excited
baryons can be also studied. In this case, the ground-state $\Sigma_{c}$
baryons are necessary since the excited baryons can decay into them. Based on
chiral symmetry, the $\Sigma_{c}$ baryons are straightforwardly included on
top of our present model Harada  _et al._ (2020).222We expect that inclusion
of a pentaquark component for the ground-state $\Sigma_{c}$ baryons is not
necessary since they are not excited states. By comparing our calculation
given by such a hybrid model and the experimental data of decays, e.g., the
$\Lambda_{c}(2765)$ decay widths, some of the model parameters are expected to
be fixed or constrained such that the uncertainty of our present model is
narrowed. Moreover, to shed light on the decay properties of the excited
baryons predicted in this paper is useful for future heavy baryon experiments
to observe them. In addition, detailed understanding of decay properties of
the excited baryons, especially from the viewpoint of the Goldberger-Treiman
relation (20) and $U(1)_{A}$ anomaly Kawakami  _et al._ (2020), leads to
further elucidation of properties of exotic constituents from symmetry aspects
of QCD. Investigation of the decays is beyond the scope of this paper and we
leave it for future publication.

##  IV Baryon masses with the restoration of chiral symmetry

The symmetry relations such as the sum rule (19) and the extended Goldberger-
Treiman relation (20) provide useful information on the chiral symmetry
properties of the baryons for future experiments and lattice simulations. In
addition, as one of the most important advantages of employing the present
chiral effective model, we can examine the properties of baryons at extreme
conditions such as finite temperature/density where chiral symmetry tends to
be restored. The mirror diquark can be regarded as an analogue of the mirror
nucleon Detar and Kunihiro (1989); Jido  _et al._ (2001); Gallas  _et al._
(2010); Yamazaki and Harada (2019), which gives rise to the mass degeneracy of
parity partners at the chiral restoration point in the nucleon sector Detar
and Kogut (1987); Zschiesche  _et al._ (2007); Motohiro  _et al._ (2015);
Aarts  _et al._ (2017); Suenaga (2018); Ishikawa  _et al._ (2019b). Thus, we
can expect that a similar mass degeneracy of the partners of the singly heavy
baryons arises.

In order to see the above mass degeneracy, we study mass modifications of the
baryons by changing the VEVs $\sigma_{q}$ and $\sigma_{s}$ in the mass
formulae (7). In Fig. 3 we show the VEV dependence of the baryon masses. The
left panel of this figure shows the $\sigma_{s}$ dependence of
$M(B_{+,i=3}^{L})$, $M(B_{-,i=3}^{L})$, $M(B_{+,i=3}^{H})$, and
$M(B_{-,i=3}^{H})$, while the right one shows the $\sigma_{q}$ dependence of
$M(B_{+,i=1,2}^{L})$, $M(B_{-,i=1,2}^{L})$, $M(B_{+,i=1,2}^{H})$, and
$M(B_{-,i=1,2}^{H})$. It should be noted that the results are identical for
the parameter sets (I) and (II) in Table 1. Figure 3 shows that

$\displaystyle M(B_{+,i}^{L})=M(B_{-,}^{L})\ ,$ $\displaystyle
M(B_{+,i}^{H})=M(B_{-,i}^{H})$ (21)

hold for $i=1,2$ and $i=3$, respectively, at the chiral restoration point
denoted by $\sigma_{q}=0$ and $\sigma_{s}=0$. Therefore, our model predicts
that $B_{+,i}^{L}$ and $B_{-,i}^{L}$ are the parity partners as in the nucleon
sector. Similarly, $B_{+,i}^{H}$ and $B_{-,i}^{H}$ are the partners. Even when
chiral symmetry is restored, due to the presence of $\mu_{3}$ which is chiral
invariant, still the baryons are interpreted as superpositions of a three-
quark state and a pentaquark state. The role of the mass parameter $\mu_{3}$
is similar to the one of the so-called chiral invariant mass $M_{0}$ for the
naive and mirror nucleons Detar and Kunihiro (1989); Jido  _et al._ (2001).
Namely, $\mu_{3}$ and $M_{0}$ are expected to share common fundamental
properties.

The mass degeneracy demonstrated above originating from the restoration of
chiral symmetry is expected to be realized at extreme conditions such as
finite temperature/density. These environments are provided in heavy-ion
collisions (HICs) and lattice simulations. Further studies on the properties
of the baryons at such extreme conditions from the viewpoint of the
restoration of chiral symmetry are left for future work.

##  V Conclusions

In this paper, we have proposed a new type of the light-quark degrees of
freedom named mirror diquark in addition to the conventional diquark, to
explain $\Lambda_{c}(2765)$ [$\Xi_{c}(2970)$] and $\Lambda_{c}(2286)$
[$\Xi_{c}(2470)$] in a unified way based on chiral symmetry. Accordingly, we
have obtained the masses of negative-parity as well as positive-parity baryons
in Eq. (7), and furthermore have derived the unique relations such as the sum
rule (19) and the extended Goldberger-Treiman relation (20). Moreover our
model can naturally explain the strong suppression of the direct decay process
of $\Lambda_{c}(2765)$ by the two pion emission Arifi  _et al._ (2020a),
because the $\sigma$-$\Lambda_{c}(2765)$-$\Lambda_{c}(2286)$ ($\sigma$ is the
light scalar meson) coupling disappears after the diagonalization in the
Lagrangian (5).

The mass spectrum of singly heavy baryons predicted in this paper will provide
useful information for future experiments including HICs. Moreover, our
finding that the mirror diquark plays a significant role is expected to
provide new direction for future lattice simulations on diquarks Hess  _et
al._ (1998); Alexandrou  _et al._ (2006); Babich  _et al._ (2007); Bi  _et
al._ (2016). For example, the importance of the mirror diquark would be tested
by examining the correlators of the conventional $qq$ and mirror $qq\bar{q}q$
diquarks, similarly to the simulation of light scalar mesons with $\bar{q}q$
and $\bar{q}q\bar{q}q$ states Wakayama  _et al._ (2015). The development of
the diquark chiral effective theory Hong  _et al._ (2004) including the mirror
diquark would be of interest. We also expect that the mirror diquark with
chiral symmetry provides a new aspect for the understanding of the Roper
resonance $N(1440)$ with the unique relations similar to Eqs. (19) and (20).

In our novel picture, the mirror diquark is the main constituent of
$\Lambda_{c}(2765)$ [$\Xi_{c}(2970)$], while conventionally such a baryon is
regarded as a radial excitation ($2S$ state) in the quark model Arifi  _et
al._ (2020b). Namely, an additional mixing from $2S$ state to the pentaquark
for $\Lambda_{c}(2765)$ [$\Xi_{c}(2970)$] is expected 333In addition, another
constituent provided by bi-local interpolating fields Dmitrasinovic and Chen
(2011); Chen and Dmitrašinović (2013) may enter.. Therefore, checking the
unique relations Eqs. (19) and (20), and the mass degeneracy of parity
partners or its precursory behaviors in future experiments or lattice
simulations would be desired for the better understanding of the hadrons from
chiral symmetry.

## Acknowledgement

We thank Ahmad Jafar Arifi for fruitful discussions and comments. We also
thank Veljko Dmitrasinovic for useful comments. A. H. is supported in part by
Grants-in Aid for Scientific Research, Grants No. 17K05441(C) and by Grants-in
Aid for Scientific Research on Innovative Areas (No. 18H05407).

## References

  * Manohar and Wise (2000) A. V. Manohar and M. B. Wise,  _Heavy Quark Physics_, Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology (Cambridge University Press, 2000).
  * Hatsuda and Kunihiro (1994) T. Hatsuda and T. Kunihiro, Phys. Rept. 247, 221 (1994), arXiv:hep-ph/9401310 [hep-ph] .
  * Harada and Yamawaki (2003) M. Harada and K. Yamawaki, Phys. Rept. 381, 1 (2003), arXiv:hep-ph/0302103 .
  * Copley  _et al._ (1979) L. A. Copley, N. Isgur, and G. Karl, Phys. Rev. D 20, 768 (1979), [Erratum: Phys.Rev.D 23, 817 (1981)].
  * Yoshida  _et al._ (2015) T. Yoshida, E. Hiyama, A. Hosaka, M. Oka, and K. Sadato, Phys. Rev. D 92, 114029 (2015), arXiv:1510.01067 [hep-ph] .
  * Detar and Kunihiro (1989) C. E. Detar and T. Kunihiro, Phys. Rev. D 39, 2805 (1989).
  * Jido  _et al._ (2001) D. Jido, M. Oka, and A. Hosaka, Prog. Theor. Phys. 106, 873 (2001), arXiv:hep-ph/0110005 .
  * Gallas  _et al._ (2010) S. Gallas, F. Giacosa, and D. H. Rischke, Phys. Rev. D 82, 014004 (2010), arXiv:0907.5084 [hep-ph] .
  * Yamazaki and Harada (2019) T. Yamazaki and M. Harada, Phys. Rev. D 99, 034012 (2019), arXiv:1809.02359 [hep-ph] .
  * Detar and Kogut (1987) C. E. Detar and J. B. Kogut, Phys. Rev. Lett. 59, 399 (1987).
  * Aarts  _et al._ (2017) G. Aarts, C. Allton, D. De Boni, S. Hands, B. Jäger, C. Praki, and J.-I. Skullerud, JHEP 06, 034 (2017), arXiv:1703.09246 [hep-lat] .
  * Nowak  _et al._ (1993) M. A. Nowak, M. Rho, and I. Zahed, Phys. Rev. D 48, 4370 (1993), arXiv:hep-ph/9209272 .
  * Bardeen and Hill (1994) W. A. Bardeen and C. T. Hill, Phys. Rev. D 49, 409 (1994), arXiv:hep-ph/9304265 .
  * Dmitrasinovic (2012) V. Dmitrasinovic, Phys. Rev. D 86, 016006 (2012).
  * Ma and Harada (2018) Y.-L. Ma and M. Harada, J. Phys. G 45, 075006 (2018), arXiv:1709.09746 [hep-ph] .
  * Kawakami and Harada (2019) Y. Kawakami and M. Harada, Phys. Rev. D 99, 094016 (2019), arXiv:1902.06774 [hep-ph] .
  * Dmitrašinović and Chen (2020) V. Dmitrašinović and H.-X. Chen, Phys. Rev. D 101, 114016 (2020).
  * Moon  _et al._ (2021) T. J. Moon _et al._ (Belle), Phys. Rev. D 103, L111101 (2021), arXiv:2007.14700 [hep-ex] .
  * Arifi  _et al._ (2020a) A. Arifi, H. Nagahiro, A. Hosaka, and K. Tanida, Phys. Rev. D 101, 111502 (2020a), arXiv:2004.07423 [hep-ph] .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Aaij  _et al._ (2020) R. Aaij _et al._ (LHCb), JHEP 06, 136 (2020), arXiv:2002.05112 [hep-ex] .
  * Chen  _et al._ (2016) H.-X. Chen, W. Chen, X. Liu, and S.-L. Zhu, Phys. Rept. 639, 1 (2016), arXiv:1601.02092 [hep-ph] .
  * Guo  _et al._ (2018) F.-K. Guo, C. Hanhart, U.-G. Meißner, Q. Wang, Q. Zhao, and B.-S. Zou, Rev. Mod. Phys. 90, 015004 (2018), arXiv:1705.00141 [hep-ph] .
  * Roper (1964) L. Roper, Phys. Rev. Lett. 12, 340 (1964).
  * Burkert and Roberts (2019) V. D. Burkert and C. D. Roberts, Rev. Mod. Phys. 91, 011003 (2019), arXiv:1710.02549 [nucl-ex] .
  * Note (1) In the quark model approach, the properties of this Roper-like baryons $\Lambda_{c}(2765)$ [$\Xi_{c}(2967)$] such as its mass and decay width can be reproduced by employing the $\tmspace+{.1667em}^{3}P_{0}$ model Chen  _et al._ (2017), or by including relativistic corrections Arifi  _et al._ (2021).
  * Suenaga  _et al._ (2015) D. Suenaga, B.-R. He, Y.-L. Ma, and M. Harada, Phys. Rev. D 91, 036001 (2015), arXiv:1412.2462 [hep-ph] .
  * Sasaki (2014) C. Sasaki, Phys. Rev. D 90, 114007 (2014), arXiv:1409.3420 [hep-ph] .
  * Suenaga  _et al._ (2017) D. Suenaga, S. Yasui, and M. Harada, Phys. Rev. C 96, 015204 (2017), arXiv:1703.02762 [nucl-th] .
  * Buchheim  _et al._ (2018) T. Buchheim, T. Hilger, B. Kämpfer, and S. Leupold, J. Phys. G 45, 085104 (2018), arXiv:1801.01472 [nucl-th] .
  * Ishikawa  _et al._ (2019a) T. Ishikawa, K. Nakayama, D. Suenaga, and K. Suzuki, Phys. Rev. D 100, 034016 (2019a), arXiv:1905.11164 [hep-ph] .
  * Montaña  _et al._ (2020) G. Montaña, A. Ramos, L. Tolos, and J. M. Torres-Rincon, Phys. Lett. B 806, 135464 (2020), arXiv:2001.11877 [hep-ph] .
  * Harada  _et al._ (2020) M. Harada, Y.-R. Liu, M. Oka, and K. Suzuki, Phys. Rev. D 101, 054038 (2020), arXiv:1912.09659 [hep-ph] .
  * Bi  _et al._ (2016) Y. Bi, H. Cai, Y. Chen, M. Gong, Z. Liu, H.-X. Qiao, and Y.-B. Yang, Chin. Phys. C 40, 073106 (2016), arXiv:1510.07354 [hep-ph] .
  * Note (2) We expect that inclusion of a pentaquark component for the ground-state $\Sigma_{c}$ baryons is not necessary since they are not excited states.
  * Kawakami  _et al._ (2020) Y. Kawakami, M. Harada, M. Oka, and K. Suzuki, Phys. Rev. D 102, 114004 (2020), arXiv:2009.06243 [hep-ph] .
  * Zschiesche  _et al._ (2007) D. Zschiesche, L. Tolos, J. Schaffner-Bielich, and R. D. Pisarski, Phys. Rev. C 75, 055202 (2007), arXiv:nucl-th/0608044 .
  * Motohiro  _et al._ (2015) Y. Motohiro, Y. Kim, and M. Harada, Phys. Rev. C 92, 025201 (2015), [Erratum: Phys.Rev.C 95, 059903 (2017)], arXiv:1505.00988 [nucl-th] .
  * Suenaga (2018) D. Suenaga, Phys. Rev. C 97, 045203 (2018), arXiv:1704.03630 [nucl-th] .
  * Ishikawa  _et al._ (2019b) T. Ishikawa, K. Nakayama, and K. Suzuki, Phys. Rev. D 99, 054010 (2019b), arXiv:1812.10964 [hep-ph] .
  * Hess  _et al._ (1998) M. Hess, F. Karsch, E. Laermann, and I. Wetzorke, Phys. Rev. D 58, 111502 (1998), arXiv:hep-lat/9804023 .
  * Alexandrou  _et al._ (2006) C. Alexandrou, P. de Forcrand, and B. Lucini, Phys. Rev. Lett. 97, 222002 (2006), arXiv:hep-lat/0609004 .
  * Babich  _et al._ (2007) R. Babich, N. Garron, C. Hoelbling, J. Howard, L. Lellouch, and C. Rebbi, Phys. Rev. D 76, 074021 (2007), arXiv:hep-lat/0701023 .
  * Wakayama  _et al._ (2015) M. Wakayama, T. Kunihiro, S. Muroya, A. Nakamura, C. Nonaka, M. Sekiguchi, and H. Wada, Phys. Rev. D 91, 094508 (2015), arXiv:1412.3909 [hep-lat] .
  * Hong  _et al._ (2004) D. K. Hong, Y. J. Sohn, and I. Zahed, Phys. Lett. B 596, 191 (2004), arXiv:hep-ph/0403205 .
  * Arifi  _et al._ (2020b) A. Arifi, H. Nagahiro, A. Hosaka, and K. Tanida, Phys. Rev. D 101, 094023 (2020b), arXiv:2003.08202 [hep-ph] .
  * Note (3) In addition, another constituent provided by bi-local interpolating fields Dmitrasinovic and Chen (2011); Chen and Dmitrašinović (2013) may enter.
  * Chen  _et al._ (2017) B. Chen, K.-W. Wei, X. Liu, and T. Matsuki, Eur. Phys. J. C 77, 154 (2017), arXiv:1609.07967 [hep-ph] .
  * Arifi  _et al._ (2021) A. J. Arifi, D. Suenaga, and A. Hosaka, Phys. Rev. D 103, 094003 (2021), arXiv:2102.03754 [hep-ph] .
  * Dmitrasinovic and Chen (2011) V. Dmitrasinovic and H.-X. Chen, Eur. Phys. J. C 71, 1543 (2011), arXiv:1101.5906 [hep-ph] .
  * Chen and Dmitrašinović (2013) H.-X. Chen and V. Dmitrašinović, Phys. Rev. D 88, 036013 (2013), arXiv:1309.0387 [hep-ph] .


11institutetext:  Nuclear Science Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA 22institutetext: Physics Department,
University of California, Berkeley, CA 94720, USA 33institutetext: Center for
Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA 44institutetext: Theory Center, Jefferson Lab, Newport News,
Virginia 23606, USA 55institutetext: Department of Physics, The College of
William & Mary, Williamsburg, Virginia 23187, USA 66institutetext: C.N. Yang
Institute for Theoretical Physics and Department of Physics and Astronomy,
Stony Brook University, Stony Brook, NY 11794, USA 77institutetext: SUBATECH
UMR 6457 (IMT Atlantique, Université de Nantes, IN2P3/CNRS), 4 rue Alfred
Kastler, 44307 Nantes, France 88institutetext: Faculty of Science and
Technology, Seikei University, Musashino, Tokyo 180-8633, Japan

# Subleading power corrections to heavy quarkonium production in QCD
factorization approach

Kyle Lee 112233 Jian-Wei Qiu 4455 George Sterman 66 Kazuhiro Watanabe 7788
<EMAIL_ADDRESS>
###### Abstract

We report the current understanding of heavy quarkonium production at high
transverse momentum ($p_{T}$) in hadronic collisions in terms of QCD
factorization. In this presentation, we highlight the role of subleading power
corrections to heavy quarkonium production, which are essential to describe
the $p_{T}$ spectrum of quarkonium at a relatively lower $p_{T}$. We also
introduce prescription to match QCD factorization to fixed-order NRQCD
factorization calculations for quarkonium production at low $p_{T}$.

##  1 Introduction

Heavy quarkonium production allows us to study fundamental QCD dynamics about
how a quark and antiquark pair is produced and then fragments into a color
singlet hadron, one of the central issues discussed in this conference. Thus
far, many theoretical efforts have been made to push the accuracy of
calculations in NRQCD factorization Bodwin:1994jh  toward the NLO level and
beyond, and to apply the NRQCD (and pNRQCD) factorization approach to help pin
down the quarkonium production mechanism Brambilla:2022ayc . We have argued
that the lack of a systematic treatment of significantly enhanced
$\ln(p_{T}^{2}/m_{Q}^{2})$-type corrections to many NRQCD calculations could
affect data fitting of long-distance-matrix-elements (LDMEs) to describe the
shape of the $p_{T}$ spectrum and polarization of heavy quarkonium
simultaneously.

We studied the $p_{T}$ spectrum of heavy quarkonium production employing the
QCD factorization approach, which expands the production cross sections in
powers of $1/p_{T}$ first and factorizes both the leading power (LP) and next-
to-leading power (NLP) contributions in terms of perturbatively calculable
partonic scattering, convoluted with universal parton distribution functions
(PDFs) and fragmentation functions (FFs), respectively Lee:2021oqr . The
factorized LP cross-section has successfully described LHC data on $J/\psi$
production at high $p_{T}\gg m_{J/\psi}$ Lee:2021oqr . Here we show that the
NLP contribution is required to describe $J/\psi$ production data for
$p_{T}\gtrsim\mathcal{O}(m_{J/\psi})$, whose shape differs from the LP
distribution.

In this presentation, we highlight the role of the NLP corrections to hadronic
$J/\psi$ production by comparing our calculations with the LHC and Tevatron
data to demonstrate that the NLP contribution starts to dominate over the LP
one when $p_{T}\lesssim 4m_{J/\psi}$. We also will provide the justification
of an assumption, described below, that makes the QCD evolution equations of
the twist-4 double parton FFs more manageable. Finally, we briefly discuss how
to match the NLP cross-section to the fixed order calculation in NRQCD
factorization.

##  2 QCD factorization formalism for inclusive quarkonium production

Inclusive production of heavy quarkonium $H$ at high $p_{T}$ ($\gg m_{H}$) in
hadronic collisions ($A+B$) can be studied in QCD factorization formalism
Nayak:2005rt ; Kang:2014tta :

$\displaystyle E_{p}\frac{d\sigma_{A+B\to H(p)+X}}{d^{3}p}$
$\displaystyle\approx E_{p}\frac{d\sigma_{A+B\to H(p)+X}}{d^{3}p}\bigg{|}_{\rm
LP}+E_{p}\frac{d\sigma_{A+B\to H(p)+X}}{d^{3}p}\bigg{|}_{\rm NLP},$ (1)
$\displaystyle E_{p}\frac{d\sigma_{A+B\to H(p)+X}}{d^{3}p}\bigg{|}_{\rm LP}$
$\displaystyle=\sum\int\frac{dz}{z^{2}}D_{f\to
H}(z,\mu^{2})E_{c}\frac{d\hat{\sigma}_{A+B\to
f(p_{c})+X}}{d^{3}p_{c}}\left(p_{c}=\frac{p}{z},\mu^{2}\right),$ (2)
$\displaystyle E_{P}\frac{d\sigma_{A+B\to H(p)+X}}{d^{3}p}\bigg{|}_{\rm NLP}$
$\displaystyle=\sum_{\kappa}\int\frac{dz}{z^{2}}D_{[Q\bar{Q}(\kappa)]\to
H}(z,\mu^{2})E_{c}\frac{d\hat{\sigma}_{A+B\to[Q\bar{Q}(\kappa)](p_{c})+X}}{d^{3}p_{c}}\left(p_{c}=\frac{p}{z},\mu^{2}\right),$
(3)

where the second and third lines, respectively, represent the factorization
formalism at LP and NLP. $d\hat{\sigma}_{A+B\to f(p_{c})+X}$ denotes the
partonic cross section for producing the fragmenting parton of flavor $f$ and
momentum $p_{c}$ with all collinear sensitivities around $p_{c}\sim p/z$
absorbed into the twist-2 $f$-to-$H$ FFs, $D_{f\to H}$ with momentum fraction
$z$. Similarly, $d\hat{\sigma}_{A+B\to[Q\bar{Q}(\kappa)](p_{c})+X}$ denotes
the partonic cross section to produce the fragmenting $Q\bar{Q}$ pair of spin-
color state $\kappa$ and momentum
$p_{c}=P_{Q}+P_{\bar{Q}}=P^{\prime}_{Q}+P^{\prime}_{\bar{Q}}$, where
$P^{\prime}_{Q}$ and $P^{\prime}_{\bar{Q}}$ are momenta in the conjugated
production amplitude. All collinear sensitivities around $p_{c}$ in the cross-
section are absorbed into the twist-4 $Q\bar{Q}(\kappa)$-to-$H$ FFs,
$D_{[Q\bar{Q}(\kappa)]\to H}$.

Since physical cross sections should not depend on the factorization scale,
the twist-2 and twist-4 FFs satisfy the following coupled evolution equations
Kang:2014tta ,

$\displaystyle\frac{\partial}{\partial\ln\mu^{2}}D_{f\to H}(z,\mu^{2})$
$\displaystyle=\frac{\alpha_{s}(\mu)}{2\pi}\sum_{f^{\prime}}\int_{z}^{1}\frac{dz^{\prime}}{z^{\prime}}P_{f\to
f^{\prime}}\left(\frac{z}{z^{\prime}}\right)D_{f^{\prime}\to
H}(z^{\prime},\mu^{2})$
$\displaystyle+\frac{\alpha_{s}^{2}(\mu)}{\mu^{2}}\sum_{\kappa}\int_{z}^{1}\frac{dz^{\prime}}{z^{\prime}}P_{f\to[Q\bar{Q}(\kappa)]}\left(\frac{z}{z^{\prime}}\right)D_{[Q\bar{Q}(\kappa)]\to
H}\left(z^{\prime},\mu^{2}\right),$ (4)
$\displaystyle\frac{\partial}{\partial\ln\mu^{2}}D_{[Q\bar{Q}(\kappa)]\to
H}(z,\mu^{2})$
$\displaystyle=\frac{\alpha_{s}(\mu)}{2\pi}\sum_{n}\int^{1}_{z}\frac{dz^{\prime}}{z^{\prime}}P_{[Q\bar{Q}(n)]\to[Q\bar{Q}(\kappa)]}\left(\frac{z}{z^{\prime}}\right)\,D_{[Q\bar{Q}(n)]\to
H}(z^{\prime},\mu^{2}),$ (5)

allowing one to use the renormalization group improved QCD factorization at
the NLP accuracy. The first line of Eq.(4) is the standard DGLAP evolution of
the twist-2 FFs and the second line of Eq.(4) represents corrections to the
evolution from the situation when the fragmenting parton fragments to a
$Q\bar{Q}$ pair at the scale $\mu$ and then the pair fragments to the observed
quarkonium $H$. Eq.(5) represents the linear DGLAP-like evolution equations
for twist-4 $Q\bar{Q}(\kappa)$-to-$H$ FFs. When combined with the NLP partonic
cross sections, the solution of Eq.(5) gives the NLP contribution and is very
significant for quarkonium production at low $p_{T}$. The solution of Eq.(4),
with the twist-4 $Q\bar{Q}(\kappa)$-to-$H$ FFs from the solution of Eq.(5),
effectively resum logarithmically enhanced contributions to the cross section
when the fragmenting parton fragments into a heavy quark pair at a scale
between $[\mu_{0}\sim 2m_{Q},\mu\sim p_{T}]$.

The factorization formalism at LP in Eq.(2) should dominate the quarkonium
production when $p_{T}\gg m_{H}\approx 2m_{Q}$ when the LP
$\ln(p_{T}^{2}/m_{Q}^{2})$-enhanced contributions are significant, and the
power corrections in $1/p_{T}$ are negligible Nayak:2005rt . Although it is
formally suppressed by $1/p_{T}^{2}$ for producing a $Q\bar{Q}$ pair, compared
to the LP formalism, contributions from the NLP factorization formalism are
enhanced by the fact that a heavy $Q\bar{Q}$ pair is more likely than a single
parton to become a quarkonium $H$. The NLP contribution plays a crucial role
in describing quarkonium production cross section at relatively lower
$p_{T}>\mathcal{O}(m_{H})$.

In order to obtain Eqs.(4) and (5) from corresponding equations and functions
in Ref.Kang:2014tta , we introduced the modified double parton FFs as
$D_{[Q\bar{Q}(\kappa)]\to H}(z,\mu^{2})\equiv\int
dudv\,\mathcal{D}_{[Q\bar{Q}(\kappa)]\to H}(z,u,v,\mu^{2})$ where $u(v)$
denotes a heavy quark’s relative longitudinal momentum fraction in the
amplitude (conjugated amplitude), and modified the evolution kernels involving
the heavy quark pair as $P_{[Q\bar{Q}(n)]\to[Q\bar{Q}(\kappa)]}(z)\equiv\int
du\int
dv\,\mathcal{P}_{[Q\bar{Q}(n)]\to[Q\bar{Q}(\kappa)]}\left(z,u,v,u^{\prime}=\frac{1}{2},v^{\prime}=\frac{1}{2}\right)$,
$P_{f\to[Q\bar{Q}(\kappa)]}(z)\equiv\mathcal{P}_{f\to[Q\bar{Q}(\kappa)]}\left(z,u^{\prime}=\frac{1}{2},v^{\prime}=\frac{1}{2}\right)$,
where $u^{\prime}(v^{\prime})$ is the momentum fraction of a heavy quark in
the amplitude (conjugated amplitude) on the lower virtuality side. The above
modification assumes that the same momentum fractions of the heavy quark pair
($u=v=1/2$) are preferred for forming a physical quarkonium.

Next, we explore whether our assumption $u=v=1/2$ can be justified by
considering the derivative of double parton FFs with respect to $\mu^{2}$ (or
the slope of their $\mu^{2}$-evolution):

$\displaystyle\mathcal{D}^{\prime}_{\kappa\to n}(z,u,v,\mu^{2})$
$\displaystyle\equiv\frac{2\pi}{\alpha_{s}}\frac{d\mathcal{D}_{\kappa\to
n}(z,u,v,\mu^{2})}{d\ln\mu^{2}}$
$\displaystyle=\int_{z}^{1}\frac{dz^{\prime}}{z^{\prime}}\int_{0}^{1}du^{\prime}\int_{0}^{1}dv^{\prime}\mathcal{P}_{\kappa\to
n}\left(\frac{z}{z^{\prime}},u,v,u^{\prime},v^{\prime}\right)\mathcal{D}(z^{\prime},u^{\prime},v^{\prime},\mu^{2}),$
(6)

where the shape of the double parton FFs in $u,v,z$-space is modeled with a
simple test function $\mathcal{D}(z,u,v)$, assuming $\mathcal{D}(z,u,v)\to
D_{z}(z;\alpha,\beta)D_{u}(u;\gamma)D_{v}(v;\gamma)$ with

$\displaystyle
D_{z}(z;\alpha,\beta)=\frac{z^{\alpha}(1-z)^{\beta}}{B[1+\alpha,1+\beta]},\quad
D_{u,v}(x;\gamma)=\frac{x^{\gamma}(1-x)^{\gamma}}{B[1+\gamma,1+\gamma]},$ (7)

satisfying $\int_{0}^{1}dz\int_{0}^{1}du\int_{0}^{1}dv\mathcal{D}(z,u,v)=1$.
The input parameters $\alpha,\beta,\gamma$ in Eq.(7) control the shape of the
test function. Compared to light hadron production, such as $\pi$ meson, heavy
quarkonium production is more sensitive to the information of FFs at large
$z$, like parton-to-open heavy flavor FFs. We imitate such FFs by combining a
large $\alpha$ and a small $\beta$ for $D_{z}$. We fix $\alpha=30$,
$\beta=0.5$ as implemented in Ref.Lee:2021oqr . In our analysis, we take
$\gamma=10$ to have a relatively narrow peak around $u=v=1/2$ as input.

Figure 1:  Contour plots of $\mathcal{D}^{\prime}_{v1\to v1}$ (Left) and
$\mathcal{D}^{\prime}_{a8\to a8}$ (Right) for $z=0.85$ in $(u,v)$-space. The
sidebar on each panel represents the numerical value of the derivative in
arbitrary units; the lighter the sidebar’s color, the larger the derivative’s
value.

Figure 1 shows $\mathcal{D}^{\prime}_{\kappa\to n}(z,u,v)$ in Eq.(6) for the
diagonal channels, $\mathcal{D}^{\prime}_{v1\to v1}$ and
$\mathcal{D}^{\prime}_{v8\to v8}$. For the color-singlet to color-singlet
(S-S) diagonal transition, the derivative is negative around $u,v\sim 1/2$ and
positive at $u,v\ll 1/2$ and $u,v\gg 1/2$, indicating that the double parton
FFs can get broadened in $(u,v)$ space after much evolution. Consequently, our
assumption $u^{\prime}=v^{\prime}=1/2$ in the evolution kernels is not a good
approximation in this case. Nevertheless, since the production rate of
quarkonium from a heavy quark pair in the color singlet is much lower than the
color octet contribution at high $p_{T}$, we expect that setting
$u^{\prime}=v^{\prime}=1/2$ in the evolution kernels would be a good
approximation for phenomenological study of the NLP effects.

In contrast, for the color-octet to color-octet (O-O) diagonal channel, the FF
can become narrower with a prominent peak around $u=v=1/2$ after much
evolution because the derivative at large $z$ is positive over the whole range
in $(u,v)$ space, as shown in Fig.1. Therefore, our assumption
$u^{\prime}=v^{\prime}=1/2$ in the evolution kernels is a good approximation
in the large $z$ region. The off-diagonal O-S and S-O channels have the same
features as the O-O diagonal channels, so that we can safely approximate
$P_{Q}=P_{\bar{Q}}=P^{\prime}_{Q}=P^{\prime}_{\bar{Q}}=p/(2z)$ in the NLP
partonic cross section Kang:2014tta , which has been assumed and used in all
perturbative NRQCD calculations. More details of our analysis will be reported
in Ref.LQSW:2022 .

Figure 2: The single and double parton FFs for producing a heavy quark pair in
the intermediate state ${{}^{1}{S}_{0}^{[8]}}$ as functions of $z$, evolved
from $\mu_{0}^{2}=(4m_{c})^{2}$ with $m_{c}=1.5\,{\rm GeV}$ up to given values
of $\mu^{2}$: $\mu^{2}=10^{2}\,{\rm GeV}^{2}$ (left column),
$\mu^{2}=10^{4}\,{\rm GeV}^{2}$ (middle column), $\mu^{2}=10^{6}\,{\rm
GeV}^{2}$ (right column).
$\hat{D}_{\kappa}=D_{\kappa}/\langle\mathcal{O}_{\kappa}\rangle$ with
$\langle\mathcal{O}_{\kappa}\rangle$ the LDME for the production of a quark
pair in $\kappa$. The single parton FFs with the linear (thin) and nonlinear
(thick) evolution differ in line width in each panel. Black solid curves:
light quarks ($u,d,s$), Blue dashed curves: heavy quarks ($c$), Red dotted
curves: gluon.

##  3 Numerical results for hadronic $J/\psi$ production

We shall show numerical results of hadronic $J/\psi$ production cross section
at high $p_{T}$ in the QCD factorization approach with Eqs.(2)-(5). We use
CT18NLO set for colinear PDFs in the proton and antiproton Hou:2019efy .
Significant uncertainty in our approach is from the non-perturbative
information of the twist-2 and twist-4 FFs. As studied in Ref.Lee:2021oqr , we
could approximate these FFs by applying NRQCD factorization to obtain analytic
forms of the twist-2 and twist-4 FFs at the input scale,
$\mu_{0}=\mathcal{O}(2m_{Q})$, in terms of a finite set of NRQCD LDMEs. The
$z$-dependence of the input FFs can be calculated perturbatively in NRQCD in
an expansion of $\alpha_{s}$ and heavy quark velocity $v$ in the pair’s rest
frame and is available at the first non-trivial order in the powers of
$\alpha_{s}$ Ma:2013yla ; Ma:2014eja .

Although the non-perturbative input FFs at
$\mu_{0}=\mathcal{O}(2m_{Q})\gg\Lambda_{\rm QCD}$ can be approximated by non-
perturbative LDMEs with perturbatively calculated coefficients in terms of
NRQCD factorization, the perturbative short-distance coefficients inevitably
involve very large, even negative, distributions near $z=1$, such as
$\delta(1-z)$, $f(z)\ln(1-z)$, $f(z)/[1-z]_{+}$, and
$f(z)[\ln(1-z)/(1-z)]_{+}$ with $f$ arbitrary finite function of $z$, which
are naturally from the perturbative cancellation of IR divergences between
real and virtual diagrams. Since heavy quarkonium FFs have peaks at large $z$,
these negative terms could have a much more important impact on heavy
quarkonium production and drive the perturbatively evaluated cross-section to
be negative. Resummation near the threshold could help us cure this problem of
perturbative calculations. At the same time, we must keep in mind that the
uncertainty from the expansion of input FFs in terms of $\alpha_{s}$ and
relative velocity $v$, not to be confused with the relative momentum fraction
"$v$" above, remains a source of systematic error. We have then taken one
particular approach Lee:2021oqr , in which all the perturbatively calculated
input FFs are converted into
$N_{\kappa}\,z^{\alpha}(1-z)^{\beta}/B[1+\alpha,1+\beta]$, where $\alpha$ and
$\beta$ are free parameters, $B$ is the Euler Beta-function, and $N_{\kappa}$
is equal to the first moment of the related term that takes into account the
relative size of different terms from perturbative calculations Ma:2013yla ;
Ma:2014eja . In our numerical calculations, we have set $\alpha=30$,
$\beta=0.5$ for the input FFs at $\mu_{0}=6\,{\rm GeV}$. We show in Fig.2 the
$z$-dependence of the single parton and double parton FFs for producing a
heavy quark pair of ${{}^{1}S_{0}^{[8]}}$ at different values of $\mu^{2}$. As
we explained in Ref.Lee:2021oqr , the subleading quark pair corrections to the
twist-2 evolution equations enhance the twist-2 single parton FFs for a wide
range of $z$, as shown in Fig.2, which lead to 10–30% enhancement to the cross
sections even at a sizeable probing scale of $p_{T}=\mathcal{O}(100\,{\rm
GeV})$. The same feature of the nonlinear corrections was found when studying
nonlinear corrections to the DGLAP evolution of collinear PDFs Mueller:1985wy
.

Figure 3: (Upper): Theoretical calculations for the differential cross section
of prompt $J/\psi$ production at mid rapidity in hadronic collisions at
collider energies. Experimental data are taken from Refs.CDF:1997ykw ;
CDF:2004jtw ; ALICE:2012vpz ; ATLAS:2015zdw ; ATLAS:2017prf ; CMS:2015lbl ;
CMS:2017dju ; ALICE:2021edd . (Lower): Ratios of mean points of collider data
to theoretical predictions in the same $p_{T}$ bin for prompt $J/\psi$
production, obtained from the upper panel.

In Ref.Lee:2021oqr , we have compared our numerical results of $J/\psi$
production from the QCD factorization formalism with CMS data on prompt
$J/\psi$ production at $p_{T}\geq 60\,{\rm GeV}$ in $\sqrt{s}=7,\,13\,{\rm
TeV}$ pp collisions in the rapidity bin $|y|<1.2$ Khachatryan:2015rra ;
Sirunyan:2017qdw , and fixed the overall normalization factor of the theory
curve embedded in the input FFs. Our particular setup there leads to the
production rate of $J/\psi$ being dominated by the ${{}^{1}S_{0}^{[8]}}$
channel, along with some kind of cancelation between other channels. However,
it is important to note that our results are only sensitive to the input FFs
from the perturbative QCD factorized single-parton and
double-$Q\bar{Q}[\kappa]$ states with $\kappa=v^{[1,8]},a^{[1,8]},t^{[1,8]}$,
and not sensitive to the details of how each NRQCD channel, such as
${{}^{1}S_{0}^{[8]}}$, ${{}^{3}S_{1}^{[8]}}$, …, along with the LDMEs,
contributes to these input FFs. In this presentation, we extend our analysis
to other LHC data and those from Tevatron experiments, for which many precise
data sets on hadronic $J/\psi$ production at low $p_{T}$ down to
$p_{T}=\mathcal{O}(2m_{Q})$ are available. Without adjusting any normalization
and parameters, in Fig.3, we compare our calculations to all published data
from the LHC and Tevatron on hadronic prompt $J/\psi$ production in the mid
rapidity. The comparison in Fig.3 clearly demonstrates that the QCD
factorization approach with both LP and NLP contributions can describe the
data on prompt $J/\psi$ production at collider energies. With the same choice
of $\alpha=30$, $\beta=0.5$ and $K_{\rm NLP}=2$, the K-factor for the NLP
contribution that was evaluated at LO in $\alpha_{s}$ only, the ratio of data
over theory calculations is of the order of unity for a wide range of $p_{T}$.
Clearly, a global fit of the input FFs can further improve the theoretical
calculations.

In Fig.4, we focus on the comparison between our calculations with the
Tevatron data at relatively low $p_{T}$ to explore the role of LP and NLP
contributions. The LP curves (dashed lines) include the quark pair corrections
to the evolution equations. The NLP contributions (dotted lines) are evaluated
with the evolved $Q\bar{Q}$-FFs and the NLP partonic cross sections at the LO
along with a $K$-factor, $K_{\rm NLP}=2$, to mimic NLO contributions. The NLP
contribution becomes predominant over the LP contribution around
$p_{T}=15\,{\rm GeV}$ and below at Tevatron energies.

Figure 4: Theoretical calculations for the differential cross section of
prompt $J/\psi$ production at mid rapidity in hadronic collisions at Tevatron
energies. Experimental data from CDF are taken from Refs. CDF:1997ykw ;
CDF:2004jtw .

We expect the LP factorization formula in Eq.(2) to work for $J/\psi$
production when $p_{T}\gtrsim 10m_{J/\psi}$ where the
$\ln(p_{T}^{2}/m_{c}^{2})$-type contributions start to dominate. Although the
factorized LP contribution provides a good description of the published LHC
data at high $p_{T}$, its extrapolation to lower $p_{T}$ is below the data at
collider energies. The underestimation by the LP contribution emphasizes the
importance of the NLP contribution because it is more likely to get the
quarkonium from a fragmenting $Q\bar{Q}$-pair than a single parton
Kang:2014tta . Including the NLP contribution provides a better description of
data at collider energies even at $p_{T}\sim 3m_{J/\psi}$. The LP and NLP
contributions distinctly show different $p_{T}$ spectra; the high $p_{T}$ data
favor the LP curve, while the NLP curve is essential to reproduce the lower
$p_{T}$ data.

We notice that the NLP contribution starts to overshoot the data points when
$p_{T}\lesssim 2m_{J/\psi}$, below which the $\ln(p_{T}^{2}/(2m_{Q})^{2})$ is
not large, and its contribution cannot dominate the production cross sections.
We need to match our calculations to fixed-order calculations in the NRQCD
factorization or pNRQCD factorization Brambilla:2022ayc . We can implement the
matching between the QCD factorization and the NRQCD factorization by defining
Kang:2014pya

$\displaystyle d\sigma_{A+B\to H+X}$ $\displaystyle=d\sigma_{A+B\to
H+X}^{\textrm{QCD-Res}}+d\sigma_{A+B\to H+X}^{\textrm{NRQCD-
Fixed}}-d\sigma_{A+B\to H+X}^{\textrm{QCD-Asym}}$
$\displaystyle\Rightarrow\begin{dcases}d\sigma_{A+B\to H+X}^{\textrm{QCD-
Res}}&\textrm{when}\leavevmode\nobreak\ \quad p_{T}\gg
m_{H};\leavevmode\nobreak\ d\sigma^{\textrm{NRQCD-Fixed}}\approx
d\sigma^{\textrm{QCD-Asym}}\\\ d\sigma_{A+B\to H+X}^{\textrm{NRQCD-
Fixed}}&\textrm{when}\leavevmode\nobreak\ \quad p_{T}\to
m_{H};\leavevmode\nobreak\ d\sigma^{\textrm{QCD-Res}}\approx
d\sigma^{\textrm{QCD-Asym}}\end{dcases},$ (8)

where $d\sigma_{A+B\to H+X}^{\textrm{QCD-Res}}$ in Eq.(8) represents the QCD
factorization formalism in Eq.(1) including the resummation of
$\ln(p_{T}^{2}/m_{Q}^{2})$-type contributions, $d\sigma_{A+B\to
H+X}^{\textrm{NRQCD-Fixed}}$ represents the fixed order NRQCD (or pNRQCD)
factorization formalism that does not include the
$\ln(p_{T}^{2}/m_{Q}^{2})$-type resummation, and $d\sigma_{A+B\to
H+X}^{\textrm{QCD-Asym}}$ is equal to the fixed-order expansion of
$d\sigma_{A+B\to H+X}^{\textrm{QCD-Res}}$. The $d\sigma_{A+B\to
H+X}^{\textrm{QCD-Asy}}$ term is necessary for avoiding double-counting and
for a smooth transition between the high and low $p_{T}$ regimes. This
matching has to happen around $p_{T}\lesssim 3m_{J/\psi}$ for $J/\psi$
production at Tevatron and the LHC energies.

##  4 Summary

We discussed the calculations of $p_{T}$ distribution of prompt $J/\psi$
production in hadronic collisions at Tevatron, and the LHC energies in the
renormalization group improved QCD factorization approach. While we focused on
the LP contribution to $J/\psi$ production at high $p_{T}\geq 60$ GeV where
the $\ln(p_{T}^{2}/m_{Q}^{2})$-type contribution dominates in Ref.Lee:2021oqr
, we highlighted in this presentation the importance of the NLP contribution
in describing the data at relatively lower $p_{T}$. The NLP partonic cross
section with the twist-4 double parton FFs was calculated with the assumption
of $u=v=1/2$, which was justified for the leading contributions with O-O and
O-S color transitions. We also presented a numerical comparison between our
calculations and data from the LHC and Tevatron, showing the consistency for
$p_{T}\gtrsim 3m_{J/\psi}$. We also introduced a matching formalism for
extending theory calculations to $p_{T}\lesssim 3m_{J/\psi}$, and leave the
details of this matching for future work LQSW:2022 .

We thank Hee Sok Chung and Xiang-Peng Wang for their valuable discussions.
K.W. would also like to thank Jefferson Lab for the computational resources
essential to perform this project. This work is supported by Jefferson Science
Associates, LLC, under U.S. DOE Contract No. DE-AC05-06OR23177. The work of
K.L. and G.S. was supported in part by the National Science Foundation, award
1915093.

## References

  * (1) G.T. Bodwin, E. Braaten, G.P. Lepage, Phys.Rev. D51, 1125 (1995), hep-ph/9407339
  * (2) N. Brambilla, H.S. Chung, A. Vairo, X.P. Wang (2022), 2210.17345
  * (3) K. Lee, J.W. Qiu, G. Sterman, K. Watanabe, SciPost Phys. Proc. 8, 143 (2022), 2108.00305
  * (4) G.C. Nayak, J.W. Qiu, G.F. Sterman, Phys. Rev. D 72, 114012 (2005), hep-ph/0509021
  * (5) Z.B. Kang, Y.Q. Ma, J.W. Qiu, G. Sterman, Phys.Rev. D90, 034006 (2014), 1401.0923
  * (6) K. Lee, J.W. Qiu, G. Sterman, K. Watanabe, in preparation (2022)
  * (7) T.J. Hou et al., Phys. Rev. D 103, 014013 (2021), 1912.10053
  * (8) Y.Q. Ma, J.W. Qiu, H. Zhang, Phys.Rev. D89, 094029 (2014), 1311.7078
  * (9) Y.Q. Ma, J.W. Qiu, H. Zhang, Phys.Rev. D89, 094030 (2014), 1401.0524
  * (10) A.H. Mueller, J.W. Qiu, Nucl. Phys. B 268, 427 (1986)
  * (11) F. Abe et al. (CDF), Phys. Rev. Lett. 79, 572 (1997)
  * (12) D. Acosta et al. (CDF), Phys. Rev. D 71, 032001 (2005), hep-ex/0412071
  * (13) B. Abelev et al. (ALICE), JHEP 11, 065 (2012), 1205.5880
  * (14) G. Aad et al. (ATLAS), Eur. Phys. J. C 76, 283 (2016), 1512.03657
  * (15) M. Aaboud et al. (ATLAS), Eur. Phys. J. C 78, 171 (2018), 1709.03089
  * (16) V. Khachatryan et al. (CMS), Phys. Rev. Lett. 114, 191802 (2015), 1502.04155
  * (17) A.M. Sirunyan et al. (CMS), Phys. Lett. B 780, 251 (2018), 1710.11002
  * (18) S. Acharya et al. (ALICE), JHEP 03, 190 (2022), 2108.02523
  * (19) V. Khachatryan et al. (CMS), Phys. Rev. Lett. 114, 191802 (2015), 1502.04155
  * (20) A.M. Sirunyan et al. (CMS), Phys. Lett. B 780, 251 (2018), 1710.11002
  * (21) Z.B. Kang, Y.Q. Ma, J.W. Qiu, G. Sterman, Phys.Rev. D91, 014030 (2015), 1411.2456