11institutetext: Nora Brambilla 22institutetext: Physik-Department, Technische
Universität München, James-Franck-Str. 1, 85748 Garching, Germany; Institute
for Advanced Study, Technische Universität München, Lichtenbergstrasse 2 a,
85748 Garching, Germany; Munich Data Science Institute, Technische Universität
München, Walther-von-Dyck-Strasse 10, 85748 Garching, Germany 22email:
<EMAIL_ADDRESS>
# Quark Nuclear Physics with Heavy Quarks

Nora Brambilla corresponding author

###### Abstract

Heavy quarks have been instrumental for progress in our exploration of strong
interactions. Quarkonium in particular, a heavy quark-antiquark
nonrelativistic bound state, has been at the root of several revolutions.
Quarkonium is endowed with a pattern of separated energy scales qualifying it
as special probe of complex environments. Its multiscale nature has made a
description in Quantum Field Theory particularly difficult up to the advent of
nonrelativistic effective field theories. We will focus on systems made by two
or more heavy quarks. After considering some historical approaches based on
the potential models and the Wilson loop approach, we will introduce the
contemporary nonrelativistic effective field theory descriptions, in
particular potential nonrelativistic QCD which entails the Schödinger equation
as zero order problem, define the potentials as matching coefficients and
allows systematic calculations of the physical properties. The effective field
theory allows us to explore quarkonium properties in the realm of QCD. In
particular it allows us to make calculations with unprecedented precision when
high order perturbative calculations are possible and to systematically
factorize short from long range contributions where observables are sensitive
to the nonperturbative dynamics of QCD. Such effective field theory treatment
can be extended at finite temperature and in presence of gluonic and light
quark excitations. We will show that in this novel theoretical framework,
quarkonium can play a crucial role for a number of problems at the frontier of
our research, from the investigation of the confinement dynamics in strong
interactions to the study of deconfinement and the phase diagram of nuclear
matter, to the precise determination of Standard Model parameters up to the
emergence of exotics X Y Z states of an unprecedented nature.

##  1 The role of heavy quarks

Heavy quarks and bound states of heavy quarks, primarily quarkonium, a bound
state of a heavy quark and a heavy antiquark, have been historically
instrumental to construct the theory of strong interactions and continue today
to be at the forefront of our research as a golden probe of the strong
dynamics.

The discovery of heavy quarks drastically changed the Standard Model (SM) of
particle physics. This happened in the November revolution of 1974 when two
labs on opposite sides of USA announced discovery of a new particle, a fact
that helped the acceptance of the Standard Model of particle physics
E598:1974sol ; SLAC-SP-017:1974ind . The new particle was the first example of
quarkonium: the $J/\psi$, the lowest excitation made by a charm and an
anticharm. The $J/\psi$ appeared as an unprecedented sharp peak, tall and
narrow, 3 GeV in mass and 90 KeV in width, at variance with the typical width
of several tens and hundreds of MeV of the hadrons discovered up to that time,
i.e. strongly interacting light quark composite particles. The $J/\psi$
discovery represented the confirmation of the quark model, the discovery of
the charm quark, the confirmation of the GIM mechanism Glashow:1970gm  (the
mechanism through which flavour-changing neutral currents are suppressed in
loop diagrams) and the first discovery of a quark of large mass moving
nonrelativistically. It triggered additional searches and in few years the
higher excitations of charmonium were discovered as well has bottomonia
(1977), i.e. bound states of bottom and antibottom, $B_{c}$ (1998), and the
top (1995). States made by top and antitop decay weakly before forming a
proper bound state, however still leaving their signature in the form of an
enhancement of the cross section at threshold QuarkoniumWorkingGroup:2004kpm .

The $J/\psi$ discovery was the confirmation of QuantumChromoDynamics (QCD)
Fritzsch:1973pi , the Quantum Field Theory describing strong interactions. QCD
has a well defined behaviour in the ultraviolet (UV) region at large energy
and a fundamental coupling constant $\alpha_{s}=g^{2}/4\pi$ running from small
values at large energy to large values at small energy. This encodes the
properties of asymptotic freedom (quarks are free at high momentum transfer)
and confinement (quarks are confined in color singlet hadrons at low energy)
Wilson:1974sk . As we will see, confinement becomes manifest in the case of
heavy quarks, where one can write the color singlet quark-antiquark
interaction potential in terms of a so called Wilson loop Brambilla:1999ja .
Confinement emerges in an area law of the Wilson loop and correspondingly in a
linear potential growing with the distance between the quarks Creutz:1976ch .
An emergent scale $\Lambda_{QCD}$ parametrizes the importance of the
nonperturbative corrections, i.e. the contributions that cannot be calculated
in an expansion in $\alpha_{s}$ and is mirrored in the hadron spectrum, the
mass of the proton being proportional to $\Lambda_{QCD}$. We consider a quark
to be heavy when its mass $m$ is larger than the scale $\Lambda_{QCD}$: this
qualifies as heavy the charm, bottom and top quarks.

These fundamental features of QCD find the best realization in the $J/\psi$
and in quarkonium in general. The small width can be explained by the fact
that $J/\psi$ is the lowest energy level and can decay only via annihilation,
which makes available in the process a large energy, of order of two times the
mass of the charm (about 2 GeV). The annihilation width is then proportional
to $\alpha_{s}^{2}(2m_{c})$ which is small due to asymptotic freedom, since
$m_{c}$ is bigger than $\Lambda_{QCD}$. On the other hand, when theorists set
up to investigate the structure of the energy levels of charmonium and
bottomonium, they noticed that it can be reproduced by using in the
Schrödinger equation a static potential superposition of an attractive Coulomb
contribution (with the appropriate $SU(3)$ color factor for a singlet
$Q\bar{Q}$) and a term linear in the distance: the famous Cornell potential
Eichten:1978tg ; Eichten:1979ms . Such form of the potential has been later
confirmed by nonperturbative calculations performed using computational
lattice QCD Creutz:1976ch ; Brambilla:1999ja ; Rothe:1992nt ; Bali:2000gf ;
Campostrini:1987ht . In 1977 the static potential was calculated at two loops
Fischler:1977yf  which enabled to reobtain the QCD $\beta$ function at order
$O(g^{5})$ and the solution of the Callan-Symanzik equation to order $g^{5}$,
which gave a hint on the existence of a color confining potential.

In the following we will address the importance of heavy quarks for quark
nuclear physics. We will focus on systems made by two or more heavy quarks and
discuss how they play a crucial role for a number of problems at the frontier
of our research, from the investigation of the confinement dynamics in QCD to
the study of deconfinement and the phase diagram of nuclear matter, to the
precise determination of Standard Model parameters up to the emergence of
exotics X Y Z states of an unprecedented nature Brambilla:2010cs ;
QuarkoniumWorkingGroup:2004kpm ; Andronic:2015wma ; Chapon:2020heu ;
Brambilla:2014jmp . In particular, we will conclude that our progress in these
strong interactions topics is connected to a broad sweep of physical problems
in settings ranging from astrophysics and cosmology to strongly correlated
systems in particle and condensed matter physics as well as to search of
physics beyond the Standard model.

Heavy quarks had and have a key role to address the phenomenon of CP (Charge
Conjugation and Parity) violation, which is crucial to understand the
asymmetry between matter and antimatter that exists in the Universe, by
testing and constraining the Cabibbo-Kobayashi-Maskawa matrix entries entering
$B$ and $D$ mesons decays and mixing and by identifying new physics
contributions, see e.g. the reviews Buchalla:2008jp ; Neubert:1996qg ;
Chang:2017wpl ; Gershon:2016fda . The Babar and Belle experiments at the B
factories Belle-II:2018jsg ; BaBar:2014omp  have been constructed to this aim
but turned out to be also formidable heavy mesons machines giving a great
boost to out knowledge of heavy quark systems and their strong interaction
dynamics.

The review is organized as follows. First, we will discuss the physics
characteristics of systems made by heavy quarks. Then, we will summarize what
have been historically the pioneering approaches and the phenomenological
models. After that, we will explain that to address in quantum field theory a
nonrelativistic multiscale system like quarkonium, it is necessary resort to
nonrelativistic effective field theories (NREFTs). We will show how to
construct NREFTs up to arriving at the simplest possible version, called
potential nonrelativistic QCD (pNRQCD) which implements the Schrödinger
equation as the zero order problem. Then, we will show how combining NREFTs
and lattice we can get systematic and under control predictions on a number of
physical processes and observables like the spectrum, decays, transitions and
production. In this framework quarkonium becomes a unique laboratory for the
study of strong interactions from the high energy to the low energy scales.
The NREFT allows us to explore quarkonium properties in the realm of QCD. In
particular it allows us to make calculations with unprecedented precision when
high order perturbative calculations are possible and to systematically
factorize short from range contributions where observables are sensitive to
the nonperturbative dynamics of QCD. Finally, we will address some of the most
interesting open problems in relation to quarkonium, i.e. the non-equilibrium
propagation of quarkonium in medium which calls for introducing open quantum
system on top of NREFTs and the new exotic stares X Y Z discovered at the
accelerator experiments.

##  2 Heavy-light mesons, quarkonia , baryons with two or more heavy quarks

Heavy quarkonia are systems composed by a heavy quark and a heavy antiquark
with mass $m$ larger than the “QCD confinement scale” $\Lambda_{\rm QCD}$, so
that $\alpha_{\rm s}(m)\ll 1$ holds. We have that $m_{c}\sim 1.5$ GeV and
$m_{b}\sim 5$ GeV. Of course the quark masses are scheme dependent objects and
their actual value depends on the considered scheme and scale, only the pole
mass being scale independent, see e.g. the review Brambilla:2004jw . From the
quarkonia spectra, see Fig. 1 (and 2 and 3), it is evident that the difference
in the orbital energy levels is much smaller than the quark mass. It scales
like $mv^{2}$, while fine and hyperfine separations scale like $mv^{4}$. Here
$v$ is the heavy quark velocity ( $v=|\vec{v}|$) in the rest frame of the
meson in units of $c$, and $v^{2}\sim 0.1$ for the $b\bar{b}$, $v^{2}\sim 0.3$
for $c\bar{c}$ systems. This is the same scaling of the hydrogen atom if one
identifies $v$ with the fine structure constant $\alpha_{\rm em}$. Therefore
quarkonia are nonrelativistic systems. Being nonrelativistic, quarkonia are
characterized by a hierarchy of energy scales: the mass $m$ of the heavy quark
(hard scale), the typical relative momentum $p\sim mv$ (in the meson rest
frame) corresponding to the inverse Bohr radius $r\sim 1/(mv)$ (soft scale),
and the typical binding energy $E\sim mv^{2}$ (ultrasoft scale). Of course,
for quarkonium there is another scale that can never be switched off in QCD,
i.e. $\Lambda_{\rm QCD}$, the scale at which non perturbative effects become
dominant. A similar pattern of scales emerge in the case of baryons composed
of two or three heavy quarks Brambilla:2005yk  and for the just discovered
state $X(6900)$ LHCb:2020bwg  made by two charm and two anticharm quarks. The
pattern of nonrelativistic scales makes all the difference between heavy
quarkonia and heavy-light mesons, which are characterized by just two scales:
$m$ and $\Lambda_{\rm QCD}$.

Being a multiscale system, heavy quarkonium is probing different energy
regimes of the strong interactions, from the hard region, where an expansion
in the coupling constant is possible and precision studies may be done, to the
low-energy region, dominated by confinement and the many manifestations of
nonperturbative dynamics. In addition, the properties of production and
absorption of quarkonium in a nuclear and hot medium are crucial inputs for
the study of QCD at high density and temperature, reaching out to cosmology.
On the experimental side the diversity, quantity and accuracy of the data
collected in the last decades is impressive and includes clean and precise
samples of quarkonia spectra decay, transition and production processes,
including the discovery of exotics X Y Z states at tau-charm (BES experiment)
and B factories (Babar and Belle experiments), hadroproduction at Fermilab
Tevatron and the Large Hadron Collider (LHC) experiments at CERN, production
in photon-gluon fusion at DESY, photoproduction at Jlab, heavy ions production
and suppression at RHIC, NA60, and LHC. New data are coming also from the
upgraded experiments BELLE II and BESIII Brambilla:2010cs ;
QuarkoniumWorkingGroup:2004kpm ; Andronic:2015wma ; Chapon:2020heu ;
Yuan:2019zfo  and more will come in future from Panda at FAIR and the Electron
Ion Collider (EIC) Yuan:2019zfo ; PANDA:2021ozp . On the theoretical side the
last few decades have seen the construction of new nonrelativistic effective
field theories and new developmentd in computational lattice QCD which supply
us with a systematic calculational framework in quantum field theory. All this
make quarkonium a golden probe of strong interactions.

From now on we will concentrate on the study of systems with two or more heavy
quarks. We will outline the rich interplay of theoretical advancement and
experimental success and its implication on our control of strong interactions
inside the Standard Model of Particle Physics. Before, however, I will
summarize the models that have been used in the past to describe the
quarkonium properties.

Figure 1: The experimental quarkonium energy levels ($b\bar{b}$ and
$c\bar{c}$) as known in the eighties (plotted as relative to the spin-average
of the $\chi_{b}(1P)$ and $\chi_{c}(1P)$ states to be able to present the two
sectors in the same plot). The states are identified by names and mass in
parenthesis, $\Upsilon$ being the vector states (with $L=0$) in bottomonium,
$J/\psi$ and $\psi$ the vector states in charmonium sector, $\eta$ the
pseudoscalar states and $\chi$ and $h$ states with $L=1$ and different total
angular momentum. The first strong decay thresholds are also shown. Figure 2:
The spectrum of states in the $c\bar{c}$ sector as of 2019, taken from ref.
Brambilla:2019esw . Thin solid lines represent the states established
experimentally and dashed lines are for those that are claimed but not (yet)
established (a state is regarded as established if it is seen in different
modes). States whose quantum numbers are undetermined are not shown. States in
the plot are labeled according to the PDG primary naming scheme that
superseded the X Y Z notation, for further details see Brambilla:2019esw .
Dashed lines show some relevant thresholds that open in the considered mass
range; here $D_{1}$ stands for $D_{1}(2420)$ and $D_{2}^{*}$ for
$D_{2}^{*}(2460)$. Thresholds with hidden strangeness or involving broad
states are not shown. The states shown in the two columns to the right are
isovectors containing a $\bar{c}c$ pair; they are necessarily exotic. The just
discovered state $X(6900)$ LHCb:2020bwg  made by two charm and two anticharm
quarks would be out of the scale of the figure. More states have been observed
after 2019.  Figure 3: The spectrum of states in the $\bar{c}c$ sector as of
2019, taken from ref. Brambilla:2019esw . Thin solid lines represent the
states established experimentally and dashed lines are for those that are
claimed but not (yet) established (a state is regarded as established if it is
seen in different modes). States whose quantum numbers are undetermined are
not shown. States in the plot are labeled according to the PDG primary naming
scheme that superseded the X Y Z notation, for further details see
Brambilla:2019esw . Dashed lines show some relevant thresholds that open in
the considered mass range; here $D_{1}$ stands for $D_{1}(2420)$ and
$D_{2}^{*}$ for $D_{2}^{*}(2460)$. Thresholds with hidden strangeness or
involving broad states are not shown. The states shown in the two columns to
the right are isovectors containing a $\bar{c}c$ pair; they are necessarily
exotic.

###  2.1  The potential and the phenomenology of quarkonium

You find in Figs. (2) and (3) our present knowledge about the states in the
charmonium and in the bottomonium sector. We will consider in the following
quarks of equal mass for simplicity (in the case of $B_{c}$ one should take
into account that the masses are different). The quark and the antiquark spins
combine to give the total spin ${\bf S}={\bf S}_{1}+{\bf S}_{2}$ which
combines with the orbital angular momentum ${\bf L}$ to give the total angular
momentum ${\bf J}$. The resulting state solution of a Schrödinger equation is
denoted by $n^{2S+1}L_{J}$ where $n-1$ is the number of radial nodes of the
wave function. As usual, to $L=0$ is given the name $S$, to $L=1$ the name
$P$, to $L=2$ the name D and so on. The experimental resonances (see Figs. (2)
and (3)) are classified via the $J^{PC}$ quantum numbers, $P=(-1)^{L+1}$ being
the parity number and $C=(-1)^{L+S}$ the C-parity. Strong decay thresholds are
marked by horizontal dashed lines that represent the energy necessary to decay
in a couple of heavy-light mesons. In the charmonium sector only 10 states of
the states presented in Fig. (2) have been discovered before 1980 and no one
between 1980 and 2002. In 2003 the new revolution started with the discovery
of the $X(3872)$ Belle:2003nnu  and a number of new states above and below the
strong decay threshold Brambilla:2010cs . Many of the states discovered at or
above threshold presented exotic features and have been initially termed X Y Z
states Brambilla:2019esw .

The first tool used to describe quarkonium in the eighties has been the quark
model with some notions of QCD incorporated. In particular, constituent heavy
quark masses have been considered and a static potential called Cornell
potential Eichten:1978tg ; Eichten:1979ms  has been used in a Schrödinger
equation, successfully reproducing the quarkonia levels measured at that time.
The Cornell potential is flavor independent and has the form

$V_{0}(r)=-{\kappa\over r}+\sigma r+{\rm const},$ (1)

$r$ being the modulus of the quark-antiquark distance. The parameters $\kappa$
should be identified with ${4\over 3}\alpha_{\rm s}$, corresponding to the one
gluon exchange that should dominate at small distances due to asymptotic
freedom, and the string tension $\sigma$ corresponds to a constant energy
density related to confinement and originating a potential growing with the
interquark distance at large distances. A fit to the states gave $\kappa=0.52$
and $\sigma=0.182$ GeV2. Since then, several different phenomenological forms
of the static potential have been exploited, see Lucha:1991vn ;
Brambilla:1999ja  for a review. In order to explain the fine and hyperfine
structure of the quarkonium spectrum, however, spin dependent relativistic
corrections to the static potential have to be considered. Moreover, for
$v^{2}\sim 0.1$ for the $b\bar{b}$, $v^{2}\sim 0.3$ for $c\bar{c}$ systems,
one expects relativistic corrections of order $20\div 30\%$ for the charmonium
spectrum and up to $10\%$ for the bottomonium spectrum and also spin
independent but momentum dependent corrections have to be considered, arriving
at a phenomenological hamiltonian of the type

$H=\sum_{j=1,2}\left(m_{j}+{p_{j}^{2}\over 2m_{j}}-{p_{j}^{4}\over
8m_{j}^{3}}\right)+V_{0}+V_{\rm SD}+V_{\rm VD}.$ (2)

The $1/m^{2}$ spin-dependent $V_{\rm SD}$ and velocity-dependent $V_{\rm VD}$
potentials were initially derived in the eighties from the semirelativistic
reduction of a Bethe–Salpeter (BS) Salpeter:1951sz  equation for the quark-
antiquark connected amputated Green functions or, equivalently at this level,
from the semirelativistic reduction of the quark-antiquark scattering
amplitude with an effective exchange equal to the BS kernel. Several
ambiguities are involved in this procedure, due on one hand to the fact that
we do not know the relevant confining nonperturbative Bethe–Salpeter kernel,
on the other hand due to the fact that we have to get rid of the temporal (or
energy $Q_{0}$, $Q=p_{1}-p_{1}^{\prime}$ being the momentum transfer)
dependence of the kernel to recover a potential (instantaneous) description.
It turned out that, at the level of the approximation involved, the spin-
independent relativistic corrections at the order $1/m^{2}$ depend on the way
in which $Q_{0}$ is fixed together with the gauge choice of the kernel
Brambilla:1992fx ; Lucha:1991vn ; Brambilla:1999ja . The Lorentz structure of
the kernel was also not known. On a phenomenological basis, the following
ansatz for the kernel was intensively studied Brambilla:1999ja

$\displaystyle
I(Q^{2})=(2\pi)^{3}\left[\gamma_{1}^{\mu}\gamma_{2}^{\nu}P_{\mu\nu}J_{v}(Q)+J_{s}(Q)\right]$
(3)

in the instantaneous approximation $Q_{0}=0$.

Notice that the effective kernel above was taken with a pure dependence on the
momentum transfer $Q$. But, of course, the dependence on the quark and
antiquark momenta could have been more complicated. The vector kernel
$J_{v}(Q)$ above would correspond to the one gluon exchange (with $P_{\mu\nu}$
depending on the adopted gauge) while the scalar kernel $J_{s}(Q)$ would
account for the nonperturbative interaction. Taking
$J_{v}\\!=\\!\displaystyle-{1\over 2\pi^{2}}{4\over 3}{1\over{\bf Q}^{2}}$ and
$J_{s}\\!=\\!\displaystyle-{\sigma\over\pi^{2}}{1\over{\bf Q}^{4}}$, reproduce
the Cornell potential and the corrections in the nonrelativistic reduction of
such kernel in the instantaneous approximation would give a form for the
$V_{\rm SD}$ and the $V_{\rm VD}$. The confining part of the kernel was
usually chosen to be a Lorentz scalar in order to match the data on the fine
separation on the $P$ states Schnitzer:1978gq . If one however takes a kernel
which is not a pure convolution kernel, as it is to be expected in
interactions generated at higher orders, and deals more appropriately with the
instantaneous approximation, one can get quite different relativistic
corrections, especially for the velocity dependent part Brambilla:1992fx . We
conclude that this type of description is model dependent and does not allow
for further progress. In particular, the parameters of the model cannot be
related to the underlying field theory and the systematics of the model cannot
be estimated.

A more systematic procedure to relate the potential to QCD has been based on
the so called Wilson loop approach. In order to obtain information on the
structure of nonperturbative corrections to the potential it is useful to
define the potential in terms of gauge invariant objects suitable for a direct
lattice QCD evaluation. Indeed lattice QCD is one of the best methods to
extract nonperturbative information from QCD, for some reviews see
Kronfeld:2012uk ; Rothe:1992nt ; Kogut:1982ds ; Bali:2000gf ; Karsch:2001vs .
Let us see how this works in the case of the static potential. Let us consider
a locally gauge invariant quark-antiquark color singlet state (for more
details see Brambilla:1999ja ; Lucha:1991vn :

$|\phi_{\alpha\beta}^{lj}\rangle\equiv{\delta_{lj}\over\sqrt{3}}\bar{\psi}^{i}_{\alpha}(x)U^{ik}(x,y,C)\psi^{k}_{\beta}(y)|0\rangle$
(4)

where $i,j,k,l$ are colour indices (that will be suppressed in the following),
$|0\rangle$ denotes the QCD ground state and the Schwinger string line has the
form

$U(x,y;C)=P\exp\left\\{ig\int_{y}^{x}A_{\mu}(z)\,dz^{\mu}\right\\}\,,$ (5)

where $A_{\mu}=A_{\mu}^{a},a=1,8$ is the gluon vector potential of QCD in the
fundamental representation, $g$ the QCD coupling constant, and the integral is
extended along the path $C$. The operator $P$ denotes the path-ordering
prescription which is necessary due to the fact that $A_{\mu}$ are non-
commuting matrices. Let us see how the quark-antiquark potential can be
extracted from the quark-antiquark Green function constructed with such color
singlet states:

$G(T)=\langle\phi({\bf x},0)|\phi({\bf y},T)\rangle=\langle\phi({\bf
x},0)|\exp{(-iHT)}|\phi({\bf y},0)\rangle.$ (6)

Inserting a complete set of energy eigenstates $\psi_{n}$ with eigenvalues
$E_{n}$ and making a Wick rotation we find

$\displaystyle G(-iT)$ $\displaystyle=$ $\displaystyle\sum_{n}\langle\phi({\bf
x},0)|\psi_{n}\rangle\langle\psi_{n}|\phi({\bf y},0)\rangle\exp{(-E_{n}T)}$
(7) $\displaystyle\to$ $\displaystyle\langle\phi({\bf
x},0)|\psi_{0}\rangle\langle\psi_{0}|\phi({\bf
y},0)\rangle\exp{(-E_{0}T)}\quad{\rm for}\quad T\to\infty$

which gives the Feynman–Kac formula for the ground state energy

$E_{0}=-\lim_{T\to\infty}{\log G(-iT)\over T}.$ (8)

The only condition for the validity of Eq. (8) is that the $\phi$ states have
a non-vanishing component over the ground state. This is precisely the way in
which hadron masses are computed on the lattice. Of course, to maximize the
overlap with the ground state in consideration appropriate operators may be
used. If the $\phi$ state denotes a state of two exactly static particles
interacting at a distance $r$, then the ground state energy is a function of
the particle separation, $E_{0}\equiv E_{0}(r)$, and gives the potential of
the first adiabatic surface. It is possible to obtain an explicit analytic
form of the quark-antiquark Green function for infinitely heavy quarks
($m\to\infty$, static limit) and for large temporal intervals ($T\to\infty$),
see Brambilla:1999ja . We consider that at a time $t=0$ a quark and an
antiquark pair is created and that they interact while propagating for a time
$t=T$ at which they are annihilated. Then ($x_{j}=({\bf x}_{j},T),y_{j}=({\bf
y}_{j},0)$, see Fig. 4) we obtain

$\displaystyle G_{\beta_{1}\beta_{2}\alpha_{1}\alpha_{2}}(T)$
$\displaystyle{\buildrel m\to\infty\over{\longrightarrow}}$
$\displaystyle\,\delta^{3}({\bf x}_{1}-{\bf y}_{1})\delta^{3}({\bf x}_{2}-{\bf
y}_{2})(P_{+})_{\beta_{1}\alpha_{1}}(P_{-})_{\alpha_{2}\beta_{2}}$ (9)
$\displaystyle\times e^{-2mT}\langle{\rm Tr}\,{\rm
P}e^{ig\oint_{\Gamma_{0}}dz_{\mu}A_{\mu}(z)}\rangle$

with $P_{\pm}\equiv(1\pm\gamma_{4})/2$. The integral in Eq. (9) extends over
the circuit $\Gamma_{0}$ which is a closed rectangular path with spatial and
temporal extension $r=|{\bf x}_{1}-{\bf x_{2}}|$ and $T$ respectively, and has
been formed by the combination of the path-ordered exponentials along the
horizontal (=time fixed) lines, coming from the Schwinger strings, and those
along the vertical lines coming from the static propagators (see Fig. 4). The
brackets in (9) denote the QCD vacuum expectation value, which in Euclidean
space is

$\langle f[A]\rangle\equiv{1\over Z}\int{\cal D}Af[A]e^{-\int
d^{4}xL^{E}_{YM}}.$ (10)

From Eq. (9) it is clear that the dynamics of the quark-antiquark interaction
is contained in

$W(\Gamma_{0})={\rm Tr\,P}e^{\displaystyle
ig\oint_{\Gamma_{0}}dz_{\mu}A_{\mu}(z)}.$ (11)

This is the famous static Wilson loop Wilson:1974sk . In the limit of infinite
quark mass considered, the kinetic energies of the quarks drop out of the
theory, the quark Hamiltonian becomes identical with the potential while the
full Hamiltonian contains also all types of gluonic excitations. According to
the Feynman–Kac formula the limit $T\to\infty$ projects out the lowest state
i.e. the one with the “glue” in the ground state. This has the role of the
quark-antiquark potential for pure mesonic states. We will see in the section
over exotics that the excited states have the role of the potential in the
case of hybrid states, i.e. states with gluons contributing to the binding.
Now, comparing Eq. (9) with Eq. (7) and considering that the exponential
factor $\exp{(-2mT)}$ just accounts for the the fact that the energy of the
quark-antiquark system includes the rest mass of the pair, we obtain

$V_{0}(r)\equiv E_{0}(r)=-\lim_{T\to\infty}{1\over T}\log\langle
W(\Gamma_{0})\rangle.$ (12)

The heavy quark degrees of freedom have now completely disappeared and the
expectation value in (12) can be evaluated in the pure Yang–Mills theory or in
unquenched QCD (i.e. taking into account sea light quarks) Wilson:1974sk ;
Brown:1979ya . Notice that the potential is given purely in terms of a gauge
invariant quantity (the Wilson loop precisely). In this way we have reduced
the calculation of the static potential to a well posed problem in field
theory: to obtain the actual form of $V_{0}$ we need to calculate the QCD
expectation value of the static Wilson loop. The static Wilson loop is the low
energy domain of QCD is dominated by an area law, i.e. can be approximated at
leading order at large $r$ as $\langle W(\Gamma_{0})\rangle\sim exp\\{-\sigma
rT\\}$ (i.e. the area of the loop in the exponent multiplied by the string
tension), as shown in the strong coupling expansion Rothe:1992nt ;
Creutz:1976ch ; Bali:2000gf ; Brambilla:1999ja ; Creutz:1980zw . In the high
energy domain the first perturbative contribution comes from the one gluon
exchange and therefore Eq. (12) would reproduce in these limits the Cornell
potential of Eq. (1), which appears to be a simple superposition of these two
behaviours. We will come back to the area law behaviour of the Wilson loop in
the section on studies of confinement. The Wilson loop formalism have been
used to obtain the form of spin dependent $V_{SD}$ potential in the famous
papers of Eichten, Feinberg Eichten:1980mw  and Gromes Gromes:1984ma  and
later to obtain the form of spin independent, momentum dependent $V_{VD}$
potential in the classical works Barchielli:1986zs ; Barchielli:1988zp . These
works express the relativistic corrections to the potential in terms of
generalized static Wilson loops containing field strengths (i.e.
chromoelectric and chromomagnetic fields) insertion in the quark lines. In
Brambilla:1993zw  a nonstatic Wilson loop has been used to obtain the
relativistic corrections to the static three-quarks potential.

Up to the establishment of the NREFT called potential nonrelativistic QCD in
the nineties Brambilla:1999qa ; Brambilla:1999xf ; Brambilla:2004jw , the
Wilson loop approach has been the best theory founded method and was used to
calculate the potentials on the lattice (see e.g. the review Bali:2000gf ).

However, this method was providing only part of the QCD result. The Wilson
loop relativistic potential corrections were missing all the non analytic term
in the mass, $\ln{m/\mu}$ terms, that were instead obtained in a direct purely
perturbative QCD one loop calculation Gupta:1982kp . Moreover, QCD in the
static quark antiquark configuration can still change colour by emission of
gluons. This introduces a new dynamical scale in the evaluation of the
potential that should be taken into account as it became evident in the
seminal paper of Appelquist Dine and Muzinich Appelquist:1977tw  that
identified infrared divergences in the three loops perturbative fixed order
calculation of the static Wilson. Moreover, already in Lucha:1991vn  the
question was raised about why potential models where working so well. It took
a while and the development of non relativistic effective field theories to
answer this question and to be able to calculate systematically the potentials
from QCD.

All these problems are addressed and solved in pNRQCD.

Figure 4: The static Wilson loop. It contains the interaction of a static
quark-antiquark pair created at a time $t=0$ and annihilated at a subsequent
large time $T$. Initial and final states are made gauge invariant by the
present of Schwinger line given in eq. (5).

##  3 Nonrelativistic effective field theories

As we mentioned, the study of quarkonium in the last few decades has witnessed
two major developments: the establishment of nonrelativistic effective field
theories (NREFTs) and progress in lattice QCD calculations of excited states
and resonances, with calculations at physical light-quark masses. Both allow
for precise and systematically improvable computations that are largely model-
independent. It is precisely this advancement in the understanding of
quarkonium and quarkonium-like systems inside QCD that makes today quarkonium
exotics as particularly valuable (see the section on Exotics). In fact, today
that we are confronted with a huge amount of high-quality data, which have
provided for the first time uncontroversial evidence for the existence of
exotic hadrons, by using modern theoretical tools that allow us to explore in
a controlled way these new forms of matter we can get a unique insight into
the low-energy dynamics of QCD.

The appearance of a hierarchy of scales calls for the application of effective
field theory (EFT) methods. However, “heavy quark effective theory” (HQET)
Neubert:2005mu ; Georgi:1991mr , the EFT description of heavy-light mesons,
where only an ultraviolet mass scale $m$ and an infrared mass scale
$\Lambda_{\rm QCD}$ appear, is not suitable for the description of heavy
quarkonia, since HQET is unable to describe the dynamics of the binding. The
existence of several physical scales makes the theoretical description of
quarkonium physics more complicated. All scales get entangled in a typical
amplitude involving a quarkonium observable. For example, quarkonium
annihilation and production take place at the scale $m$, quarkonium binding
takes place at the scale $mv$, which is the typical momentum exchanged inside
the bound state, while very low-energy gluons and light quarks (also called
ultrasoft (US) degrees of freedom) live long enough that a bound state has
time to form and, therefore, are sensitive to the scale $mv^{2}$. Ultrasoft
gluons are responsible for phenomena similar to the the Lamb shift in hydrogen
atom. This pattern of scales has made the description in Quantum Field Theory
particularly difficult. The solution is to take advantage of the existence of
the different energy scales to substitute QCD with simpler but equivalent
NREFTs. A hierarchy of NREFTs, see Fig.5, may be constructed by systematically
integrating out modes associated with high-energy scales not relevant for the
quarkonium system. Such integration is made in a matching procedure that
enforces the equivalence between QCD and the EFT at a given order of the
expansion in $v$. The EFT Lagrangian is factorized in matching coefficients,
encoding the high energy degrees of freedom and low energy operators. The
relativistic invariance is realized via exact relations among the matching
coefficients Manohar:1997qy ; Brambilla:2003nt ; Heinonen:2012km ;
Berwein:2018fos  The EFT displays a power counting in the small parameter $v$,
therefore we are able to attach a definite power of $v$ to the contribution of
each EFT operators to the physical observables. We will introduce in the next
sections nonrelativistic QCD (NRQCD) Bodwin:1994jh ; Caswell:1985ui ;
Lepage:1992tx  and potential nonrelativistic QCD (pNRQCD) Pineda:1997bj ;
Brambilla:1999xf , see the review Brambilla:2004jw .

Figure 5: The hard, soft and ultrasoft scales of quarkonium and the
corresponding NREFTs that can be constructed. NRQCD follows from the
integration of the gluons and quarks at the hard scale, pNRQCD follows from
integration of gluons at the soft scale. If the soft scale is bigger than
$\Lambda_{\rm QCD}$ then the matching from NRQCD to pNRQCD is perturbative,
otherwise it is nonperturbative.

It turns to be instrumental to combine NREFTs and lattice. In fact, on one
hand, the NREFT is pulling out scales from observables in a controlled way,
separating them, and delivering results for calculations at large scales
including the logs resummation. On the other hand, the low energy
contributions that the EFT has factorized are often nonperturbative and can be
evaluated on the lattice. This greatly reduces the complexity of the problem
allowing the lattice to target the nonperturbative part directly,
appropriately defined in the EFT in terms of gauge invariant, purely gluonic
objects, a big simplification with respect to an ab initio lattice calculation
of an observable. The latter is much more difficult because still contains all
the physical scales of the problem. In this framework also quenched lattice
calculations can be pretty useful because at least at higher energy the flavor
dependence is accounted for by the EFT matching coefficients. We founded the
TUMQCD lattice collaboration precisely to address these calculations directly
inside the EFT tumqcd .

##  4 Nonrelativistic QCD

Nonrelativistic QCD (NRQCD) Caswell:1985ui ; Bodwin:1994jh , follows from QCD
integrating out the scale $m$. As a consequence, the effective Lagrangian is
organized as an expansion in $1/m$ and $\alpha_{\rm s}(m)$:

${\cal L}_{\rm NRQCD}=\sum_{n}\frac{c_{n}(\alpha_{\rm s}(m),\mu)}{m^{n}}\times
O_{n}(\mu,mv,mv^{2},...),$ (13)

where $c_{n}$ are Wilson coefficients that contain the contributions from the
scale $m$. They can be calculated via a well defined procedure call matching,
see e.g. Brambilla:2004jw , in which one calculates in perturbation theory
Green Functions or amplitudes in QCD and in NRQCD, expand in $1/m$ and insert
the difference between the two calculations (i.e. the non analytic terms in
the scale that is been integrated out, the mass) inisde the matching
coefficients of the EFT. In fact since in NRQCD we expand in the mass no
contributions of the form $\log(m/\mu)$ can be generated in calculations in
the EFT.

The $O_{n}$ are local operators of NRQCD; the matrix elements of these
operators contain the physics of scales below $m$, in particular of the scales
$mv$ and $mv^{2}$ and also of the nonperturbative scale $\Lambda_{\rm QCD}$.
Finally, the parameter $\mu$ is the NRQCD factorization scale. The low-energy
operators $O_{n}$ are constructed out of two or four heavy quark/antiquark
fields plus gluons. At the lowest order in the $1/m$ expansion the NRQCD
Lagrangian density is

$\displaystyle{\cal L}_{\rm NRQCD}\\!$ $\displaystyle=$
$\displaystyle\\!\psi^{\dagger}\\!\left\\{iD^{0}+{{\bf D}^{2}\over
2m}\right\\}\\!\psi+\chi^{\dagger}\\!\left\\{iD^{0}-{{\bf D}^{2}\over
2m}\right\\}\\!\chi$ (14) $\displaystyle-{1\over
4}F_{\mu\nu}^{a}F^{\mu\nu\,a}\,,$

where $\psi$ is the Pauli spinor field that annihilates the fermion and $\chi$
is the Pauli spinor field that creates the antifermion;
$iD^{0}=i\partial_{0}-gA^{0}$ and $i{\bf D}=i{\mbox{\boldmath$\nabla$}}+g{\bf
A}$. At order $1/m^{2}$ bilinear terms in the quark (antiquark) field
containing the chhromoelectric and chromomagnetic fields, as well covariant
derivatives and soin start to appear. These are of the same form that one
would obtain from a Foldy-Wouthuysen transformation but in addition are
multiplied by matching coefficients that contain the UV behaviour of QCD. In
addition, four fermion terms start to appear at order $1/m^{2}$ with matching
coefficients containing both real and imaginary parts. The presence of
imaginary parts describe the decays of quarkonium via annihilation of hard
gluons that have been integrated out from the theory.

The quarkonium state $|H\rangle$ in NRQCD is expanded in the number of partons

$|H\rangle=|\bar{Q}Q\rangle+|\bar{Q}Qg\rangle+|\bar{Q}Q\bar{q}q\rangle+\cdots$
(15)

where the states including one or more light parton are shown to be suppressed
by powers of $v$. In the $|\bar{Q}Qg\rangle$ for example the quark-antiquark
are in a color octet state. The NRQCD lagrangian has been extensively used on
the lattice to calculate quarkonium spectra and decays. On the other hand,
NRQCD has been deeply impactful on the study of quarkonium production at the
LHC putting forward a factorization formula for the inclusive cross section
for the direct production of the quarkonium $H$ at large transverse momentum
written as a sum of products of NRQCD matrix elements and short-distance
coefficients:

$\sigma[H]=\sum_{n}\sigma_{n}\langle{\cal K}_{n}^{4{\rm fermions}}\rangle$
(16)

where the $\sigma_{n}$ are short-distance coefficients, and the matrix
elements $\langle{\cal K}_{n}^{4{\rm fermions}}\rangle$ are vacuum-expectation
values containing four-fermions operators in color singlet and color octet
configurations identified by $Q\bar{Q}$angular momentum state quantum numbers.
They contain in the middle a projector over the $Q\bar{Q}$ pair plus anything
and some Schwinger lines with a particular path to make them gauge invariant.
The factorization for production is proved only at order NNLO and in some
cases Nayak:2005rt ; Bodwin:2019bpf . It is different from the NRQCD
factorization of inclusive decays that is proven at all orders. In such case
the low energy part contains quarkonium expectation values of color singlet
and color octet four quark operators without the projector in the middle. The
matrix elements (LDMEs, long distance matrix elements) $\langle{\cal
K}_{n}^{4{\rm fermions}}\rangle$ contain all of the nonperturbative physics
associated with the evolution of the $Q\bar{Q}$ pair into a quarkonium state
and they are universal. Their NRQCD definition does not allow a lattice
calculation and their extractions from collider data is still challenging
Brambilla:2010cs ; Chapon:2020heu ; Andronic:2015wma ; Bodwin:2013nua ;
Chung:2020gyb . Differently from HQET, the power counting of NRQCD is not
unique, due to the fact that many physical scales are still dynamical
($mv,mv^{2},\Lambda_{\rm QCD}$). This still complicates bound state
calculations as the soft and US scale can get entangled.

##  5  potential Nonrelativistic QCD

Quarkonium formation happens at the scale $mv$. The suitable NREFT is
potential Nonrelativistic QCD, pNRQCD Pineda:1997bj ; Brambilla:1999xf ;
Brambilla:2004jw , which follows from NRQCD by integrating out the scale
$mv\sim r^{-1}$. The pNRQCD description directly addresses the bound state
dynamics, implements the Schödinger equation as zero order problem, properly
defines the potentials as matching coefficients, and allows to systematically
calculate relativistic and retardation corrections. In this way even quantum
mechanics can also be reinterpreted as a pNREFT. Moreover, since the scale
$mv$ has been integrated out, the power counting of pNRQCD is less ambiguous
than the one of NRQCD.

The soft scale, proportional to the inverse quarkonium radius $r$, may be
either perturbative ($mv\gg\Lambda_{\rm QCD}$) or nonperturbative
($mv\sim\Lambda_{\rm QCD}$) depending on the physical system under
consideration. We do not have any direct information on the radius of the
quarkonia, and thus the attribution of some of the lowest bottomonia and
charmonia states to the perturbative or the nonperturbative soft regime may be
ambiguous, but it is likely that the lowest bottomonium and possibly also the
lowest charmonium states have small enough radii that the scale $mv$ is in
fact still perturbative. In a case with $\Lambda_{\rm QCD}\ll mv^{2}$ both
scales are still perturbative and the system would be similar to a Coulombic
system, For such a quarkonium we would have - similar as for the hydrogen atom
- $v\sim\alpha_{s}(mv)$. However, none of the $b\bar{b}$ and $c\bar{c}$ states
satisfy this condition. In all realistic quarkonia ($b\bar{b}$ and $c\bar{c}$)
the ultrasoft scale is nonperturbative. Only for $t\bar{t}$ threshold states
the ultrasoft scale may be considered still perturbative. In Fig. 5 we
schematically show the various scales. The short-range quarkonia are small
enough to allow for a perturbative treatment of the scale $mv$, while for the
long range quarkonia already this scale requires a non perturbative treatment.

In the next section we will present weakly coupled and strongly coupled
pNRQCD, discuss the form of the potential and briefly address the countless
phenomenological applications of this that has become the standard treatment
for quarkonium. On one hand, pNRQCD allows us systematic and precise
determinations of processes at high energy collider experiments and the
definition and computation of quantities of large phenomenological interest,
on the other hand the systematic factorization enables us for studies of
confinement. The low energy nonperturbative factorized effects depend on the
size of the physical system: when $\Lambda_{\rm QCD}$ is comparable or smaller
than $mv$ they are carried by correlators of chromoelectric or chromomagnetic
fields nonocal or local in time, when $\Lambda_{\rm QCD}$ is of order $mv$
they are carried by generalized static Wilson loops, non local in space, with
insertion on chromoelectric and chromomagnetic fields. The EFT allows us to
make model independent predictions and we can use the power counting to attach
an error to the theoretical predictions. The nonperturbative physics in pNRQCD
is encoded in few low energy correlators that depend only on the glue and are
gauge invariant: these are objects in principle ideal for lattice
calculations. Even more interesting is that this EFT description allows
modifications to be used to describe exotics states (BOEFT) and the
nonequilibrium evolution of quarkonium in medium (pNRQCD at finite temperature
with open quantum system) as we will address in the last two sections.

###  5.1 Weakly coupled pNRQCD

When $mv\gg\Lambda_{\rm QCD}$, we speak about weakly-coupled pNRQCD because
the soft scale is perturbative and the matching from NRQCD to pNRQCD may be
performed in perturbation theory. The lowest levels of quarkonium, like
$J/\psi$, $\Upsilon(1S),\Upsilon(2S)\dots$, may be described by weakly coupled
pNRQCD, while the radii of the excited states are larger and presumably need
to be described by strongly coupled pNRQCD. All this is valid for states away
from strong-decay threshold, i.e. the threshold for a decay into two heavy-
light hadrons. The effective Lagrangian is organized as an expansion in $1/m$
and $\alpha_{\rm s}(m)$, inherited from NRQCD, and an expansion in $r$
(multipole expansion) Brambilla:1999xf ; Brambilla:2004jw

$L_{\rm pNRQCD}^{\rm weak}=\int d^{3}R\,\int
d^{3}r\,\sum_{n}\sum_{k}\frac{c_{n}(\alpha_{\rm
s}(m),\mu)}{m^{n}}V_{n,k}(r,\mu^{\prime},\mu)\;r^{k}\times
O_{k}(\mu^{\prime},mv^{2},...),$ (17)

where ${\bf R}$ is the center of mass position and $O_{k}$ are the operators
of pNRQCD. We do not show explicitly the part of the lagrangian involving
gluons and light quarks, as this part coincides with the QCD one. The matrix
elements of the operators above depend on the low-energy scale $mv^{2}$ and
$\mu^{\prime}$, where $\mu^{\prime}$ is the pNRQCD factorization scale. The
$V_{n,k}$ are the Wilson coefficients of pNRQCD that encode the contributions
from the scale $r$ and are nonanalytic in $r$. The $c_{n}$ are the NRQCD
matching coefficients as given in (13).

The degrees of freedom, which are relevant below the soft scale, and which
appear in the operators $O_{k}$, are $Q\overline{Q}$ states (a color-singlet
$S$ and a color-octet $O=O_{a}T^{a}$ state, depending on ${\bf r}$ and ${\bf
R}$) and (ultrasoft) gluon fields, which are expanded in $r$ as well
(multipole expanded and depending only on ${\bf R}$). pNRQCD makes apparent
that the correct zero order problem is the Schrödinger equation. Looking at
the equations of motion of pNRQCD, we may identify $V_{n,0}=V_{n}$ with the
$1/m^{n}$ potentials that enter the Schrödinger equation and $V_{n,k\neq 0}$
with the couplings of the ultrasoft degrees of freedom, which provide
corrections to the Schrödinger equation. These nonpotential interactions,
associated with the propagation of low-energy degrees of freedom start to
contribute at NLO in the multipole expansion. They are typically related to
nonperturbative effects and are carried by purely gluonic correlator local or
nonlocal in time: they need to be calculated on the lattice.

In particular at the NLO in the multipole expansion at leading order in the
expansion in $1/m$, the pNRQCD Lagrangian density is

$\displaystyle{\cal L}_{\rm pNRQCD}^{\rm weak}={\rm Tr}\Biggl{\\{}{\rm
S}^{\dagger}\left(i\partial_{0}-{{\bf p}^{2}\over m}-V_{s}(r)+\dots\right){\rm
S}$ $\displaystyle\qquad+{\rm O}^{\dagger}\left(iD_{0}-{{\bf p}^{2}\over
m}-V_{o}(r)+\dots\right){\rm O}\Biggr{\\}}$ (18)
$\displaystyle\qquad+gV_{A}(r){\rm Tr}\left\\{{\rm O}^{\dagger}{\bf
r}\cdot{\bf E}\,{\rm S}+{\rm S}^{\dagger}{\bf r}\cdot{\bf E}\,{\rm
O}\right\\}+g{V_{B}(r)\over 2}{\rm Tr}\left\\{{\rm O}^{\dagger}{\bf
r}\cdot{\bf E}\,{\rm O}+{\rm O}^{\dagger}{\rm O}{\bf r}\cdot{\bf E}\right\\},$

where ${\bf R}\equiv({\bf x}_{1}+{\bf x}_{2})/2$, ${\rm S}={\rm S}({\bf
r},{\bf R},t)$ and ${\rm O}={\rm O}({\bf r},{\bf R},t)$ are the singlet and
octet quark-antiquark composite nonrelativistic fields respectively. All the
gauge fields in Eq. (5.1) are evaluated in ${\bf R}$ and $t$. In particular
${\bf E}\equiv{\bf E}({\bf R},t)$ and $iD_{0}{\rm O}\equiv i\partial_{0}{\rm
O}-g[A_{0}({\bf R},t),{\rm O}]$. $V_{s}^{0}$ and $V_{o}^{0}$ are the singlet
and octet heavy $Q$-$\bar{Q}$ static potentials Higher-order potentials in the
$1/m$ expansion and the centre-of-mass kinetic term are not shown here and are
indicated by the dots. At higher order in the multipole expansion and in the
$1/m$ expansion more operators appears containing US chromomagnetic and
chromoelectric fields but they are suppressed in the power counting in $v$. We
define

$\displaystyle V_{s}^{0}(r)$ $\displaystyle\equiv$ $\displaystyle-
C_{F}{\alpha_{V_{s}}(r)\over r},$ (19) $\displaystyle V_{o}^{0}(r)$
$\displaystyle\equiv$ $\displaystyle\left({C_{A}\over
2}-C_{F}\right){\alpha_{V_{o}}(r)\over r}$

which shows that the singlet potential is attractive and the octet one is
repulsive. $V_{A}$ and $V_{B}$ are the matching coefficients associated in the
Lagrangian (5.1) to the leading corrections in the multipole expansion. Both
the potentials and the coefficients $V_{A}$ and $V_{B}$ have to be determined
by matching pNRQCD with NRQCD at a scale $\mu$ smaller than $mv$ and larger
than the US scales. Since, in particular, $\mu$ is larger than $\Lambda_{\rm
QCD}$ the matching can be done perturbatively. At the lowest order in the
coupling constant we get $\alpha_{V_{s}}=\alpha_{V_{o}}=\alpha_{\rm s}$,
$V_{A}=V_{B}=1$. In order to have the proper free-field normalization in the
colour space we define

${\rm S}\equiv{1\\!\\!{\rm l}_{c}\over\sqrt{N_{c}}}S\quad\quad{\rm
O}\equiv{T^{a}\over\sqrt{T_{F}}}O^{a},$ (20)

where $T_{F}=1/2,N_{c}=3$.

The perturbative QCD potential

For system with a perturbative soft scale the potentials can be calculated in
perturbation theory, i.e. no nonperturbative quantities enter the potential
Brambilla:2004jw . The potentials can be calculated at all order in the
perturbative expansion in $\alpha_{\rm s}$ via a well defined matching
procedure that entails to compare appropriate Green functions in NRQCD and in
pNRQCD calculated in perturbation theory up to the desired order of the
perturbative expansion, order by order in $1/m$ and order by order in the
multipole expansion. The difference between the two is encoded in the matching
coefficients. On the pNRCD side lower energy scales may be expanded in the
loop integral and give no contribution in the matching in dimensional
regularization. The singlet static potential is calculated by matching the
relevant NRQCD green function (the static Wilson loop) and the pNRQCD singlet
Green function see Fig. 6). At three loops, an ultrasoft (US) divergence is
emerging at fixed order calculation that can be regularized resumming series
of diagrams containing the US scale which is the difference between the
singlet and the octet potential, see Brambilla:1999qa ; Brambilla:2004jw .

Figure 6: The matching of the static potential. On the right side are the
pNRQCD fields: simple lines are singlet propagator, double lines are octet
propagators, circled-crosses are the singlet-octet vertices of Eq. (5.1) and
the wavy line is the US gluon propagator.

This solves the problem raised by Appelquist Dine and Muzinich
Appelquist:1977tw , it explains how this divergence cancels between the
potential ($\log r\mu$) and the contribution of the US chomoelectric
correlator $\log((V_{o}-V_{S})/\mu)$ (second contribution in the r.h. s. of
Fig. (6)) leaving a term in $\log\alpha_{s}$ in the static energy and
qualifies the potential as a matching coefficient of an EFT which is depending
on the US scale. Being matching coefficients of the effective field theory,
the potentials undergo renormalization, develop a scale dependence and satisfy
renormalization group equations, which allow to resum large US logarithms see
e.g. Pineda:2000gza ; Brambilla:2006wp ; Brambilla:2009bi .

In particular one obtains the complete expression for the static potential, at
3 loops, up to the relative order $\alpha_{\rm s}^{3}\ln\mu r$ in coordinate
space ($\alpha_{\rm s}$ is in the $\overline{\rm MS}$ scheme) as:

$\displaystyle{\alpha}_{V_{s}}(r,\mu)=\alpha_{\rm
s}(r)\left\\{1+\left(a_{1}+2{\gamma_{E}\beta_{0}}\right){\alpha_{\rm
s}(r)\over 4\pi}\right.$ (21)
$\displaystyle+\left[\gamma_{E}\left(4a_{1}\beta_{0}+2{\beta_{1}}\right)+\left({\pi^{2}\over
3}+4\gamma_{E}^{2}\right){\beta_{0}^{2}}+a_{2}\right]{\alpha_{\rm
s}^{2}(r)\over 16\,\pi^{2}}\left.+{C_{A}^{3}\over 12}{\alpha_{\rm
s}^{3}(r)\over\pi}\ln{r\mu}\right\\},$

where $\beta_{n}$ are the coefficients of the beta function, $a_{1}$ and
$a_{2}$ are the constants at 1 and 2 loops respectively. We emphasize that
this US contribution to the static potential would be zero in QED. From Eq.
(5.1) is clear that $\alpha_{V_{s}}$ depends on the US scale and as such is
not a short distance quantity.

The evaluated terms clarify the long-standing issue of how the perturbative
static potential should be defined at higher order in the perturbative series.
In perturbation theory he static potential does not coincide with the static
Wilson loop starting from three loops. It should be emphasized that the
separation between soft and US contributions is not an artificial trick but a
necessary procedure if one wants to use the static potential in a Schrödinger-
like equation in order to study the dynamics of $Q$-$\bar{Q}$ states of large
but finite mass. In that equation the kinetic term of the $Q$-$\bar{Q}$ system
is US and so is the energy. Since the US gluons interact with the
$Q$-$\bar{Q}$ system, their dynamics is sensitive to the energies of the (non-
static) system and hence it is not correct to include them in the static
potential. When calculating a physical observable the $\mu$ dependence must
cancel against $\mu$-dependent contributions coming from the US gluons.
Finally it is worth mentioning that the static potential suffers from IR
renormalons ambiguities with the following structure

$\delta V_{s}\sim C+C_{2}r^{2}+.\,.\,.$

The constant $C$ is known to be cancelled by the IR pole mass renormalon
Brambilla:2004jw . The static singlet potential is currently know at NNNLL
Brambilla:2009bi , the only unknown at four loops is the constant term.

The static energy in the perturbative regime has the form

$E_{0}(r)=V_{s}^{0}(r;\mu)+\delta_{\rm US}(r,V_{s}^{0},V_{o}^{0},\dots;\mu),$
(22)

where $\delta_{\rm US}(r,V_{s},V_{o},V_{A},V_{B},\dots;\mu)$ ($\delta_{\rm
US}$ for short) contains contributions from the ultrasoft gluons. As
mentioned, $V_{s}(r;\mu)$ and $V_{o}(r;\mu)$ do not depend on $\mu$ up to N2LO
Brambilla:1999xf . The former coincides with $E_{0}(r)$ at this order and the
latter may be found in Kniehl:2004rk . The fact that the $\mu$ dependence of
$\delta_{\rm US}$ must cancel the one in $V_{s}(r;\mu)$ is the key observation
that leads to a drastic simplification in the calculation of the
$\log\alpha_{\rm s}$ terms in $E_{0}(r)$. So, for instance, the logarithmic
contribution at N3LO, which is part of the three-loop contributions to
$V_{s}(r;\mu)$, may be extracted from a one-loop calculation of $\delta_{\rm
US}$ Brambilla:1999qa ; Brambilla:1999xf  and the single logarithmic
contribution at N4LO, which is part of the four-loop contributions to
$V_{s}(r;\mu)$, may be extracted from a two-loop calculation of $\delta_{\rm
US}$. Notice that we denote NnLO, contributions to the potential of order
$\alpha_{\rm s}^{n+1}$ and NnLL, contributions of order $\alpha_{\rm
s}^{n+2}\log^{n-1}\alpha_{\rm s}$.

The general form of the relativistic corrections to the singlet potential in
the center of mass is (we drop the index $s$ for simplicity):

$\displaystyle V(r)$ $\displaystyle=$ $\displaystyle
V^{(0)}(r)+{V^{(1)}(r)\over m}+{V^{(2)}\over m^{2}}+\cdots,$ (23)
$\displaystyle V^{(2)}$ $\displaystyle=$ $\displaystyle
V^{(2)}_{SD}+V^{(2)}_{SI},$ $\displaystyle V^{(2)}_{SI}$ $\displaystyle=$
$\displaystyle{1\over 2}\left\\{{\bf p}^{2},V_{{\bf
p}^{2}}^{(2)}(r)\right\\}+{V_{{\bf L}^{2}}^{(2)}(r)\over r^{2}}{\bf
L}^{2}+V_{r}^{(2)}(r),$ (24) $\displaystyle V^{(2)}_{SD}$ $\displaystyle=$
$\displaystyle V_{LS}^{(2)}(r){\bf L}\cdot{\bf S}+V_{S^{2}}^{(2)}(r){\bf
S}^{2}+V_{{\bf S}_{12}}^{(2)}(r){\bf S}_{12}({\hat{\bf r}}),$

${\bf S}={\bf S}_{1}+{\bf S}_{2}$, ${\bf L}={\bf r}\times{\bf p}$, ${\bf
S}_{1}=\mbox{\boldmath$\sigma$}_{1}/2$, ${\bf
S}_{2}=\mbox{\boldmath$\sigma$}_{2}/2$, and ${\bf S}_{12}({\hat{\bf r}})\equiv
3{\hat{\bf r}}\cdot\mbox{\boldmath$\sigma$}_{1}\,{\hat{\bf
r}}\cdot\mbox{\boldmath$\sigma$}_{2}-\mbox{\boldmath$\sigma$}_{1}\cdot\mbox{\boldmath$\sigma$}_{2}$.

We see that differently from what discussed previously the corrections to the
potential start at order $1/m$. The potential proportional to $V_{LS}^{(2)}$
may be identified with the spin-orbit potential, the potential proportional to
$V_{S^{2}}^{(2)}$ with the spin-spin potential and the potential proportional
to $V_{S_{12}}^{(2)}$ with the spin tensor potential. The above potentials
read at leading (non-vanishing) order in perturbation theory (see, e.g., Ref.
Brambilla:2004jw ):

$\displaystyle V^{(1)}(r)=-\frac{2\alpha_{s}^{2}}{r^{2}}\,,$ (25)
$\displaystyle V_{r}^{(2)}(r)=\frac{4\pi}{3}\alpha_{s}\delta^{(3)}({\bf
r}\,)\,,\hskip 25.6073ptV_{p^{2}}^{(2)}(r)=-\frac{4\alpha_{s}}{3r}\,,\hskip
25.6073ptV_{L^{2}}^{(2)}(r)=\frac{2\alpha_{s}}{3r}\,,$ (26) $\displaystyle
V_{LS}^{(2)}(r)=\frac{2\alpha_{s}}{r^{3}}\,,\hskip
28.45274ptV_{S^{2}}^{(2)}(r)=\frac{16\pi\alpha_{s}}{9}\delta^{(3)}({\bf
r}\,)\,,\hskip 28.45274ptV_{S_{12}}^{(2)}(r)=\frac{\alpha_{s}}{3r^{3}}\,.$
(27)

These potentials and the static octet potential have been calculated at higher
order in perturbation theory see e.g. Brambilla:2004jw ; Pineda:2011dg ;
Kniehl:2004rk . Once one has obtained the potentials as matching coefficients
of pNRQCD one can start doing calculations inside the EFT obtaining for
example the energy levels of quarkonium Brambilla:1999xj ; Kniehl:2002br .
Nonperturbative effects in the form of local or time nonlocal chromoeletric
correlators start to appear at order $\alpha_{\rm s}^{5}$ i.e. at NNNLO
Brambilla:2004jw ; Pineda:2011dg : their contribution is suppressed.

The case of baryons made by three heavy quarks or two heavy quarks and a light
quark has also been addressed in pNRQCD Brambilla:2005yk  and the static three
quark potential has been calculated at one loop accuracy, showing some
interesting features Brambilla:2009cd .

###  5.2 Strongly coupled pNRQCD

When $mv\sim\Lambda_{\rm QCD}$, we speak about strongly-coupled pNRQCD because
the soft scale is nonperturbative and the matching from NRQCD to pNRQCD is
nonperturbative i.e. it cannot be performed within an expansion in
$\alpha_{s}$. In this case, since we work at a nonperturbative scale, only
color singlet degrees of freedom remain dynamical and they include $Q\bar{Q}$
states, hybrids $Q\bar{Q}g$ states and glueballs (and if we consider light
quarks as part of the binding, color singlet states formed by a heavy quark
and heavy antiquark and light quarks, as in the case of tetraquarks). Since
the physics is nonperturbative we need to use input coming from the lattice to
construct the EFT. In particular we can use the lattice evaluation of the
gluonic static energies of a $Q\bar{Q}$ pair: they have been calculated since
long Juge:2002br ; Bali:2003jq  and recently updated Muller:2019joq . Such
lattice calculations use insertions of gluonic fields at the initial and final
Schwinger strings, inside a generalized static Wilson loop, to select some
given symmetries. The gluonic static energies, $E_{\Gamma}$ in Fig. 8, are
classified according to representations of the symmetry group $D_{\infty\,h}$,
typical of diatomic molecules, and labeled by $\Lambda_{\eta}^{\sigma}$ (see
Fig. 8): $\Lambda$ is the rotational quantum number $|\hat{\bf r}\cdot{\bf
K}|=0,1,2,\dots$, with ${\bf K}$ the angular momentum of the gluons, that
corresponds to $\Lambda=\Sigma,\Pi,\Delta,\dots$; $\eta$ is the CP eigenvalue
($+1\equiv g$ (gerade) and $-1\equiv u$ (ungerade)); $\sigma$ is the
eigenvalue of reflection with respect to a plane passing through the
$Q\bar{Q}$ axis. The quantum number $\sigma$ is relevant only for $\Sigma$
states. In general there can be more than one state for each irreducible
representation of $D_{\infty\,h}$: higher states are denoted by primes, e.g.,
$\Pi_{u}$, $\Pi_{u}^{\prime}$, $\Pi_{u}^{\prime\prime},\dots$ Notice that this
set of static energies will be fundamental later to address the exotics.

Figure 7: Quarkonium hybrid symmetries.

Figure 8: Hybrid static energies, $E_{\Gamma}$, in lattice units, from
Juge:2002br .

Since at this moment we are dealing with $Q\bar{Q}$ states below threshold we
are interested merely in the $\Sigma_{g}^{+}$ static energy (that in this case
coincides with the static singlet $Q\bar{Q}$ potential) and in the information
that such curve develops a gap of order $\Lambda_{\rm QCD}$ at a distance
$r\sim\Lambda_{\rm QCD}^{-1}$: therefore all the hybrid static energies can be
integrated out as pNRQCD follows from integrating out all degrees of freedom
with energy up to $mv^{2}$. The quarkonium singlet field $S$ is now the only
low-energy dynamical degree of freedom in the pNRQCD Lagrangian (up to US
pions) which reads Brambilla:2000gk ; Pineda:2000sz ; Brambilla:2004jw :

$\quad L_{\rm pNRQCD}=\int d^{3}R\,\int
d^{3}r\,{S}^{\dagger}\left(i\partial_{0}-\frac{{\bf
p}^{2}}{2m}-V_{S}(r)\right){S}\,$ (28)

and lends support to potentials models in this particular regime. Indeed Eq.
(28) originates a Schödinger equation governed by the singlet potential.

The nonperturbative QCD potential

The difference with the phenomenological potential models is that now the
singlet potential $V_{s}(r)=V_{0}+V^{(1)}/m+V^{(2)}/m^{2}$, is the QCD
potential calculated in the matching obtained in Brambilla:2000gk ;
Pineda:2000sz  and display novel characteristics. These potentials are
nonperturbative and are given in terms of generalized Wilson loops. However
relevant differences emerge also with respect to the Wilson loop approach. The
general decomposition of the potential is of the form of Eq. (23). It features
a contribution already at order $1/m$ Brambilla:2000gk  given in terms of a
novel expectation value of a static Wilson loop with insertion of a
chromoelectric field. The potential at order $1/m^{2}$ appears factorized in
the product of NRQCD matching coefficients, carrying contribution in the
$\log(m/\mu)$ and generalized Wilson loops with chromomagnetic and
chromoelectric insertions. This solves the problem of the incompatibility of
the Wilson loop potentials and the perturbative calculation. Moreover,
additional generalized Wilson loops contribution of novel type emerge see
Pineda:2000sz ; Brambilla:2004jw  for all the explicit expressions. Spin
effects in quarkonium are appearing only at order $1/m^{2}$ and are therefore
suppressed in the spectrum and in the transitions: we will see that things are
different for exotics. Some of these generalized Wilson loops have been
calculated on the lattice (only quenched) Bali:1997am ; Bali:2000gf ;
Koma:2006si ; Koma:2007jq  but some contributions at order $1/m^{2}$ are still
to be calculated. The EFT then lends a clean definition and an interpretation
of the static Wilson loops measured on the lattice as actual potentials in
this regime, together with a prescription to use them to calculate
observables. Once the nonperturbative potentials are given in terms of
generalized Wilson loops one can use a model of low energy QCD to evaluate
them. A minimal area law for example gives for the potentials of
Eq.(LABEL:pott) a linear term for $V^{(0)}(r)=\sigma r$ and
$V_{L^{2}}^{(2)}(r)=-{\sigma\over 6}r$ and nonzero contributions to
$V_{r}^{(2)}(r)$ and $V_{LS}^{(2)}(r)$. The velocity dependent relativistic
nonperturbative contribution corresponds to the angular momentum and the
energy of the flux tube between the quark and the antiquark while the one in
the spin dependent part is relevant to obtain agreement with the fine
separations of the quarkonia multiplets (as it was pursued previously using a
scalar confinement kernel in the Bethe-Salpeter equation). These findings are
in agreement with what is obtained from the lattice calculations of the Wilson
loops. Using a QCD effective string model add several further subleading
corrections see Brambilla:2014eaa ; Perez-Nadal:2008wtr .

Using these potentials, all the masses for heavy quarkonia below threshold can
be obtained by solving the Schrödinger equation with such potentials. Lorentz
invariance is still there in the form of exact relations among potentials and
it has been observed on the lattice. Decays are described by calculating the
imaginary parts of the potentials Brambilla:2002nu  where nonperturbative
contributions enter in the form of gauge invariant time non local
chromoelectric and chromomagnetic correlators that have still to be calculated
on the lattice. Summarizing, strongly coupled pNRQCD factorizes low energy
nonperturbative contributions in terms of generalised gauge invariant Wilson
loops, opening the way to a systematic study of the confinement mechanism and
systematic applications to quarkonium spectrum and decay as we will see in the
dedicated section.

The three heavy quark static potential has been studied on the lattice using
Wilson loops Takahashi:2000te ; Takahashi:2002bw .

###  5.3 Spectra, transitions, decays and production, SM parameters
extractions

Together with lattice, pNRQCD is the theoretical framework that is nowadays
mostly used for calculations and predictions of quarkonium properties. The
power counting of the EFT allows to attach an error to each prediction. In the
regime in which the soft scale is perturbative, pNRQCD enables precise and
systematic higher order calculations on bound state allowing the extraction of
precise determinations of standard model parameters like the quark masses and
$\alpha_{s}$. For example, based on Eq. (22), it has been possible to use
lattice calculations of the static energy with 2+1 flavour and the NNNLO
pNRQCD calculation of the static energy, including the US log resummation, to
extract a precise determination of $\alpha_{\rm s}$ at rather low energy and
run it at the mass of the Z, obtaining $\alpha_{\rm
s}(M_{Z})=0.11660^{+0.00110}_{-0.00056}$, which is a competitive extraction
made at a pretty high orders of the perturbative expansion Bazavov:2019qoo ;
Brambilla:2009bi ; Bazavov:2012ka . This method of $\alpha_{\rm s}$ extraction
is now used by several groups, see e.g. Ayala:2020odx ; Takaura:2018vcy  and
the force, defined in terms of a single chromolectric insertion in the Wilson
loop could be used as well Brambilla:2021wqs ; Vairo:2016pxb .

In the same way one can extract precise determinations of the bottom and charm
masses using the experimental measurements of the mass of the lowest states
and comparing it to the formula for the energies in pNRQCD at NNNLO, which
depends on the mass in a given scheme. The renormalon ambiguity cancels
between the mass and the static potential and a pretty good determination is
possible see e.g. Peset:2018ria ; Kiyo:2014uca  and references therein.
Moreover, the energy levels of some of the lowest quarkonia states have been
obtained at high orders in perturbation theory using weakly coupled pNRQCD
showing that the constant separation on the energy levels may be generated in
this way, see e.g. Brambilla:2001qk ; Brambilla:2001qk ; Brambilla:2000db ;
Peset:2015vvi  Electromagnetic M1 and E1 transitions have been calculated in
pNRQCD, see e.g. Brambilla:2005zw ; Brambilla:2012be ; Pineda:2013lta ;
Segovia:2018qzb . There are so many results that it is impossible to discuss
all of them here and we refer you to some reviews Pineda:2011dg ;
Brambilla:2004jw ; Brambilla:2010cs ; QuarkoniumWorkingGroup:2004kpm ;
Brambilla:2014jmp .

For what concern production it is very promising that by factorizing the
quarkonium productions cross section at lower energy in pNRQCD
Brambilla:2021abf ; Brambilla:2020ojz ; Brambilla:2020xod , one can rewrite
the octet NRQCD LDMEs, which are the nonperturbative unknowns, in terms of
product of wave functions and gauge invariant low energy correlators that
depend only on the glue and not the on flavor quantum numbers. This allows us
to reduce by half the number of nonperturbative unknowns and promises to have
a great impact in the progress of the field.

Figure 9: The origin of the linear potential between the static quark and
antiquark may be traced back to a flux tube: a string of gluon energy between
the quark pair. Here we present the hystorical picture of the action density
distribution between a static quark antiquark couple in $SU(2)$ at a physical
distance of 1.2 fm, from Bali:1994de .

##  6 Confinement and low energy QCD models

Strongly coupled pNRQCD realizes a scale factorization encoding the low energy
physics in generalized gauge invariant Wilson loops, i.e. Wilson loops with
insertions of chromoelectric and chromomagnetic fields. Such objects no longer
depends on the heavy quark degrees of freedom and on the quark flavor. It
turns out that this is a systematic method to study the QCD confinement
properties and put them directly in relation to the quarkonium phenomenology.
Indeed lattice QCD seems more suitable to ask ’what’ instead of ’why’ and to
understand the mechanism underlying confinement it may be useful to build
models of low energy QCD and compare to the lattice results: the interface is
offered by the EFT. We have seen that the area law emerging in the static
Wilson loop at large distance is responsible of confinement, this in turn
corresponds to the formation of a chromolectric flux tube between the quark
and the antiquark that sweeps the area of the Wilson loop, see Fig. 9. This
effects is originated by the nonperturbative QCD vacuum that could be imagined
as a disordered medium with whirlpools of colour on different scales, thus
densely populated by fluctuating fields whose amplitude is so large that they
cannot be described by perturbation theory. A QCD vacuum model can be
established by making an assumption on the behaviour of the Wilson loop that
gives the static potential. The relativistic corrections that involve
insertions of gluonic fields in the Wilson loop follow then via functional
derivative with respect to the quark path see Brambilla:1993zw ;
Brambilla:1999ja . Then, one can notice for example that the part proportional
to the square of the angular momentum in the of the $V_{vd}$ potential
obtained in strongly coupled pNRQCD takes into account the energy and the
angular momentum of the flux tube, which is something that could not be
obtained in any Bethe Salpeter approach with a confinement scalar convolution
kernel. Lattice simulations of the action density or the energy density
between the static quark and the static antiquark show indeed a chromoelectrix
flux tube configuration see Fig. 9 from ref. Bali:1994de  and more recent
calculations in Bicudo:2018jbb ; Yanagihara:2019foh ; Baker:2021jnr . The
mechanism underlying confinement and flux tube formation has been investigated
since long on the lattice Greensite:2003bk  using the Wilson loops and the ’t
Hooft abelian projection, to identify the roles of magnetic monopoles
Amemiya:1998jz ; Sasaki:1994sa  and center vortices, see e.g. the review
Brambilla:2014jmp .

In the continuum several models of low energy QCD has been used to explain the
flux tube formation ranging from the dual Meissner effect and a dual abelian
Higgs model picture, dual QCD Baker:1991bc , the stochastic vacuum
Dosch:1988ha , the flux tube model Isgur:1984bm  and an effective QCD string
description. Each of these models can be used to obtain analytic estimates of
the behaviour of the generalized Wilson loops for large distance, which in
turn give the static potential and the relativistic corrections
$V_{1},V_{sd},V_{vd}$ as function of $r$, see Baker:1996mk ; Brambilla:1996aq
; Baker:1998jw ; Brambilla:2014eaa ; Perez-Nadal:2008wtr . Similar
nonperturbative configurations leading to confinement can be studied analyzing
the Wilson loop in case of baryons with three or two heavy quarks Nawa:2006gv
; Soto:2021cgk .

##  7 BOEFT and X Y Z exotics

Exotic states, i.e. states with a composition different from a quark-antiquark
or three quarks in a color singlet, have been predicted before and after the
advent of QCD. In the last decades a large number of states, either with a
manifest different composition (with isospin and electric charge different
from zero) or with other exotic characteristics have been observed in the
sector with two heavy quarks $Q\bar{Q}$, at or above the quarkonium strong
decay threshold at the B-Factories, tau-charm and LHC and Tevatron collider
experiments Brambilla:2019esw ; Brambilla:2010cs ;
QuarkoniumWorkingGroup:2004kpm ; Brambilla:2021mpo . These states have been
termed $X,Y,Z$ in the discovery publications, without any special criterion,
apart from $Y$ being used for exotics with vector quantum numbers, i.e.,
$J^{PC}=1^{--}$ . Meanwhile, the Particle Data Group (PDG) has proposed a new
naming scheme ParticleDataGroup:2018ovx , that extends the scheme used for
ordinary quarkonia, in which the new names carry information on the $J^{PC}$
quantum numbers, see Brambilla:2019esw  for more details. Some of these
exotics have quantum numbers that cannot be obtained with ordinary hadrons. In
this case, the identification of these states as exotic is straightforward. In
the other cases, the distinction requires a careful analysis of experimental
observations and theoretical predictions. Of course all hadrons should be
color singlets but allowing combinations beyond the $Q\bar{Q}$ and $QQQ$ in
the two heavy quarks sector calls for tetraquarks like $Q\bar{Q}q\bar{q}$,
$QQ\bar{q}\bar{q}$, pentaquarks like $Q\bar{Q}qqq$, hybrids $Q\bar{Q}g$ and so
on Ali:2017jda  . We have observed these exotics up to now only in the sector
with two heavy quarks likely due to the fact that the presence of the two
heavy quarks stabilizes them. Some of the discovered states have an
unprecedentely small width even if they are at or above the strong decay
threshold. $XYZ$ states offer us unique possibilities for the investigation of
the dynamical properties of strongly correlated systems in QCD: we should
develop the tools to gain a solid interpretation from the underlying field
theory, QCD. This is a very significant problem with trade off to other fields
featuring strong correlations and a pretty interesting connections to heavy
ion physics, as propagation of these states in medium may help us to
scrutinize their properties.

Since the new quarkonium revolution i.e. the discovery of the first exotic
state, the $X(3872)$ at BELLE in 2003 Belle:2003nnu , a wealth of theoretical
papers appeared to supply interpretation and understanding of the
characteristics of the exotics. Many models are based on the choice of some
dominant degrees of freedom and an assumption on the related interaction
hamiltonian. An effective field theory molecular description of some of these
states particularly close to threshold was also put forward, see e.g.
AlFiky:2005jd ; Guo:2017jvc ; Braaten:2007dw ; Braaten:2003he ; Fleming:2011xa
. A priori the simplest system consisting of only two quarks and two
antiquarks (generically called tetraquarks) is already a very complicated
object and it is unclear whether or not any kind of clustering occurs in it.
However, to simplify the problem it is common to focus on certain
substructures and investigate their implications: in hadroquarkonia the heavy
quark and antiquark form a compact core surrounded by a light-quark cloud; in
compact tetraquarks the relevant degrees of freedom are compact diquarks and
antidiquarks; in the molecular picture two color singlet mesons are
interacting at some typical distance, for a review see Brambilla:2019esw .
Discussions about exotics therefore often concentrate on the ’pictures’ of the
states, like for example the tetraquark interpretation against the molecular
one (of which both several different realizations exist). However, as a matter
of fact all the light degrees of freedom (light quarks, glue, in different
configurations) should be there in QCD close or above the strong decay
threshold: it is a result of the strong dynamics which one sets in and which
configuration dominates in a given regime.

Even in an ordinary quarkonium or in a heavy baryon, which has a dominant
$Q\bar{Q}$ or $QQQ$ configuration, subleading contributions of the quantum
field theoretical Fock space may contribute, with have additional quark-
antiquark pairs and active gluons. However, in the most interesting region,
close or above the strong decay threshold, where the X Y Z have been
discovered, the situation is much more complicated: there is no mass gap
between quarkonium and the creation of a heavy-light mesons couple, nor to
gluon excitations, therefore many additional states built on the light quark
quantum numbers may appear Brambilla:2021mpo . Still, $m$ is a large scale and
a first scale factorization is applicable so that Nonrelativistic QCD is still
valid. Then, if we want to introduce a description of the bound state similar
to pNRQCD, making apparent that the zero order problem is the Schrödinger
equation, we can still count on another scale separation. Let us consider
bound states of two nonrelativistic particles and some light d.o.f., e.g.
molecules in QED or quarkonium hybrids ($Q\bar{Q}g$ states) or tetraquarks
($Q\bar{Q}q\bar{q}$ states) in QCD: electron/gluon fields/light quarks change
adiabatically in the presence of heavy quarks/nuclei. The heavy quarks/nuclei
interaction may be described at leading order in the nonrelativistic expansion
by an effective static energy (or potential) $E_{\kappa}$ between the static
sources where $\kappa$ labels different excitations of the light degrees of
freedom. A plethora of states can be built on each on the static energies
$E_{\kappa}$ by solving the corresponding Schrödinger equation. This picture
corresponds to the Born-Oppenheimer (BO) approximation Juge:2002br ;
Braaten:2014qka ; Braaten:2014ita . Starting from pNRQED/pNRQCD the BO
approximation can be made rigorous and cast into a suitable EFT called Born-
Oppenheimer EFT (BOEFT) Berwein:2015vca ; Brambilla:2017uyf ; Oncala:2017hop ;
Soto:2020xpm ; Brambilla:2019jfi ; Brambilla:2018pyn ; Brambilla:2019esw
which exploits the hierarchy of scales $\Lambda_{\rm QCD}\gg mv^{2}$, $v$
being the velocity of the heavy quark.

Figure 10: Pictorial view of electronic static energies in QED, labelled by a
collective quantum number $\kappa$.

Figure 11: Pictorial view of the gluonic (or hybrid) static energies,
$E_{\Gamma}$, in QCD. The collective quantum number $\kappa$ has been detailed
in $\Lambda^{\sigma}_{\eta}$ as explained in the section on strongly coupled
pNRQCD. It corresponds to the actual lattice results in Fig. 8.

In Berwein:2015vca  we have obtained the BOEFT that describes hybrids. In
particular we have obtained the static potentials and the set of coupled
Schrödinger equations, solved them and produced all the hybrids multiplets,
see Fig. 12, for the case of the two first static energies $\Sigma_{u}^{-}$
and $\Pi_{u}$. Such static energies are degenerated at short distance where
the cylindrical symmetry gets subdue to a $O(3)$ symmetry and are then
labelled by the quantum number of a gluonic operator $1^{+-}$ called a
gluelump. The hybrid static energies are described by a repulsive octet
potential plus the gluelump mass in the short distance limit and the $O(3)$
symmetry is broken at order $r^{2}$ of the multipole expansion. In the long
distance regime the static energies display a linear r behaviour. The gluelump
correlator can be calculated on the lattice to determine the gluelump mass. It
is depending on the scheme used (the scheme dependence cancels against the
analogous dependence in the quark mass and in the octet static potential) but
it is of the order of 800 MeV. The hybrid multiplets $H_{i}$ are constructed
from the solution of the Schrödinger equations in correspondence of their
$J^{PC}$ quantum numbers. The coupling between the different Schrödinger
equations is induced by a nonadiabatic coupling, known in the Born Oppenheimer
description of diatomic molecules, induced by the noncommutation between the
kinetic term and a projector of the cylindrical symmetry in the BOEFT
lagrangian.

The degeneracy of the static energies at small distance induces a phenomenon
called $\Lambda$ doubling, removing the degenerations between multiplets of
opposite parity. This phenomenon is known in molecular physics but with
smaller size. This and the structure of the multiplets differ from what is
obtained in models for the hybrids, cf. Berwein:2015vca . We used lattice
input on the hybrid static energies and on the gluelump mass.

Figure 12: Mass spectrum of neutral exotic charmonium states obtained by
solving the BOEFT coupled Schrödinger equations. The neutral experimental
states that have matching quantum numbers are plotted in solid blue lines. In
the figure $T$ stay for the total angular momentum. $H_{1}^{\prime}$ is the
first $H_{1}$ radial excitation of $H_{1}$. The multiplets have been plotted
with error bands corresponding to a gluelump mass uncertainty of 0.15 GeV.
Figure taken from Brambilla:2019esw .

In Brambilla:2019jfi ; Brambilla:2018pyn  we obtained the spin dependent
potentials at order $1/m$ and $1/m^{2}$ in the quark mass expansion and we
could calculate all the hybrids spin multiplets. Notice that on one hand it is
seldom to find the spin interaction considered in models, on the other hand it
would be different. In fact the $1/m$ contributions couple the angular
momentum of the gluonic excitation with the total spin of the heavy-quark-
antiquark pair. These operators are characteristic of the hybrid states and
are absent for standard quarkonia. Among the $1/m^{2}$ operators besides the
standard spin-orbit, total spin squared, and tensor spin operators
respectively, which appear for standard quarkonia, three novel operators
appear. So interestingly, differently from the quarkonium case, the hybrid
potential gets a first contribution already at order $\Lambda^{2}_{\rm
QCD}/m$. Hence, spin splittings are remarkably less suppressed in heavy
quarkonium hybrids than in heavy quarkonia: this will have a notable impact on
the phenomenology of exotics. We extracted the nonperturbative low energy
correlators appearing in the factorization fixing them on lattice data on the
masses of charmonium hybrids and we could then predict all the bottomonium
hybrids spin multiplets, that are more difficult to evaluate on the lattice.
In this same framework it is also possible to calculate hybrids decays and
quarkonium/hybrids mixing Oncala:2017hop ; Abi2021. The BOEFT may be used to
describe also tetraquarks, double heavy baryons and pentaquarks
Brambilla:2017uyf ; Soto:2020xpm . In the case of tetraquarks, a necessary
input of the theory is the calculation of the lattice generalized Wilson loops
with appropriate symmetry and light quark operators, so that besides the
quantum number $\kappa$ also the isospin quantum numbers $I=0,1$ have to be
considered. The BOEFT approach reconciles the different pictures of exotics
based on tetraquarks, molecules, hadroquarkonium …In fact in the plot of a
static energy as a function of $r$ for a state with $Q\bar{Q}q\bar{q}$ or
$Q\bar{Q}g$ we will have different regions: for short distance a
hadroquarkonium picture would emerge, then a tetraquark (or hybrid) one and
when passing the the heavy-light mesons line, avoided cross level effects
should have to be taken into account and a molecular picture would emerge.
However QCD would dictate, through the lattice correlators and the BOEFT
characteristics and power counting, which structure would dominate and in
which precise way. In addition production and suppression in medium may be
described in the same approach Brambilla:2021abf ; Brambilla:2020qwo .
Production of hybrid states have been studied in Chiladze:1998ti ;
Petrov:2005tp .

##  8 pNRQCD at finite T, open quantum system and quarkonium in medium

Quarkonium is a special probe also for deconfinement, besides confinement. A
prediction of QCD is that at a certain value of temperature (or energy
density) hadronic matter undergoes a transition to a deconfined state of
quarks and gluons called the quark-gluon plasma (QGP). Lattice QCD have shown
that such transition takes place at a critical temperature around 150 MeV
Karsch:2001vs ; Borsanyi:2020fev . Experimentally heavy ion collisions make it
possible to study strongly interacting matter under extreme conditions in the
laboratory, recreating the QGP that should have primordially existed
microseconds after the Big Bang. These experiments have therefore a great
importance also for cosmology and allows to investigate the nuclear matter
phase diagram, i.e. how nuclear matter change by varying the temperature and
the chemical potential, see Brambilla:2014jmp ; Andronic:2015wma . In heavy
ions experiments at the LHC at CERN and at the RHIC at BNL the produced matter
is characterized by small net baryon densities and high temperature.
Karsch:2001vs . Heavy quarks are good probes of the QGP. They are produced at
the beginning of the collision and remain up to the end. The heavy quark mass
$m$ introduces a large scale, whose contribution may be factorized and
computed in perturbation theory. Low-energy scales are sensitive to the
temperature $T$ and even if nonperturbative they may be accessible via lattice
calculations (for what concerns equilibrium physics). As we discussed,
quarkonia are special hard probes as they are multi-scale systems. In the hot
QCD medium also the thermal scales of the Quark Gluon Plasma (QGP) are
emerging: the scale related to the temperature $(\pi)T$, the Debye mass
$m_{D}\sim gT$ related to the (chromo) electric screening and the scale
$g^{2}T$ related to the (chromo)magnetic screening. In a weakly coupled plasma
the scales are separated and hierarchically ordered, in a strongly coupled
plasma $m_{D}\sim T$. To address QCD at finite T calculations, EFTs have been
developed too to resum contributions related to the thermal scales and to
address IR sensitivities. In real time the Hard Thermal Loop EFT (HTL)
Braaten:1989mz ; Braaten:1991gm  is taking care of integrating out the
temperature scale.

Heavy quarkonium dissociation has been proposed long time ago as a clear probe
of the quark-gluon plasma formation in colliders through the measurement of
the dilepton decay-rate signal- In Matsui:1986dk  this was related to the
screening of the quark-antiquark interaction due to Debye mass and it was
suggested that this would have manifested in an exponential screening term
$\exp(-m_{D}r)$ in the static potential. One of the key quantities measured in
experiments is the nuclear modification factor
$R_{AA}=Y(PbPb)/(N_{coll}Y(pp))$ where $Y(PbPb)$ and $Y(pp)$ are the
quarkonium yield in PbPb and in pp collisions respectively. $R_{AA}$ is a
measure for the difference in particle production in pp and nucleus-nucleus
collisions. Since higher excited quarkonium states have larger radius then the
expectation was that, as the temperature increases, quarkonium would
dissociate subsequently from the higher to the lower states giving origin to
sequential melting. In order to study quarkonium properties in a thermal bath
at a temperature $T$, the quantity to be determined is the quarkonium
potential $V$ which dictates, through the Schrödinger equation the real-time
evolution of the wave function of a $Q\bar{Q}$ pair in the medium. This has
been investigated for years using many phenomenological assumptions, spanning
from the internal energy to the free energy, either the average free energy or
the singlet one which is gauge dependent. pNRQCD has given us the possibility
to define in QCD what is this potential: it is the matching coefficient of the
EFT that results from the integration of all the scales above the scale of the
binding energy. In a series of papers Brambilla:2008cx ; Brambilla:2010vq ;
Brambilla:2013dpa ; Brambilla:2011sg ; Escobedo:2008sy ; Brambilla:2020esl  a
pNRQCD at finite $T$ description has been constructed. One has to proceed
integrating out all scales up to the binding energy and if the temperature is
higher than the binding energy then also the temperature has to be integrated
out using HTL EFT. When $T$ is bigger than the energy, the potential depends
on the temperature, otherwise not. Thermal effects appear in any case in the
nonpotential contributions to the energy levels. We assumed that the bound
state exist for $T\ll m$ and $1/r\geq m_{D}$, we worked in the weak coupling
limit and we consider all possible scales hierarchies Brambilla:2008cx . We
found that the thermal part of the potential has a real part (roughly
described by the free energy) and an imaginary part. The imaginary part comes
from two effects: the Landau damping Laine:2006ns ; Escobedo:2008sy ;
Brambilla:2008cx , an effect existing also in QED, and the singlet to octet
transition, existing only in QCD Brambilla:2008cx . Which one dominates
depends on the ratio between $m_{D}$ and $E$. In the EFT one could show that
the imaginary part of the potential related to the Landau damping comes from
inelastic parton scattering Brambilla:2013dpa  and the singlet to octet
transitition from gluon dissociation Brambilla:2011sg . The existence of the
imaginary part, first realized in Laine:2006ns  changed our paradigm for
quarkonium suppression as the state dissociates for this reason well before
that the conventional screening becomes active Laine:2006ns ; Escobedo:2008sy
. The pattern of thermal corrections is pretty interesting Brambilla:2008cx :
when $T<E$ thermal corrections are only in the energy; for $T>1/r,1/r>m_{D}$
or $1/r>T>E$ there is no exponential screening and $T$ dependent power like
corrections appear; if $T>1/r,1/r\sim m_{D}$ we have exponential screening but
the imaginary part of the static potential is already bigger than the real one
and dissociation already happened. Once the potential has been calculated, the
EFT gives the systematic framework to obtain the thermal energies: in
Brambilla:2010vq  it was performed the first QCD calculation of the thermal
contributions to the $\Upsilon(1S)$ mass and width at order $m\alpha_{s}^{5}$
at LHC below the dissociation temperature of about 500 MeV. This calculation
is very important because it gives the parametric $T$ dependence of this
observables. The width goes linear in $T$ in the dominant term and this has
been confirmed by lattice calculations of the spectrum Aarts:2011sm . These
findings in the EFT in perturbation theory have inspired many subsequent
nonperturbative calculations of the static potential at finite $T$. In
particular the Wilson loop at finite $T$ has been calculated on the lattice
Rothkopf:2011db ; Rothkopf:2019ipj  finding hints of a large imaginary parts.
These calculations are pretty challenging and refining of the extraction
methods are currently in elaboration Bala:2021fkm .

Free energies defined by Polyakov loop correlators have been always very
prominent in QCD at finite $T,$ see e.g. the reviews Bazavov:2020teh ;
Ghiglieri:2020dpq . The Polyakov loop correlators of a single heavy quark and
of a quark-antiquark pair have been calculated both in perturbation theory
using pNRQCD to resum scales contributions in Berwein:2015ayt ;
Berwein:2013xza; Berwein:2017thy ; Brambilla:2010xn  and on the lattice to
obtain these quantities fully nonperturbatively in Bazavov:2016uvm ;
Bazavov:2018wmo . In particular: the Polyakov loop has been computed up to
order $g^{6}$, the (subtracted) $Q\bar{Q}$ free energy has been computed at
short distances up to corrections of order $g^{7}(rT)^{4}$, $g^{8}$, the
(subtracted) $Q\bar{Q}$ free energy has been computed at screening distances
up to corrections of order $g^{8}$; the singlet free energy has been computed
at short distances up to corrections of order $g^{4}(rT)^{5}$, $g^{6}$; the
singlet free energy has been computed at screening distances up to corrections
of order $g^{5}$ Berwein:2015ayt ; Berwein:2017thy ; Brambilla:2010xn . From
the lattice simulations and from comparison to the perturbative results some
interesting information can be obtained Bazavov:2016uvm ; Bazavov:2018wmo :
lattice calculations are consistent with weak-coupling expectations in the
regime of application of the weakly coupled resummed perturbation theory which
confirms the predictive power of the EFT; the crossover temperature to the
quark-gluon plasma is $153+6.5-5$ MeV as extracted from the entropy of the
Polyakov loop; the screening sets in at $rT\sim 0.3-0.4$ (observable
dependent), consistent with a screening length of about $1/m_{D}$; asymptotic
screening masses are about $2m_{D}$ (observable dependent); the first
determination of the color octet $Q\bar{Q}$ free energy has been obtained
Bazavov:2018wmo . The free energies turn out not to be the objects to be used
as a potential in the Schrödinger equation even if the singlet free energy may
provide a good approximation of the real part of the static potential.

These are all results in thermal equilibrium. However, the evolution of
quarkonium in the QGP is an out of equilibrium process in which many effects
enter: the hydrodynamical evolution of the plasma and the production,
dissociation and regeneration of quarkonium in the medium, to quote the most
prominent ones. It is necessary therefore to introduce an appropriate
framework to describe the real time nonequilibrium evolution of quarkonium in
the QGP medium. This has been realized in Brambilla:2016wgg ;
Brambilla:2017zei ; Brambilla:2019tpt  where an open quantum system (OQS)
framework rooted in pNRQCD at finite $T$ has been developed that is fully
quantum, conserves the number of heavy quarks and considers both color singlet
and color octet quarkonium degrees of freedom. For a review of open quantum
system approach see Akamatsu:2020ypb ; Akamatsu:2014qsa ; Sharma:2019xum ;
Yao:2021lus  and references therein.

We consider the density matrix associated to our system. We distinguish
between the environment (QGP), characterized by the scale $T$, and the
subsystem of system (quarkonium) characterized by the scale $E$. We identify
the inverse of E with the intrinsic time scale of the subsystem: $\tau_{S}\sim
1/E$ and the inverse of $\pi T$ with the correlation time of the environment:
$\tau_{E}\sim 1/(\pi T)$. If the medium is in thermal equilibrium, or locally
in thermal equilibrium, we may understand $T$ as a temperature, otherwise is
just a parameter. The medium can be strongly coupled. Then we trace the
density matrix over the environment and we are left with a color singlet and
color octet density matrix that can be written in terns of the pNRQCD fields,
working in the close time path formalism. The evolution of the system is
characterized by a relaxation time $\tau_{R}$ that is estimated by the inverse
of the color singlet self-energy diagram in pNRQCD at finite $T$. We select
quarkonia states with a small radius (Coulombic) for which $1/r\gg\pi
T,\Lambda_{QCD}$ and we consider $\pi T\gg E$. In this framework ,in
Brambilla:2017zei , a set of master equations governing the time evolution of
heavy quarkonium in a medium can be derived. The equations follow from
assuming the inverse Bohr radius of the quarkonium to be greater than the
energy scale of the medium, and model the quarkonium as evolving in the vacuum
up to a time $t=t_{0}$ at which point interactions with the medium begin. The
equations express the time evolution of the density matrices of the heavy
quark-antiquark color singlet, $\rho_{s}$, and octet states, $\rho_{o}$, in
terms of the color singlet and octet Hamiltonians, $h_{s}={\bf
p}^{2}/M-C_{F}\alpha_{s}/r+...$ and $h_{o}={\bf
p}^{2}/M+\alpha_{s}/(2N_{c}r)+...$, and interaction terms with the medium,
which, at order $r^{2}$ in the multipole expansion, are encoded in the self-
energy diagrams of the EFTs. These interactions account for the mass shift of
the heavy $Q\bar{Q}$ pair induced by the medium, its decay width induced by
the medium, the generation of $Q\bar{Q}$ color singlet states from $Q\bar{Q}$
color octet states interacting with the medium and the generation of
$Q\bar{Q}$ color octet states from $Q\bar{Q}$ (color singlet or octet) states
interacting with the medium (recombination terms). Both Landau damping and
singlet to octet effects are contained in the self energy. The leading order
interaction between a heavy $Q\bar{Q}$ field and the medium is encoded in
pNRQCD in a chromoelectric dipole interaction, which appears at order
$r/m^{0}$ in the EFT Lagrangian. The approach gives us master equations, in
general non Markovian, for the out of equilibrium evolution of the color
singlet and color octet matrix densities. The system is in non-equilibrium
because through interaction with the environment (quark gluon plasma) singlet
and octet quark-antiquark states continuously transform in each other although
the total number of heavy quarks is conserved.

Assuming that any energy scale in the medium is larger than the heavy
$Q\bar{Q}$ binding energy $E$, in particular that $\tau_{R}\gg\tau_{E}$, we
obtain a Markovian evolution while the chosen hierarchy of scales implies
$\tau_{s}\gg\tau_{E}$ qualifying the regime as quantum Brownian motion. In
this situation we can reduce the general master equation to a Linblad form. In
this case the properties of the QGP are encoded in two transport coefficients:
the heavy quark momentum diffusion coefficient, $\kappa$, and its dispersive
counterpart $\gamma$ which are given by time integrals of appropriate gauge
invariant correlators at finite $T$ given by the integral of gauge invariant
finite T QCD correlators of chromoelectric fields:

$\displaystyle\kappa$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{6N_{c}}\int^{\infty}_{0}\mathrm{d}s~{}\Big{\langle}\left\\{E^{a,i}(s,{\bf
0}),E^{a,i}(0,{\bf 0})\right\\}\Big{\rangle},$ (29) $\displaystyle\gamma$
$\displaystyle=-i$
$\displaystyle\frac{g^{2}}{6N_{c}}\int^{\infty}_{0}\mathrm{d}s~{}\Big{\langle}\left[E^{a,i}(s,{\bf
0}),E^{a,i}(0,{\bf 0})\right]\Big{\rangle}.$ (30)

They come from the pNRQCD self energies that in this regime could be
factorized between this contribution and the bound state dependent part. In
the case of a nonperturbative QGP, these objects are nonperturbative and
should be evaluated on the lattice at a given temperature in an extended
window of temperatures Brambilla:2020siz . In this way one could use lattice
to give input to a nonequilibrium calcultation. In a series of papers
Brambilla:2021wkt ; Brambilla:2020qwo  we have solved the Linblad equation
using the highly efficient quantum trajectory method and realistically
implementing a 3+1D dissipative hydrodynamics code Omar:2021kra , and we have
obtained the bottomonium ($\Upsilon(1S),\Upsilon(2S),\Upsilon(3S)$) nuclear
modification factors and the anisotropic flows, see Fig. 13 for a comparison
between LHC data and our results on $R_{AA}$. The computation does not rely on
any free parameter, as it depends only on the two transport coefficients that
have been evaluated independently in lattice QCD. Our final results, which
include late-time feed down of excited states, agree well with the available
data from LHC 5.02 TeV PbPb collisions. Notice that the bands in the Fig. 13
depend on the undetermination with which the nonperturbative transport
coefficients $\kappa$ and $\gamma$, that are properties of the QGP, are
presently known. More precise experimental data can help to narrow down their
range. In this way we can use quarkonia with small radii as a diagnostic and
investigation tool of the characteristics of the strongly coupled QGP. To
apply the same description to charmonium would require to modify the master
equations considering terms beyond leading order in the quarkonium density.
Further work in a similar approach at the level of the Boltzmann equation has
been done in Yao:2020kqy ; Yao:2018nmy .

Figure 13: The nuclear modification factor $R_{AA}$ of the $\Upsilon(1S)$,
$\Upsilon(2S)$, and $\Upsilon(3S)$ as a function of $N_{{part}}$ compared to
experimental measurements from the ALICE Acharya:2020kls , ATLAS ATLAS5TeV ,
and CMS Sirunyan:2018nsz  collaborations. The bands in the theoretical curves
indicate variation with respect to $\hat{\kappa}(T)$ (left) and $\hat{\gamma}$
(right). The central curves represent the central values of $\hat{\kappa}(T)$
and $\hat{\gamma}$, and the dashed and dot-dashed lines represent the lower
and upper values, respectively, of $\hat{\kappa}(T)\equiv\kappa/T^{3}$ and
$\hat{\gamma}\equiv\gamma/T^{3}$ , see text for the definition of these
parameters. Figure taken from Brambilla:2021wkt  .

##  9 Outlook

The great progress of the last few decades in the construction of
nonrelativistic effective quantum field theories and the progress in lattice
QCD calculation allows to treat heavy quark bound states systematically in
QCD. In this way quarkonium becomes a golden probe of strong interactions at
zero and finite temperature.

Moreover, quarkonium is the prototype of a nonrelativistic multiscale system.
Systems of such kind are ubiquitous in matter and in fields of physics, from
particle to nuclear physics, to condensed matter and to astro and cosmological
applications and play a key role in several open challenges at the frontier of
particle physics. Therefore, all the progress obtained in this framework holds
the promise to have an impact in a number of relevant contemporary problems.

For example. combining the techniques of nonrelativistic effective field
theories and the open quantum framework, we can address the nonequilibrium
evolution of heavy dark matter pairs in the early universe. On the other hand
the XYZ exotics bear similarities to atomic and molecular physics, and results
obtained in BOEFT can be used in those fields Brambilla:2017uyf ;
Brambilla:2017ffe  as well as maybe extended to the consideration of strongly
correlated system in condensed matter at finite temperature.

In summary, the physics of heavy quarks not only has been and is extremely
important for our advancement in nuclear and particle physics, but has the
promise to impact also many other nearby fields.

## References

  * (1) J. J. Aubert et al. [E598], Phys. Rev. Lett. 33, 1404-1406 (1974) doi:10.1103/PhysRevLett.33.1404
  * (2) J. E. Augustin et al. [SLAC-SP-017], Phys. Rev. Lett. 33, 1406-1408 (1974) doi:10.1103/PhysRevLett.33.1406
  * (3) H. Fritzsch, M. Gell-Mann and H. Leutwyler, Phys. Lett. B 47, 365-368 (1973) doi:10.1016/0370-2693(73)90625-4
  * (4) S. L. Glashow, J. Iliopoulos and L. Maiani, Phys. Rev. D 2, 1285-1292 (1970) doi:10.1103/PhysRevD.2.1285
  * (5) W. Fischler, Nucl. Phys. B 129 (1977), 157-174 doi:10.1016/0550-3213(77)90026-8
  * (6) K. G. Wilson, Phys. Rev. D 10, 2445-2459 (1974) doi:10.1103/PhysRevD.10.2445
  * (7) L. S. Brown and W. I. Weisberger, Phys. Rev. D 20, 3239 (1979) doi:10.1103/PhysRevD.20.3239
  * (8) M. Creutz, Phys. Rev. D 15, 1128 (1977) doi:10.1103/PhysRevD.15.1128
  * (9) M. Creutz, Phys. Rev. D 21, 2308-2315 (1980) doi:10.1103/PhysRevD.21.2308
  * (10) M. Creutz and K. J. M. Moriarty, Phys. Rev. D 26, 2166 (1982) doi:10.1103/PhysRevD.26.2166
  * (11) S. N. Gupta, S. F. Radford and W. W. Repko, Phys. Rev. D 26, 3305 (1982) doi:10.1103/PhysRevD.26.3305
  * (12) T. Appelquist, M. Dine and I. J. Muzinich, Phys. Lett. B 69, 231-236 (1977) doi:10.1016/0370-2693(77)90651-7
  * (13) N. Brambilla and A. Vairo, [arXiv:hep-ph/9904330 [hep-ph]].
  * (14) W. Lucha, F. F. Schoberl and D. Gromes, Phys. Rept. 200, 127-240 (1991) doi:10.1
  * (15) N. Brambilla, Nuovo Cim. A 105, 949-958 (1992) doi:10.1007/BF02730836
  * (16) E. E. Salpeter and H. A. Bethe, Phys. Rev. 84, 1232-1242 (1951) doi:10.1103/PhysRev.84.1232
  * (17) N. Brambilla, Nuovo Cim. A 105, 949-958 (1992) doi:10.1007/BF02730836
  * (18) H. J. Schnitzer, Phys. Rev. D 18, 3482 (1978) doi:10.1103/PhysRevD.18.3482
  * (19) G. Buchalla, T. K. Komatsubara, F. Muheim, L. Silvestrini, M. Artuso, D. M. Asner, P. Ball, E. Baracchini, G. Bell and M. Beneke, et al. Eur. Phys. J. C 57, 309-492 (2008) doi:10.1140/epjc/s10052-008-0716-1 [arXiv:0801.1833 [hep-ph]].
  * (20) M. Neubert, Int. J. Mod. Phys. A 11, 4173-4240 (1996) doi:10.1142/S0217751X96001966 [arXiv:hep-ph/9604412 [hep-ph]].
  * (21) P. Chang, K. F. Chen and W. S. Hou, Prog. Part. Nucl. Phys. 97, 261-311 (2017) doi:10.1016/j.ppnp.2017.07.001 [arXiv:1708.03793 [hep-ph]].
  * (22) T. Gershon and V. V. Gligorov, Rept. Prog. Phys. 80, no.4, 046201 (2017) doi:10.1088/1361-6633/aa5514 [arXiv:1607.06746 [hep-ex]].
  * (23) E. Kou et al. [Belle-II], PTEP 2019, no.12, 123C01 (2019) [erratum: PTEP 2020, no.2, 029201 (2020)] doi:10.1093/ptep/ptz106 [arXiv:1808.10567 [hep-ex]].
  * (24) A. J. Bevan et al. [BaBar and Belle], Eur. Phys. J. C 74, 3026 (2014) doi:10.1140/epjc/s10052-014-3026-9 [arXiv:1406.6311 [hep-ex]].
  * (25) E. J. Eichten and C. Quigg, Phys. Rev. D 49, 5845-5856 (1994) doi:10.1103/PhysRevD.49.5845 [arXiv:hep-ph/9402210 [hep-ph]].
  * (26) 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) doi:10.1140/epjc/s10052-010-1534-9 [arXiv:1010.5827 [hep-ph]].
  * (27) N. Brambilla, S. Eidelman, P. Foka, S. Gardner, A. S. Kronfeld, M. G. Alford, R. Alkofer, M. Butenschoen, T. D. Cohen and J. Erdmenger, et al. Eur. Phys. J. C 74, no.10, 2981 (2014) doi:10.1140/epjc/s10052-014-2981-5 [arXiv:1404.3723 [hep-ph]].
  * (28) N. Brambilla et al. [Quarkonium Working Group], doi:10.5170/CERN-2005-005 [arXiv:hep-ph/0412158 [hep-ph]].
  * (29) 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) doi:10.1140/epjc/s10052-015-3819-5 [arXiv:1506.03981 [nucl-ex]].
  * (30) E. Chapon, D. d’Enterria, B. Ducloue, M. G. Echevarria, P. B. Gossiaux, V. Kartvelishvili, T. Kasemets, J. P. Lansberg, R. McNulty and D. D. Price, et al. Prog. Part. Nucl. Phys. 122, 103906 (2022) doi:10.1016/j.ppnp.2021.103906 [arXiv:2012.14161 [hep-ph]].
  * (31) G. T. Bodwin, E. Braaten, E. Eichten, S. L. Olsen, T. K. Pedlar and J. Russ, [arXiv:1307.7425 [hep-ph]].
  * (32) H. S. Chung, doi:10.1142/9789811219313_0064
  * (33) C. Z. Yuan and S. L. Olsen, Nature Rev. Phys. 1, no.8, 480-494 (2019) doi:10.1038/s42254-019-0082-y [arXiv:2001.01164 [hep-ex]].
  * (34) S. Joosten and Z. E. Meziani, PoS QCDEV2017, 017 (2018) doi:10.22323/1.308.0017 [arXiv:1802.02616 [hep-ex]].
  * (35) G. Barucca et al. [PANDA], Eur. Phys. J. A 57, no.6, 184 (2021) doi:10.1140/epja/s10050-021-00475-y [arXiv:2101.11877 [hep-ex]].
  * (36) E. Eichten, K. Gottfried, T. Kinoshita, K. D. Lane and T. M. Yan, Phys. Rev. D 17, 3090 (1978) [erratum: Phys. Rev. D 21, 313 (1980)] doi:10.1103/PhysRevD.17.3090
  * (37) E. Eichten, K. Gottfried, T. Kinoshita, K. D. Lane and T. M. Yan, Phys. Rev. D 21, 203 (1980) doi:10.1103/PhysRevD.21.203
  * (38) R. Aaij et al. [LHCb], Sci. Bull. 65, no.23, 1983-1993 (2020) doi:10.1016/j.scib.2020.08.032 [arXiv:2006.16957 [hep-ex]].
  * (39) A. S. Kronfeld, Ann. Rev. Nucl. Part. Sci. 62, 265-284 (2012) doi:10.1146/annurev-nucl-102711-094942 [arXiv:1203.1204 [hep-lat]].
  * (40) H. J. Rothe, World Sci. Lect. Notes Phys. 43, 1-381 (1992) doi:10.1142/8229
  * (41) J. B. Kogut, Rev. Mod. Phys. 55, 775 (1983) doi:10.1103/RevModPhys.55.775
  * (42) G. S. Bali, Phys. Rept. 343, 1-136 (2001) doi:10.1016/S0370-1573(00)00079-X [arXiv:hep-ph/0001312 [hep-ph]].
  * (43) F. Karsch, Nucl. Phys. A 698, 199-208 (2002) doi:10.1016/S0375-9474(01)01365-3 [arXiv:hep-ph/0103314 [hep-ph]].
  * (44) M. Campostrini, K. J. M. Moriarty, J. Potvin and C. Rebbi, Phys. Lett. B 193, 78-84 (1987) doi:10.1016/0370-2693(87)90459-X
  * (45) M. Neubert, doi:10.1142/9789812773579_0004 [arXiv:hep-ph/0512222 [hep-ph]].
  * (46) H. Georgi, HUTP-91-A039.
  * (47) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Rev. Mod. Phys. 77, 1423 (2005) [arXiv:hep-ph/0410047 [hep-ph]].
  * (48) E. Eichten and F. Feinberg, Phys. Rev. D 23, 2724 (1981) doi:10.1103/PhysRevD.23.2724
  * (49) D. Gromes, Z. Phys. C 26, 401 (1984) doi:10.1007/BF01452566
  * (50) A. Barchielli, E. Montaldi and G. M. Prosperi, Nucl. Phys. B 296, 625 (1988) [erratum: Nucl. Phys. B 303, 752 (1988)] doi:10.1016/0550-3213(88)90036-3
  * (51) A. Barchielli, N. Brambilla and G. M. Prosperi, Nuovo Cim. A 103, 59 (1990) doi:10.1007/BF02902620
  * (52) N. Brambilla, P. Consoli and G. M. Prosperi, Phys. Rev. D 50, 5878-5892 (1994) doi:10.1103/PhysRevD.50.5878 [arXiv:hep-th/9401051 [hep-th]].
  * (53) G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51, 1125-1171 (1995) [erratum: Phys. Rev. D 55, 5853 (1997)] doi:10.1103/PhysRevD.55.5853 [arXiv:hep-ph/9407339 [hep-ph]].
  * (54) W. E. Caswell and G. P. Lepage, Phys. Lett. B 167, 437-442 (1986) doi:10.1016/0370-2693(86)91297-9
  * (55) G. P. Lepage, L. Magnea, C. Nakhleh, U. Magnea and K. Hornbostel, Phys. Rev. D 46, 4052-4067 (1992) doi:10.1103/PhysRevD.46.4052 [arXiv:hep-lat/9205007 [hep-lat]].
  * (56) G. C. Nayak, J. W. Qiu and G. F. Sterman, Phys. Rev. D 72, 114012 (2005) doi:10.1103/PhysRevD.72.114012 [arXiv:hep-ph/0509021 [hep-ph]].
  * (57) G. T. Bodwin, H. S. Chung, J. H. Ee, U. R. Kim and J. Lee, Phys. Rev. D 101, no.9, 096011 (2020) doi:10.1103/PhysRevD.101.096011 [arXiv:1910.05497 [hep-ph]].
  * (58) A. Pineda and J. Soto, Nucl. Phys. B Proc. Suppl. 64, 428-432 (1998) doi:10.1016/S0920-5632(97)01102-X [arXiv:hep-ph/9707481 [hep-ph]].
  * (59) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Phys. Rev. D 60, 091502 (1999) doi:10.1103/PhysRevD.60.091502 [arXiv:hep-ph/9903355 [hep-ph]].
  * (60) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Nucl. Phys. B 566, 275 (2000) doi:10.1016/S0550-3213(99)00693-8 [arXiv:hep-ph/9907240 [hep-ph]].
  * (61) https://einrichtungen.ph.tum.de/T30f/tumqcd/index.html
  * (62) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Phys. Rev. D 63, 014023 (2001) doi:10.1103/PhysRevD.63.014023 [arXiv:hep-ph/0002250 [hep-ph]].
  * (63) A. Pineda and A. Vairo, Phys. Rev. D 63, 054007 (2001) [erratum: Phys. Rev. D 64, 039902 (2001)] doi:10.1103/PhysRevD.64.039902 [arXiv:hep-ph/0009145 [hep-ph]].
  * (64) A. V. Manohar, Phys. Rev. D 56, 230-237 (1997) doi:10.1103/PhysRevD.56.230 [arXiv:hep-ph/9701294 [hep-ph]].
  * (65) N. Brambilla, D. Gromes and A. Vairo, Phys. Lett. B 576, 314-327 (2003) doi:10.1016/j.physletb.2003.09.100 [arXiv:hep-ph/0306107 [hep-ph]].
  * (66) J. Heinonen, R. J. Hill and M. P. Solon, Phys. Rev. D 86, 094020 (2012) doi:10.1103/PhysRevD.86.094020 [arXiv:1208.0601 [hep-ph]].
  * (67) M. Berwein, N. Brambilla, S. Hwang and A. Vairo, Phys. Rev. D 99, no.9, 094008 (2019) doi:10.1103/PhysRevD.99.094008 [arXiv:1811.05184 [hep-th]].
  * (68) A. Pineda and J. Soto, Phys. Lett. B 495, 323-328 (2000) doi:10.1016/S0370-2693(00)01261-2 [arXiv:hep-ph/0007197 [hep-ph]].
  * (69) N. Brambilla, X. Garcia Tormo, i, J. Soto and A. Vairo, Phys. Lett. B 647, 185-193 (2007) doi:10.1016/j.physletb.2007.02.015 [arXiv:hep-ph/0610143 [hep-ph]].
  * (70) N. Brambilla, A. Vairo, X. Garcia Tormo, i and J. Soto, Phys. Rev. D 80, 034016 (2009) doi:10.1103/PhysRevD.80.034016 [arXiv:0906.1390 [hep-ph]].
  * (71) A. Bazavov, N. Brambilla, X. Garcia Tormo, i, P. Petreczky, J. Soto and A. Vairo, Phys. Rev. D 86, 114031 (2012) doi:10.1103/PhysRevD.86.114031 [arXiv:1205.6155 [hep-ph]].
  * (72) C. Ayala, X. Lobregat and A. Pineda, JHEP 09, 016 (2020) doi:10.1007/JHEP09(2020)016 [arXiv:2005.12301 [hep-ph]].
  * (73) C. Peset, A. Pineda and J. Segovia, JHEP 09, 167 (2018) doi:10.1007/JHEP09(2018)167 [arXiv:1806.05197 [hep-ph]].
  * (74) Y. Kiyo and Y. Sumino, Nucl. Phys. B 889, 156-191 (2014) doi:10.1016/j.nuclphysb.2014.10.010 [arXiv:1408.5590 [hep-ph]].
  * (75) H. Takaura, T. Kaneko, Y. Kiyo and Y. Sumino, JHEP 04, 155 (2019) doi:10.1007/JHEP04(2019)155 [arXiv:1808.01643 [hep-ph]].
  * (76) A. Pineda, Prog. Part. Nucl. Phys. 67, 735-785 (2012) doi:10.1016/j.ppnp.2012.01.038 [arXiv:1111.0165 [hep-ph]].
  * (77) B. A. Kniehl, A. A. Penin, Y. Schroder, V. A. Smirnov and M. Steinhauser, Phys. Lett. B 607, 96-100 (2005) doi:10.1016/j.physletb.2004.12.024 [arXiv:hep-ph/0412083 [hep-ph]].
  * (78) N. Brambilla, A. Pineda, J. Soto and A. Vairo, Phys. Lett. B 470, 215 (1999) doi:10.1016/S0370-2693(99)01301-5 [arXiv:hep-ph/9910238 [hep-ph]].
  * (79) B. A. Kniehl, A. A. Penin, V. A. Smirnov and M. Steinhauser, Nucl. Phys. B 635, 357-383 (2002) doi:10.1016/S0550-3213(02)00403-0 [arXiv:hep-ph/0203166 [hep-ph]].
  * (80) K. J. Juge, J. Kuti and C. Morningstar, Phys. Rev. Lett. 90, 161601 (2003) doi:10.1103/PhysRevLett.90.161601 [arXiv:hep-lat/0207004 [hep-lat]].
  * (81) G. S. Bali and A. Pineda, Phys. Rev. D 69, 094001 (2004) doi:10.1103/PhysRevD.69.094001 [arXiv:hep-ph/0310130 [hep-ph]].
  * (82) G. S. Bali, K. Schilling and A. Wachter, Phys. Rev. D 56, 2566-2589 (1997) doi:10.1103/PhysRevD.56.2566 [arXiv:hep-lat/9703019 [hep-lat]].
  * (83) Y. Koma, M. Koma and H. Wittig, Phys. Rev. Lett. 97, 122003 (2006) doi:10.1103/PhysRevLett.97.122003 [arXiv:hep-lat/0607009 [hep-lat]].
  * (84) Y. Koma, M. Koma and H. Wittig, PoS LATTICE2007, 111 (2007) doi:10.22323/1.042.0111 [arXiv:0711.2322 [hep-lat]].
  * (85) L. Müller, O. Philipsen, C. Reisinger and M. Wagner, Phys. Rev. D 100, no.5, 054503 (2019) doi:10.1103/PhysRevD.100.054503 [arXiv:1907.01482 [hep-lat]].
  * (86) A. Bazavov et al. [TUMQCD], Phys. Rev. D 100, no.11, 114511 (2019) doi:10.1103/PhysRevD.100.114511 [arXiv:1907.11747 [hep-lat]].
  * (87) N. Brambilla, Y. Sumino and A. Vairo, Phys. Rev. D 65, 034001 (2002) doi:10.1103/PhysRevD.65.034001 [arXiv:hep-ph/0108084 [hep-ph]].
  * (88) N. Brambilla, Y. Jia and A. Vairo, Phys. Rev. D 73, 054005 (2006) doi:10.1103/PhysRevD.73.054005 [arXiv:hep-ph/0512369 [hep-ph]].
  * (89) N. Brambilla, P. Pietrulewicz and A. Vairo, Phys. Rev. D 85, 094005 (2012) doi:10.1103/PhysRevD.85.094005 [arXiv:1203.3020 [hep-ph]].
  * (90) A. Pineda and J. Segovia, Phys. Rev. D 87, no.7, 074024 (2013) doi:10.1103/PhysRevD.87.074024 [arXiv:1302.3528 [hep-ph]].
  * (91) J. Segovia, S. Steinbeißer and A. Vairo, Phys. Rev. D 99, no.7, 074011 (2019) doi:10.1103/PhysRevD.99.074011 [arXiv:1811.07590 [hep-ph]].
  * (92) N. Brambilla and A. Vairo, Phys. Rev. D 71, 034020 (2005) doi:10.1103/PhysRevD.71.034020 [arXiv:hep-ph/0411156 [hep-ph]].
  * (93) N. Brambilla, D. Eiras, A. Pineda, J. Soto and A. Vairo, Phys. Rev. D 67, 034018 (2003) doi:10.1103/PhysRevD.67.034018 [arXiv:hep-ph/0208019 [hep-ph]].
  * (94) N. Brambilla and A. Vairo, Phys. Rev. D 62, 094019 (2000) doi:10.1103/PhysRevD.62.094019 [arXiv:hep-ph/0002075 [hep-ph]].
  * (95) C. Peset, A. Pineda and M. Stahlhofen, JHEP 05, 017 (2016) doi:10.1007/JHEP05(2016)017 [arXiv:1511.08210 [hep-ph]].
  * (96) N. Brambilla, A. Vairo and T. Rosch, Phys. Rev. D 72, 034021 (2005) doi:10.1103/PhysRevD.72.034021 [arXiv:hep-ph/0506065 [hep-ph]].
  * (97) N. Brambilla, J. Ghiglieri and A. Vairo, Phys. Rev. D 81, 054031 (2010) doi:10.1103/PhysRevD.81.054031 [arXiv:0911.3541 [hep-ph]].
  * (98) T. T. Takahashi, H. Matsufuru, Y. Nemoto and H. Suganuma, Phys. Rev. Lett. 86, 18-21 (2001) doi:10.1103/PhysRevLett.86.18 [arXiv:hep-lat/0006005 [hep-lat]].
  * (99) T. T. Takahashi, H. Suganuma, Y. Nemoto and H. Matsufuru, Phys. Rev. D 65, 114509 (2002) doi:10.1103/PhysRevD.65.114509 [arXiv:hep-lat/0204011 [hep-lat]].
  * (100) G. S. Bali, K. Schilling and C. Schlichter, Phys. Rev. D 51, 5165-5198 (1995) doi:10.1103/PhysRevD.51.5165 [arXiv:hep-lat/9409005 [hep-lat]].
  * (101) P. Bicudo, N. Cardoso and M. Cardoso, Phys. Rev. D 98, no.11, 114507 (2018) doi:10.1103/PhysRevD.98.114507 [arXiv:1808.08815 [hep-lat]].
  * (102) R. Yanagihara and M. Kitazawa, PTEP 2019, no.9, 093B02 (2019) [erratum: PTEP 2020, no.7, 079201 (2020)] doi:10.1093/ptep/ptz093 [arXiv:1905.10056 [hep-ph]].
  * (103) M. Baker, V. Chelnokov, L. Cosmai, F. Cuteri and A. Papa, [arXiv:2111.06618 [hep-lat]].
  * (104) J. Greensite, Prog. Part. Nucl. Phys. 51, 1 (2003) doi:10.1016/S0146-6410(03)90012-3 [arXiv:hep-lat/0301023 [hep-lat]].
  * (105) M. Baker, J. S. Ball and F. Zachariasen, Phys. Rept. 209, 73-127 (1991) doi:10.1016/0370-1573(91)90123-4
  * (106) N. Isgur and J. E. Paton, Phys. Rev. D 31, 2910 (1985) doi:10.1103/PhysRevD.31.2910
  * (107) H. G. Dosch and Y. A. Simonov, Phys. Lett. B 205, 339-344 (1988) doi:10.1016/0370-2693(88)91675-9
  * (108) M. Baker, J. S. Ball, N. Brambilla, G. M. Prosperi and F. Zachariasen, Phys. Rev. D 54, 2829-2844 (1996) [erratum: Phys. Rev. D 56, 2475 (1997)] doi:10.1103/PhysRevD.56.2475 [arXiv:hep-ph/9602419 [hep-ph]].
  * (109) M. Baker, N. Brambilla, H. G. Dosch and A. Vairo, Phys. Rev. D 58, 034010 (1998) doi:10.1103/PhysRevD.58.034010 [arXiv:hep-ph/9802273 [hep-ph]].
  * (110) N. Brambilla, M. Groher, H. E. Martinez and A. Vairo, Phys. Rev. D 90, no.11, 114032 (2014) doi:10.1103/PhysRevD.90.114032 [arXiv:1407.7761 [hep-ph]].
  * (111) G. Perez-Nadal and J. Soto, Phys. Rev. D 79, 114002 (2009) doi:10.1103/PhysRevD.79.114002 [arXiv:0811.2762 [hep-ph]].
  * (112) N. Brambilla and A. Vairo, Phys. Rev. D 55, 3974-3986 (1997) doi:10.1103/PhysRevD.55.3974 [arXiv:hep-ph/9606344 [hep-ph]].
  * (113) K. Nawa, H. Suganuma and T. Kojo, Phys. Rev. D 75, 086003 (2007) doi:10.1103/PhysRevD.75.086003 [arXiv:hep-th/0612187 [hep-th]].
  * (114) K. Amemiya and H. Suganuma, Phys. Rev. D 60, 114509 (1999) doi:10.1103/PhysRevD.60.114509 [arXiv:hep-lat/9811035 [hep-lat]].
  * (115) S. Sasaki, H. Suganuma and H. Toki, Prog. Theor. Phys. 94, 373-384 (1995) doi:10.1143/PTP.94.373
  * (116) J. Soto and J. Tarrús Castellà, Phys. Rev. D 104, 074027 (2021) doi:10.1103/PhysRevD.104.074027 [arXiv:2108.00496 [hep-ph]].
  * (117) N. Brambilla, V. Leino, O. Philipsen, C. Reisinger, A. Vairo and M. Wagner, [arXiv:2106.01794 [hep-lat]].
  * (118) A. Vairo, Mod. Phys. Lett. A 31, no.34, 1630039 (2016) doi:10.1142/S0217732316300391
  * (119) N. Brambilla, H. S. Chung and A. Vairo, JHEP 09, 032 (2021) [arXiv:2106.09417 [hep-ph]].
  * (120) N. Brambilla, H. S. Chung and A. Vairo, Phys. Rev. Lett. 126, no.8, 082003 (2021) doi:10.1103/PhysRevLett.126.082003 [arXiv:2007.07613 [hep-ph]].
  * (121) N. Brambilla, H. S. Chung, D. Müller and A. Vairo, JHEP 04, 095 (2020) doi:10.1007/JHEP04(2020)095 [arXiv:2002.07462 [hep-ph]]. Copy to ClipboardDownload
  * (122) N. Brambilla, [arXiv:2111.10788 [hep-lat]].
  * (123) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan, Phys. Rept. 873, 1-154 (2020) doi:10.1016/j.physrep.2020.05.001 [arXiv:1907.07583 [hep-ex]].
  * (124) M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98 (2018) no.3, 030001.
  * (125) S. K. Choi et al. [Belle], Phys. Rev. Lett. 91 (2003), 262001.
  * (126) A. Ali, J. S. Lange and S. Stone, Prog. Part. Nucl. Phys. 97 (2017), 123-198.
  * (127) F. K. Guo, C. Hanhart, U. G. Meißner, Q. Wang, Q. Zhao and B. S. Zou, Rev. Mod. Phys. 90 (2018) no.1, 015004.
  * (128) M. T. AlFiky, F. Gabbiani and A. A. Petrov, Phys. Lett. B 640 (2006), 238-245 doi:10.1016/j.physletb.2006.07.069 [arXiv:hep-ph/0506141 [hep-ph]].
  * (129) E. Braaten and M. Lu, Phys. Rev. D 76 (2007), 094028.
  * (130) E. Braaten and M. Kusunoki, Phys. Rev. D 69 (2004), 074005.
  * (131) S. Fleming and T. Mehen, Phys. Rev. D 85 (2012), 014016.
  * (132) E. Braaten, C. Langmack and D. H. Smith, Phys. Rev. D 90 (2014) no.1, 014044.
  * (133) E. Braaten, C. Langmack and D. H. Smith, Phys. Rev. Lett. 112 (2014), 222001.
  * (134) G. Chiladze, A. F. Falk and A. A. Petrov, Phys. Rev. D 58 (1998), 034013 doi:10.1103/PhysRevD.58.034013 [arXiv:hep-ph/9804248 [hep-ph]].
  * (135) A. A. Petrov, J. Phys. Conf. Ser. 9 (2005), 83-86 doi:10.1088/1742-6596/9/1/013
  * (136) M. Berwein, N. Brambilla, J. Tarrús Castellà and A. Vairo, Phys. Rev. D 92 (2015) no.11, 114019.
  * (137) N. Brambilla, G. Krein, J. Tarrús Castellà and A. Vairo, Phys. Rev. D 97 (2018) no.1, 016016.
  * (138) N. Brambilla, W. K. Lai, J. Segovia and J. Tarrús Castellà, Phys. Rev. D 101 (2020) no.5, 054040.
  * (139) N. Brambilla, W. K. Lai, J. Segovia, J. Tarrús Castellà and A. Vairo, Phys. Rev. D 99 (2019) no.1, 014017.
  * (140) J. Soto and J. Tarrús Castellà, Phys. Rev. D 102 (2020) no.1, 014012.
  * (141) R. Oncala and J. Soto, Phys. Rev. D 96 (2017) no.1, 014004.
  * (142) S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti and K. K. Szabo, Phys. Rev. Lett. 125, no.5, 052001 (2020) doi:10.1103/PhysRevLett.125.052001 [arXiv:2002.02821 [hep-lat]].
  * (143) T. Matsui and H. Satz, Phys. Lett. B 178, 416-422 (1986) doi:10.1016/0370-2693(86)91404-8
  * (144) E. Braaten and R. D. Pisarski, Nucl. Phys. B 337, 569-634 (1990) doi:10.1016/0550-3213(90)90508-B
  * (145) E. Braaten and R. D. Pisarski, Phys. Rev. D 45, no.6, R1827 (1992) doi:10.1103/PhysRevD.45.R1827
  * (146) M. Laine, O. Philipsen, P. Romatschke and M. Tassler, JHEP 03, 054 (2007) doi:10.1088/1126-6708/2007/03/054 [arXiv:hep-ph/0611300 [hep-ph]].
  * (147) N. Brambilla, PoS HardProbes2020 (2021), 017.
  * (148) N. Brambilla, J. Ghiglieri, A. Vairo and P. Petreczky, Phys. Rev. D 78, 014017 (2008) doi:10.1103/PhysRevD.78.014017 [arXiv:0804.0993 [hep-ph]].
  * (149) N. Brambilla, M. A. Escobedo, J. Ghiglieri, J. Soto and A. Vairo, JHEP 09, 038 (2010) doi:10.1007/JHEP09(2010)038 [arXiv:1007.4156 [hep-ph]].
  * (150) N. Brambilla, M. A. Escobedo, J. Ghiglieri and A. Vairo, JHEP 12, 116 (2011) doi:10.1007/JHEP12(2011)116 [arXiv:1109.5826 [hep-ph]].
  * (151) N. Brambilla, M. A. Escobedo, J. Ghiglieri and A. Vairo, JHEP 05, 130 (2013) doi:10.1007/JHEP05(2013)130 [arXiv:1303.6097 [hep-ph]].
  * (152) M. A. Escobedo and J. Soto, Phys. Rev. A 78, 032520 (2008) doi:10.1103/PhysRevA.78.032520 [arXiv:0804.0691 [hep-ph]].
  * (153) G. Aarts, C. Allton, S. Kim, M. P. Lombardo, M. B. Oktay, S. M. Ryan, D. K. Sinclair and J. I. Skullerud, JHEP 11, 103 (2011) doi:10.1007/JHEP11(2011)103 [arXiv:1109.4496 [hep-lat]].
  * (154) A. Rothkopf, T. Hatsuda and S. Sasaki, Phys. Rev. Lett. 108, 162001 (2012) doi:10.1103/PhysRevLett.108.162001 [arXiv:1108.1579 [hep-lat]].
  * (155) A. Rothkopf, Phys. Rept. 858, 1-117 (2020) doi:10.1016/j.physrep.2020.02.006 [arXiv:1912.02253 [hep-ph]].
  * (156) D. Bala, O. Kaczmarek, R. Larsen, S. Mukherjee, G. Parkar, P. Petreczky, A. Rothkopf and J. H. Weber, [arXiv:2110.11659 [hep-lat]].
  * (157) J. Ghiglieri, A. Kurkela, M. Strickland and A. Vuorinen, Phys. Rept. 880 (2020), 1-73 [arXiv:2002.10188 [hep-ph]].
  * (158) M. Berwein, N. Brambilla, P. Petreczky and A. Vairo, Phys. Rev. D 96 (2017) no.1, 014025 [arXiv:1704.07266 [hep-ph]].
  * (159) M. Berwein, N. Brambilla, P. Petreczky and A. Vairo, Phys. Rev. D 93 (2016) no.3, 034010 [arXiv:1512.08443 [hep-ph]].
  * (160) N. Brambilla, J. Ghiglieri, P. Petreczky and A. Vairo, Phys. Rev. D 82 (2010), 074019 [arXiv:1007.5172 [hep-ph]].
  * (161) A. Bazavov et al. [TUMQCD], Phys. Rev. D 98 (2018) no.5, 054511 [arXiv:1804.10600 [hep-lat]].
  * (162) A. Bazavov, N. Brambilla, H. T. Ding, P. Petreczky, H. P. Schadler, A. Vairo and J. H. Weber, Phys. Rev. D 93 (2016) no.11, 114502 [arXiv:1603.06637 [hep-lat]].
  * (163) A. Bazavov and J. H. Weber, Prog. Part. Nucl. Phys. 116, 103823 (2021) doi:10.1016/j.ppnp.2020.103823 [arXiv:2010.01873 [hep-lat]].
  * (164) N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, Phys. Rev. D 97 (2018) no.7, 074009 [arXiv:1711.04515 [hep-ph]].
  * (165) N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, Phys. Rev. D 96 (2017) no.3, 034021 [arXiv:1612.07248 [hep-ph]].
  * (166) N. Brambilla, M. A. Escobedo, A. Vairo and P. Vander Griend, Phys. Rev. D 100 (2019) no.5, 054025 [arXiv:1903.08063 [hep-ph]].
  * (167) N. Brambilla, M. Á. Escobedo, M. Strickland, A. Vairo, P. Vander Griend and J. H. Weber, Phys. Rev. D 104, no.9, 094049 (2021) [arXiv:2107.06222 [hep-ph]].
  * (168) N. Brambilla, M. Á. Escobedo, M. Strickland, A. Vairo, P. Vander Griend and J. H. Weber, JHEP 05, 136 (2021) [arXiv:2012.01240 [hep-ph]].
  * (169) H. B. Omar, M. Á. Escobedo, A. Islam, M. Strickland, S. Thapa, P. Vander Griend and J. H. Weber, Comput. Phys. Commun. 273, 108266 (2022) doi:10.1016/j.cpc.2021.108266 [arXiv:2107.06147 [physics.comp-ph]].
  * (170) M. Á. Escobedo, Phys. Rev. D 103, no.3, 034010 (2021) doi:10.1103/PhysRevD.103.034010 [arXiv:2010.10424 [hep-ph]].
  * (171) Y. Akamatsu, Prog. Part. Nucl. Phys. 123, 103932 (2022) doi:10.1016/j.ppnp.2021.103932 [arXiv:2009.10559 [nucl-th]].
  * (172) Y. Akamatsu, Phys. Rev. D 91, no.5, 056002 (2015) doi:10.1103/PhysRevD.91.056002 [arXiv:1403.5783 [hep-ph]].
  * (173) X. Yao, W. Ke, Y. Xu, S. A. Bass, T. Mehen and B. Müller, [arXiv:2002.04079 [hep-ph]]; X. Yao and T. Mehen, [arXiv:2009.02408 [hep-ph]].
  * (174) X. Yao and T. Mehen, Phys. Rev. D 99 (2019) no.9, 096028 [arXiv:1811.07027 [hep-ph]].
  * (175) X. Yao, Int. J. Mod. Phys. A 36, no.20, 2130010 (2021) doi:10.1142/S0217751X21300106 [arXiv:2102.01736 [hep-ph]].
  * (176) R. Sharma and A. Tiwari, Phys. Rev. D 101 (2020) no.7, 074004 [arXiv:1912.07036 [hep-ph]].
  * (177) N. Brambilla, V. Leino, P. Petreczky and A. Vairo, Phys. Rev. D 102, no.7, 074503 (2020) doi:10.1103/PhysRevD.102.074503 [arXiv:2007.10078 [hep-lat]].
  * (178) S. Acharya et al. [ALICE], Phys. Lett. B 822, 136579 (2021) doi:10.1016/j.physletb.2021.136579 [arXiv:2011.05758 [nucl-ex]].
  * (179) A. M. Sirunyan et al. [CMS], Phys. Lett. B 790, 270-293 (2019) doi:10.1016/j.physletb.2019.01.006 [arXiv:1805.09215 [hep-ex]].
  * (180) Songkyo Lee (ATLAS Collaboration) Quark Matter 2020 https://indico.cern.ch/event/792436/contributions/3535775/
  * (181) N. Brambilla, V. Shtabovenko, J. Tarrús Castellà and A. Vairo, Phys. Rev. D 95, no.11, 116004 (2017) doi:10.1103/PhysRevD.95.116004 [arXiv:1704.03476 [hep-ph]].