# Second-Order Dissociation and Transition of Heavy Quarkonia in the Quark-
Gluon Plasma

Shouxing Zhao  Department of Applied Physics, Nanjing University of Science
and Technology, Nanjing 210094, China  Min He  Department of Applied Physics,
Nanjing University of Science and Technology, Nanjing 210094, China  Shanghai
Research Center for Theoretical Nuclear Physics, NSFC and Fudan University,
Shanghai 200438, China

###### Abstract

We revisit the dissociation of heavy quarkonia by thermal partons at the next-
to-leading order (NLO, also known as inelastic parton scattering dissociation)
in the Quark-Gluon Plasma (QGP). Utilizing the chromo-electric dipole coupling
from QCD multipole expansion as an effective Hamiltonian, this has been
conducted in the approach of second-order quantum mechanical perturbation
theory, which allows us to systematically incorporate the bound state wave
functions. Employing the quarkonium wave functions and binding energies
obtained from an in-medium potential model, we then numerically evaluate the
dissociation cross sections and rates for various charmonia and bottomonia,
where the infrared and collinear divergences are regularized by the thermal
masses of medium partons. We demonstrate that distinct from the leading order
(LO, also known as gluo-dissociation) counterparts peaking at relatively low
gluon energy and falling off thereafter, the NLO cross sections first grow and
then nearly saturate as the incident parton energy increases, as a result of
the outgoing parton carrying away the excess energy. The resulting NLO
dissociation rates increase with temperature and take over from the LO
counterparts toward high temperatures, similar to pertinent findings from
previous studies. We also evaluate the in-medium second-order transition
between different bound states, which may contribute to the total thermal
decay widths of heavy quarkonia in the QGP.

Heavy quarkonium, Quark-Gluon Plasma, Heavy Ion collision

###### pacs:

25.75.Dw, 12.38.Mh, 25.75.Nq

##  I Introduction

Heavy quarkonia, the bound states of charm or bottom quark-antiquark pair,
represent an versatile laboratory for testing properties of strong
interactions not only in vacuum but also at finite temperatures and densities
such as the environment of the deconfined medium (known as Quark-Gluon Plasma,
QGP) created in relativistic heavy-ion collisions Brambilla:2010cs ;
Rapp:2008tf ; Braun-Munzinger:2009dzl ; Kluberg:2009wc ; Mocsy:2013syh ;
Andronic:2015wma ; Rothkopf:2019ipj ; Zhao:2020jqu ; Andronic:2024oxz . It was
first advocated by Matsui and Satz Matsui:1986dk  that the screening of color
charges in the deconfined medium should weaken the binding between the heavy
quark and its antiquark, leading eventually to the melting of their bound
states and thus serving as a signal of the formation of QGP. This picture
further suggests that different heavy quarkonium should melt sequentially at
temperatures in the order of their vacuum binding energies Digal:2001ue .

However, the complexity of in-medium quarkonium dynamics prevents one from
using them as a straightforward thermometer of the medium based on the above
purely static screening scenario. On the one hand, mechanisms involving
collisions with medium constituents can lead to dynamical dissociation of the
bound states Kharzeev:1994pz ; Grandchamp:2001pf . On the other hand, in
situations where heavy quarks are abundantly produced, e.g., charm production
in Pb-Pb collisions at the LHC energy, individual heavy quarks and antiquarks
diffusing and thermalizing in the QGP may recombine and regenerate the bound
states Braun-Munzinger:2000csl ; Thews:2000rj ; Grandchamp:2003uw . The
dissociation rates/collisional widths generated from partonic inelastic
scatterings with thermal partons serve as a pivotal input for both
semiclassical Zhao:2010nk ; Song:2011nu ; Strickland:2011mw ; Zhou:2014kka
and quantum Brambilla:2016wgg ; Brambilla:2020qwo ; Yao:2021lus  transport
descriptions of the heavy quarkonium production in the QGP. These reaction
rates, as facilitated by the principle of detailed balance and implemented in
the kinetic rate equations Zhao:2010nk ; Wu:2024gil , also govern the
regeneration processes.

From the theoretical point of view, the dissociation of heavy quarkonium
through collisions with thermal partons can be classified into leading order
(LO) and next-to-leading order (NLO) processes. The LO process refers to the
gluo-dissociation, $g+\Psi\rightarrow Q+\bar{Q}$ ($\Psi$ denotes a heavy
quarkonium, and $Q$ ($\bar{Q}$) the heavy quark (antiquark)), in which the
bound state absorbs a thermal gluon from the medium and breaks up into an
unbound color octet, $(Q\bar{Q})_{8}$, by overcoming the binding energy. The
analysis of this process was initiated by Peskin within the operator-product-
expansion approach, who identified the coupling of the tightly bound heavy
quarkonium with external soft gluons as a gauge invariant color-electric
dipole type Peskin:1979va ; Bhanot:1979vb . Using the same coupling or its
variants (e.g., nonrelativistic approximation of the associated Bethe-Salpeter
amplitude Oh:2001rm , the color-magnetic dipole transition Chen:2017jje ),
many authors have thereafter investigated the same process Kharzeev:1994pz ;
Oh:2001rm ; Wong:2004zr ; Brezinski:2011ju ; Liu:2013kkg ; Chen:2017jje . The
cross section of this process typically shows a peak at the incident gluon
energy when the gluon wavelength matches the bound state size, but quickly
drops off toward higher gluon energies Chen:2017jje .

Going to the NLO, $p+\Psi\rightarrow p+Q+\bar{Q}$ (the thermal parton $p=g,q$
or $\bar{q}$), also known as dissociation by inelastic parton scattering, the
outgoing parton in the final state carries away the excess energy, which thus
may help maintain the efficiency for the break-up of the bound state even for
incident gluons with higher energies. The phenomenological relevance of this
process was first put forward from the “quasifree” perspective
Grandchamp:2001pf , by noting that in the limit of small binding (or large
bound state size) at high temperatures, the incident thermal parton with
energy $\sim T$ (temperature of the medium) “sees” the individual constituents
within the bound state, rather than the quarkonium as a whole, and therefore
the scattering effectively occurs on a single $Q$ with the other $\bar{Q}$
being a spectator (or vice versa) Riek:2010py ; the interference effect
between the parton scattering off the $Q$ and $\bar{Q}$ as a result of the
residual binding was originally ignored and later corrected for Du:2017qkv  in
the “quasifree” treatment. The cross section of the NLO process was also
computed in the perturbative QCD approach Song:2005yd  (and recently updated
in Hong:2018vgp ), by identifying an effective vertex from the non-
relativistic Bethe-Salpeter amplitude used earlier in the LO calculation
Oh:2001rm .

The identification of an imaginary part of the heavy quark potential in the
effective field theory (EFT) approach Laine:2006ns ; Beraudo:2007ky ;
Brambilla:2008cx ; Brambilla:2013dpa , as arising from the Landau damping of
the space-like gluon exchanged between $Q$ and $\bar{Q}$ in the environment of
QGP, represents a paradigm shift for the NLO process of heavy quarkonium
dissociation. Microscopically, the Landau damping underlying the imaginary
potential originates from the energy transfer from the low-frequency gluons
mediating the interaction between $Q$ and $\bar{Q}$ to the “hard” particles
with energy of order $\sim T$ in the thermal bath Laine:2006ns ;
Beraudo:2007ky . This can be clearly seen from “cutting” the (one-loop) self-
energy diagram of the exchanged space-like gluon, thereby making a close
connection to aforementioned “quasifree” process as well as the elastic
collisional broadening/thermalization of a single heavy quark. The imaginary
static potential has been now verified from first-principle lattice QCD
analysis of the Wilson loop spectral functions Burnier:2014ssa ; Bala:2021fkm
; Bazavov:2023dci .

It has been established that the LO process dominates the heavy quarkonium
dissociation at low temperatures when the associated binding energy is still
large compared to the screening mass, whereas the NLO takes over as the
binding becomes weak toward high temperatures Zhao:2010nk ; Wu:2024gil ;
Du:2017qkv ; Hong:2018vgp ; Brambilla:2013dpa . However, the dissociation
rates currently used in different heavy quarkonium transport models suffer
from large discrepancies Andronic:2024oxz , calling further for reliable
evaluations of the heavy quarkonium dissociation processes. The purpose of the
present work is to revisit the cross sections and rates for the NLO heavy
quarkonium dissociation by thermal partons in the QGP from the dynamical
scattering point of view, in the similar spirit of Hong:2018vgp  but from
different approach. We calculate the NLO inelastic scattering cross sections
within the framework of the second-order quantum mechanical perturbation
theory, utilizing the color-electric dipole coupling of the heavy quarkonium
with external gluons Peskin:1979va ; Yan:1980uh ; Brambilla:2004jw ;
Voloshin:2007dx ; Sumino:2014qpa  as the effective interaction Hamiltonian.
This approach was first introduced in Chen:2018dqg  by one of the present
authors as an attempt to study the NLO break-up of heavy quarkonium by thermal
gluons. The present work aims to present a more systematic and clear treatment
of the approach and generalize the calculation to include thermal quarks (and
antiquarks). Another theme tackled in the present work is the transitions
between different eigenstates of quarkonium bound states, which also occur at
the second order involving two gluons, since the bound states in both of the
initial and final state are color singlets. We compute the pertinent
transition rates, which may contribute to the total thermal spectral width of
heavy quarkonium as evaluated in lattice QCD Larsen:2019zqv . The current
approach allows us to systematically incorporate the bound state wave
functions, thus going beyond the “quasifree” approximation. The direct
calculation of the cross sections from dynamical scattering processes also
makes the current approach distinct from the method in Brambilla:2013dpa
where the NLO cross section was extracted from an formula of the thermal decay
width expressed in terms of the imaginary potential Laine:2006ns . The color-
electric dipole coupling was also employed in Yao:2018sgn  to compute the
heavy quarkonium scattering amplitudes in the QGP from the potential non-
relativistic QCD (pNRQCD) perspective, but no derivations of the cross
sections were given.

The article is organized as follows. In Sec. II, the QCD system of a heavy
quarkonium and external light quarks and gluons is formulated in terms of an
effective Hamiltonian. In Sec. III, we derive the cross sections for the NLO
dissociation of heavy quarkonium by thermal gluons and quarks within the
second-order quantum mechanical perturbation framework. In Sec. IV, we employ
the bound state wave functions and binding energies computed from an in-medium
potential model and numerically evaluate the NLO cross sections and
dissociation rates for various kinds of charmonia and bottomonia. In Sec. V,
we study the second-order transition processes between different eigenstates
of quarkonium singlets and compare the pertinent cross sections and rates to
the NLO dissociation counterparts. We finally summarize in Sec. VI.

##  II The effective Hamiltonian

In order to compute the NLO processes in the second-order quantum mechanical
perturbation approach, we first formulate the QCD system of a heavy quarkonium
and external light quarks and gluons in terms of an effective Hamiltonian. The
non-relativistic heavy quarkonium system possesses a Hamiltonian

$\displaystyle H_{0}=\frac{{\vec{p}}^{2}}{m_{Q}}+V_{1}(r),$ (1)

where $V_{1}$ is the binding potential for the color singlet $Q{\bar{Q}}$ (we
ignore the weak repulsive potential in the color octet configuration). The
coupling of the compact heavy quarkonium with external soft gluons can be
summarized into an effective color-electric dipole type Peskin:1979va ;
Yan:1980uh ; Brambilla:2004jw ; Voloshin:2007dx ; Sumino:2014qpa  with two
vertices

$\displaystyle H_{Q\bar{Q}g}=V_{SO}+V_{OO},$ (2)

where $V_{SO}$ represents the transition of a $Q\bar{Q}$ singlet state
$|S\rangle$ to an octet $|O,a\rangle$ ($a=1,...,8$ denotes its color index)
through interacting with an external gluon, which can be expressed in terms of
the matrix element

$\displaystyle\langle O,a|V_{SO}|S\rangle$ $\displaystyle=\langle
O,a|\frac{1}{2}g_{s}\vec{r}(\frac{\lambda^{b}}{2}-\frac{\bar{\lambda}^{b}}{2})\cdot\vec{E}^{b}|S\rangle$
$\displaystyle=\frac{g_{s}}{\sqrt{2N_{c}}}\vec{E}^{a}(t,\vec{x})\cdot\langle
O|\vec{r}|S\rangle.$ (3)

$V_{OO}$ in Eq. (2) represents the transition from an octet of color $b$ to
another octet of color $a$ with the matrix element

$\displaystyle\langle
O,a|V_{OO}|O,b\rangle=\frac{ig_{s}}{2}d^{abc}\vec{E}^{c}(t,\vec{x})\cdot\langle
O|\vec{r}|O\rangle.$ (4)

In Eqs.(II) and (4), $N_{c}=3$ is the number of colors in the fundamental
representation of $SU(3)_{c}$, $\frac{\lambda^{b}}{2}$
($\frac{\bar{\lambda}^{b}}{2}$) the color matrix of $Q$ ($\bar{Q}$), and
$d^{abc}=2tr[\lambda^{a}/2\\{\lambda^{b}/2,\lambda^{c}/2\\}]$ a totally
symmetric $SU(3)_{c}$ group invariant.

The Hamiltonian of the external light quarks, antiquarks and gluons reads
Peskin:1995

$\displaystyle H_{q\bar{q}g}$ $\displaystyle=\int
d^{3}x{\psi^{i}}^{\dagger}(x)\Big{[}\beta
m_{q}-i\vec{\alpha}\cdot\Big{(}\vec{\partial}-ig_{s}\vec{A}^{a}(\frac{\lambda^{a}}{2})^{ij}\Big{)}\Big{]}\psi^{j}(x)$
$\displaystyle+\frac{1}{2}\int
d^{3}x\Big{[}{\vec{E^{a}}}^{2}(x)+\vec{B^{a}}^{2}(x)\Big{]},$ (5)

where the $\vec{\alpha}$ and $\beta$ are Dirac matrices, $i,j=1,2,3$ the color
indices of quarks, and $m_{q}$ the light quark mass. The self-interaction of
gluons is implicit in the color-electric and color-magnetic field energy
densities. We will work with the Weyl gauge $A_{0}^{a}(x)=0$ Sakurai-AQM:2008
, such that the color-electric field $\vec{E}^{a}(x)=\nabla
A_{0}^{a}-\partial_{0}\vec{A}^{a}+g_{s}f^{abc}\vec{A}^{b}A_{0}^{c}$ appearing
in Eqs.(II) and (4) reduces to

$\displaystyle\vec{E}^{a}(t,\vec{x})=-\frac{\partial}{\partial
t}\vec{A}^{a}(t,\vec{x}).$ (6)

Consequently, the self-interaction of gluons now arises only from the energy-
density $1/2\vec{B}^{a}\cdot\vec{B}^{a}$ of the color-magnetic field
$\vec{B}^{a}(x)=\nabla\times\vec{A}^{a}(x)-\frac{1}{2}g_{s}f^{abc}\vec{A}^{b}(x)\times\vec{A}^{c}(x)$.
In particular, the three-gluon vertex (the four-gluon vertex appears at higher
order in $g_{s}$ and is thus ignored here) reads

$\displaystyle V_{3g}$ $\displaystyle=-\frac{1}{4}g_{s}f^{abc}\int
d^{3}x\Big{[}(\nabla\times\vec{A}^{a})\cdot(\vec{A}^{b}\times\vec{A}^{c})$
$\displaystyle+(\vec{A}^{b}\times\vec{A}^{c})\cdot(\nabla\times\vec{A}^{a})\Big{]},$
(7)

where $f^{abc}$ is the totally antisymmetric $SU(3)_{c}$ structure constant
and the transverse gluon field is quantized in a box of volume $V$ Sakurai-
AQM:2008

$\displaystyle\vec{A}^{a}(t,\vec{x})=\sum_{\vec{k},\lambda}\sqrt{\frac{1}{2V\omega_{k}}}\vec{\epsilon}_{\vec{k},\lambda}\Big{[}a^{a}_{\vec{k},\lambda}e^{i\vec{k}\cdot\vec{x}-i\omega_{k}t}+h.c.\Big{]},$
(8)

with $\omega_{\vec{k}}$ being the gluon energy with momentum $\vec{k}$ and
$\vec{\epsilon}_{\vec{k},\lambda}$ the polarization vector ($\lambda$=1, 2
denotes two transverse polarizations) that satisfies

$\displaystyle\sum_{\lambda=1,2}\epsilon^{(i)}_{\vec{k}\lambda}\epsilon^{(j)}_{\vec{k}\lambda}=\delta^{ij}-\frac{k^{i}k^{j}}{\vec{k}^{2}}.$
(9)

The creation and annihilation operators in Eq.(8) satisfies the commutation
relation
$[a_{\vec{k},\lambda}^{a}(t),{a_{\vec{k}^{\prime},\lambda^{\prime}}^{a^{\prime}\dagger}}(t)]=\delta_{\vec{k}\vec{k}^{\prime}}\delta_{\lambda\lambda^{\prime}}\delta^{aa^{\prime}}$,
with one-gluon state
$|g(\vec{k},\lambda,a)\rangle=a_{\vec{k},\lambda}^{a\dagger}|0\rangle$.

Figure 1: The four coupling vertices. The upper left, upper right, lower left
and lower right refer to the $V_{SO}$, $V_{OO}$, $V_{3g}$ and $V_{q\bar{q}g}$
vertex, respectively. The single solid line, the double solid line, the single
dashed line and the wavy line denote the heavy quarkonium singlet, the octet,
the light quark/antiquark and the gluon, respectively.

The coupling of the light quarks with gluons from Eq.(II) reads

$\displaystyle V_{q\bar{q}g}=g_{s}\int
d^{3}x{\psi^{i}}^{\dagger}(x)\vec{\alpha}\cdot\vec{A}^{a}(x)(\frac{\lambda^{a}}{2})^{ij}\psi^{j}(x).$
(10)

The light quark fields are also quantized in a box of volume $V$ Sakurai-
AQM:2008

$\displaystyle{\psi^{j}}(t,\vec{x})=\sum_{\vec{k},s}\sqrt{\frac{m_{q}}{E_{\vec{k}}V}}\Big{(}$
$\displaystyle
b^{j}_{\vec{k},s}u^{s}(\vec{k})e^{i\vec{k}\cdot\vec{x}-iE_{\vec{k}}t}$
$\displaystyle+d^{j\dagger}_{\vec{k},s}v^{s}(\vec{k})e^{-i\vec{k}\cdot\vec{x}+iE_{\vec{k}}t}\Big{)},$
(11)
$\displaystyle{\psi^{i}}^{\dagger}(t,\vec{x})=\sum_{\vec{k},s}\sqrt{\frac{m_{q}}{E_{\vec{k}}V}}\Big{(}$
$\displaystyle
d^{i}_{\vec{k},s}v^{s\dagger}(\vec{k})e^{i\vec{k}\cdot\vec{x}-iE_{\vec{k}}t}$
$\displaystyle+b^{i\dagger}_{\vec{k},s}u^{s\dagger}(\vec{k})e^{-i\vec{k}\cdot\vec{x}+iE_{\vec{k}}t}\Big{)},$
(12)

where $E_{\vec{k}}=\sqrt{\vec{k}^{2}+m_{q}^{2}}$ is the energy of light quark
(or antiquark) of momentum $\vec{k}$. The light quark/antiquark creation and
annihilation operators satisfy the anti-commutation relation
$\\{b_{\vec{k},s}^{i}(t),{b_{\vec{k}^{\prime},s^{\prime}}^{i^{\prime}\dagger}}(t)\\}=\\{d_{\vec{k},s}^{i}(t),{d_{\vec{k}^{\prime},s^{\prime}}^{i^{\prime}\dagger}}(t)\\}=\delta_{\vec{k}\vec{k}^{\prime}}\delta_{ss^{\prime}}\delta^{ii^{\prime}}$,
with one-quark state
$|q(\vec{k},s,i)\rangle=b_{\vec{k},s}^{i\dagger}|0\rangle$ ($s,s^{\prime}=1,2$
denotes the two spin states of the light quark or antiquark) and one-antiquark
state $|\bar{q}(\vec{k},s,i)\rangle=d_{\vec{k},s}^{i\dagger}|0\rangle$.
Following the conventions in Sakurai-AQM:2008 , the orthogonal and
normalization relation for spinors reads
$u^{s\dagger}(\vec{k})u^{s^{\prime}}(\vec{k})=v^{s\dagger}(\vec{k})v^{s^{\prime}}(\vec{k})=\frac{E_{\vec{k}}}{m_{q}}\delta^{ss^{\prime}}$,
and the completeness relations are

$\displaystyle\sum_{s=1}^{2}u_{\alpha}^{s}(\vec{k})\bar{u}_{\beta}^{s}(\vec{k})=\Big{(}\frac{-i\gamma\cdot
k+m}{2m}\Big{)}_{\alpha\beta},$ (13)
$\displaystyle\sum_{s=1}^{2}v_{\alpha}^{s}(\vec{k})\bar{v}_{\beta}^{s}(\vec{k})=-\Big{(}\frac{i\gamma\cdot
k+m}{2m}\Big{)}_{\alpha\beta},$ (14)

where $\bar{u}=u^{\dagger}\gamma_{4}$, $\bar{v}=v^{\dagger}\gamma_{4}$ and
$\gamma\cdot k=\gamma_{\mu}k_{\mu}$ ($\,k_{\mu}=(\vec{k},iE_{\vec{k}}))$ with
the Hermitian matries $\gamma_{\mu}$ satisfying
$\alpha_{k}=i\gamma_{4}\gamma_{k}$, $\gamma_{4}=\beta$ and
$\\{\gamma_{\mu},\,\gamma_{\nu}\\}=2\delta_{\mu\nu}$.

To sum up, four coupling vertices come out of the QCD system of a heavy
quarkonium and external light quarks and gluons, as depicted in Fig. 1, which
will be used in Sec. III to construct the pertinent Feynman diagrams for the
heavy quarkonium NLO dissociation processes.

##  III Deriving the NLO dissociation cross sections

In this section we derive the NLO dissociation cross sections of heavy
quarkonium due to inelastic scatterings by medium partons. We use the $s$-wave
ground state charmonium $J/\psi$ as the example to illustrate the derivations,
but in Sec. IV, we will demonstrate numerical results for all heavy quarkonia
including also the $p$-wave states. The derivations are performed within the
approach of the second-order old-fashioned quantum mechanical perturbation
theory Sakurai-AQM:2008 ; Schwartz-QFT:2014 ; Greiner:1998 ; Gottfried-
Yan:2003 , where all states are on-shell at all times and energy is not
conserved at each vertex. The key quantity is the second-order transition
amplitude

$\displaystyle T_{fi}=\sum_{m}\frac{\langle f|V|m\rangle\langle
m|V|i\rangle}{E_{i}-E_{m}+i\epsilon},$ (15)

where the initial state $|i\rangle$, final state $|f\rangle$ and intermediate
states $|m\rangle$ are all eigenstates of the unperturbed (i.e., the free
heavy quarkonium plus free light quarks/antiquarks and gluons system)
Hamiltonian with eigenenergies $E_{i}=E_{f}$ and $E_{m}$, respectively, and
$V$ denotes the coupling vertices (interaction Hamiltonian).

###  III.1 Dissociation involving $V_{SO}$ and $V_{OO}$

Using vertices $V_{SO}$ and $V_{OO}$, we can construct the Feynman diagrams
(a) and (b) for the NLO dissociation of $J/\psi$ by an external gluon,
$g+J/\psi\to g+c+\bar{c}$, as shown in Fig. 2. These two diagrams are the
counterparts of the $s$\- and $u$-channel diagrams in relativistically
covariant perturbation calculations for the photon-electron Compton scattering
Peskin:1995 . In the following, we work in the rest frame of the $J/\psi$ and
neglect the three-momentum transfer from the incident gluon to the $J/\psi$.
The latter approximation is justified by the fact that the mass of the heavy
quarkonium is much larger than the typical momentum of thermal partons in the
QGP. As a result, the center-of-mass momentum of the final state unbound octet
$(c\bar{c})_{8}$ is also neglected and one deals only with the internal
relative momentum of $(c\bar{c})_{8}$. The initial state for diagrams (a) and
(b) consists of a $J/\psi$ at rest and an incident gluon of momentum
$\vec{k}$, polarization $\lambda$ and color $a$,
$|i\rangle=|J/\psi,g(\vec{k},\lambda,a)\rangle$, whereas the final state
involves an unbound octet $(c\bar{c})_{8}$ of internal relative momentum
$\vec{p}$ and color $b$, plus an outgoing gluon of momentum $\vec{\kappa}$,
polarization $\sigma$ and color $c$,
$|f\rangle=|(c\bar{c})_{8}(\vec{p},b),g(\vec{\kappa},\sigma,c)\rangle$.

Figure 2: Feynman diagrams constructed from $V_{SO}$ and $V_{OO}$ for
$g+J/\psi\to g+c+\bar{c}$. The time direction goes upward. Between the initial
time $t_{i}$ and final time $t_{f}$ is the intermediate state.

For diagram (a), the intermediate state is an octet of internal relative
momentum $\vec{q}$ and color $d$,
$|m\rangle=|(c\bar{c})_{8}(\vec{q},d)\rangle$. Then the transition matrix
element due to $V_{SO}$ is calculated making use of Eqs. (II), (6), (8)

$\displaystyle\langle
m|V_{SO}|i\rangle=\langle(c\bar{c})_{8}(\vec{q},d)|V_{SO}|J/\psi,g(\vec{k},\lambda,a)\rangle$
$\displaystyle=\frac{g_{s}}{\sqrt{6}}\langle
0|\vec{E}^{d}(t,\vec{x})|g(\vec{k},\lambda,a)\rangle\cdot\langle(c\bar{c})_{8}(\vec{q})|\vec{r}|J/\psi\rangle$
$\displaystyle=i\delta^{da}\frac{g_{s}}{\sqrt{6}}\sqrt{\frac{\omega_{\vec{k}}}{2V}}e^{i\vec{k}\cdot\vec{x}}\vec{\epsilon}_{\vec{k}\lambda}\cdot\langle(c\bar{c})_{8}(\vec{q})|\vec{r}|J/\psi\rangle$
$\displaystyle=g_{s}\frac{\delta^{da}}{V}\sqrt{\frac{\pi\omega_{\vec{k}}}{3}}e^{i\vec{k}\cdot\vec{x}}(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{q})\frac{1}{q}\int
r^{3}drj_{1}(qr)R_{10}(r),$ (16)

where $R_{10}(r)$ is the normalized radial wave function for $J/\psi$. We have
neglected the weak repulsive potential for the unbound octet $(c\bar{c})_{8}$,
so that the wave function for the internal relative motion is a plane wave
$e^{i\vec{q}\cdot\vec{r}}/\sqrt{V}$. Upon a spherical wave expansion for the
plane wave,
$e^{i\vec{q}\cdot\vec{r}}=4\pi\sum_{l,m}i^{l}j_{l}(qr)Y_{lm}(\theta,\phi)Y_{lm}^{*}(\theta^{\prime},\phi^{\prime})$
(primed angles for $\vec{q}$, and unprimed for $\vec{r}$), and use of the
orthogonality relation for the spherical harmonics $\int\int\sin\theta d\theta
d\phi
Y_{lm}^{*}(\theta,\phi)Y_{l^{\prime}m^{\prime}}(\theta,\phi)=\delta_{ll^{\prime}}\delta_{mm^{\prime}}$,
only the first-order spherical Bessel function survives. Similarly, the other
matrix element due to $V_{OO}$ is calculated using Eqs. (4), (6), (8)

$\displaystyle\langle
f|V_{OO}|m\rangle=\langle(c\bar{c})_{8}(\vec{p},b),g(\vec{\kappa},\sigma,c)|V_{OO}|(c\bar{c})_{8}(\vec{q},d)\rangle$
$\displaystyle=\frac{ig_{s}}{2}d^{bde}\langle
g(\vec{\kappa},\sigma,c)|\vec{E}^{e}(t,\vec{x})|0\rangle\cdot\langle(c\bar{c})_{8}(\vec{p})|\vec{r}|(c\bar{c})_{8}(\vec{q})\rangle$
$\displaystyle=-d^{bdc}\frac{ig_{s}}{2V}(2\pi)^{3}\sqrt{\frac{\omega_{\vec{\kappa}}}{2V}}e^{-i\vec{\kappa}\cdot\vec{x}}\vec{\epsilon}_{\vec{\kappa}\sigma}\cdot\nabla_{\vec{q}}\delta^{3}(\vec{q}-\vec{p}),$
(17)

where we have used $\int
d^{3}re^{i(\vec{q}-\vec{p})\cdot\vec{r}}=(2\pi)^{3}\delta^{3}(\vec{q}-\vec{p})$
for the continuous momenta.

Now substituting these two matrix elements into Eq. (15) and carrying out the
summation over the momentum and color of the intermediate state
$\sum_{m}=\sum_{\vec{q}}\sum_{d}=V/(2\pi)^{3}\int d^{3}\vec{q}\sum_{d}$, one
arrives at the transition amplitude for diagram (a) after an integration by
part to integrate out the $\delta^{3}(\vec{q}-\vec{p})$

$\displaystyle T_{fi}^{(a)}=$ $\displaystyle
d^{abc}\frac{ig_{s}^{2}}{2V}\sqrt{\frac{\pi\omega_{\vec{k}}\omega_{\vec{\kappa}}}{6V}}e^{i\vec{k}\cdot\vec{x}-i\vec{\kappa}\cdot\vec{x}}$
$\displaystyle\times\vec{\epsilon}_{\vec{\kappa}\sigma}\cdot\Big{[}A(p,k)\vec{\epsilon}_{\vec{k}\lambda}+(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{p})\frac{\vec{p}}{p^{2}}B(p,k)\Big{]},$
(18)

where

$\displaystyle A(p,k)=\frac{\int
r^{3}drj_{1}(pr)R_{10}(r)}{p(-\epsilon_{B}+\omega_{\vec{k}}-\frac{p^{2}}{m_{Q}}+i\epsilon)},$
$\displaystyle B(p,k)=\frac{\frac{2p}{m_{Q}}\int
r^{3}drj_{1}(pr)R_{10}(r)}{(-\epsilon_{B}+\omega_{\vec{k}}-\frac{p^{2}}{m_{Q}}+i\epsilon)^{2}}-\frac{\int
r^{4}drj_{2}(pr)R_{10}(r)}{(-\epsilon_{B}+\omega_{\vec{k}}-\frac{p^{2}}{m_{Q}}+i\epsilon)}.$
(19)

contain the bound state wave function and binding energy $\epsilon_{B}$. Here
all energies are measured relative to the threshold of the $c\bar{c}$ pair.

For diagram (b), the intermediate
$|m\rangle=|(c\bar{c})_{8}(\vec{q},d),\,g(\vec{k},\lambda,a),\,g(\vec{\kappa},\sigma,c)\rangle$
involves two gluons that have the same quantum numbers as the incident and
outgoing gluons, respectively. Similar manipulations as above lead to the
matrix elements due to $V_{SO}$ and $V_{OO}$, respectively

$\displaystyle\langle
m|V_{SO}|i\rangle=\langle(c\bar{c})_{8}(\vec{q},d),g(\vec{\kappa},\sigma,c)|V_{SO}|J/\psi\rangle$
$\displaystyle=-g_{s}\frac{\delta^{dc}}{V}\sqrt{\frac{\pi\omega_{\vec{\kappa}}}{3}}e^{-i\vec{\kappa}\cdot\vec{x}}(\vec{\epsilon}_{\vec{\kappa}\sigma}\cdot\vec{q})\frac{1}{q}\int
r^{3}drj_{1}(qr)R_{10}(r),$ (20) $\displaystyle\langle
f|V_{OO}|m\rangle=\langle(c\bar{c})_{8}(\vec{p},b)|V_{OO}|(c\bar{c})_{8}(\vec{q},d),g(\vec{k},\lambda,a)\rangle$
$\displaystyle=d^{bda}\frac{ig_{s}}{2V}(2\pi)^{3}\sqrt{\frac{\omega_{\vec{k}}}{2V}}e^{i\vec{k}\cdot\vec{x}}\vec{\epsilon}_{\vec{k}\lambda}\cdot\nabla_{\vec{q}}\delta^{3}(\vec{q}-\vec{p}).$
(21)

Finally the transition amplitude for diagram (b) reads

$\displaystyle T_{fi}^{(b)}=$ $\displaystyle
d^{abc}\frac{ig_{s}^{2}}{2V}\sqrt{\frac{\pi\omega_{\vec{k}}\omega_{\vec{\kappa}}}{6V}}e^{i\vec{k}\cdot\vec{x}-i\vec{\kappa}\cdot\vec{x}}$
$\displaystyle\times\vec{\epsilon}_{\vec{k}\lambda}\cdot\Big{[}C(p,\kappa)\vec{\epsilon}_{\vec{\kappa}\sigma}+(\vec{\epsilon}_{\vec{\kappa}\sigma}\cdot\vec{p})\frac{\vec{p}}{p^{2}}D(p,\kappa)\Big{]},$
(22)

where

$\displaystyle C(p,\kappa)=\frac{\int
r^{3}drj_{1}(pr)R_{10}(r)}{p(-\epsilon_{B}-\omega_{\vec{\kappa}}-\frac{p^{2}}{m_{Q}}+i\epsilon)},$
$\displaystyle D(p,\kappa)=\frac{\frac{2p}{m_{Q}}\int
r^{3}drj_{1}(pr)R_{10}(r)}{(-\epsilon_{B}-\omega_{\vec{\kappa}}-\frac{p^{2}}{m_{Q}}+i\epsilon)^{2}}-\frac{\int
r^{4}drj_{2}(pr)R_{10}(r)}{(-\epsilon_{B}-\omega_{\vec{\kappa}}-\frac{p^{2}}{m_{Q}}+i\epsilon)}.$
(23)

again contain the information of bound state wave function and binding energy.

###  III.2 Dissociation involving $V_{SO}$ and $V_{3g}$

The NLO dissociation of $J/\psi$ by an external gluon, $g+J/\psi\to
g+c+\bar{c}$ can be also constructed using the vertices of $V_{SO}$ and
$V_{3g}$, as represented by the Feynman diagrams (c) and (d) in Fig. 3. These
two diagrams have the same initial and final state as specified for diagrams
(a) and (b).

Figure 3: Feynman diagrams constructed from $V_{SO}$ and $V_{3g}$ for
$g+J/\psi\to g+c+\bar{c}$. The time direction goes upward. Between the initial
time $t_{i}$ and final time $t_{f}$ is the intermediate state.

For diagram (c), the intermediate state contains an octet with the same
quantum numbers as the $(c\bar{c})_{8}$ in the final state, a gluon with the
same quantum numbers as the incident gluon, and another gluon whose momentum,
polarization and color are not fixed and thus denoted as $(\vec{k_{1}}$,
$\lambda_{1}$ and $d_{1})$, respectively, such that
$|m\rangle=|(c\bar{c})_{8}(\vec{p},b),g(\vec{k},\lambda,a),g(\vec{k_{1}},\lambda_{1},d_{1})\rangle$.
Similar to Eq. (III.1), the transition matrix element due to $V_{SO}$ reads

$\displaystyle\langle
m|V_{SO}|i\rangle=\langle(c\bar{c})_{8}(\vec{p},b),\,g(\vec{k}_{1},\lambda_{1},d_{1})|V_{SO}|J/\psi\rangle$
$\displaystyle=-g_{s}\frac{\delta^{d_{1}b}}{V}\sqrt{\frac{\pi\omega_{\vec{k}_{1}}}{3}}e^{-i\vec{k}_{1}\cdot\vec{x}}(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\cdot\vec{p})\frac{1}{p}\int
r^{3}drj_{1}(pr)R_{10}(r).$ (24)

To calculate the transition matrix element due to $V_{3g}$, we substitute Eq.
(8) into Eq. (II) and obtain the relevant term involving the combination of
two annihilation operators and one creation operator

$\displaystyle
V_{3g}(a^{\dagger}aa)=f^{abc}\frac{ig_{s}}{4}\sum_{\vec{k}_{1}\lambda_{1}}\sum_{\vec{k}_{2}\lambda_{2}}\sum_{\vec{k}_{3}\lambda_{3}}\sqrt{\frac{1}{2V\omega_{\vec{k}_{1}}\omega_{\vec{k}_{2}}\omega_{\vec{k}_{3}}}}$
$\displaystyle(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\times\vec{k}_{1})\cdot(\vec{\epsilon}_{\vec{k}_{2}\lambda_{2}}\times\vec{\epsilon}_{\vec{k}_{3}\lambda_{3}})\Big{(}a_{\vec{k}_{2}\lambda_{2}}^{b\dagger}a_{\vec{k}_{3}\lambda_{3}}^{c}a_{\vec{k}_{1}\lambda_{1}}^{a}\delta_{\vec{k}_{2},\vec{k}_{3}+\vec{k}_{1}}$
$\displaystyle+a_{\vec{k}_{3}\lambda_{3}}^{c\dagger}a_{\vec{k}_{1}\lambda_{1}}^{a}a_{\vec{k}_{2}\lambda_{2}}^{b}\delta_{\vec{k}_{3},\vec{k}_{1}+\vec{k}_{2}}-a_{\vec{k}_{1}\lambda_{1}}^{a\dagger}a_{\vec{k}_{2}\lambda_{2}}^{b}a_{\vec{k}_{3}\lambda_{3}}^{c}\delta_{\vec{k}_{1},\vec{k}_{2}+\vec{k}_{3}}\Big{)},$
(25)

which yields

$\displaystyle\langle f|V_{3g}|m\rangle=\langle
g(\vec{\kappa},\sigma,c)|V_{3g}(a^{\dagger}aa)|g(\vec{k},\lambda,a),\,g(\vec{k}_{1},\lambda_{1},d_{1})\rangle$
$\displaystyle=f^{abc}\frac{ig_{s}}{2}\sqrt{\frac{1}{2V\omega_{\vec{k}}\omega_{\vec{\kappa}}\omega_{\vec{k}_{1}}}}\delta_{\vec{\kappa},\vec{k}+\vec{k}_{1}}\Big{[}(\vec{\epsilon}_{\vec{k}\lambda}\times\vec{k})\cdot(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\times\vec{\epsilon}_{\vec{\kappa}\sigma})$
$\displaystyle-(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\times\vec{k}_{1})\cdot(\vec{\epsilon}_{\vec{k}\lambda}\times\vec{\epsilon}_{\vec{\kappa}\sigma})+(\vec{\epsilon}_{\vec{\kappa}\sigma}\times\vec{\kappa})\cdot(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\times\vec{\epsilon}_{\vec{k}\lambda})\Big{]}$
$\displaystyle=f^{abc}\frac{ig_{s}}{2}\sqrt{\frac{1}{2V\omega_{\vec{k}}\omega_{\vec{\kappa}}\omega_{\vec{k}_{1}}}}\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\cdot\Big{[}2\vec{\epsilon}_{\vec{k}\lambda}(\vec{k}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})$
$\displaystyle+2\vec{\epsilon}_{\vec{\kappa}\sigma}(\vec{\kappa}\cdot\vec{\epsilon}_{\vec{k}\lambda})-(\vec{k}+\vec{\kappa})(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})\Big{]}\delta_{\vec{\kappa},\vec{k}+\vec{k}_{1}}.$
(26)

Substituting these two matrix elements into Eq. (15) and performing the
summation over the intermediate state
$\sum_{m}=\sum_{\vec{k}_{1}}\sum_{\lambda_{1}}\sum_{d_{1}}$, we obtain the
transition amplitude for diagram (c)

$\displaystyle T_{fi}^{(c)}$
$\displaystyle=f^{abc}\frac{-ig_{s}^{2}}{V}\sqrt{\frac{\pi}{6V\omega_{\vec{k}}\omega_{\vec{\kappa}}}}e^{-i(\vec{\kappa}-\vec{k})\cdot\vec{x}}$
$\displaystyle\times\Big{\\{}(\vec{p}\cdot\vec{\epsilon}_{\vec{k}\lambda})(\vec{k}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})+(\vec{p}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})(\vec{\kappa}\cdot\vec{\epsilon}_{\vec{k}\lambda})$
$\displaystyle-\frac{(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})}{(\vec{\kappa}-\vec{k})^{2}}\Big{[}(\vec{p}\cdot\vec{k})(\vec{\kappa}^{2}-\vec{\kappa}\cdot\vec{k})+(\vec{p}\cdot\vec{\kappa})(\vec{k}^{2}-\vec{\kappa}\cdot\vec{k})\Big{]}\Big{\\}}$
$\displaystyle\times\frac{\int
r^{3}drj_{1}(pr)R_{10}(r)}{p(-\epsilon_{B}-\frac{p^{2}}{m_{Q}}-\omega_{(\vec{\kappa}-\vec{k})}+i\epsilon)}.$
(27)

where $\omega_{(\vec{\kappa}-\vec{k})}$ is the energy of the intermediate
gluon with $\vec{k}_{1}=\vec{\kappa}-\vec{k}$. To arrive at Eq. (III.2), the
completeness relation for the gluon polarization vector Eq. (9) has been used
for $\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}$.

For diagram (d), the intermediate state
$|m\rangle=|J/\psi,g(\vec{\kappa},\sigma,c),g(\vec{k_{1}},\lambda_{1},d_{1})\rangle$.
The term proportional to the combination of one annihilation operator and two
creation operators in the expansion of the $V_{3g}$ upon using Eq. (8), i.e.,

$\displaystyle
V_{3g}(a^{\dagger}a^{\dagger}a)=f^{abc}\frac{ig_{s}}{4}\sum_{\vec{k}_{1},\lambda_{1}}\sum_{\vec{k}_{2},\lambda_{2}}\sum_{\vec{k}_{3},\lambda_{3}}\sqrt{\frac{1}{2V\omega_{\vec{k}_{1}}\omega_{\vec{k}_{2}}\omega_{\vec{k}_{3}}}}$
$\displaystyle(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\times\vec{k}_{1})\cdot(\vec{\epsilon}_{\vec{k}_{2}\lambda_{2}}\times\vec{\epsilon}_{\vec{k}_{3}\lambda_{3}})\Big{(}-a_{\vec{k}_{3}\lambda_{3}}^{c\dagger}a_{\vec{k}_{1}\lambda_{1}}^{a\dagger}a_{\vec{k}_{2}\lambda_{2}}^{b}\delta_{\vec{k}_{2},\vec{k}_{3}+\vec{k}_{1}}$
$\displaystyle-
a_{\vec{k}_{1}\lambda_{1}}^{a\dagger}a_{\vec{k}_{2}\lambda_{2}}^{b\dagger}a_{\vec{k}_{3}\lambda_{3}}^{c}\delta_{\vec{k}_{3},\vec{k}_{1}+\vec{k}_{2}}+a_{\vec{k}_{2}\lambda_{2}}^{b\dagger}a_{\vec{k}_{3}\lambda_{3}}^{c\dagger}a_{\vec{k}_{1}\lambda_{1}}^{a}\delta_{\vec{k}_{1},\vec{k}_{2}+\vec{k}_{3}}\Big{)},$
(28)

is responsible for the transition matrix element due to $V_{3g}$, which reads

$\displaystyle\langle m|V_{3g}|i\rangle=\langle
g(\vec{\kappa},\sigma,c),g(\vec{k_{1}},\lambda_{1},d_{1})|V_{3g}(a^{\dagger}a^{\dagger}a)|g(\vec{k},\lambda,a)\rangle$
$\displaystyle=f^{abc}\frac{ig_{s}}{2}\sqrt{\frac{1}{2V\omega_{\vec{k}}\omega_{\vec{\kappa}}\omega_{\vec{k}_{1}}}}\delta_{\vec{k},\vec{k}_{1}+\vec{\kappa}}\Big{[}(\vec{\epsilon}_{\vec{k}\lambda}\times\vec{k})\cdot(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\times\vec{\epsilon}_{\vec{\kappa}\sigma})$
$\displaystyle+(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\times\vec{k}_{1})\cdot(\vec{\epsilon}_{\vec{k}\lambda}\times\vec{\epsilon}_{\vec{\kappa}\sigma})+(\vec{\epsilon}_{\vec{\kappa}\sigma}\times\vec{\kappa})\cdot(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\times\vec{\epsilon}_{\vec{k}\lambda})\Big{]}$
$\displaystyle=f^{abc}\frac{ig_{s}}{2}\sqrt{\frac{1}{2V\omega_{\vec{k}}\omega_{\vec{\kappa}}\omega_{\vec{k}_{1}}}}\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\cdot\Big{[}2\vec{\epsilon}_{\vec{k}\lambda}(\vec{k}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})$
$\displaystyle+2\vec{\epsilon}_{\vec{\kappa}\sigma}(\vec{\kappa}\cdot\vec{\epsilon}_{\vec{k}\lambda})-(\vec{k}+\vec{\kappa})(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})\Big{]}\delta_{\vec{k},\vec{k}_{1}+\vec{\kappa}}.$
(29)

The other transition matrix element due to $V_{SO}$ is

$\displaystyle\langle
f|V_{SO}|m\rangle=\langle(c\bar{c})_{8}(\vec{p},b)|V_{SO}|J/\psi,\,g(\vec{k}_{1},\lambda_{1},d_{1})\rangle$
$\displaystyle=g_{s}\frac{\delta^{d_{1}b}}{V}\sqrt{\frac{\pi\omega_{\vec{k}_{1}}}{3}}e^{i\vec{k}_{1}\cdot\vec{x}}(\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\cdot\vec{p})\frac{1}{p}\int
r^{3}drj_{1}(pr)R_{10}(r),$ (30)

Similar treatments as taken to arrive at Eq. (III.2) lead to the transition
amplitude for diagram (d)

$\displaystyle T_{fi}^{(d)}$
$\displaystyle=f^{abc}\frac{ig_{s}^{2}}{V}\sqrt{\frac{\pi}{6V\omega_{\vec{k}}\omega_{\vec{\kappa}}}}e^{i(\vec{k}-\vec{\kappa})\cdot\vec{x}}$
$\displaystyle\times\Big{\\{}(\vec{p}\cdot\vec{\epsilon}_{\vec{k}\lambda})(\vec{k}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})+(\vec{p}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})(\vec{\kappa}\cdot\vec{\epsilon}_{\vec{k}\lambda})$
$\displaystyle-\frac{(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{\epsilon}_{\vec{\kappa}\sigma})}{(\vec{\kappa}-\vec{k})^{2}}\Big{[}(\vec{p}\cdot\vec{k})(\vec{\kappa}^{2}-\vec{\kappa}\cdot\vec{k})+(\vec{p}\cdot\vec{\kappa})(\vec{k}^{2}-\vec{\kappa}\cdot\vec{k})\Big{]}\Big{\\}}$
$\displaystyle\times\frac{\int
r^{3}drj_{1}(pr)R_{10}(r)}{p(\omega_{\vec{k}}-\omega_{\vec{\kappa}}-\omega_{(\vec{\kappa}-\vec{k})}+i\epsilon)},$
(31)

which is of opposite sign compared to $T_{fi}^{(c)}$ in Eq. (III.2).

###  III.3 Dissociation involving $V_{SO}$ and $V_{q\bar{q}g}$

The Feynman diagrams for NLO dissociation of heavy quarkonium by thermal
quarks/antiquarks, $q/\bar{q}+J/\psi\to q/\bar{q}+c+\bar{c}$, as shown in Fig.
4, are constructed using the vertices $V_{SO}$ and $V_{q\bar{q}g}$. Here the
initial state $|i\rangle=|J/\psi,q/\bar{q}(\vec{k},r,i)\rangle$ and the final
state
$|f\rangle=|(c\bar{c})_{8}(\vec{p},b),q/\bar{q}(\vec{\kappa},s,j)\rangle$,
with $i,j=1,2,3$ labeling the quark/antiquark colors and $r,s=1,2$ their
spins. In the following, we focus on the light quark of a given flavor and the
terms in $V_{q\bar{q}g}$ (upon using expansions Eqs. (8), (II), (II)) that are
relevant for diagram (e) and (f) are

Figure 4: Feynman diagrams constructed from $V_{SO}$ and $V_{q\bar{q}g}$ for
$q/\bar{q}+J/\psi\to q/\bar{q}+c+\bar{c}$. The time direction goes upward.
Between the initial time $t_{i}$ and final time $t_{f}$ is the intermediate
state.

$\displaystyle
V_{q\bar{q}g}(b^{\dagger}ba)=g_{s}m_{q}\sum_{\vec{k},\lambda}\sum_{\vec{k}_{1},s_{1}}\sum_{\vec{k}_{2},s_{2}}\sqrt{\frac{1}{2V\omega_{\vec{k}}E_{\vec{k}_{1}}E_{\vec{k}_{2}}}}\delta_{\vec{k}_{2},\vec{k}+\vec{k}_{1}}$
$\displaystyle\times
b_{\vec{k}_{2},s_{2}}^{i_{2}\dagger}(t)u^{s_{2}\dagger}(\vec{k}_{2})\vec{\alpha}\cdot\vec{\epsilon}_{\vec{k}\lambda}a_{\vec{k}\lambda}^{a}(t)\Big{(}\frac{\lambda^{a}}{2}\Big{)}^{i_{2}i_{1}}b_{\vec{k}_{1},s_{1}}^{i_{1}}(t)u^{s_{1}}(\vec{k}_{1}),$
(32) $\displaystyle
V_{q\bar{q}g}(b^{\dagger}ba^{\dagger})=g_{s}m_{q}\sum_{\vec{k},\lambda}\sum_{\vec{k}_{1},s_{1}}\sum_{\vec{k}_{2},s_{2}}\sqrt{\frac{1}{2V\omega_{\vec{k}}E_{\vec{k}_{1}}E_{\vec{k}_{2}}}}\delta_{\vec{k}_{1},\vec{k}+\vec{k}_{2}}$
$\displaystyle\times
b_{\vec{k}_{2},s_{2}}^{i_{2}\dagger}(t)u^{s_{2}\dagger}(\vec{k}_{2})\vec{\alpha}\cdot\vec{\epsilon}_{\vec{k}\lambda}a_{\vec{k}\lambda}^{a\dagger}(t)\Big{(}\frac{\lambda^{a}}{2}\Big{)}^{i_{2}i_{1}}b_{\vec{k}_{1},s_{1}}^{i_{1}}(t)u^{s_{1}}(\vec{k}_{1}).$
(33)

For diagram (e), the intermediate state
$|m\rangle=|(c\bar{c})_{8}(\vec{p},b),q(\vec{k},r,i),g(\vec{k_{1}},\lambda_{1},d_{1})\rangle$.
The transition matrix element $\langle m|V_{SO}|i\rangle$ is same as Eq.
(III.2). The other transition matrix element due to $V_{q\bar{q}g}$ reads

$\displaystyle\langle f|V_{q\bar{q}g}|m\rangle=\langle
q(\vec{\kappa},s,j)|V_{q\bar{q}g}(b^{\dagger}ba)|g(\vec{k}_{1},\lambda_{1},d_{1}),q(\vec{k},r,i)\rangle$
$\displaystyle=g_{s}m_{q}\sqrt{\frac{1}{2V\omega_{\vec{k}_{1}}E_{\vec{k}}E_{\vec{\kappa}}}}u^{s\dagger}(\vec{\kappa})\vec{\alpha}\cdot\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\Big{(}\frac{\lambda^{d_{1}}}{2}\Big{)}^{ji}u^{r}(\vec{k})\delta_{\vec{\kappa},\vec{k}+\vec{k}_{1}}.$
(34)

The transition amplitude for diagram (e) is then computed from Eq. (15) with
$\sum_{m}=\sum_{\vec{k}_{1}}\sum_{\lambda_{1}}\sum_{d_{1}}$:

$\displaystyle T_{fi}^{(e)}$ $\displaystyle=\sum_{m}\frac{\langle
f|V_{q\bar{q}g}|m\rangle\langle m|V_{SO}|i\rangle}{E_{i}-E_{m}+i\epsilon}$
$\displaystyle=-\frac{g_{s}^{2}m_{q}}{V}\sqrt{\frac{\pi}{6VE_{\vec{k}}E_{\vec{\kappa}}}}e^{-i(\vec{\kappa}-\vec{k})\cdot\vec{x}}\Big{[}u^{s\dagger}(\vec{\kappa})\vec{\alpha}\cdot\vec{p}u^{r}(\vec{k})$
$\displaystyle-\frac{[u^{s\dagger}(\vec{\kappa})(\vec{\kappa}-\vec{k})\cdot\vec{\alpha}u^{r}(\vec{k})][(\vec{\kappa}-\vec{k})\cdot\vec{p}]}{|\vec{\kappa}-\vec{k}|^{2}}\Big{]}\Big{(}\frac{\lambda^{b}}{2}\Big{)}^{ji}$
$\displaystyle\times\frac{\int
r^{3}drj_{1}(pr)R_{10}(r)}{p(-\epsilon_{B}-\frac{p^{2}}{m_{Q}}-\omega_{(\vec{\kappa}-\vec{k})}+i\epsilon)}.$
(35)

To arrive at Eq. (III.3), the completeness relation for the gluon polarization
vector Eq. (9) has been again used for
$\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}$.

Finally for diagram (f), the intermediate state
$|m\rangle=|J/\psi,q(\vec{\kappa},s,j),g(\vec{k_{1}},\lambda_{1},d_{1})\rangle$.
The transition matrix element due to $V_{q\bar{q}g}$ reads

$\displaystyle\langle m|V_{q\bar{q}g}|i\rangle=\langle
q(\vec{\kappa},s,j),\,g(\vec{k_{1}},\lambda_{1},d_{1})|V_{q\bar{q}g}(b^{\dagger}ba^{\dagger})|q(\vec{k},r,i)\rangle$
$\displaystyle=g_{s}m_{q}\sqrt{\frac{1}{2V\omega_{\vec{k}_{1}}E_{\vec{k}}E_{\vec{\kappa}}}}u^{s\dagger}(\vec{\kappa})\vec{\alpha}\cdot\vec{\epsilon}_{\vec{k}_{1}\lambda_{1}}\Big{(}\frac{\lambda^{d_{1}}}{2}\Big{)}^{ji}u^{r}(\vec{k})\delta_{\vec{k},\vec{k}_{1}+\vec{\kappa}}.$
(36)

The other matrix element due to $V_{SO}$ ($\langle f|V_{SO}|m\rangle$) is same
as Eq. (III.2). The transition amplitude for diagram (f) is obtained by
combining these two matrix elements into Eq. (15) and carrying out the
summation over intermediate states
$\sum_{m}=\sum_{\vec{k}_{1}}\sum_{\lambda_{1}}\sum_{d_{1}}$, such that

$\displaystyle T_{fi}^{(f)}$
$\displaystyle=+\frac{g_{s}^{2}m_{q}}{V}\sqrt{\frac{\pi}{6VE_{\vec{k}}E_{\vec{\kappa}}}}e^{i(\vec{k}-\vec{\kappa})\cdot\vec{x}}\Big{[}u^{s\dagger}(\vec{\kappa})\vec{\alpha}\cdot\vec{p}u^{r}(\vec{k})$
$\displaystyle-\frac{[u^{s\dagger}(\vec{\kappa})(\vec{k}-\vec{\kappa})\cdot\vec{\alpha}u^{r}(\vec{k})][(\vec{k}-\vec{\kappa})\cdot\vec{p}]}{|\vec{k}-\vec{\kappa}|^{2}}\Big{]}\Big{(}\frac{\lambda^{b}}{2}\Big{)}^{ji}$
$\displaystyle\times\frac{\int
r^{3}drj_{1}(pr)R_{10}(r)}{p(E_{\vec{k}}-E_{\vec{\kappa}}-\omega_{(\vec{\kappa}-\vec{k})}+i\epsilon)},$
(37)

which is of opposite sign relative to $T_{fi}^{(e)}$ in Eq. (III.3).

The transition amplitudes for $\bar{q}+J/\psi\to\bar{q}+c+\bar{c}$ are
obtained upon replacement of spinors in Eqs. (III.3) and (III.3) via
$u_{s}^{\dagger}(\vec{\kappa})\rightarrow v_{r}^{\dagger}(\vec{k})$,
$u_{r}(\vec{k})\rightarrow v_{s}(\vec{\kappa})$.

###  III.4 Dissociation cross sections

Using the second-order transition amplitudes calculated above, the transition
rates for the NLO dissociation processes are calculated from Fermi’s golden
rule Sakurai-AQM:2008 ; Greiner:1998 ; Gottfried-Yan:2003

$\displaystyle W_{i\to f}=2\pi|T_{fi}|^{2}\delta(E_{i}-E_{f}).$ (38)

The cross section is obtained by dividing the transition rate by the flux of
incident parton $v_{rel}/V$ ($v_{rel}$ being the relative velocity between the
incident parton and the $J/\psi$ target) and averaging (summing) over the
initial (final) state degeneracies.

We note that both $T_{fi}^{(a)}$, $T_{fi}^{(b)}\propto d^{abc}$, whereas both
$T_{fi}^{(c)}$, $T_{fi}^{(d)}\propto f^{abc}$. The totally symmetric and
totally antisymmetry property of $d^{abc}$ and $f^{abc}$, respectively,
implies $\sum_{abc}d^{abc}f^{abc}=0$, so that there’s no interference between
$T_{fi}^{(a)}+T_{fi}^{(b)}$ and $T_{fi}^{(c)}+T_{fi}^{(d)}$ and they have
independent cross sections. For diagrams (a) and (b), the summation and
averaging procedure gives the cross section

$\displaystyle\sigma^{(a+b)}(E_{g})=2\pi V\frac{V}{(2\pi)^{3}}\int
d^{3}p\sum_{b}\frac{V}{(2\pi)^{3}}\int d^{3}\kappa\sum_{\sigma}\sum_{c}$
$\displaystyle\times\frac{1}{4\pi}\int
d\Omega_{\vec{k}}\frac{1}{2}\sum_{\lambda}\frac{1}{8}\sum_{a}|T_{fi}^{(a)}+T_{fi}^{(b)}|^{2}\delta(E_{i}-E_{f})$
$\displaystyle=\frac{5}{216\pi^{2}}g_{s}^{4}E_{g}\int
p^{2}dp\int\kappa^{2}\omega_{\vec{\kappa}}d\kappa\\{\cdots\\}$
$\displaystyle\times\delta(-\epsilon_{B}+\omega_{\vec{k}}-\frac{p^{2}}{m_{Q}}-\omega_{\vec{\kappa}}),$
(39)

where $E_{g}=\omega_{\vec{k}}$ is the incident gluon energy and

$\displaystyle\\{\cdots\\}=$ $\displaystyle
A^{2}(p,k)+C^{2}(p,\kappa)+2A(p,k)C(p,\kappa)$
$\displaystyle+\frac{1}{3}\Big{\\{}B^{2}(p,k)+D^{2}(p,\kappa)+2\Big{[}A(p,k)B(p,k)$
$\displaystyle+A(p,k)D(p,\kappa)+B(p,k)C(p,\kappa)$
$\displaystyle+B(p,k)D(p,\kappa)+C(p,\kappa)D(p,\kappa)\Big{]}\Big{\\}}.$ (40)

To arrive at second equality of Eq. (III.4), two identities regarding the
totally symmetric $SU(3)_{c}$ group invariant $\sum_{abc}d^{abc}d^{abc}=40/3$
and the polarization vector ($\vec{\rho}$ being an arbitrary vector
independent of $\Omega_{\vec{k}}$)

$\displaystyle\frac{1}{4\pi}\int
d\Omega_{\vec{k}}\frac{1}{2}\sum_{\lambda=1,2}|\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{\rho}|^{2}=\frac{1}{3}|\vec{\rho}|^{2}$
(41)

have been used. Note that upon using the energy-conserving $\delta$ function
(the last line in Eq. (III.4)), the denominators of $A(p,k)$, $B(p,k)$,
$C(p,\kappa)$ and $D(p,\kappa)$ defined in Eqs. (III.1) and (III.1) would be
$\propto\omega_{\vec{\kappa}},~{}\omega_{\vec{k}}$, respectively, which are
simply the corresponding momentum magnitudes if the gluons are massless and
thus lead to infrared divergences. In practice, these divergences are
regularized and removed by the finite thermal gluon mass in the QGP. For
finite gluon mass $m_{g}$, we define
$p_{c}=\sqrt{(\omega_{\vec{k}}-\epsilon_{B}-m_{g})m_{Q}}$. Apparently if
$p>p_{c}$, one has
$-\epsilon_{B}+\omega_{\vec{k}}-p^{2}/m_{Q}<m_{g}\leq\omega_{\vec{\kappa}}$,
such that the zero point of the argument of the $\delta$ function can never be
reached for any $\kappa>0$, and the corresponding integrand in Eq. (III.4)
vanishes. Therefore, the $p_{c}$ acts as an cutoff for the integration over
$p$. On the other hand, for $p<p_{c}$, after integrating out the $\delta$
function via $\int d\kappa$, $\int_{0}^{p_{c}}dp$ yields a finite result for
the cross section.

For diagrams (c) and (d), similar procedure gives the cross section

$\displaystyle\sigma^{(c+d)}(E_{g})=2\pi V\frac{V}{(2\pi)^{3}}\int
d^{3}p\sum_{b}\frac{V}{(2\pi)^{3}}\int d^{3}\kappa\sum_{\sigma}\sum_{c}$
$\displaystyle\times\frac{1}{4\pi}\int
d\Omega_{\vec{k}}\frac{1}{2}\sum_{\lambda}\frac{1}{8}\sum_{a}|T_{fi}^{(c)}+T_{fi}^{(d)}|^{2}\delta(E_{i}-E_{f})$
$\displaystyle=\frac{g_{s}^{4}}{32\pi^{4}}\frac{1}{E_{g}}\int
dp\int\kappa^{2}\frac{1}{\omega_{\vec{\kappa}}}d\kappa\int
d\Omega_{\vec{k}}\int d\Omega_{\vec{\kappa}}$ $\displaystyle\Big{[}\int
r^{3}drj_{1}(pr)R_{10}(r)\Big{]}^{2}g(\vec{\kappa},\vec{p},\vec{k})\Big{[}\frac{\omega_{\vec{\kappa}}-\omega_{\vec{k}}}{\omega_{(\vec{\kappa}-\vec{k})}^{2}-(\omega_{\vec{\kappa}}-\omega_{\vec{k}})^{2}}\Big{]}^{2}$
$\displaystyle\times\delta(-\epsilon_{B}+\omega_{\vec{k}}-\frac{p^{2}}{m_{Q}}-\omega_{\vec{\kappa}}),$
(42)

where $\sum_{abc}f^{abc}f^{abc}=24$ is used and the polynomial function
$g(\vec{\kappa},\vec{p},\vec{k})$ arises from handling the summation over
gluon polarizations using the completeness relation Eq. (9):

$\displaystyle
g(\vec{\kappa},\vec{p},\vec{k})=\sum_{\sigma}\sum_{\lambda}[(\vec{p}\cdot{\vec{\epsilon}}_{\vec{k}\lambda})(\vec{k}\cdot{\vec{\epsilon}}_{\vec{\kappa}\sigma})+(\vec{p}\cdot{\vec{\epsilon}}_{\vec{\kappa}\sigma})(\vec{\kappa}\cdot{\vec{\epsilon}}_{\vec{k}\lambda})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-({\vec{\epsilon}}_{\vec{k}\lambda}\cdot{\vec{\epsilon}}_{\vec{\kappa}\sigma})X(\vec{\kappa},\vec{p},\vec{k})]^{2}$
$\displaystyle=({\vec{p}}^{2}-\frac{{(\vec{p}\cdot\vec{k})}^{2}}{{\vec{k}}^{2}})({\vec{k}}^{2}-\frac{{(\vec{\kappa}\cdot\vec{k})}^{2}}{{\vec{\kappa}}^{2}})$
$\displaystyle+({\vec{p}}^{2}-\frac{{(\vec{p}\cdot\vec{\kappa})}^{2}}{{\vec{\kappa}}^{2}})({\vec{\kappa}}^{2}-\frac{{(\vec{\kappa}\cdot\vec{k})}^{2}}{{\vec{k}}^{2}})$
$\displaystyle+(1+\frac{{(\vec{\kappa}\cdot\vec{k})}^{2}}{{\vec{\kappa}}^{2}{\vec{k}}^{2}})X^{2}(\vec{\kappa},\vec{p},\vec{k})$
$\displaystyle+2(\vec{p}\cdot\vec{\kappa}-\frac{(\vec{p}\cdot\vec{k})(\vec{\kappa}\cdot\vec{k})}{{\vec{k}}^{2}})(\vec{p}\cdot\vec{k}-\frac{(\vec{p}\cdot\vec{\kappa})(\vec{\kappa}\cdot\vec{k})}{{\vec{\kappa}}^{2}})$
$\displaystyle+2\frac{\vec{\kappa}\cdot\vec{k}}{{\vec{\kappa}}^{2}}(\vec{p}\cdot\vec{\kappa}-\frac{(\vec{p}\cdot\vec{k})(\vec{\kappa}\cdot\vec{k})}{{\vec{k}}^{2}})X(\vec{\kappa},\vec{p},\vec{k})$
$\displaystyle+2\frac{\vec{\kappa}\cdot\vec{k}}{{\vec{k}}^{2}}(\vec{p}\cdot\vec{k}-\frac{(\vec{p}\cdot\vec{\kappa})(\vec{k}\cdot\vec{\kappa})}{{\vec{\kappa}}^{2}})X(\vec{\kappa},\vec{p},\vec{k})$
(43)

with

$\displaystyle
X(\vec{\kappa},\vec{p},\vec{k})=\frac{(\vec{p}\cdot\vec{k})({\vec{\kappa}}^{2}-\vec{\kappa}\cdot\vec{k})+(\vec{p}\cdot\vec{\kappa})({\vec{k}}^{2}-\vec{\kappa}\cdot\vec{k})}{(\vec{\kappa}-\vec{k})^{2}}.$
(44)

When the incident and outgoing gluons are collinear, the denominator in the
integrand of Eq. (III.4)
$\omega_{(\vec{\kappa}-\vec{k})}^{2}-(\omega_{\vec{\kappa}}-\omega_{\vec{k}})^{2}=-(\omega_{\vec{k}}-\omega_{\vec{\kappa}}-\omega_{(\vec{\kappa}-\vec{k})})(\omega_{\vec{k}}-\omega_{\vec{\kappa}}+\omega_{(\vec{\kappa}-\vec{k})})$
would be vanishing for massless gluons
($|\vec{k}|-|\vec{\kappa}|-|\vec{k}-\vec{\kappa}|=0$ for
$\vec{k}//\vec{\kappa}$) and thus a divergence occurs, which will be
practically removed by finite gluon masses. In evaluating Eq. (III.4) with Eq.
(III.4), we have chosen the momentum $\vec{p}$ along the third axis in the
momentum space spanned by $\vec{k}$ and $\vec{\kappa}$.

Using the completeness relation for spinors Eq. (13), the summation and
averaging procedure can be performed for the combination of scattering
amplitudes from diagrams (e) and (f), to obtain the cross section for
$q+J/\psi\to q+c+\bar{c}$

$\displaystyle\sigma^{(e+f)}(E_{q})=2\pi V\frac{V}{(2\pi)^{3}}\int
d^{3}p\sum_{b}\frac{V}{(2\pi)^{3}}\int d^{3}\kappa\sum_{s}\sum_{j}$
$\displaystyle\times\frac{1}{4\pi}\int
d\Omega_{\vec{k}}\frac{1}{2}\sum_{r}\frac{1}{3}\sum_{i}|T_{fi}^{(e)}+T_{fi}^{(f)}|^{2}\delta(E_{i}-E_{f})$
$\displaystyle=\frac{g_{s}^{4}}{72\pi^{4}}\frac{1}{E_{q}}\int
dp\int\kappa^{2}\frac{1}{E_{\vec{\kappa}}}d\kappa\int d\Omega_{\vec{k}}\int
d\Omega_{\vec{\kappa}}$ $\displaystyle\Big{[}\int
r^{3}drj_{1}(pr)R_{10}(r)\Big{]}^{2}f(\vec{\kappa},\vec{p},\vec{k})\Big{[}\frac{E_{\vec{\kappa}}-E_{\vec{k}}}{\omega_{(\vec{\kappa}-\vec{k})}^{2}-(E_{\vec{\kappa}}-E_{\vec{k}})^{2}}\Big{]}^{2}$
$\displaystyle\times\delta(-\epsilon_{B}+E_{\vec{k}}-\frac{p^{2}}{m_{Q}}-E_{\vec{\kappa}}),$
(45)

where $E_{q}=E_{\vec{k}}$ and the polynomial function

$\displaystyle
f(\vec{\kappa},\vec{p},\vec{k})=m_{q}^{2}\sum_{r}\sum_{s}|u^{s\dagger}(\vec{\kappa})\vec{\alpha}\cdot\vec{p}u^{r}(\vec{k})$
$\displaystyle-\frac{[u^{s\dagger}(\vec{\kappa})(\vec{\kappa}-\vec{k})\cdot\vec{\alpha}u^{r}(\vec{k})][(\vec{\kappa}-\vec{k})\cdot\vec{p}]}{(\vec{\kappa}-\vec{k})^{2}}|^{2}$
$\displaystyle=2(\vec{p}\cdot\vec{k})(\vec{p}\cdot\vec{\kappa})+\Big{[}\frac{(\vec{\kappa}-\vec{k})\cdot\vec{p}}{(\vec{\kappa}-\vec{k})^{2}}\Big{]}^{2}2(\vec{\kappa}-\vec{k})\cdot\vec{k}(\vec{\kappa}-\vec{k})\cdot\vec{\kappa}$
$\displaystyle-2\frac{(\vec{\kappa}-\vec{k})\cdot\vec{p}}{(\vec{\kappa}-\vec{k})^{2}}\Big{[}(\vec{\kappa}-\vec{k})\cdot\vec{k}(\vec{p}\cdot\vec{\kappa})+(\vec{\kappa}-\vec{k})\cdot\vec{\kappa}(\vec{p}\cdot\vec{k})\Big{]}$
$\displaystyle+(\vec{k}\cdot\vec{\kappa}-E_{\vec{k}}E_{\vec{\kappa}}+m_{q}^{2})\Big{[}-\vec{p}^{2}+\frac{[(\vec{\kappa}-\vec{k})\cdot\vec{p}]^{2}}{(\vec{\kappa}-\vec{k})^{2}}\Big{]},$
(46)

When the outgoing quark moves in parallel to the incident quark,
$\vec{\kappa}//\vec{k}$, and all partons are massless, Eq. (III.4) would have
a collinear divergence similar to that identified for Eq. (III.4), which will
be practically cured by finite thermal parton masses. The cross section for
$\bar{q}+J/\psi\to\bar{q}+c+\bar{c}$ reads the same.

##  IV Dissociation of various heavy quarkonia in an in-medium potential model

###  IV.1 NLO dissociation cross sections

Medium effects enter the dissociation cross sections derived in Sec. III.4
through the temperature dependence of bound state wave functions, binding
energies and thermal masses of light partons involved. For the in-medium
radial wave functions and binding energies of various heavy quarkonium bound
states, we have solved the radial Schr$\ddot{o}$digner equation for the
$Q\bar{Q}$ system Karsch:1987pv ; Chen:2017jje

$\displaystyle\frac{1}{r^{2}}\frac{d}{dr}(r^{2}\frac{dR_{nl}}{dr})+$
$\displaystyle[m_{Q}(E_{n,l}-2m_{Q}-V_{1}(r,T))-\frac{l(l+1)}{r^{2}}]R_{nl}(r)=0,$
(47)

with a temperature dependent singlet potential Karsch:1987pv

$V_{1}(r,T)=-\frac{\alpha}{r}e^{-\mu(T)r}+\frac{\sigma}{\mu(T)}(1-e^{-\mu(T)r}),$
(48)

where $\mu/T=-4.058+6.32\cdot(T/T_{c}-0.885)^{0.1035}$ (critical temperature
$T_{c}$=172.5MeV) parameterizes the screening of the potential Burnier:2015tda
. The Eq. (48) represents a modification of the vacuum Cornell potential
$V_{1}(r,0)=-\alpha/r+\sigma r$ by color screening, which reproduces well the
vacuum masses of various charmonia and bottomonia below threshold with
coupling strength $\alpha=4/3\alpha_{s}=0.471$, string tension $\sigma=0.192$
GeV2 and charm and bottom quark masses $m_{c}=1.320$ GeV and $m_{b}=4.746$
GeV, respectively Karsch:1987pv . The temperature-dependent binding energy of
a bound state is then obtained via
$\epsilon_{B}(T)=2m_{Q}+\sigma/\mu(T)-E_{n,l}(T)$, whose zero point defines
the melting temperature for the bound state Karsch:1987pv ; Chen:2017jje . For
example, the melting temperatures $T_{d}\simeq 1.6T_{c},~{}1.2T_{c}$ for
$J/\psi$, $\chi_{c}$ (and $\psi(2S)$), respectively, and $T_{d}\simeq
2.6T_{c},~{}1.5T_{c},~{}1.3T_{c}$ for $\Upsilon(1S)$, $\Upsilon(2S)$ (and
$\chi_{b}(1P)$) and $\Upsilon(3S)$, respectively.

For the thermal masses of gluons and light quarks/antiquarks that help to
remove the infrared and collinear divergences as discussed in Sec. III.4, we
use their perturbative values $m_{g}(T)=\sqrt{3/4}g_{s}T$ and
$m_{q}(T)=\sqrt{1/3}g_{s}T$ with fixed $g_{s}=2.3$ for three active light
flavors Cao:2018ews . However, it was noted Gossiaux:2008jv ; Zhao:2023ucp
that to mimic the same partonic energy loss from some nonperturbative
calculations Braaten:1991jj  with a Born diagram, the screening on the
exchanged gluon propagator in the $t$-channel scattering should be much
reduced. In this spirit, we use a significantly reduced effective mass
$m_{g_{\rm ex}}(T)=rg_{s}T$, with $r=0.2$ Zhao:2023ucp , for the intermediate
gluon associated with momentum $\vec{k}-\vec{\kappa}$ in diagrams (c), (d),
(e) and (f), that would correspond to the exchanged gluon in the $t$-channel
Born diagram in relativistically covariant perturbation theory, such that
$\omega_{(\vec{k}-\vec{\kappa})}=\sqrt{(\vec{k}-\vec{\kappa})^{2}+m_{g_{\rm
ex}}^{2}}$ in Eqs. (III.4) and (III.4). In this way, we expect to partially
account for the pertinent nonperturabtive effects in the current perturbative
calculations.

Figure 5: NLO dissociation cross sections for charmonia at varying
temperatures. For (a) $J/\psi$, $\sigma^{(a+b)}$, $\sigma^{(c+d)}$ and
$\sigma^{(e+f)}$ are displayed separately, and LO (gluo-dissociation) cross
sections taken from Chen:2017jje  are also shown for comparison. For (b)
$\psi(2S)$ and (c) $\chi_{c}$, NLO dissociation cross sections by gluons
($\sigma^{(a+b)}+\sigma^{(c+d)}$) and quarks ($\sigma^{(e+f)}$) are separately
shown.

Figure 6: Same as Fig. 5 but for bottomonia: (a) $\Upsilon(1S)$, (b)
$\Upsilon(2S)$ and (c) $\chi_{b}(1P)$.

The numerical results of NLO dissociation cross sections for various charmonia
and bottomonia are displayed in Fig. 5 and Fig. 6, respectively, as a function
of the incident parton energy at varying temperatures. A general observation
is that the NLO cross sections first grow with the incident parton energy and
then saturate toward large energies, in marked contrast to the LO (gluo-
dissocation) cross sections which exhibit a pronounced peak at energy
(wavelength) matching the bound state binding energy (size) and quickly drop
off thereafter, as shown in Fig. 5(a) and Fig. 6(a) for the $1S$ states
$J/\psi$ and $\Upsilon(1S)$, respectively. This is presumably attributed to
the outgoing parton in the final state of the NLO processes carrying away the
excess energy, thereby overcoming the mismatch between the incident parton
wavelength and the bound state size in the LO process. Among the three types
of NLO dissociation cross sections as derived in Sec. III.4, $\sigma^{(a+b)}$
turns out much smaller than $\sigma^{(c+d)}$ and $\sigma^{(e+f)}$, as
demonstrated for $J/\psi$ and $\Upsilon(1S)$. This may be expected from the
fact that what is exchanged in $s$\- and $u$-channel-like diagrams (a) and (b)
(Fig. 2) is a very massive $(Q\bar{Q})_{8}$ octet, whereas underlying
$\sigma^{(c+d)}$ and $\sigma^{(e+f)}$ are the $t$-channel-like one-gluon-
exchange diagrams (c) and (d) (Fig. 3), and (e) and (f) (Fig. 4). As
temperature increases, the bound state wave functions (dipole size) expand and
the binding energies decrease. While the latter is responsible for the shift
of the onset point of cross sections toward lower incident parton energy, the
former is expected to enhance the magnitude of cross sections. However, this
naive expectation is quantitatively offset by the suppression from the
increasing effective mass of the exchanged gluon $m_{g_{\rm ex}}$ (as well as
the decreasing binding energy parameter entering the cross section formulas,
as we’ve numerically checked), leading to an overall reduction of cross
sections toward higher temperatures.

Comparing different heavy quarkonium states, the NLO cross sections are
roughly ordered by their vacuum binding energies: more loosely bound excited
states (a dipole of larger size) are more easily dissociated by thermal
partons and thus possess larger dissociation cross sections. The NLO
dissociation cross sections for Bottomonia are consistently smaller than those
for charmonia. The $\Upsilon(1S)$ is the most tightly bound state, for which
the color dipole coupling mechanism underlying our calculations may be most
applicable. Furthermore, the technical approximation made in our calculations
that the rest frame of the heavy quarkonium is also considered to be the rest
frame of the unbound octet $(Q\bar{Q})_{8}$ in the final state, i.e., the
recoil effect on the bound state by the incident parton is neglected, should
be more reasonable for the much more massive bottomonia. Therefore, our
results of the NLO dissociation cross sections may be deemed most reliable for
the $\Upsilon(1S)$. To compare with the corresponding result shown in
Hong:2018vgp , where the NLO cross sections for $\Upsilon(1S)$ were calculated
from another perturbative approach using an effective vertex constructed from
the non-relativistic Bethe-Salpeter amplitude, we choose a comparison point of
$T\sim 300$ MeV and incident parton energy $\sim 1.5$ GeV, at which our value
of $\sigma(g+\Upsilon(1S)\to g+b+\bar{b})\sim 0.2$ mb is a factor of $\sim 2$
larger, which might be in part due to the use of a reduced effective mass
($m_{g_{\rm ex}}$) for the exchanged gluon in our calculations. However, the
cross section for $\Upsilon(1S)$ shown in Hong:2018vgp  keeps increasing
rapidly with the incident parton energy, differing remarkably from the near-
saturation of the NLO cross section toward higher energies as found in our
calculations.

###  IV.2 Dissociation rates

The dissociation rates, serving as inputs for the phenomenological modelling
of heavy quarknonium transport in the QGP Zhao:2010nk ; Song:2011nu ;
Strickland:2011mw ; Zhou:2014kka ; Du:2017qkv ; Brambilla:2020qwo ; Wu:2024gil
, are obtained by folding the dissociation cross section with the distribution
function of the incident parton. For a bound state sitting at rest in the QGP,
the dissociation rate reads

$\displaystyle\Gamma(T)=d_{g,q}\int\frac{d^{3}k}{(2\pi)^{3}}f_{g,q}(E_{g,q}(\vec{k}),T)v_{\rm
rel}\sigma(E_{g,q},T),$ (49)

where $d_{g}=2\cdot 8=16$ for thermal gluons and $d_{q}=3\cdot 2\cdot 3\cdot
2=36$ for thermal ligth quarks/antiquarks is the incident parton degeneracy,
$v_{\rm rel}$ the relative velocity between the incident parton and the heavy
quarkonium at rest, and $f_{g,q}(E_{g,q},T)=1/(e^{(E_{g,q}(\vec{k})/T)}\mp 1)$
the Bose ($-$)/Fermi ($+$) distribution with
$E_{g,q}=\sqrt{\vec{k}^{2}+m_{g,q}^{2}(T)}$ for gluons and light
quarks/antiquarks, respectively. We note that for typical temperatures (a
couple of hundreds of MeV) reached in the QGP created in current relativistic
heavy-ion collisions, the thermal parton energy is on the order of $\leq 1$
GeV, which is much smaller than the mass of charmonia and bottomonia.
Therefore, in calculating the dissociation rate, thermal partons are generally
not probing the NLO cross sections at very high energies, which also renders
the approximation of neglecting the recoil effect on the heavy quarkonium
relatively safe.

Figure 7: (a) NLO dissociation rates for $J/\psi$ at finite temperatures;
contributions from different NLO processes and from LO (gluo-dissociation)
Chen:2017jje  process are separately displayed for comparison. (b) Total NLO
dissociation rates for various charmonia.

Figure 8: (a) NLO dissociation rates for $\Upsilon(1S)$ at finite
temperatures; contributions from different NLO processes and from LO (gluo-
dissociation) process are separately displayed for comparison. (b) Total NLO
dissociation rates for various bottomonia.

The calculated dissociation rates for charmonia and bottomonia are presented
in Figs. 7 and 8 up to their respective melting temperatures. The dissociation
rates from different NLO processes are displayed separately and compared to
the counterpart from LO (gluo-dissociation) process (taken from Chen:2017jje )
for the $1S$ groud states $J/\psi$ (Fig. 7(a)) and $\Upsilon(1S)$ (Fig. 8(a)).
At low temperatures (up to $\sim 1.2T_{c}$ and $\sim 1.6T_{c}$ for $J/\psi$
and $\Upsilon(1S)$, respectively), the ground states are still sufficiently
tightly bound and the incident parton of long enough wavelength does not
resolve the substructure of the bound state, so that the LO dissociation
remains more effective and its dissociation rates dominate over the total NLO
results. As temperature increases, the peak of the LO dissociation cross
section shifts toward lower energies (cf. Fig. 5(a) and Fig. 6(a)), where the
incident parton has a smaller and smaller phase space Rapp:2008tf , leading to
decreasing dissociation rates. In contrast, the NLO dissociation rates keep
growing with temperature and take over from the LO counterparts toward high
temperatures. The total NLO dissociation rates for $\Upsilon(1S)$ from our
calculation turn out to be quite comparable to corresponding results in
Hong:2018vgp  as computed from another approach. However, we note that in
Hong:2018vgp , the dissociation rates were calculated using massless incident
parton (therefore of significantly higher density) with a smaller NLO cross
section, whereas we have taken into account the thermal masses of the partons
colliding with the heavy quarkonium.

The total NLO dissociation rates for the excited charmonia and bottomonia as
shown in Fig. 7(b) and Fig. 8(b), respectively, are generally larger than
their ground state counterparts at the same temperatures, in accord with the
larger dissociation cross sections for these states. Here, we note that for
$\psi(2S)$ and $\Upsilon(2S)$, the NLO dissociation rates turn out to drop off
when approaching their respective melting temperatures, in contrast to the
monotonously increasing behavior for the $1S$ and $1P$ states. This is due to
the fast decline of their NLO cross sections with temperature (cf. Fig. 5(b)
and Fig. 6(b)), which, in turn, can be technically attributed to the peculiar
behavior of their radial functions having a node, as we’ve numerically
checked.

##  V Second-order transition between different eigen-bound states

###  V.1 Deriving the transition cross sections

In addition to the NLO dissociation, the heavy quarkonium may experience
quantum transition to another eigen-bound state in the QGP (including
excitation and elastic scattering), which also occurs at the second order as
required by color neutrality, by first absorbing and then emitting a gluon.
This process may also contribute to the total thermal decay rates for heavy
quarknoium as evaluated in lattice QCD Larsen:2019zqv . We compute the cross
sections and rates for these second-order transition processes in the same
framework as for the NLO dissociation. Using the vertex $V_{SO}$ twice, the
Feynamn diagrams for the transition $g+\Psi\to g+\Psi^{\prime}$ ($\Psi$ and
$\Psi^{\prime}$ denote two eigen-bound states of heavy quarkonium) can be
constructed, as illustrated in Fig. 9 for $g+J/\psi\to g+\psi(2S)$, which
looks rather similar to the dissociation diagrams in Fig. 2, except that now
the final state involves another bound state instead of the unbound octet. We
note that a $s$-wave bound state can only transition into another $s$-wave or
$d$-wave bound state, i.e., the change of orbital angular momentum is
constrained to be even, as dictated by the selection rules ($\Delta l=1,\Delta
s=0$) for the color-electric dipole coupling Chen:2017jje .

Figure 9: Feynman diagrams for $g+J/\psi\to g+\psi(2S)$. The time direction
goes upward. Between the initial time $t_{i}$ and final time $t_{f}$ is the
intermediate state.

For the transition represented in Fig. 9, the initial state
$|i\rangle=|J/\psi,g(\vec{k},\lambda,a)\rangle$ and final state
$|f\rangle=|\psi(2S),g(\vec{\kappa},\sigma,c)\rangle$. For diagram (g), the
intermediate state reads $|m\rangle=|(c\bar{c})_{8}(\vec{q},d)\rangle$. The
transition matrix element $\langle m|V_{SO}|i\rangle$ is the same as the Eq.
(III.1), and the other one

$\displaystyle\langle
f|V_{SO}|m\rangle=\langle\psi(2S),g(\vec{\kappa},\sigma,c)|V_{SO}|(c\bar{c})_{8}(\vec{q},d)\rangle$
$\displaystyle=g_{s}\frac{\delta^{dc}}{V}\sqrt{\frac{\pi\omega_{\vec{\kappa}}}{3}}e^{-i\vec{\kappa}\cdot\vec{x}}(\vec{\epsilon}_{\vec{\kappa}\sigma}\cdot\vec{q})\frac{1}{q}\int
r^{3}drj_{1}(qr)R_{20}(r).$ (50)

By carrying out the summation over intermediate state
$\sum_{m}=\sum_{\vec{q}}\sum_{d}=\frac{V}{(2\pi)^{3}}\int d^{3}q\sum_{d}$, the
transition amplitude is obtained

$\displaystyle T_{fi}^{(g)}$ $\displaystyle=\sum_{m}\frac{\langle
f|V_{SO}|m\rangle\langle m|V_{SO}|i\rangle}{E_{i}-E_{m}+i\epsilon}$
$\displaystyle=\delta^{ac}\frac{g_{s}^{2}\pi}{3V}\sqrt{\omega_{\vec{\kappa}}\omega_{\vec{k}}}e^{i(\vec{k}-\vec{\kappa})\cdot\vec{x}}\int\frac{d^{3}q}{(2\pi)^{3}}\Big{[}\frac{(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{q})(\epsilon_{\vec{\kappa}\sigma}\cdot\vec{q})}{q^{2}}$
$\displaystyle\times\frac{\int r^{3}drj_{1}(qr)R_{20}(r)\times\int
r^{3}drj_{1}(qr)R_{10}(r)}{-\epsilon_{B}^{J/\psi}+\omega_{\vec{k}}-\frac{q^{2}}{m_{Q}}+i\epsilon}\Big{]}.$
(51)

For diagram (h), the intermediate state is
$|m\rangle=|(c\bar{c})_{8}(\vec{q},d),\,g(\vec{k},\lambda,a),\,g(\vec{\kappa},\sigma,c)\rangle$.
The transition matrix element $\langle m|V_{SO}|i\rangle$ is the same as Eq.
(III.1), and the other one reads

$\displaystyle\langle
f|V_{SO}|m\rangle=\langle\psi^{\prime}|V_{SO}|(c\bar{c})_{8}(\vec{q},d),\,g(\vec{k},\lambda,a)\rangle$
$\displaystyle=(-)g_{s}\frac{\delta^{da}}{V}\sqrt{\frac{\pi\omega_{\vec{k}}}{3}}e^{i\vec{k}\cdot\vec{x}}(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{q})\frac{1}{q}\int
r^{3}drj_{1}(qr)R_{20}(r).$ (52)

Upon summing over the intermediate states, the transition amplitude reads

$\displaystyle T_{fi}^{(h)}$
$\displaystyle=\delta^{ac}\frac{g_{s}^{2}\pi}{3V}\sqrt{\omega_{\vec{\kappa}}\omega_{\vec{k}}}e^{i(\vec{k}-\vec{\kappa})\cdot\vec{x}}\int\frac{d^{3}q}{(2\pi)^{3}}\Big{[}\frac{(\vec{\epsilon}_{\vec{k}\lambda}\cdot\vec{q})(\vec{\epsilon}_{\vec{\kappa}\sigma}\cdot\vec{q})}{q^{2}}$
$\displaystyle\times\frac{\int r^{3}drj_{1}(qr)R_{20}(r)\times\int
r^{3}drj_{1}(qr)R_{10}(r)}{-\epsilon_{B}^{J/\psi}-\frac{q^{2}}{m_{Q}}-\omega_{\vec{\kappa}}+i\epsilon}\Big{]}.$
(53)

Combining these two transition amplitudes, the cross section for $g+J/\psi\to
g+\psi(2S)$ is obtained via the standard averaging and summation procedure

$\displaystyle\sigma^{(g+h)}(E_{g})=2\pi V\frac{V}{(2\pi)^{3}}\int
d^{3}\kappa\sum_{\sigma}\sum_{c}\frac{1}{4\pi}\int
d\Omega_{\vec{k}}\frac{1}{2}\sum_{\lambda}\frac{1}{8}\sum_{a}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times|T_{fi}^{(g)}+T_{fi}^{(h)}|^{2}\delta(E_{i}-E_{f})$
$\displaystyle=\frac{{g_{s}}^{4}}{2^{5}3^{4}{\pi}^{5}}\omega_{\vec{k}}\int_{0}^{\infty}d\kappa{\kappa}^{2}{\omega}_{\vec{\kappa}}Tr[\mathcal{M}^{2}]\delta(-{\epsilon}_{B}^{J/\psi}+{\omega}_{\vec{k}}+{\epsilon}_{B}^{\psi(2S)}-{\omega}_{\vec{\kappa}}),$
(54)

with the matrix $\mathcal{M}$ being

$\displaystyle\mathcal{M}_{ij}=\int{d^{3}q}\Big{[}\frac{q_{i}q_{j}}{q^{2}}{\int
r^{3}drj_{1}(qr)R_{10}(r)\int r^{3}drj_{1}(qr)R_{20}(r)}$
$\displaystyle\times(\frac{1}{-\epsilon_{B}^{J/\psi}+\omega_{\vec{k}}-\frac{q^{2}}{m_{Q}}+i\epsilon}+\frac{1}{-\epsilon_{B}^{J/\psi}-\frac{q^{2}}{m_{Q}}-\omega_{\vec{\kappa}}+i\epsilon})\Big{]},$
(55)

where $i,j$=1, 2, 3 label the Cartesian momentum components. It is not hard to
prove that all off-diagonal elements vanish and one is left with only the
diagonal elements

$\displaystyle\mathcal{M}_{11}=\mathcal{M}_{22}=\mathcal{M}_{33}$
$\displaystyle=\frac{4\pi}{3}\int_{0}^{\infty}dq\Big{[}q^{2}{\int
r^{3}drj_{1}(qr)R_{10}(r)\int r^{3}drj_{1}(qr)R_{20}(r)}$
$\displaystyle\times(\frac{1}{-\epsilon_{B}^{J/\psi}+\omega_{\vec{k}}-\frac{q^{2}}{m_{Q}}+i\epsilon}+\frac{1}{-\epsilon_{B}^{J/\psi}-\frac{q^{2}}{m_{Q}}-\omega_{\vec{\kappa}}+i\epsilon})\Big{]},$
(56)

where, $q=|\vec{q}|$ is the magnitude of the internal relative momentum of
intermediate octet $(c\bar{c})_{8}$. We use the analytical Coulomb radial wave
function to mimic the numerically solved radial function with full potential
(Eq. (48)) by tuning the Bohr radius parameter. This way, the two integrations
over radial wave functions in Eq. (V.1) can be analytically performed and
what’s left is two polynomial functions of $q$

$\displaystyle\frac{2^{4}a^{7/2}q}{(1+a^{2}q^{2})^{3}},~{}~{}~{}~{}~{}~{}\frac{2^{19/2}a^{7/2}q(2a^{2}q^{2}-1)}{(1+4a^{2}q^{2})^{4}},$
(57)

with $a$ being the Bohr radius parameter. These two polynomials contain higher
order poles with respect to $q$, which, together with the pole
$q=\sqrt{m_{Q}(\omega_{\vec{k}}-\epsilon_{B}^{J/\psi})}$ in the first term of
the third line of Eq. (V.1) when $\omega_{\vec{k}}>\epsilon_{B}^{J/\psi}$, can
be handled by extending the integration to the upper half plane of complex $q$
(noting that the integrand is an even function with respect to $q$) and then
making use of the residual theorem.

We’ve also derived the transition cross sections between two $p$-wave bound
states. These cross sections, including Eq. (V.1), suffer from no divergence.

###  V.2 Numerical results

Figure 10: Cross sections for the second-order transitions between different
charmonium bound states at finite temperatures. (a) $g+J/\psi\to g+J/\psi)$
(elastic collisions) and $g+J/\psi\to g+\psi(2S)$ (excitation). (b)
$g+\chi_{c}\to g+\chi_{c}$ (elastic collision).

Figure 11: Cross sections for the second-order transitions between different
bottomonium bound states at finite temperatures. (a) $g+\Upsilon(1S)\to
g+\Upsilon(1S)$ (elastic collisions) and $g+\Upsilon(1S)\to g+\Upsilon(2S)$
(excitation). (b) $g+\chi_{b}(1P)\to g+\chi_{b}(1P)$ (elastic collision) and
$g+\chi_{b}(1P)\to g+\chi_{b}(2P)$ (excitation).

The computed second-order transition (elastic scattering or excitation) cross
sections are presented in Fig. 10 and Fig. 11 for charmonia and bottomonia,
respectively. These cross sections share the same feature that they rise fast
with incident gluon energy and peak near the threshold, which is followed by a
decline, and gradually level off at large energies. In terms of magnitude,
these elastic or excitation cross sections are much smaller than the total NLO
dissociation cross sections, but are comparable to the component
$\sigma^{(a+b)}$ (cf. Fig. 5(a) and Fig. 6(a)) arising from diagrams (a) and
(b) in Fig. 2. The latter is not surprising since underlying all these
diagrams is the same exchange of a color octet state. The cross sections for
transition between $p$-wave bound states are generally larger than those
involving $s$-wave bound states.

The second-order transition rates are computed using Eq. (49). The numerical
results for these rates are displayed in Fig. 12 and Fig. 13 for charmonia and
bottomonia, respectively. Contrary to the NLO dissociation rates growing with
temperature (cf. Fig. 7 and Fig. 8), these rates decrease with temperature
(except for the case of $\Upsilon(1S)$ elastic scattering which exhibits a
moderate increase toward high temperatures). At rather low temperatures near
$T_{c}$, these rates turn out to be comparable to the corresponding NLO
dissociation rates, implying that it might be necessary to take these rates
into account when interpreting the thermal spectral width of heavy quarkonium
at low temperatures as computed in lattice QCD Larsen:2019zqv .

Figure 12: Second-order transition rates between different charmonium bound
states at finite temperatures. (a) $g+J/\psi\to g+J/\psi)$ (elastic
collisions) and $g+J/\psi\to g+\psi(2S)$ (excitation). (b) $g+\chi_{c}\to
g+\chi_{c}$ (elastic collision).

Figure 13: Second-order transition rates between different bottomonium bound
states at finite temperatures. (a) $g+\Upsilon(1S)\to g+\Upsilon(1S)$ (elastic
collisions) and $g+\Upsilon(1S)\to g+\Upsilon(2S)$ (excitation). (b)
$g+\chi_{b}(1P)\to g+\chi_{b}(1P)$ (elastic collision) and $g+\chi_{b}(1P)\to
g+\chi_{b}(2P)$ (excitation).

##  VI Summary

In this work, we have presented a systematic study of the parton-induced
second-order dissociation and transition processes for heavy quarkonia in the
QGP within the approach of quantum mechanical perturbation theory, utilizing
the effective chromo-electric dipole coupling of the heavy quarkonium with
external gluons and a Hamiltonian formulation of QCD vertices. In this
approach, the NLO dissociation and transition cross sections have been derived
directly from the dynamical scattering point of view, which systematically
incorporate the bound state wave functions and binding energies. This is in
contrast to the indirect extraction of the cross sections from the imaginary
heavy quark potential as done in Brambilla:2013dpa  and also goes beyond the
quasifree treatment Grandchamp:2001pf ; Rapp:2008tf . Our approach also
differs from the non-relativistic Bethe-Salpeter amplitude method Hong:2018vgp
which considers only barely bound quarkonia with near-threshold energy,
whereas our approach highlights the binding effects. Employing the bound state
wave functions and binding energies calculated from an in-medium potential
model, we have numerically evaluated the NLO dissociation and transition cross
sections and rates for various charmonia and bottomonia at finite
temperatures. Medium effects come into play through the temperature dependence
of bound state wave functions (dipole size), binding energies and the thermal
masses of light partons (which also cure the otherwise infrared and collinear
divergences in evaluating cross sections). The NLO dissociation cross sections
exhibit a rise and then a peculiar saturation as incident parton energy
increases, as a result of the outgoing parton carrying away the excess energy,
making them much more efficient in breaking up the bound states than the LO
processes for energetic partons. Consequently, the NLO dissociation rates
quickly take over from the LO counterparts at increasing temperatures. Finally
we have also revealed that the second-order elastic scattering or excitation
of bound state may induce significant widths for in-medium heavy quarkonium at
temperatures not far from $T_{c}$.

The results presented here may be useful for phenomenological studies of heavy
quarkonia transport in the relativistic heavy-ion collision systems. In
particular, the parton-induced heavy quarkonium scattering amplitudes derived
here might be of important use for more differential (beyond the simple rate
equation approach Rapp:2008tf ) and real-time simulations of heavy quarkonium
dissociation and regeneration in the QGP in combination with the single heavy
quark transport Dong:2019byy ; He:2021zej ; He:2022ywp .

###### Acknowledgements.

M. H. thanks Profs. Y.-G. Ma, G.-L. Ma and J.-H. Chen for the hospitality
during his stay as a visiting scholar at Shanghai Research Center for
Theoretical Nuclear Physics of Fudan University where the work was finalized.
This work was supported by the National Natural Science Foundation of China
(NSFC) under Grants No.12075122 (M. H.) and No.12147101 (via Shanghai Research
Center for Theoretical Nuclear Physics).

## References

  * (1) N. Brambilla, S. Eidelman, B. K. Heltsley, R. Vogt, G. T. Bodwin, E. Eichten, A. D. Frawley, A. B. Meyer, R. E. Mitchell and V. Papadimitriou, et al. Eur. Phys. J. C 71, 1534 (2011).
  * (2) R. Rapp, D. Blaschke and P. Crochet, Prog. Part. Nucl. Phys. 65, 209-266 (2010).
  * (3) P. Braun-Munzinger and J. Stachel, Landolt-Bornstein 23, 6-53 (2010).
  * (4) L. Kluberg and H. Satz, Landolt-Bornstein 23, 6-1 (2010).
  * (5) A. Mocsy, P. Petreczky and M. Strickland, Int. J. Mod. Phys. A 28, 1340012 (2013).
  * (6) A. Andronic, F. Arleo, R. Arnaldi, A. Beraudo, E. Bruna, D. Caffarri, Z. C. del Valle, J. G. Contreras, T. Dahms and A. Dainese, et al. Eur. Phys. J. C 76, no.3, 107 (2016).
  * (7) A. Rothkopf, Phys. Rept. 858, 1-117 (2020).
  * (8) J. Zhao, K. Zhou, S. Chen and P. Zhuang, Prog. Part. Nucl. Phys. 114, 103801 (2020).
  * (9) A. Andronic, P. B. Gossiaux, P. Petreczky, R. Rapp, M. Strickland, J. P. Blaizot, N. Brambilla, P. Braun-Munzinger, B. Chen and S. Delorme, et al. arXiv:2402.04366 [nucl-th].
  * (10) T. Matsui and H. Satz, Phys. Lett. B 178, 416-422 (1986).
  * (11) S. Digal, P. Petreczky and H. Satz, Phys. Rev. D 64, 094015 (2001).
  * (12) D. Kharzeev and H. Satz, Phys. Lett. B 334, 155-162 (1994).
  * (13) L. Grandchamp and R. Rapp, Phys. Lett. B 523, 60-66 (2001).
  * (14) P. Braun-Munzinger and J. Stachel, Phys. Lett. B 490, 196-202 (2000).
  * (15) R. L. Thews, M. Schroedter and J. Rafelski, Phys. Rev. C 63, 054905 (2001).
  * (16) L. Grandchamp, R. Rapp and G. E. Brown, Phys. Rev. Lett. 92, 212301 (2004).
  * (17) X. Zhao and R. Rapp, Phys. Rev. C 82, 064905 (2010).
  * (18) T. Song, K. C. Han and C. M. Ko, Phys. Rev. C 85, 014902 (2012).
  * (19) M. Strickland, Phys. Rev. Lett. 107, 132301 (2011).
  * (20) K. Zhou, N. Xu, Z. Xu and P. Zhuang, Phys. Rev. C 89, no.5, 054911 (2014).
  * (21) N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, Phys. Rev. D 96, no.3, 034021 (2017).
  * (22) N. Brambilla, M. Á. Escobedo, M. Strickland, A. Vairo, P. Vander Griend and J. H. Weber, JHEP 05, 136 (2021).
  * (23) X. Yao, Int. J. Mod. Phys. A 36, no.20, 2130010 (2021).
  * (24) B. Wu and R. Rapp, [arXiv:2404.09881 [nucl-th]].
  * (25) M. E. Peskin, Nucl. Phys. B 156, 365-390 (1979).
  * (26) G. Bhanot and M. E. Peskin, Nucl. Phys. B 156, 391-416 (1979).
  * (27) Y. s. Oh, S. Kim and S. H. Lee, Phys. Rev. C 65, 067901 (2002).
  * (28) S. Chen and M. He, Phys. Rev. C 96, no.3, 034901 (2017).
  * (29) C. Y. Wong, Phys. Rev. C 72, 034906 (2005).
  * (30) F. Brezinski and G. Wolschin, Phys. Lett. B 707, 534-538 (2012).
  * (31) Y. Liu, C. M. Ko and T. Song, Phys. Rev. C 88, no.6, 064902 (2013).
  * (32) F. Riek and R. Rapp, New J. Phys. 13, 045007 (2011).
  * (33) X. Du, M. He and R. Rapp, Phys. Rev. C 96, no.5, 054901 (2017).
  * (34) T. Song and S. H. Lee, Phys. Rev. D 72, 034002 (2005).
  * (35) J. Hong and S. H. Lee, Phys. Rev. C 99, no.3, 034905 (2019).
  * (36) M. Laine, O. Philipsen, P. Romatschke and M. Tassler, JHEP 03, 054 (2007).
  * (37) A. Beraudo, J. P. Blaizot and C. Ratti, Nucl. Phys. A 806, 312-338 (2008).
  * (38) N. Brambilla, J. Ghiglieri, A. Vairo and P. Petreczky, Phys. Rev. D 78, 014017 (2008).
  * (39) N. Brambilla, M. A. Escobedo, J. Ghiglieri and A. Vairo, JHEP 05, 130 (2013).
  * (40) Y. Burnier, O. Kaczmarek and A. Rothkopf, Phys. Rev. Lett. 114, no.8, 082001 (2015).
  * (41) D. Bala et al. [HotQCD], Phys. Rev. D 105, no.5, 054513 (2022).
  * (42) A. Bazavov et al. [HotQCD], Phys. Rev. D 109, no.7, 074504 (2024).
  * (43) T. M. Yan, Phys. Rev. D 22, 1652 (1980).
  * (44) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Rev. Mod. Phys. 77, 1423 (2005).
  * (45) M. B. Voloshin, Prog. Part. Nucl. Phys. 61, 455-511 (2008).
  * (46) Y. Sumino, arXiv:1411.7853 [hep-ph].
  * (47) S. Chen and M. He, Phys. Lett. B 786, 260-267 (2018).
  * (48) R. Larsen, S. Meinel, S. Mukherjee and P. Petreczky, Phys. Lett. B 800, 135119 (2020).
  * (49) X. Yao and B. Müller, Phys. Rev. D 100, no.1, 014008 (2019).
  * (50) M. E. Peskin and D. V. Schroeder, An Introduction to Quantum Field Theory, Addison-Wesley Publishing Company, 1995.
  * (51) J. J. Sakurai, Advanced Quantum Mechanics, World Scientific, 2008.
  * (52) M. D. Schwartz, Quantum Field Theory and the Standard Model, Cambridge University Press, 2014.
  * (53) W. Greiner, Quantum Mechanics: Special Chapters, Springer, 1998.
  * (54) K. Gottfried and T.-M. Yan, Quantum Mechanics: Fundamentals, Springer, 2003.
  * (55) F. Karsch, M. T. Mehr and H. Satz, Z. Phys. C 37, 617 (1988).
  * (56) Y. Burnier, O. Kaczmarek and A. Rothkopf, JHEP 12, 101 (2015).
  * (57) S. Cao, G. Coci, S. K. Das, W. Ke, S. Y. F. Liu, S. Plumari, T. Song, Y. Xu, J. Aichelin and S. Bass, et al. Phys. Rev. C 99, no.5, 054907 (2019).
  * (58) P. B. Gossiaux and J. Aichelin, Phys. Rev. C 78, 014904 (2008).
  * (59) J. Zhao, J. Aichelin, P. B. Gossiaux and K. Werner, Phys. Rev. D 109, no.5, 054011 (2024).
  * (60) E. Braaten and M. H. Thoma, Phys. Rev. D 44, 1298-1310 (1991).
  * (61) X. Dong, Y. J. Lee and R. Rapp, Ann. Rev. Nucl. Part. Sci. 69, 417-445 (2019).
  * (62) M. He, B. Wu and R. Rapp, Phys. Rev. Lett. 128, no.16, 162301 (2022).
  * (63) M. He, H. van Hees and R. Rapp, Prog. Part. Nucl. Phys. 130, 104020 (2023).