*[enumerate,1]label=0) tee

# Navigating the Noise: A CBF Approach for Nonlinear Control with Integral
Constraints

Idris Seidu and Roberto Tron This work was supported by a grant from NASA.
Idris Seidu and Roberto Tron are with the Department of Mechanical
Engineering, Boston University, 110 Cummington Mall, MA 02215, United States
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Many physical phenomena involving mobile agents involve time-varying scalar
fields, e.g., quadrotors that emit noise. As a consequence, agents can
influence and can be influenced by various environmental factors such as
noise. This paper delves into the challenges of controlling such agents,
focusing on scenarios where we would like to prevent excessive accumulation of
some quantity over select regions and extended trajectories. We use quadrotors
that emit noise as a primary example, to regulate the trajectory of such
agents in the presence of obstacles and noise emitted by the aerial vehicles
themselves. First, we consider constraints that are defined over accumulated
quantities, i.e functionals of the entire trajectory, as opposed to those that
depend solely on the current state as in traditional Higher order Control
Barrier Functions (HOCBF). Second, we propose a method to extend constraints
from individual points to lines and sets by using efficient over-
approximations. The efficacy of the implemented strategies is verified using
simulations.Although we use quadrotors as an example, the same principles can
equally apply to other scenarios, such as light emission microscopy or vehicle
pollution dispersion. The technical contribution of this paper is twofold.

##  I INTRODUCTION

Quadrotors have become a transformative force in multiple sectors, including
transportation, surveillance, and aerial mapping, revolutionizing the way we
approach these fields. However, one critical challenge facing quadrotors, and
serving as a potential barrier to their broader acceptance in urban
environments, is the noise they emit.

In order to overcome these challenges, this work explores the application of
Control Barrier Functions (CBFs) for the navigation and control of agents
(quadrotor). CBFs have proven to be an effective tool in the design of
controllers for real-time collision avoidance with safety guarantees in
nonlinear systems [1].

The foundational work by [1] demonstrates how CBFs can be combined with
Control Lyapunov Functions (CLFs) using Quadratic Programs (QPs), providing a
robust framework for ensuring the forward-invariance of a set with conditions
that are linear in the inputs, hence facilitating their use as QPs. CBFs find
common application in safety-critical systems, particularly when used
alongside CLFs (e.g., adaptive cruise control scenarios [2, 1]). Beyond this,
CBFs have been applied to multi-robot systems, illustrating their versatility.
For instance, [3] and subsequent extensions in [8, 9] have showcased the
application of safety barrier certificates in ensuring collision-free
interactions among robots.

In another line of work for safety-critical control of quadrotors, CBFs have
been integrated with geometric control [11, 10]. These works, as the majority
of the literature, have considered CBFs that are defined on the value of a
single, current state of the system. In the type of applications considered in
this paper, however, we are interested in constraints that depend on the
_entire_ trajectory of the system (including past states).

Another significant body of work has focused on the application of CBFs over
time-varying sets. For instance, [7] introduces a time-varying CBF for
nonautonomous control-affine systems, addressing time-dependent safety
constraint problems and introducing a control law assisted by humans. Again,
however, the applicability of such approaches is limited in scenarios where
the definition of the set does not depend on past states of the system.

In parallel, [6, 13] proposed frameworks and planners based on Signal Temporal
Logic and time-varying CBFs. These contributions are crucial for the
development of computationally efficient control methods under temporal logic
tasks, especially in multi-robot systems. However, while these works discuss
spatio-temporal constraints over entire trajectories, they do not address the
type of constraints based on the accumulation of a scalar field at specific
locations.

With respect to the state of the art reviewed above, our paper provides the
following main contributions.

  * •

We introduce the use of constraints based on integral cost functionals that
track an accumulated cost $J(x,t)$ at specific locations $x$ over time $t$.
While we rely on existing HOCBF theory, its application to constraints of the
form $J(x,t)\leq J_{\text{limit}}$, where $J$ is a trajectory-dependent
integral, is novel.

  * •

We introduce a method that saves on computational effort by over-approximating
constraints across lines and dense sets, rather than individual points. By
setting a Control Barrier Function (CBF) on an upper bound
$\bar{J}>J_{\text{max}}$, where $J_{\max}$ is the theoretical maximum value of
$J$ over a set, we avoid the need to constantly track $J_{\text{max}}$,
significantly reducing computational demands.

Together, these contribution allow us to control mobile agents such as
quadrotors on paths that are not only void of collisions, but also enforce
limits on the accumulation of a scalar field such as noise. It’s important to
note that while this paper primarily focuses on quadrotors as an example of
mobile agents emitting noise, the underlying principles of our approach have
broader applications. These principles can be equally applied to scenarios
such as managing light emission in microscopy to minimize sample bleaching, or
controlling vehicle pollution dispersion in an urban environments.

##  II PRELIMINARIES

Consider the following nonlinear control-affine system:

$\dot{x}=f(x)+g(x)u$ (1)

where x $\in\mathbb{R}^{n}$ is the system state and u $\in\mathbb{R}^{m}$ is
the control input while $f(x)$ and $g(x)$ are the smooth vector fields. The
$\nabla$ operator represents the gradient for scalar-valued functions that are
differentiable with respect to $x$. Additionally, the time derivative of a
function $h\bigl{(}x(t)\bigr{)}$ with respect to time $t$ is denoted as
$\dot{h}(x)=\frac{d}{dt}h(x(t))$. The Lie derivative $L_{f}h=\nabla h^{T}f$
measure the changes in a function $h$ along the vector field(s) $f$.

######  Definition 1 (Time-Varying Control Barrier Functions)

A function $h(x,t)$ is a Time-Varying Control Barrier Function (TV-CBF) for
the system (1) if, for all $x\in\mathbb{R}^{n}$ and $t\geq 0$, it satisfies:

$\sup_{u\in\mathbb{R}^{m}}\left[L_{f}h(x,t)+L_{g}h(x,t)u+\frac{\partial
h}{\partial t}(x,t)\right]\geq-\alpha(h(x,t))$ (2)

A TV-CBF reduces to a regular CBF if $h$ is constant with respect to time.

The concept of Higher Order Control Barrier Functions (HOCBFs) extends
traditional barrier functions to accommodate systems where safety constraints
are dependent on higher-order derivatives of the system state. This extension
is necessary for systems where the control does not directly influence lower-
order derivatives. This concept is made more rigorous with the following.

######  Definition 2 (Relative degree)

The relative degree of a sufficiently differentiable function
$h:\mathbb{R}^{n}\to\mathbb{R}$ with respect to the dynamics (1) is defined as
the number of Lie derivatives needed until the control input $u$ explicitly
appears.

######  Definition 3 (Higher Order Barrier Functions)

For a function $h:\mathbb{R}^{n}\times[t_{0},\infty)\rightarrow\mathbb{R}$
that is differentiable up to order $m$, we construct a sequence of functionals
$\Phi_{0},\Phi_{1},\ldots,\Phi_{m}$, where
$\Phi_{i}:\mathbb{R}^{n}\times[t_{0},\infty)\rightarrow\mathbb{R}$ for
$i\in\\{0,\ldots,m\\}$, specified by:

$\displaystyle\Phi_{0}(x,t)$ $\displaystyle:=h(x,t),$ (3)
$\displaystyle\Phi_{i}(x,t)$
$\displaystyle:={\dot{\Phi}}_{i-1}(x,t)+\alpha_{i}(\Phi_{i-1}(x,t)),\quad\text{for
}i=1,\ldots,m,$ (4)

with each $\alpha_{i}$ being a class $\mathcal{K}$ function.

Correspondingly, we associate a collection of sets $B_{i}(t)$ for
$i\in\\{1,\ldots,m\\}$, defined as:

$B_{i}(t):=\left\\{x\in\mathbb{R}^{n}:\Phi_{i-1}(x,t)\geq 0\right\\}.$ (5)

######  Definition 4 (Higher Order Control Barrier Functions)

Consider the sets $B_{1}(t),\ldots,B_{m}(t)$ established by the preceding
definition, along with their associated functions $\Phi_{0},\ldots,\Phi_{m}$.
A function $h:\mathbb{R}^{n}\times[t_{0},\infty)\rightarrow\mathbb{R}$ is a
High-Order Control Barrier Function (HOCBF) with relative degree $m$ for the
given system 1 if there exist differentiable class $K$ functions
$\alpha_{1},\ldots,\alpha_{m}$ and control input $u$ that satisfy:

$\displaystyle
L_{f}^{m}h(x,t)+L_{g}L_{f}^{m-1}h(x,t){u}+\frac{\partial^{m}h(x,t)}{\partial
t^{m}}$ (6) $\displaystyle+G(h(x,t))+\alpha_{m}(\Phi_{m-1}(x,t))\geq 0,$
$\text{for all }(x,t)\in B_{1}(t)\cap B_{2}(t)\cap\ldots\cap
B_{m}(t)\times[t_{0},\infty)$

. The equation above involves $G(h(x,t))$ which represents the remaining Lie
derivatives along $f$ and partial derivatives with respect to $t$ that have a
degree equal to or lower than $m-1$ as discussed in [12].

###### Theorem 1

From the HOCBF given in Definition 4 with its related sets $B_{1}(t)$,
$B_{2}(t)$, …, $B_{m}(t)$ defined in (5), if $x(t_{0})\in B_{1}(t_{0})\cap
B_{2}(t_{0})\cap\ldots\cap B_{m}(t_{0})$, then any Lipschitz continuous
controller $u(t)\in\mathcal{U}$ that satisfies (6) for all $t\geq t_{0}$
renders the sets $B_{1}(t)$, $B_{2}(t)$, …, $B_{m}(t)$ forward invariant for
system (1).

##  III CBF OVER TIME INTEGRALS

In this section, we extend the traditional CBF formulation to encompass
constraints that are based on the accumulated effect over time, rather than
relying solely on instantaneous state values. Let $x(t)$ be the trajectory of
an agent, and $p(q;x)$ be a scalar field function $d\to_{+}$ representing the
effect of the agent at each point $q$ when the agent is at location $x$. We
define an integral cost functional of the form:

$J(q,t)=\int_{0}^{t}p\bigl{(}q,x(\tau)\bigr{)}\mathrm{d}\tau,$ (7)

which is a functional representing the cumulative effect at a given location
$q$.

###### Problem 1

Given a set $\mathcal{Q}\subseteq\mathbb{R}^{d}$, let $J_{\text{limit}}$ be
the maximum permissible value of $J(q,t)$. Find constraints on the control
input $u$ for the dynamical system (1) such that the resulting trajectory
$x(t)$ satisfies $J(q,t)\leq J_{\text{limit}}$ for every time $t>0$ and every
location $q\in\mathcal{Q}$.

We propose to tackle Problem 1 using the following function as a CBF

$h_{J}(x,t)=\min_{q\in\mathcal{Q}}\bigl{(}J_{\textrm{limit}}-J(q,t)\bigr{)}.$
(8)

In the following we derive constraints for different types of sets
$\mathcal{Q}$: single points, line segments, and general polygons. Our
analysis primarily focuses on two-dimensional (2-D) spaces, but the same
theory could be extended to three-dimensional (3-D) spaces, although we leave
this generalization for future research.

###  III-A Constraint for a single point

If the set $\mathcal{Q}$ contains a single point, $\mathcal{Q}=\\{q_{0}\\}$,
the Control Barrier Function (8) reduces to:

$h_{J}(x,t;q_{0})=J_{\text{{limit}}}-J(q_{0},t).$ (9)

Taking derivatives until the control $u$ appears explicitly we obtain:

$\displaystyle\dot{h}_{J}(x,t;q_{0})$
$\displaystyle=-p\bigl{(}q_{0},x\bigr{)},$ (10)
$\displaystyle\ddot{h}_{J}(x,t;q_{0})$
$\displaystyle=-\nabla_{x}p\bigl{(}q_{0},x\bigr{)}^{\mathrm{T}}(f(x)+g(x)u);$
(11)

here and for the remainder of the paper, we assume that
$\nabla_{x}p\bigl{(}q_{0},x\bigr{)}^{\mathrm{T}}g(x)\neq 0$, which implies
that $h_{J}$ has relative degree $m=2$. Applying the HOCBF framework reviewed
above, we have:

$\displaystyle\Phi_{J,1}(x,t;q_{0})$
$\displaystyle:=-p\bigl{(}q,x(t)\bigr{)}+\alpha_{1}h_{J}(x,t)$ (12)
$\displaystyle\Phi_{J,2}(x,t;q_{0})$
$\displaystyle:=\dot{\Phi}_{J,1}(x,t)+\alpha_{2}\Phi_{J,1}(x,t)$ (13)

The constraint then becomes:

$\Phi_{J,2}(x,t;q_{0})\geq 0$ (14)

This constraint is linear in $u$, and will be incorporated in a Quadratic
Program (QP). Overall, this case is a relatively straightfoward application of
the HOCBF framework, with the only consideration being that the integral in
the cumulative cost (7) makes the use of high-order CBFs necessary (i.e., the
relative degree is $m\geq 2$) even when the dynamics (1) is first-order.

###  III-B Constraint for a line segment using upper bounds

This section considers the case where the set $\mathcal{Q}$ is a line segment
with endpoints $q_{0},q_{K}$. We use the following parametrization of the set,
with $s\in[0,1]$:

$q(s)=(1-s)q_{0}+sq_{K};$ (15)

Ideally, we would like to satisfy the constraint on $\Phi_{J2}$ for all the
points in the set (i.e., all the values of $s$). However, this leads to having
an infinite number of constraints because the parameter $s$, which dictates
the point’s location on the line segment, is a continuous value. An
alternative strategy is to keep track of $s^{*}$ and $J_{\max}$ defined by the
equation:

$\displaystyle s^{*}(t)$
$\displaystyle=\operatorname*{\arg\\!\max}_{s}J(q(s),t)$ (16) $\displaystyle
J_{\max}(t;q_{0},q_{1})$
$\displaystyle=J\Bigl{(}q\bigl{(}s^{*}(t)\bigr{)},t\Bigr{)}$ (17)

Note that we employ the $\max$ operator (instead of the minimum operator) due
to the presence of a negative sign in the minimization of (8).

However, this strategy would require updating $s^{*}$, which also requires
keeping track of the full function $J(q(s),t)$ (which is infinite-dimensional,
since it is a continuous function).

In practice, it is necessary to introduce some form of appoximation. We
consider two strategies: a naïve approximation of the segment with points, and
an approximation of the CBF constraint with an upper bound that can be easily
updated.

####  III-B1 Naïve solution using discretization

Let $\bar{\mathcal{Q}}=\\{q_{k}\\}_{k=0}^{K}$ be a discretization of the
original set $\mathcal{Q}$ with $K$ points. The most naïve way to approximate
(8) is to transform it into a series of constraints $h(x,t;q_{k})\geq 0$ for
every $k$ in $\\{0,\ldots,K\\}$, where each point $q_{k}$ is an element of
$\bar{\mathcal{Q}}$. We then define a HOCBF constraint

$\Phi_{J,2}(x,t;q_{k})\geq 0$ (18)

for every $k$.

Effectively, this replaces $\min_{q\in\mathcal{Q}}J(q,t)$ in (8) in with
$\min_{q\in\bar{\mathcal{Q}}}J(q,t)$. This, however, causes two problems: in
order to obtain a good approximation, we might need a large number of points
(thus incurring in the curse of dimensionality). More importantly, there might
always be a point $q\in\mathcal{Q}$ for which $J(q,t)>J_{\textrm{limit}}$,
independently from how fine the discretization is (i.e., independently from
the value of $K$). This inherent deficiency prompts the next approach.

####  III-B2 Approximation without discretization

Rather than attempting to track the maximum cumulative effect over $s$
directly, we propose to compute a bound. Let $q_{0}$, $q_{K}$ be two
consecutive points on the boundary of the region we want to protect from the
cumulative noise. We would like to enforce a CBF constraint on $J_{\max}$ in
(17), as that would protect all the points on the line between $q_{0}$ and
$q_{K}$. However, keeping track of $J_{\max}$ is computationally costly
because we would need to keep track of $s^{*}$, which in turn would require a
maximization of $J(q(s),t)$ (which would require evaluating integrals for the
sequence of $s$ decided by the solver, for every time step).

Instead, we define an upper bound $\bar{J}$ on $J_{\max}$ as:

$\displaystyle\bar{J}(0;q_{0},q_{K})$ $\displaystyle=J_{\max}(0)$ (19)
$\displaystyle\dot{\bar{J}}(t;q_{0},q_{K})$
$\displaystyle=\max_{s}p(q(s),x(t)).$ (20)

Note that the argument of the maximum in (20) is the kernel $p$, not the
integral over time. Intuitively, we propose to use $\bar{J}$ to define a lower
bound on (8), and then use this lower bound as a CBF we can enforce safety for
every point.

###### Proposition 1

Define the _lower-bound_ CBF

$h_{\bar{J}}(x,t;q_{0},q_{K})=\min_{q\in\mathcal{Q}}\bigl{(}J_{limit}-\bar{J}(t;q_{0},q_{K})\bigr{)},$
(21)

and the corresponding HOCBFs functions

$\displaystyle\bar{\Phi}_{1}(x,t;q_{0},q_{K})$
$\displaystyle:=-p\bigl{(}q,x(t)\bigr{)}+\alpha_{1}h_{\bar{J}}(x,t),$ (22)
$\displaystyle\bar{\Phi}_{2}(x,t;q_{0},q_{K})$
$\displaystyle:=\dot{\bar{\Phi}}_{1}(x,t)+\alpha_{2}\bar{\Phi}_{1}(x,t).$ (23)

If the control $u(t)$ satisfies the HOCBF constraint

$\bar{\Phi}_{2}(x,t;q_{0},q_{K})\geq 0$ (24)

then $J_{\max}(t)\leq J_{limit}$ for all $t\geq 0$.

###### Proof:

First, from the definition of $J_{\max}$ in (17) and the fundamental theorem
of calculus,

$\dot{J}_{\max}(t)=p(q(s^{*}),x(t)).$ (25)

From the definition of $\bar{J}$ in (20), we have

$\dot{\bar{J}}(t)=\max_{s}\dot{J}(t)=\max_{s}p(q(s),x(t))\geq\dot{J}_{\max}(t),$
(26)

where the last inequality is given by the fact that the location $s^{*}$ where
$J$ is maximum is not necessarily also the one that is currently increasing
the most, i.e., where $\dot{J}$ is maximum.

From 19 and (26), and applying Gronwall’s comparison lemma [4] we have (24).

$J_{\max}(t)\leq\bar{J}(t)$ (27)

for all $t\geq 0$. Then, if (24) holds, from Theorem 1 we have

$J_{\max}(t)\leq\bar{J}(t)\leq J_{limit}.$ (28)

The claim follows. ∎

###  III-C Constraints for a general 2-D polygon using upper bounds

Given a polygon $\mathcal{P}$ defined by vertices
$\\{v_{1},v_{2},\ldots,v_{N}\\}$ in $\mathbb{R}^{n}$, the boundary of $P$ can
be defined as a sequence of linear segments $[v_{i},v_{i+1}]$ for
$i=1,2,\ldots,N$, with $v_{N+1}=v_{1}$ to complete the polygon. We assume the
use of a regular CBF to enforce $x\notin\mathcal{Q}$. As a result, we can
reduce the goal of protecting all segments in $\mathcal{Q}$ to the goal of
protecting the boundary of $\mathcal{P}$. we extend the equation (8), where
the CBF $h_{J}$ is adapted to address the cumulative effect constraints across
the polygon’s boundary segments. For each segment $[v_{i},v_{i+1}]$, we can
apply the methodology described in Section III-B to define the lower-bound CBF
$h_{\bar{J}}(x,t;v_{i},v_{i+1})$ for each segment $i\in\\{1,\ldots,N\\}$,
which lead to $N$ HOCBF constraints of the form

$\bar{\Phi}_{2}(x,t;v_{i},v_{i+1})\geq 0.$ (29)

Furthermore, to ensure that the agent doesn’t pass through this polygon, the
obstacle CBF constraint is applied alongside equation 29.

##  IV CASE STUDY

This section outlines a case study focusing on the application of our CBF
approach on a simplified model of a quadrotor with first-order-integrator
dynamics navigating through a noise-sensitive environment. We consider a
simple isotropic distance-based noise model as the scalar field. The goal is
then to guarantee that the accumulated noise at any point on a building
(modeled as a polygon) is below the desired limit.

###  IV-A Quadrotor model

In our representation, the quadrotor is illustrated as comprising four
propellers, represented by four disks, tangent to each other, each with its
own center. This arrangement surrounds the quadrotor’s central point, which we
denote as $x$. These individual circles are referred to as “propeller disks” ,
and $r_{\text{quad}}$ is the radius of each disk.

###  IV-B Obstacle model

The obstacle is represented as a static box and for the purpose of collision
avoidance, the quadrotor’s proximity to the obstacle is quantified by the
position of the nearest propeller disk. Within the context of the quadrotor
model, $x_{\text{quad}}$ is defined as the center of the propeller disk that
is closest to the box at any given time. Also, $x_{\text{box}}$ is defined as
the point on the surface of the box that is nearest to the aforementioned
propeller disk. These are used to define the CBF for the obstacle.

###  IV-C Application of CBF-QP for obstacles

The dynamics of the quadrotor is modeled as a single integrator $\dot{x}=u$.
We define the barrier function, $h_{obs}(x)$, as $h_{obs}(x)=\lVert
x_{\text{{quad}}}-x_{\text{{box}}}\rVert-r_{\text{{quad}}}$. The CBF-QP for
the quadrotor becomes

$\displaystyle\min_{u}$ $\displaystyle\lVert u-k_{\textrm{ref}}(x)\rVert^{2}$
(30) subject to $\displaystyle\dot{h}_{obs}(x)+\alpha_{1}h_{obs}(x)\geq 0$

where $k_{\textrm{ref}}(x)=x_{\textrm{goal}}-x(t)$.

Herein, $x(t)$ denotes the quadrotor’s position at discrete time intervals,
while $x_{\text{{goal}}}$ represents the quadrotor’s intended destination. An
example of this is shown in the simulation section. In the case of multiple
obstacles, the CBF-QP can be extended with multiple constraints, each
representing a different obstacle. For each obstacle, we can define a barrier
function $h^{(i)}_{obs}(x)$ where $i$ is the obstacle index, ensuring that the
quadrotor maintains a safe distance from all obstacles. The QP can be
augmented to accommodate all these constraints, ensuring a collision-free
trajectory even in cluttered environments.

###  IV-D Parabolic kernel noise model

We represent the state-dependent noise scalar field $p(q,x(\tau))$ for the
quadrotor using a parabolic kernel:

$p(q,x(\tau))=\begin{cases}A-\sigma\lVert q-x(\tau)\rVert^{2}&\text{if }\lVert
q-x(\tau)\rVert\leq\sqrt{\frac{A}{\sigma}}\\\ 0&\text{otherwise}\end{cases}$
(31)

In this formulation, $A$ represents the peak value of the kernel, indicative
of the maximum intensity of the noise (which it a the quadrotor’s center). The
parameter $\sigma$ controls the width of the parabola, essentially determining
the spread of the kernel. The right hand side $\sqrt{\frac{A}{\sigma}}$
represents the effective radius within which the kernel possesses non-zero
values, marking the boundary of the noise influence. Note that this is a very
simplified model, but it captures the fact that the influence of the noise on
point $q$ is nonlinear and state-dependent.

In general, the computation of the bound $\bar{J}$ requires the maximization
of $p$ (equation (26)). For the case of the parabolic kernel (31) and where
the set $\mathcal{Q}$ is a line parametrized by $s$, the optimization problem
can be solved in closed form.

###### Lemma 1

Given the parabolic kernel noise model (31) defined over the line segment
parameterized by $s$ with endpoints $q_{0}$ and $q_{K}$, the optimization
problem is to maximize $p(q(s),x(\tau))$. The solution to this optimization
problem is:

$s^{*}=\frac{(x(\tau)-q_{0})^{T}(q_{K}-q_{0})}{|q_{K}-q_{0}|^{2}}.$ (32)

The maximum corresponds to the point $q^{*}\in\ q_{0},q_{K}$:

$q^{*}=\begin{cases}q_{0}&\text{if }s^{*}<0\\\ q_{K}&\text{if }s^{*}>0\\\
q(s^{*})&\text{otherwise}\end{cases}$ (33)

###### Proof:

We first express $p(q(s),x(t))$ in terms of $s$:

$p_{s}(q(s),x(t))=A-\sigma\lVert(1-s)q_{0}+sq_{K}-x\rVert^{2}$

subject to:

$\lVert{(1-s)q_{0}+sq_{K}-x(t)}\rVert\leq\sqrt{\frac{A}{\sigma}}$

We differentiate $p_{s}$ with respect to $s$:

$\frac{dp_{s}}{ds}=-2\sigma\bigl{(}(1-s)q_{0}+sq_{K}-x(t)\bigr{)}^{\mathrm{T}}(q_{K}-q_{0})$

Setting this to zero, we get:

$\bigl{(}q(s)-x(t)\bigr{)}^{T}(q_{K}-q_{0})=0$

Solving for $s$:

$0=((1-s)q_{0}+sq_{K}-x(t))^{T}(q_{K}-q_{0})\\\
=(q_{0}-x(t)+(q_{K}-q_{0})s)^{T}(q_{K}-q_{0})\\\
=(q_{K}-q_{0})^{T}(q_{K}-q_{0})s+(q_{0}-x(t))^{T}(q_{K}-q_{0}),$ (34)

from which (32) follows. ∎

The value of $s^{*}$ is used to determine the optimal point $q^{*}$ on the
line segment as shown in (33). With $q^{*}$ identified, we can now define the
barrier function $\bar{\Phi}_{2}(x,t;q_{0},q_{1})$, which incorporates the
noise footprint constraints into the control strategy, ensuring that the
cumulative effect at $q^{*}$ remains within the specified limits.

$\bar{\Phi}_{2}(x,t;q_{0},q_{K}):=-2\sigma(q-x(t))\,u\\\
+\alpha_{2}\Bigl{(}-\bigl{(}A-\sigma||q-x(t)||^{2}\bigr{)}+\alpha_{1}\bigl{(}J_{\text{limit}}-\bar{J}(q,t)\bigr{)}\Bigr{)}\geq
0$ (35)

Consequently, from (30), the refined Quadratic Programming (QP) formulation is
given as:

$\displaystyle\min_{u}$ $\displaystyle\|u-k_{\textrm{ref}}(x)\|^{2}$ (36) s.t.
$\displaystyle\dot{h}_{obs}(x)+\alpha_{1}h_{obs}(x)\geq 0,$
$\displaystyle\bar{\Phi}_{2}(x,t;q_{0},q_{K})\geq 0.$

The solution to (36) (assumin it is feasible) will respect both the obstacle
and the noise footprint constraints.

##  V SIMULATION RESULTS

This section details simulation results demonstrating the quadrotor’s
navigation using CBFs for obstacle avoidance and noise management. We showcase
how these strategies enable the quadrotor to safely navigate and comply with
environmental constraints, highlighting the practical effectiveness of our
control approach.

###  V-A Simulation parameters

We evaluate our approach with a simulation in an environment with a single
obstacle and the simulation is based on Euler’s method for integration through
time.The focus of the paper is on local control, and a single obstacle allows
us to better demonstrate the effect of the various parameters. Navigating more
complex environments would require integration with high-level path planning
algorithms such as $\mathtt{RRT^{*}}$[5], which, however, are out of scope for
this paper. The following parameters were used for the simulation: the
class-$\mathcal{K}$ function for the obstacle CBF and noise CBF constraint
uses $\alpha_{1}=3.0$; for the noise CBF constraint, we use $\alpha_{2}=6.0$;
the size of the quadrotor is given by $r_{\text{quad}}=0.1$; the Euler
integration time step is $\delta=0.1\,\mathrm{s}$; the obstacle is a square
$\mathcal{Q}=[0,1]\times[0,1]$; finally, the initial position is
$x(0)=\left[\begin{smallmatrix}3\\\ 3\end{smallmatrix}\right]$. As the
reference controller, we use a simple proportional controller
$k_{\textrm{ref}}(x)=x_{\textrm{goal}}-x$, where
$x_{\textrm{goal}}=\left[\begin{smallmatrix}-2\\\ -1\end{smallmatrix}\right]$.

###  V-B Obstacle CBF alone

Figure 1a shows the result of the simulation using the QP with only the
obstacle CBF constraints (30). In the figure, the trajectory’s goal is
symbolized by a green star while the starting position is shown in orange, and
it is observable that the quadrotor maintains a safe distance from the
obstacle towards the goal, adhering to established safety measures. Since this
simulation does not take into account the constraint $J(t)<J_{\text{limit}}$,
while the quadrotor effectively avoids the specified physical obstacles, it
may exceed the allowable cumulative impact at certain locations. As seen in
Figure 2 the cumulative noise $J_{obs}$ exceeded the limit $J_{limit}$ for the
trajectory.

(a) System trajectories for the QP-based control of the quadcopter with box
obstacle

(b) Application of Discretization

(c) QP-based control of the quadcopter with noise and obstacles using
approximation without discretization

(d) Approximation without discretization with different parameters

Figure 1: System trajectories with just the obstacle, noise using
discretization, noise but with approximation without discretization and
approximation without discretization with different parameters  Figure 2:
$J_{\textrm{obs}}(t)$, $J_{discr}(t)$, $J_{\max}(t)$, $\bar{J}$ and
$J_{limit}$ with time

###  V-C Cumulative noise CBF with discretization

In this section, we show the simulation of the quadrotor with discretization
of sets to account for the noise around the obstacle. For this simulation, we
used $J_{\textrm{limit}}=0.2$, $A=0.5$, $\sigma=0.35$ and $K=8$. Figure 1b
shows the result of the simulation, where we use eight constraints obtained by
discretizing the boundary of the obstacle and applying individual point-wise
constraints as described in Equation 18. Incorporating the noise CBF
constraints in the QP formulation, as visualized in Figure 1b, results in a
trajectory that is more considerate of the noise impact of the quadrotor
compared to the simulation shown in Figure 1a. This is evident from Figure 2,
where the cumulative noise $J_{\textrm{discr}}(t)$, representing the
cummulative noise subject to the discretized constraints, remains below
$J_{\textrm{limit}}$. Also Figure 2 shows the plots of $J_{\max}(t)$ being
below the limit, where $J_{\text{max}}(t)=\max_{i=0}^{99}J(q_{i},t)$ and
$q_{0},q_{1},\ldots,q_{99}$ are the discretized points of the square.

###  V-D Cumulative noise CBF with bound $\bar{J}$

In this section, we share the results of simulations where we applied a
Control Barrier Function (CBF) constraint on $\bar{J}$. The parameters used
are the same as above. Fig. 1c shows the quadrotor navigating toward its goal,
detouring upon encountering an obstacle. The escalating noise around this
obstacle prompts the quadrotor to maintain distance, thereby steering clear of
the accumulating noise while persistently moving toward its target. The
difference between this method and the method in Section V-C, is that in the
discretization method depicted in Figure 1b, the noise constraints are applied
only at specific points along the boundary of the obstacle. This method risks
overlooking some areas along the boundary where noise accumulation might
exceed the limit, because it does not account for the entirety of the line
segment, while in the bound $\bar{J}$ approach, a CBF constraint is imposed on
an over-approximation of the maximum cumulative noise, ensuring that no point
along the line segment will exceed the noise threshold which can be seen in
Figure 2 as $\bar{J}<J_{\text{limit}}$. The corners of the box, are points
used to symbolize the segments forming the boundary of the obstacle $q_{0}$,
$q_{1}$, $q_{2}$, and $q_{3}$.

###  V-E Other parameters

We also ran the simulation with different parameters where $A=1.0$,
$\sigma=0.55$ and the result can be seen in Fig. 1d which shows the quadrotor
giving more safe distance away from the obstacle while heading towards the
goal. This shows that by using a larger noise radius, the quadrotor adjusts
its path to keep a wider distance from the obstacle as the noise radius
increases. Finally, we ran a simulation with
$x_{\textrm{goal}}=\left[\begin{smallmatrix}0.5\\\ 1\end{smallmatrix}\right]$
near the obstacle, and from Figure 3a, it shows the quadrotor unable to get to
the goal due to the accumulated noise at the obstacle which can be also seen
in Figure 3b which shows the position of the robot at each time step and the
plot shows the quadrotor bouncing back.

(a) Goal at the obstacle

(b) Position at each time step

Figure 3: Trajectory of the quadrotor when the goal is at the boundary of the
obstacle

##  VI CONCLUSIONS

We presented a method to consider integral cost functionals into the framework
of higher-order control barrier functions, enabling constraints over entire
trajectories rather than just instantaneous states. Additionally, we present
an efficient over-approximation method to handle constraints over large
regions. Our simulations validate the ability of the agent (a simplified
quadrotor) to navigate safely through environments while mitigating the impact
of noise in sensitive areas. While we used acoustic noise emitted from
quadrotors as an illustrative application, the theory we presented could be
applied to other scenarios, such as managing light emission in microscopy or
controlling vehicle pollution dispersion. Future research may extend this work
to multi-agent systems, consider more accurate noise propagation and quadrotor
dynamical models, and combine the our low-level control constraints into
complete path planning solutions for cluttered environments.

## References

  * [1] A. D. Ames, J. W. Grizzle, and P. Tabuada.  Control Barrier Function based quadratic programs with application to adaptive cruise control.  In 53rd IEEE Conference on Decision and Control, pages 6271–6278, 2014.
  * [2] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada.  Control barrier function based quadratic programs for safety critical systems.  IEEE Transactions on Automatic Control, 62(8):3861–3876, 2017\.
  * [3] U. Borrmann, L. Wang, A. D. Ames, and M. Egerstedt.  Control barrier certificates for safe swarm behavior.  IFAC-PapersOnLine, 48(27):68–73, 2015.
  * [4] T. H. Gronwall.  Note on the derivatives with respect to a parameter of the solutions of a system of differential equations.  Ann. Math. (2), 20:292–296, 1919.
  * [5] S. Karaman and E. Frazzoli.  Sampling-based algorithms for optimal motion planning, 2011.
  * [6] L. Lindemann and D. V. Dimarogonas.  Control Barrier Functions for signal temporal logic tasks.  IEEE Control Systems Letters, 3(1):96–101, 2019.
  * [7] I. Tezuka and H. Nakamura.  Safety assist control for nonautonomous control-affine systems via time-varying control barrier function.  IFAC-PapersOnLine, 56(1):187–192, 2023.  12th IFAC Symposium on Nonlinear Control Systems NOLCOS 2022.
  * [8] L. Wang, A. Ames, and M. Egerstedt.  Safety barrier certificates for heterogeneous multi-robot systems.  In IEEE American Control Conference, 2016.
  * [9] L. Wang, A. D. Ames, and M. Egerstedt.  Safety barrier certificates for collisions-free multirobot systems.  IEEE Transactions on Robotics, 33(3):661–674, 2017.
  * [10] G. Wu and K. Sreenath.  Safety-critical and constrained geometric control synthesis using control lyapunov and control barrier functions for systems evolving on manifolds.  In IEEE American Control Conference, pages 2038–2044. IEEE, 2015\.
  * [11] G. Wu and K. Sreenath.  Safety-critical control of a planar quadrotor.  In IEEE American Control Conference, pages 2252–2258. IEEE, 2016\.
  * [12] W. Xiao and C. Belta.  Control Barrier Functions for systems with high relative degree.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 474–479, 2019.
  * [13] G. Yang, C. Belta, and R. Tron.  Continuous-time signal temporal logic planning with control barrier functions.  In 2020 American Control Conference (ACC), pages 4612–4618, 2020\.


[a,b,c] Nora Brambilla

# Effective Field Theories and Lattice QCD for the X Y Z frontier

###### Abstract

Exotic states have been predicted before and after the advent of QCD. In the
last decades they have been observed at accelerator experiments in the sector
with two heavy quarks, at or above the quarkonium strong decay threshold and
called X Y Z states. These states offer a unique possibility for investigating
the dynamical properties of strongly correlated systems in QCD. I will show
how an alliance of nonrelativistic effective field theories and lattice can
allow us to address these states in QCD. In particular I will explain what are
the opportunities and challenges of lattice QCD in this respect and which new
tools should be developed.

##  1 Introduction

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
inception 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 [1, 2, 3, 4]. 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^{--}$ (where $J$
is the spin, $P$ is the parity and $C$ is the charge-conjugation quantum
number). Meanwhile, the Particle Data Group (PDG) has proposed a new naming
scheme [5], that extends the scheme used for ordinary quarkonia, in which the
new names carry information on the $J^{PC}$ quantum numbers, see [1] 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 [7] . 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 and surprinsingly 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 have by now measured their signatures in the
spectrum and 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 [6], 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. [8, 9, 10, 11]. 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.
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 in
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. Addressed like that, the problem turns out to be extremely
challenging. I will show how the presence of two heavy quarks gives us a way
to simplify the problem and to address it using nonrelativistic effective
field theories (NR EFTs) and lattice. The existence of a large scale, the
scale of the mass of the heavy quark $m$, allows us to take advantage of scale
factorization and to use lattice QCD to evaluate purely glue dependent
correlators at the low energy scale. Such correlators are not depending on
flavor and are universal, i.e. the can be evaluated once and used for
prediction on a large number of states. I will argue that this method is
complementary and often preferable with respect to an ab initio lattice
evaluation of the properties of each single state. To pursue this aim, i.e. to
perform lattice evaluation directly on factorized low energy objects in the
framework of EFTs, we have founded a dedicated lattice collaboration called
TUMQCD [14].

The paper is organized in the following way: in section two I summarize how
the combination of NREFTs and lattice has allowed us to address quarkonium
properties in vacuum and in medium in QCD. In the following sections I put
forward a similar framework for exotica. In particular in section three I
introduce the so called BOEFT (Born Oppenheimer Effective Field Theory) to
describe hybrids and I report the predictions of the theory and how progress
depends on more precise and/or new lattice evaluations of low energy
correlators, even quenched. In Section four I extend this description to
states with light quarks, which would give a theory of the tetraquarks and
molecular pictures. In also examine how heavy ion collisions can help us
supplying complementary information and what is the role of the EFT and
lattice in this. A list of lattice calculations that would be needed to gain
progress is given in each section and in the outlook. In this latter I also
discuss tools, methods and open challenges to get these low energy correlators
from the lattice.

##  2 NREFTs and Lattice for Quarkonium

The study of quarkonium in the last few decades has witnessed two major
developments: the establishment of NREFTs and progress in lattice QCD
calculations of excited states and resonances, with calculations at light-
quark masses. Both allow for precise and systematically improvable
computations that are (to a large extent) model-independent. It is precisely
this advancement in the understanding of quarkonium and quarkonium-like
systems inside QCD that makes quarkonium exotics particularly valuable. 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.

###  2.1 Physical scales of quarkonium

Heavy quarkonia are systems composed by two heavy quarks (charm, bottom), with
mass $m$ larger than the “QCD confinement scale” $\Lambda_{\rm QCD}$, so that
$\alpha_{\rm s}(m)\ll 1$ holds. From the quarkonia spectra it is evident that
the difference in energy levels is much smaller than the quark mass and
therefore quarkonia are nonrelativistic systems. As such they are multiscale
systems. Being nonrelativistic, quarkonia are characterized by a small
parameter, the heavy-quark relative velocity $v=|\vec{v}|$ in the rest frame
of the meson, where $v^{2}\sim 0.1$ for the $b\bar{b}$, $v^{2}\sim 0.3$ for
$c\bar{c}$ systems. This parameter induces a hierarchy of dynamically
generated energy scales: starting from the mass $m$ of the heavy quark (hard
scale), there is the typical relative momentum $p\sim mv$, corresponding to
the inverse Bohr radius $r\sim 1/(mv)$ (soft scale), and a typical binding
energy $E\sim mv^{2}$ (ultrasoft (US) scale). The description of the heavy
$Q\bar{Q}$ systems depends finally on the relation of $\Lambda_{\rm QCD}$ to
the above mentioned scales. Clearly, for energy scales close to $\Lambda_{\rm
QCD}$ there is no longer any perturbative description and one has to rely on
nonperturbative methods. Regardless of this, the nonrelativistic hierarchy
$m\gg mv\gg mv^{2}$ persists also below the $\Lambda_{\rm QCD}$ threshold, as
long as $v$ is small. While the hard scale $m$ is always assumed to be larger
than $\Lambda_{\rm QCD}$, different situations may arise for the other two
scales.

The hierarchy of nonrelativistic scales makes the difference between heavy
quarkonia and heavy-light mesons, which are characterized by just two scales:
$m$ and $\Lambda_{\rm QCD}$. This renders the theoretical description of
quarkonium much more complicated. All these scales get entangled in a typical
amplitude involving a quarkonium observable. Moreover the nonrelativistic
bound state problem is “nonperturbative”, even at small $\alpha_{s}$, i.e.
infinite series of diagram needs to be resummed, differently from an on shell
scattering calculation. 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 (US)
gluons and light quarks live long enough that a bound state has time to form
and, therefore, are sensitive to the scale $mv^{2}$. US gluons are responsible
for phenomena similar to the the Lamb shift in QED. To address this
nonrelativistic multiscale system in quantum field theory is a real challenge.
It is so even in QED with positronium and it is even more difficult here with
the nonperturbative effects. Notice that even if you were interested in a
lattice first principle simulations, still the existence of a set of widely
separated scales will make the lattice calculation very challenging, 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
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 [15, 16]. 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.

###  2.2 Physics at the scale $m$: NRQCD

Nonrelativistic QCD (NRQCD) [17, 18], 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},...),$ (1)

where $c_{n}$ are Wilson coefficients that contain the contributions from the
scale $m$; they can be perturbatively calculated by matching the full QCD
result to 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. The quarkonium states
$|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$
(2)

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 with great success by the HPQCD collaboration to calculate
quarkonium spectra and decays111Notice that the NRQCD leading piece of the
Lagrangian differently from the HQET one is not renormalizable, so that
lattice calculation cannot go to the continuum limit but have to stay in the
scaling window.. 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 momentum in the center of mass frame 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$ (3)

where the $\sigma_{n}$ are short-distance coefficients, and the matrix
elements $\langle{\cal K}_{n}^{4{\rm fermions}}\rangle$ are vacuum-expectation
values of objects similar to the four-fermion operators in decays and
containing both color singlet and color octet contributions. 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. Notice however
that it has never been possible up to now to calculate these LDMEs on the
lattice due to their cumbersome definition.

###  2.3 Physics at the scale $mv$: pNRQCD

Quarkonium formation happens at the scale $mv$. The suitable EFT is potential
non relativistic QCD, pNRQCD [20, 21, 19], which follows from NRQCD by
integrating out the scale $mv\sim r^{-1}$. The soft scale $mv$ may be either
larger or smaller than the confinement scale $\Lambda_{\rm QCD}$ depending on
the radius of the quarkonium system. 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. 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 and cannot be calculated with an expansion in $\alpha_{s}$.
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 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 prediction. The nonperturbative phyiscs 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.222We
notice that leading piece of the pNRQCD lagrangian is renormalizable.

####  2.3.1 Weakly coupled pNRQCD

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) [21]:

$L_{\rm pNRQCD}=\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},...),$ (4)

where ${\bf R}$ is the center of mass position and $O_{k}$ are the operators
of pNRQCD. The matrix elements of these operators 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 (1).

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. Nonpotential interactions, associated
with the propagation of low-energy degrees of freedom are, in general, present
as well, and 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. 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}$ from quarkonium.

###  2.4 Strongly coupled pNRQCD

When $mv\sim\Lambda_{\rm QCD}$ the soft scale is nonperturbative, the matching
cannot be performed in perturbation theory any more. Only singlet degrees of
freedom exist and they include $Q\bar{Q}$ states, hybrids $Q\bar{Q}g$ states
and glueballs. Since the physics is nonperturbative we need to use lattice
input to construct the EFT. In particular we need the lattice evaluation of
the gluonic static energies of $Q\bar{Q}$ pair: the have been calculated on
the lattice since long [22, 23, 24] and recently updated in [25, 26] and they
use generalized Wilson loops. The gluonic static energies, $E_{\Gamma}$ in
Fig. 2, are classified according to representations of the symmetry group
$D_{\infty\,h}$, typical of molecules, and labeled by
$\Lambda_{\eta}^{\sigma}$ (see Fig. 2): $\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$
This set of static energies will be fundamental to address the exotics.

Figure 1: Quarkonium hybrid symmetries.

Figure 2: Hybrid static energies, $E_{\Gamma}$, in lattice units, from [23].

Since at this moment we are dealing with $Q\bar{Q}$ states below threshold we
are interested merely in $\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=\bar{Q}Q$ is the
only low-energy dynamical degree of freedom in the pNRQCD Lagrangian (up to US
pions) which reads [27, 28, 19]:

$\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}\,$ (5)

and lends support to potentials models in this regime, with the difference
that the singlet potential $V_{s}(r)=V_{0}+V_{1}/m+V_{2}/m^{2}$ is the QCD
potential calculated in the matching [27, 28] and appear to have differences
with respect to the phenomenological potential models. The potential at order
$1/m^{2}$ contains a velocity dependent and a spin dependent part, with a spin
orbit, spin-spin and spin tensor interaction. Notice that spin effects 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.
These potentials are nonperturbative quantities in the form of expectation
values of generalized (with chromolectric and chromagnetic insertions) gauge
invariant Wilson-loop operators to be evaluated on the lattice. Actually it is
the EFT that 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. Some of
these generalized Wilson loops have been calculated on the lattice (only
quenched) [30, 31, 32, 33, 34] but some contributions at order $1/m^{2}$ are
still to be calculated. Using these potentials, all the masses for heavy
quarkonia away from threshold can be obtained by the solution of the
Schrödinger equation with such potentials. Lorentz invariance is still there
in the form if exact relations among potentials and it has been observed on
the lattice. Decays are described by calculating the imaginary parts of the
potentials [35] 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.

####  2.4.1 Quarkonium production

NRQCD is the theory used to study quarkonium production, however one of the
main problem is the proliferation of LDMEs, i.e. the nonperturbative matrix
elements of singlet and octet four fermions operators: they are flavor
dependent and increase in number when one seeks precision in the expansion in
$v$. While those containing singlet operators can be reduced to the wave
functions, the octet ones had to be extracted from the data, with
contradictory results. The lower energy factorization obtained in pNRQCD may
have a great impact also on quarkonium production: it has been shown that, for
the case of inclusive quarkonium production in P wave, the octet LDMEs can be
factorized in the product of the quarkonium wave function derivative in the
origin and a chromoelectric correlator [36, 37]

$\langle\Omega|{\cal
O}^{h_{Q}}({}^{1}S_{0}^{[8]})\Omega\rangle=3\times\frac{3N_{c}}{2\pi}|R^{(0)}{}^{\prime}(0)|^{2}\frac{1}{9N_{c}m^{2}}{\cal
E},$ (6)

where ${\cal E}$ is the dimensionless gluonic correlator

${\cal
E}=\frac{3}{N_{c}}\int_{0}^{\infty}dt\,t\;\int_{0}^{\infty}dt^{\prime}\,t^{\prime}\;\langle|\Omega|\Phi_{\ell}^{{\dagger}ab}\Phi_{0}^{{\dagger}ad}(0;t)g{E}^{d,i}(t)g{E}^{e,i}(t^{\prime})\Phi_{0}^{ec}(0;t^{\prime})\Phi_{\ell}^{bc}|\Omega\rangle.$
(7)

where $\Phi_{\ell}$ is a Wilson line along the (arbitrary) direction $\ell$
(up to infinity) in the adjoint representation and $|\Omega\rangle$ is the QCD
vacuum state. This is similar to what one gets in the decay factorization of
the NRQCD matrix elements [35]

$\langle
h_{Q}|\psi^{\dagger}T^{a}\chi\chi^{\dagger}T^{a}\psi|h_{Q}\rangle=\frac{3N_{c}}{2\pi}|R^{(0)}{}^{\prime}(0)|^{2}\frac{1}{9N_{c}m^{2}}{\cal
E}_{3},$ (8)

where the correlator ${\cal E}_{3}$ is defined by

${\cal
E}_{3}=\frac{1}{2N_{c}}\int_{0}^{\infty}t^{3}\,dt\langle\Omega|g{E}^{a,i}(t)\Phi_{0}^{ab}(t,0)g{E}^{b,i}(0)|\Omega\rangle.$
(9)

Hence the pNRQCD expressions have more predictive power than the corresponding
NRQCD ones. In particular, correlators determined with charmonium data may be
used to compute color octet matrix elements and hence observables in the
bottomonium sector.

###  2.5 Quarkonium in medium

pNRQCD has been extended to the case of finite temperature $T$ which allows to
calculate in the weak coupling the finite $T$ quarkonium interaction potential
and the thermal effects to the energies and widths [38]. In the strong
coupling the potential is obtained from a lattice calculation of the Wilson
loop at finite $T$ [39]. In particular combining pNRQCD and open quantum
system [41], it has been possible to describe the nonequilibrium evolution of
small quarkonia systems (bottomonium) inside the strongly coupled Quark Gluon
Plasma in the hierarchy $m\gg 1/r\gg T,m_{D}\gg mv^{2}$ with an evolution
equation for the singlet and octet density matrix of the Linblad type [42].
The interesting thing is that the properties of the strongly coupled QGP are
carried just by two transport coefficients that are appropriate correlators of
electric fields at finite $T$ [41]:

$\kappa=\frac{g^{2}}{6N_{c}}{\rm
Re}\int^{\infty}_{-\infty}\,\mathrm{d}s~{}\Big{\langle}{\rm
T}E^{a,i}(s)\phi^{ab}(s,0)E^{a,i}(0)\Big{\rangle},\gamma=-\frac{g^{2}}{6N_{c}}{\rm
Im}\int^{\infty}_{\infty}\,\mathrm{d}s~{}\Big{\langle}{\rm
T}E^{a,i}(s)\phi^{ab}(s,0)E^{a,i}(0)\Big{\rangle},$ (10)

where $\langle\dots\rangle$ is a shortcut to indicate the QCD average at
finite T on the QCD vacuum $|\Omega\rangle$ (and ${\rm T}$ is time ordering.)
The quantity $\kappa$ is the heavy quark momentum diffusion coefficient.
$\kappa$ has been evaluated on the lattice [43] for a medium locally in
equilibrium at a temperature $T$. In particular, it is important to obtain
their $T$ dependence in a large range of $T$ as done in [29]. Using the EFT
one can relate these transport coefficients to the thermal modification of the
energy levels and to the thermal widths of quarkonium, which allows us to use
unquenched lattice calculations [44, 45] to determine them [40].

This formalism can be extended to describe the evolution of exotic states
inside the QGP created in heavy ion collisions.

###  2.6 Role of the Lattice and Input needed from Lattice

Summarizing, for what concerns quarkonium below threshold, lattice results
have been key in order to formulate strongky coupled pNRQCD and lattice input
is still needed, precisely, we need:

  * •

Lattice calculation of chromoelectric and chromomagnetic correlators local in
time and nonlocal in time: they contribute to quarkonia spectra and decay [19,
35].

  * •

Lattice calculation of generalized static Wilson loops with insertions of
chromoelectric and chromomagnetic fields: they contain the information on the
nonperturbative singlet static potential and its relativistic corrections.

  * •

Lattice calculation of the chromolelectric correlators appearing in
production: they allow to calculate the LDMEs necessary to describe quarkonium
production. They are somehow similar to what is done for TMDs and PDFs.

  * •

Lattice calculation of the transports coefficients: finite $T$ lattice
calculations of such objects in a wide range of $T$ will enable us to use the
lattice to study nonequilibrium processes [41, 42, 29]. There is not a yet a
direct evaluation of $\gamma$.

This will be instrumental to sharpen our control of quarkonium spectrum,
decays, production and propagation in medium and to prepare for the treatment
of the exotics.

##  3 Exotics at or above the strong decay threshold

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. The threshold region
remains troublesome also for the lattice, although several ab initio
calculations of the exotic masses have been recently being pionereed. Still,
$m$ is a large scale and NRQCD is applicable. 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 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 potential $V_{\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
potentials $V_{\kappa}$ by solving the corresponding Schrödinger equation, see
Figs. 4 and 4. This picture corresponds to the Born-Oppenheimer (BO)
approximation. Starting from pNRQED/pNRQCD the BO approximation can be made
rigorous and cast into a suitable EFT called Born-Oppenheimer EFT (BOEFT) [51,
52, 59, 58, 53, 54] which exploits the hierarchy of scales $\Lambda_{\rm
QCD}\gg mv^{2}$ (the analogous in QED of the energy of the electrons being
greater than the energy of nuclei or conversely the typical time of the
electrons (fast degrees of freedom) being bigger than the typical time of the
nuclei (slow degrees of freedom)).

Figure 3: Electronic static energies in QED

Figure 4: Gluonic (or Hybrid) static energies, $E_{\Gamma}$, in QCD.

BOEFT allows to address all exotics pictures (hybrids, tetraquarks, molecules
and pentaquarks) in the same framework. Of course in QCD the scale
$\Lambda_{\rm QCD}$ is nonperturbative and we need to use the lattice to
calculate the appropriate gluonic static energies, see Fig. 2.

###  3.1 Hybrids: Born Oppenheimer Effective Field Theory

In the following we consider the example of hybrids, i.e. states made by
$Q\bar{Q}g$. We consider hybrids static energies that are labelled by the
quantum numbers $\Lambda^{\sigma}_{\eta}$ , see Sec. 2.4 and Fig. 2 and we
restrict to considering the two lowest-lying static energies ones $\Pi_{u}$
and $\Sigma_{u}^{-}$. Modes associated to higher potentials are separated from
the $\Pi_{u}$ and $\Sigma_{u}^{-}$ potentials by an energy gap of order
$\Lambda_{\rm QCD}$ and are integrated out with that scale. In the $r\to 0$
limit $\Pi_{u}$ and $\Sigma_{u}^{-}$ are degenerate and correspond to a
gluonic operator with quantum numbers $\kappa\equiv K^{PC}=1^{+-}$, ${\bf K}$
being the gluon angular momentum, because the group $D_{\infty\,h}$ becomes
the more symmetric group O(3)$\times$C [21, 19, 22]. We recall that $\lambda$
is the quantum number $\lambda\equiv{\bf K}\cdot\hat{r};|\lambda|=\Lambda$.

The hybrid static energies are defined by ($\Gamma$ stays for
$\Lambda^{\sigma}_{\eta}$)

$E_{\Gamma}(r)=\lim_{T\to\infty}{i\over T}\log\langle
X_{\Gamma},T/2|X_{\Gamma},-T/2\rangle;\quad|X_{\Gamma}\rangle=\chi({\bf
x}_{2})\phi({\bf x}_{2},{\bf R})T^{a}P_{\Gamma}^{a}({\bf R})\phi({\bf R},{\bf
x}_{1})\psi^{\dagger}({\bf x}_{1})|\Omega\rangle$ (11)

where $P^{a}_{\Gamma}$ is some gluonic operator that generates the right
quantum numbers $\Gamma$. For our case we have $P^{a}_{\Gamma}=\hat{\bf
r}\cdot{\bf B}^{a},\hat{\bf r}\times{\bf B}^{a}$ respectively for the
$\Sigma_{u}^{-}$ and the $\Pi_{u}$, ${\bf B}^{a}$ being the chromomagnetic
field. In the short-interquark-distance limit $r\rightarrow 0$, quarkonium
hybrids reduce to gluelumps, which are color-singlet combinations of a local
static octet color $Q\bar{Q}$ source coupled to a gluonic field [21, 51]. We
define the gluelump operators, ${G}^{ia}_{\kappa}$, as the Hermitian color-
octet operators that generate the eigenstates of the zero order NRQCD
Hamiltonian $H_{0}$ in the presence of a local heavy-quark-antiquark octet
source:

$\displaystyle H_{0}{G}^{ia}_{\kappa}(\bm{R},t)|0\rangle$
$\displaystyle=\Lambda_{\kappa}{G}^{ia}_{\kappa}(\bm{R},t)|0\rangle\,,$ (12)

where $a$ is the color index and $i$ the spin component (it labels the
component of an irreducible representation of $SO(3)$ with Casimir K),
$\kappa$ labels the quantum numbers $K^{PC}$ of the gluonic degrees of
freedom. For our case ${\bf G}^{a}={\bf B}^{a}$ with $K^{PC}=1^{+-}$. The
spectrum of the mass eigenvalues, $\Lambda_{\kappa}$, called the gluelump
mass, has been computed on the lattice in Refs. [22, 24, 49] and it is given
by the gluelump correlator
$\Lambda_{\kappa}=\lim_{T\to\infty}{i/T}\ln\langle\Omega|G^{a}_{\kappa}(T/2)\phi_{ab}(T/2,-T/2)G^{b}_{\kappa}(-T/2)|\Omega\rangle$.

The static energy $E_{\Gamma}$ calculated in NRQCD at small r can be matched
to the following objects in BOEFT:

$E_{\Gamma}\equiv E_{\kappa,\lambda}=V^{(0)}_{\kappa,\lambda}={\alpha_{s}\over
6r}+\Lambda_{\kappa}+a_{\kappa\lambda}r^{2}+\cdots$ (13)

where $\alpha_{s}/6r$ is the static octet potential, $a_{\kappa\lambda}$ is a
nonperturbative coefficient that has a field theoretical definition and starts
to depend on $\lambda$ (instead that only on $\kappa$), which makes clear that
at the NLO order in the multipole expansion the system is no longer
spherically symmetric but acquires instead a cylindrical symmetry around the
heavy-quark-antiquark axis. Therefore it is convenient to work with a basis of
states with good transformation properties under $D_{\infty h}$. Such states
can be constructed by projecting the gluelump operators on various directions
with respect to the heavy-quark-antiquark axis. We define the degrees of
freedom of the BOEFT as the operator
$\hat{\Psi}_{\kappa\lambda}(\bm{r},\,\bm{R},\,t)$ defined by

$P^{i{\dagger}}_{\kappa\lambda}O^{a}\left(\bm{r},\bm{R},t\right)G_{\kappa}^{ia}(\bm{R},t)=Z_{\kappa}^{1/2}\hat{\Psi}_{\kappa\lambda}(\bm{r},\bm{R},t)\,,$
(14)

the $Z$ is a field renormalization, $P^{i}_{\kappa\lambda}$ is a projector
that projects the gluelump operator to an eigenstate of
$\bm{K}\cdot\hat{\bm{r}}$ with eigenvalue $\lambda$. For our case the
projectors $P^{i}_{\lambda}$ read
$P^{i}_{0}=\hat{r}_{0}^{i}=\hat{r}^{i}\,,P^{i}_{\pm
1}=\hat{r}^{i}_{\pm}=\mp\left(\hat{\theta}^{i}\pm
i\hat{\phi}^{i}\right)/\sqrt{2}\,,$ with
$\hat{\bm{r}}=(\sin(\theta)\cos(\phi),\,\sin(\theta)\sin(\phi)\,,\cos(\theta))$,
$\hat{\bm{\theta}}=(\cos(\theta)\cos(\phi),\,\cos(\theta)\sin(\phi)\,,-\sin(\theta))$,
$\hat{\bm{\phi}}=(-\sin(\phi),\,\cos(\phi)\,,0)$. The BOEFT is obtained by
integrating out modes of scale $\Lambda_{\rm QCD}$, i.e. the gluonic
excitation and the Lagrangian reads as

$L_{\rm BOEFT}=\int
d^{3}Rd^{3}r\,\sum_{\kappa}\sum_{\lambda\lambda^{\prime}}\textrm{Tr}\\{\hat{\Psi}^{\dagger}_{\kappa\lambda}(\bm{r},\,\bm{R},\,t)[i\partial_{t}-V_{\kappa\lambda\lambda^{\prime}}(r)+P^{i{\dagger}}_{\kappa\lambda}\frac{{\bm{\nabla}}^{2}_{r}}{m}P^{i}_{\kappa\lambda^{\prime}}]\hat{\Psi}_{\kappa\lambda^{\prime}}(\bm{r},\,\bm{R},\,t)\\}+\dots\,,$
(15)

where the trace is over spin indices of the heavy quark and antiquark, and the
ellipsis stands for operators producing transitions to standard quarkonium
states and transitions between hybrid states of different $\kappa$. The latter
are suppressed at least by $1/\Lambda_{\rm QCD}$ since the static energies for
different $\kappa$ are separated by a gap $\sim\Lambda_{\rm QCD}$. The
interesting points of the Lagrangian above is that the kinetic operator is not
commuting with the projector, which induces a non adiabatic coupling, already
known from molecular physics.

The potential $V_{\kappa\lambda\lambda^{\prime}}$ can be organized into an
expansion in $1/m$ and a sum of static, spin-dependent (SD) and independent
(SI) parts to be calculated in the matching with NRQCD. In Ref. [51] the
static potential $V^{(0)}_{\kappa\lambda}(r)$ was matched to the quark-
antiquark hybrid static energies computed on the lattice and the matrix
elements of
$P^{i{\dagger}}_{\kappa\lambda}\frac{{\bm{\nabla}}^{2}_{r}}{m}P^{i}_{\kappa\lambda^{\prime}}$
were obtained for $\kappa=1^{+-}$ and shown to contain off-diagonal terms in
$\lambda$-$\lambda^{\prime}$ that lead to coupled Schrödinger equations:

$\displaystyle\left[-\frac{1}{Mr^{2}}\,\partial_{r}r^{2}\partial_{r}+\frac{1}{Mr^{2}}\begin{pmatrix}l(l+1)+2&2\sqrt{l(l+1)}\\\
2\sqrt{l(l+1)}&l(l+1)\end{pmatrix}+\begin{pmatrix}E_{\Sigma_{u}^{-}}&0\\\
0&E_{\Pi_{u}}\end{pmatrix}\right]\begin{pmatrix}\psi_{\Sigma}^{(N)}\\\
\psi_{-\Pi}^{(N)}\end{pmatrix}=\mathcal{E}_{N}\begin{pmatrix}\psi_{\Sigma}^{(N)}\\\
\psi_{-\Pi}^{(N)}\end{pmatrix}\,,$ (16)
$\displaystyle\left[-\frac{1}{Mr^{2}}\,\partial_{r}\,r^{2}\,\partial_{r}+\frac{l(l+1)}{Mr^{2}}+E_{\Pi_{u}}\right]\psi_{+\Pi}^{(N)}=\mathcal{E}_{N}\,\psi_{+\Pi}^{(N)}\,,$
(17)

where $l$ is the total orbital angular momentum quantum number sum ${\bf
L}={\bf L}_{Q\bar{Q}}+{\bf K}$, $\mathcal{E}_{N}$ are the hybrids energies and
$E_{\Sigma_{u}^{-}}$ and $E_{\Pi_{u}}$ have been written in BOEFT in term of
the static potentials plus the gluelump mass taken from [24]:
$\Lambda_{1^{+-}}=0.87(15)$ GeV333Note that the gluelump mass is a scheme
dependent quantity like the quark mass, here it is calculated in RS scheme..

The Schrödinger equations were solved numerically and the spectrum and wave
functions of hybrid states generated the static energies labeled by
$\Sigma^{-}_{u}$ and $\Pi_{u}$ were obtained, precisely:

The functions $\psi_{\Sigma}^{(N)}$ and $\psi_{\pm\Pi}^{(N)}$ are radial wave
functions; $\psi_{\Sigma}^{(N)}$ and $\psi_{-\Pi}^{(N)}$ have negative parity
and $\psi_{+\Pi}^{(N)}$ positive one. The off-diagonal terms change the
$\Sigma$ wave function to $\Pi$ and vice versa, but they do not change the
parity. Hence $\psi_{\Sigma}^{(N)}$ mixes only with $\psi_{-\Pi}^{(N)}$, and
$\psi_{+\Pi}^{(N)}$ decouples. For $l=0$ the off-diagonal terms vanish, so the
equations for $\psi_{\Sigma}^{(N)}$ and $\psi_{-\Pi}^{(N)}$ decouple; there
exists only one parity state, and its radial wave function is given by a
Schrödinger equation with the $E_{\Sigma}$ potential and an angular part
$2/Mr^{2}$. The eigenstates of the equations (16) and (17) are organized in
the multiplets shown in Table 1 (all the $J^{PC}$ quantum numbers of a given
multiplet will be degenerate up to when we will consider spin dependent
interactions in Sec. 3.1.2).

$\,l\,$ $J^{PC}\\{s=0,s=1\\}$ $E_{\Gamma}$ $H_{1}$ $1$
$\\{1^{--},(0,1,2)^{-+}\\}$ $E_{\Sigma_{u}^{-}}$, $E_{\Pi_{u}}$ $H_{2}$ $1$
$\\{1^{++},(0,1,2)^{+-}\\}$ $E_{\Pi_{u}}$ $H_{3}$ $0$ $\\{0^{++},1^{+-}\\}$
$E_{\Sigma_{u}^{-}}$ $H_{4}$ $2$ $\\{2^{++},(1,2,3)^{+-}\\}$
$E_{\Sigma_{u}^{-}}$, $E_{\Pi_{u}}$ $H_{5}$ $2$ $\\{2^{--},(1,2,3)^{-+}\\}$
$E_{\Pi_{u}}$

Table 1: $J^{PC}$ multiplets with $l\leq 2$ for the $\Sigma_{u}^{-}$ and
$\Pi_{u}$ gluonic states. The last column shows the gluonic static energies
that appear in the Schrödinger equation of the respective multiplet.

####  3.1.1 Spectrum

Figure 5: Comparison of results from the direct lattice computations of
charmonium hybrid masses of [50] with the results of [51]. The direct lattice
mass predictions are plotted in solid lines with error bars corresponding to
the mass uncertainties. Our results for the H1, H2, H3, and H4 multiplets have
been plotted in error bands corresponding to the gluelump mass uncertainty of
$\pm 0.15$ GeV.

Keeping in the equations (16) and (17) only the heavy quark-antiquark kinetic
energy and the hybrid static energies, $E_{\Gamma}$, amounts at the
Born–Oppenheimer approximation. Keeping only the diagonal terms amounts at the
adiabatic approximation. The exact leading order equations include both
diagonal and off-diagonal terms that define the so-called non-adiabatic
coupling. As it is clear from (16) and (17) both diagonal and off-diagonal
terms contribute at the same order to the energy levels in the quarkonium
hybrid case. This situation is different from the case of (electromagnetic)
molecules, where the non-adiabatic coupling is subleading [52]. A physical
consequence of the mixing is the so-called $\Lambda$ doubling, i.e., the
lifting of degeneracy between states with the same parity. We show this effect
in Fig. 5, where you see that the multiplets $H_{1}$ and $H_{2}$ are not
degenerate and this effect is also confirmed in direct lattice evaluations.
The same effect would be seen also between $H_{4}$ and $H_{5}$ [51]. The
effect is also present in molecular physics, however, $\Lambda$ doubling is a
subleading effect there.

Figure 6: Mass spectrum of neutral exotic charmonium states obtained by
solving Eqs. (16) and (17). The experimental states of Sec. 3.1 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 [1].

Several charmonium-like states have been found by the B factories in the last
decade (for a recent review see [1]). We compare some of these states with the
hybrid spectrum obtained from solving the coupled Schrödinger equations (16)
and (17) in Fig. 6. In the bottomonium sector, the state
$Y_{b}(10890)[1^{--}]$ with mass $M_{Y_{b}}=(10.8884\pm 3.0)$ GeV found by
BELLE may be a possible $H_{1}$ candidate, for which we find
$M_{H_{1}}=(10.79\pm 0.15)$ GeV. This picture allows us to describe the
spectrum of hybrids within an EFT description of QCD (which should be
equivalent to QCD in the taken window) plus lattice input. It is complementary
to direct lattice evaluations of the spectrum in many ways: first of all it is
simpler and allows to get all the flavour of hybrids on the basis of the
evaluation of few correlators only dependent on the gluonic degrees of
freedom, it takes advantage of factorization and it lends an interpretation to
the underlying physics. The multiplets that we have obtained are distinctive
of the underlying physics: if one would construct the hybrids multiplets using
the constituent gluon picture for example a different multiplet structure
would be obtained [51]. Of course in order to validate the experimental
identifications one has to add the spin structure, and especially consider
decays and production features: also this may turn out to be simpler in the
BOEFT approach.

####  3.1.2 Spin effects in Hybrids

The spin of the quarks has been considered in direct lattice calculations of
the hybrids mass, mostly in charmonium for quenched and noncontinuum
extrapolated data [55], recently in unquenched and continuum extrapolated data
[57] and recently in bottomonium [56]. It is rare to see spin effects treated
in models and, in case they are considered, it is on the basis of a gluon
constitutent picture or extrapolating from the standard quarkonium spin
structure, which we will see is not the appropriate thing to do. In BOEFT we
could obtain for the first time the spin dependent potentials (for
$\kappa=1^{+-}$, $\lambda$ takes the values $0,\pm 1$) [53, 54, 58]:

$\displaystyle V_{\lambda\lambda^{\prime}\,SD}^{(1)}(r)$
$\displaystyle=V_{SK}(r)\left(P^{i{\dagger}}_{\lambda}\bm{K}^{ij}P^{j}_{\lambda^{\prime}}\right)\cdot\bm{S}$
$\displaystyle\quad+V_{{SK}b}(r)\left[\left(\bm{r}\cdot\bm{P}^{{\dagger}}_{\lambda}\right)\left(r^{i}\bm{K}^{ij}P^{j}_{\lambda^{\prime}}\right)\cdot\bm{S}-\left(r^{i}\bm{K}^{ij}P^{j{\dagger}}_{\lambda}\right)\cdot\bm{S}\left(\bm{r}\cdot\bm{P}_{\lambda^{\prime}}\right)\right]\,,$
(18) $\displaystyle V_{\lambda\lambda^{\prime}\,SD}^{(2)}(r)$
$\displaystyle=V_{SLa}(r)\left(P^{i{\dagger}}_{\lambda}\bm{L}_{Q\bar{Q}}P^{i}_{\lambda^{\prime}}\right)\cdot\bm{S}+V_{SLb}(r)P^{i{\dagger}}_{\lambda}\left(L_{Q\bar{Q}}^{i}S^{j}+S^{i}L_{Q\bar{Q}}^{j}\right)P^{j}_{\lambda^{\prime}}$
$\displaystyle\quad+V_{SLc}(r)\left[\left(\bm{r}\cdot\bm{P}^{{\dagger}}_{\lambda}\right)\left(\bm{p}\times\bm{S}\right)\cdot\bm{P}_{\lambda^{\prime}}+\bm{P}^{{\dagger}}_{\lambda}\cdot\left(\bm{p}\times\bm{S}\right)\left(\bm{r}\cdot\bm{P}_{\lambda^{\prime}}\right)\right]$
$\displaystyle\quad+V_{S^{2}}(r)\bm{S}^{2}\delta_{\lambda\lambda^{\prime}}+V_{S_{12}a}(r)S_{12}\delta_{\lambda\lambda^{\prime}}+V_{S_{12}b}(r)P^{i{\dagger}}_{\lambda}P^{j}_{\lambda^{\prime}}\left(S^{i}_{1}S^{j}_{2}+S^{i}_{2}S^{j}_{1}\right)\,,$
(19)

where $\left(K^{ij}\right)^{k}=i\epsilon^{ikj}$ is the angular momentum
operator for the spin-$1$ gluonic excitation and $\bm{L}_{Q\bar{Q}}$ is the
orbital angular momentum of the heavy-quark-antiquark pair, ${\bf S}$ is the
total spin of the quarks, and the projectors for $\kappa=1^{+-}$ have been
previously defined. The $1/m$ operators in Eq. (18), with coefficients
$V_{SK}(r)$ and $V_{SKb}(r)$ 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 in Eq. (19), the operators with
coefficients $V_{SLa}(r)$, $V_{S^{2}}(r)$, and $V_{S_{12}a}$ are the standard
spin-orbit, total spin squared, and tensor spin operators respectively, which
appear for standard quarkonia also. In addition to them, three novel operators
appear at order $1/m^{2}$. The operators with coefficients $V_{SLb}(r)$ and
$V_{SLc}(r)$ are generalizations of the spin-orbit operator to the hybrid
states. Similarly, the operator with coefficient $V_{S_{12}b}(r)$ is
generalization of the tensor spin operator to the hybrid states. Quite
interestingly, differently from the quarkonium case, the hybrid potential gets
a first contribution already at order $\Lambda^{2}_{\text{QCD}}/m$. The
corresponding operator does not contribute at LO to matrix elements of
quarkonium states as its projection on quark-antiquark color singlet states
vanishes. 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.

The coefficients $V_{i}(r)$ on the right-hand side of Eqs. (18) and (19) have
the form $V_{i}(r)=V_{oi}(r)+V_{i}^{np}(r)$, where $V_{oi}(r)$ is the
perturbative octet potential and $V_{i}^{np}(r)$ is the nonperturbative
contribution in the case in which $1/r\ll\Lambda_{QCD}$. The octet part can be
calculated in perturbation theory while the nonperturbative part is given in
terms of purely gluonic correlators whose detailed form is given in [53] (or
they are given in term of generalized Wilson loops with hybrids as initial and
final states in the case in which $1/r\sim\Lambda_{QCD}$). There are eight of
these nonperturbative correlators that have not yet been calculated on the
lattice, therefore we fix them on direct lattice determination of the
charmonium hybrid masses. In [53] we used data from the Hadron Spectrum
Collaboration one set from Ref. [50] with a pion mass of $m_{\pi}\approx 400$
MeV and a more recent set from Ref. [55] with a pion mass of $m_{\pi}\approx
240$ MeV. For all the details of the fit we refer to [53] . We take the values
$m^{RS}_{c}(1{\rm GeV})=1.477$ GeV and $\alpha_{s}$ at $4$-loops with three
light flavors, $\alpha_{s}(2.6\textrm{~{}GeV})=0.26$.

Figure 7: Spectrum of $H_{1}$ and $H3$ charmonium hybrid multiplets (including
the spin effects). The lattice results from [55] [2+1 flavors, $m_{\pi}=240$
MeV] are plotted in purple. In green we plotted the perturbative contributions
to the spin-dependent operators in Eq. (19) added to the spin average of the
lattice results (red dashed line). In blue we show the full result of the
spin-dependent operators of Eqs. (18)-(19) including both the perturbative and
nonperturbative contributions. The unknown nonperturbative matching
coefficients are determined by comparing the charmonium hybrid spectrum
obtained from the BOEFT to the lattice data. The height of the boxes indicate
the uncertainty as detailed in [53], there you find the spectrum for all the
multiplets.

Figure 8: Spectrum of $H_{2}$ and $H_{1}$ bottomonium hybrids multiplets,
nonperturbative correlators fixed on charmonium hybrids. In green we plotted
the perturbative contributions to the spin-dependent operators in Eq. (19)
added to the spin average of the lattice results (red dashed line). In blue we
show the full result of the spin-dependent operators of Eqs. (18)-(19)
including both the perturbative and nonperturbative contributions. The height
of the boxes indicate the uncertainty as detailed in [53], there you find the
spectrum for all the multiplets.

An interesting feature is that for the spin-triplets, the value of the
perturbative contributions decreases with $J$. This trend is opposite to that
of the lattice results. This discrepancy can be reconciled thanks to the
nonperturbative contributions, in particular due to the contribution from
$V^{np\,(0)}_{SK}$, which is only suppressed by $1/m$, and has no perturbative
counterpart. A consequence of this is that a model calculation with a spin
interaction inspired by quarkonium physics would give the wrong result.

All the dependence on the heavy-quark mass of the $V^{np\,(i)}$ is encoded in
the NRQCD matching coefficients $c_{F}$ and $c_{s}$. At leading order in
$\alpha_{s}$ these coefficients are known to be equal to $1$ and the
dependence on the heavy-quark mass only appears when the next-to-leading order
is considered. At the order we are working, only the heavy-quark mass
dependence of $c_{F}$ is relevant, we use its one-loop expression, with the
renormalization scale set as the heavy-quark mass. Taking this mass dependence
into account, we can use the set of nonperturbative parameters (that are
flavor independent) to predict the spin contributions in the bottomonium
hybrid sector, for which lattice determinations are difficult to obtain. We
present the results for two multiplets in Fig. 8, you can find the other
multiplets in [53].

This is a nice example of the huge predictivity that BOEFT plus lattice input
can achieve: once the nonperturbative correlators are fixed on lattice data
for charmonium hybrids (or better when these low energy correlators will be
directly calculated on the lattice) then we can predict all the bottomonia
(and $B_{c}$) hybrids multiplets.

####  3.1.3 Tetraquarks, Decays, Production

BOEFT may be used to describe also tetraquarks. In this case, it is necessary
to compute on the lattice generalized Wilson loops of the type given in Eq.
(11) where the operator $P^{a}_{\Gamma}$ contains also appropriate light quark
operators and besides the $\Lambda^{\sigma}_{\eta}$ and $\kappa$ also the
isospin quantum numbers $I=0,1$ have to be considered. The same approach can
be used also to describe pentaquarks [52, 58]. The light quark plays a similar
role to the gluonic excitations in the hybrid case. These static energies, for
which few pioneering lattice studies exist [61, 62, 63], are a crucial
ingredient to provide for the first time a dynamical description of
tetraquarks in QCD. When these objects will be made available, it will be
interesting to study the behaviour of the static energies near the avoided
level crossing with the heavy-light meson-meson thresholds [66]. Decays of
hybrids and tetraquarks as well as and mixing of hybrids and tetraquarks with
other quarkonium states could be studied and computed in BOEFT, however
lattice input on few correlators purely depending on the glue would be needed
[59].

###  3.2 Tetraquarks, molecules, hadroquarkonium ?

In this BOEFT/Lattice picture we will have all of the above descriptions
contributing [60]. When looking at the actual 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, threshold effects should have to be taken into account and a
molecular picture would emerge. However QCD would dictate, throught the
lattice correlators and the BOEFT characteristics and power counting, which
structure would dominate and in which precise way.

###  3.3 Input needed from the lattice

Summarizing, to address the X Y Z states we need input from the lattice:

  * •

Besides the hybrids static energies we need all the tetraquark static energies
in the short and long distance, extrapolated to continuum, calculated with the
generalized Wilson loop mentioned in 3.2. Some pioneering calculations exist
[61, 62, 63] but more should be done.

  * •

We need lattice calculations of the nonperturbative correlators entering the
spin dependent hybrids and tetraquark potentials: they are gauge invariant
correlators of chromolectric and chromomagnetic fields and gluelump operators
(without the color octet quark part).

  * •

We need lattice calculations of three point functions containing the
information of the quarkonium/hybrids quarkonium/tetraquarks mixing.

  * •

Lattice results about the crosstalk of the static energies with a pair of
heavy-light mesons in the lattice appeared [64, 65, 66] but further
investigations appear to be necessary.

This will be instrumental to tackle the X Y Z states in QCD.

##  4 Outlook

We have shown that there is a vital interplay between EFTs and lattice for
what concerns Exotics. NREFTs are vital to define objects of great
phenomenological interest and give a systematic scheme to calculate physical
observables, as well as an understanding of the underlying degrees of freedom
and dynamics.

NREFTs are based on factorization which allows to encode all the
nonperturbative effects in low energy correlators, depending only on
chromoelectric and chromomagnetic fields. This enhances the predictivity of
the theory since the same correlators can be used in the charmonium and
bottomonium sectors eliminating for example the difficulties of direct lattice
calculations of exotics at the high mass of the bottom. Moreover, once the
lattice has provided the nonperturbative input on the correlators, the
calculation of the physical properties ranging from masses, decays,
transitions and production, follow inside the BOEFT, which is much simpler and
under control that repeating a lattice calculation for each of these
observables. Of course there are challenges: the calculation of correlator of
chromoelectric and chromagnetic fields display bad convergence properties on
the lattice see e.g. [46] and references therein. The gradient flow method
seems however to help tremendously here, e.g. see [47, 67] There may be still
difficulties related to change of renormalization schemes between the
continuum of the BOEFT and the lattice regularization. For physical objects
like the force there are no issues and a gradient flow calculation in
continuum can help the lattice extrapolations [48].

We regard the combination of BOEFTs and Lattice QCD as very promising to
attack the physics of the X, Y, Z on the basis of QCD. We have seen that
Lattice is crucial to gain predictivity in these NREFTs calculations, on the
other hand the NREFTs allow to use lattice results in an array of problems
that often are not in its reach like nonequilibrium evolution, production
processes and in general to develop a convenient and complementary avenue to
address exotics.

## References

  * [1] N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan, Phys. Rept. 873 (2020), 1-154.
  * [2] 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 (2011), 1534.
  * [3] N. Brambilla et al. [Quarkonium Working Group], [arXiv:hep-ph/0412158 [hep-ph]].
  * [4] 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 (2014) no.10, 2981.
  * [5] M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98 (2018) no.3, 030001.
  * [6] S. K. Choi et al. [Belle], Phys. Rev. Lett. 91 (2003), 262001.
  * [7] A. Ali, J. S. Lange and S. Stone, Prog. Part. Nucl. Phys. 97 (2017), 123-198.
  * [8] 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.
  * [9] E. Braaten and M. Lu, Phys. Rev. D 76 (2007), 094028.
  * [10] E. Braaten and M. Kusunoki, Phys. Rev. D 69 (2004), 074005.
  * [11] S. Fleming and T. Mehen, Phys. Rev. D 85 (2012), 014016.
  * [12] E. Braaten, C. Langmack and D. H. Smith, Phys. Rev. D 90 (2014) no.1, 014044.
  * [13] E. Braaten, C. Langmack and D. H. Smith, Phys. Rev. Lett. 112 (2014), 222001.
  * [14] https://einrichtungen.ph.tum.de/T30f/tumqcd/index.html
  * [15] A. V. Manohar, Phys. Rev. D 56 (1997), 230-237.
  * [16] N. Brambilla, D. Gromes and A. Vairo, Phys. Lett. B 576 (2003), 314-327.
  * [17] W. E. Caswell and G. P. Lepage, Phys. Lett. B 167 (1986), 437-442.
  * [18] G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51 (1995), 1125-1171.
  * [19] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Rev. Mod. Phys. 77 (2005), 1423.
  * [20] A. Pineda and J. Soto, Nucl. Phys. B Proc. Suppl. 64 (1998), 428-432.
  * [21] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Nucl. Phys. B 566 (2000), 275.
  * [22] M. Foster et al. [UKQCD], Phys. Rev. D 59 (1999), 094509.
  * [23] K. J. Juge, J. Kuti and C. Morningstar, Phys. Rev. Lett. 90 (2003), 161601.
  * [24] G. S. Bali and A. Pineda, Phys. Rev. D 69 (2004), 094001.
  * [25] L. Müller, O. Philipsen, C. Reisinger and M. Wagner, Phys. Rev. D 100 (2019) no.5, 054503.
  * [26] C. Schlosser and M. Wagner, [arXiv:2111.00741 [hep-lat]].
  * [27] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Phys. Rev. D 63 (2001), 014023.
  * [28] A. Pineda and A. Vairo, Phys. Rev. D 63 (2001), 054007.
  * [29] N. Brambilla, V. Leino, P. Petreczky and A. Vairo, Phys. Rev. D 102 (2020) no.7, 074503.
  * [30] G. S. Bali, K. Schilling and A. Wachter, Phys. Rev. D 56 (1997), 2566-2589.
  * [31] G. S. Bali, Phys. Rept. 343 (2001), 1-136.
  * [32] Y. Koma, M. Koma and H. Wittig, Phys. Rev. Lett. 97 (2006), 122003.
  * [33] Y. Koma and M. Koma, PoS LATTICE2012 (2012), 140.
  * [34] Y. Koma and M. Koma, AIP Conf. Proc. 1322 (2010) no.1, 298-306.
  * [35] N. Brambilla, D. Eiras, A. Pineda, J. Soto and A. Vairo, Phys. Rev. D 67 (2003), 034018. [36]
  * [36] N. Brambilla, H. S. Chung and A. Vairo, Phys. Rev. Lett. 126 (2021) no.8, 082003.
  * [37] N. Brambilla, H. S. Chung and A. Vairo, JHEP 09 (2021), 032.
  * [38] N. Brambilla, PoS HardProbes2020 (2021), 017.
  * [39] D. Bala, O. Kaczmarek, R. Larsen, S. Mukherjee, G. Parkar, P. Petreczky, A. Rothkopf and J. H. Weber, [arXiv:2110.11659 [hep-lat]].
  * [40] N. Brambilla, M. A. Escobedo, A. Vairo and P. Vander Griend, Phys. Rev. D 100 (2019) no.5, 054025.
  * [41] N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, Phys. Rev. D 97 (2018) no.7, 074009.
  * [42] N. Brambilla, M. Á. Escobedo, M. Strickland, A. Vairo, P. Vander Griend and J. H. Weber, JHEP 05 (2021), 136.
  * [43] L. Altenkort, A. M. Eller, O. Kaczmarek, L. Mazur, G. D. Moore and H. T. Shu, Phys. Rev. D 103 (2021) no.1, 014511.
  * [44] G. Aarts, C. Allton, S. Kim, M. P. Lombardo, M. B. Oktay, S. M. Ryan, D. K. Sinclair and J. I. Skullerud, JHEP 11 (2011), 103.
  * [45] S. Kim, P. Petreczky and A. Rothkopf, JHEP 11 (2018), 088.
  * [46] N. Brambilla, V. Leino, O. Philipsen, C. Reisinger, A. Vairo and M. Wagner, [arXiv:2106.01794 [hep-lat]].
  * [47] V. Leino, N. Brambilla, O. Philipsen, C. Reisinger, A. Vairo and M. Wagner, [arXiv:2111.07916 [hep-lat]].
  * [48] N. Brambilla, H. S. Chung, A. Vairo and X. P. Wang, [arXiv:2111.07811 [hep-ph]].
  * [49] K. Marsh and R. Lewis, Phys. Rev. D 89 (2014) no.1, 014502 [arXiv:1309.1627 [hep-lat]].
  * [50] L. Liu et al. [Hadron Spectrum], JHEP 07 (2012), 126 [arXiv:1204.5425 [hep-ph]].
  * [51] M. Berwein, N. Brambilla, J. Tarrús Castellà and A. Vairo, Phys. Rev. D 92 (2015) no.11, 114019.
  * [52] N. Brambilla, G. Krein, J. Tarrús Castellà and A. Vairo, Phys. Rev. D 97 (2018) no.1, 016016.
  * [53] N. Brambilla, W. K. Lai, J. Segovia and J. Tarrús Castellà, Phys. Rev. D 101 (2020) no.5, 054040.
  * [54] N. Brambilla, W. K. Lai, J. Segovia, J. Tarrús Castellà and A. Vairo, Phys. Rev. D 99 (2019) no.1, 014017.
  * [55] G. K. C. Cheung et al. [Hadron Spectrum], JHEP 12 (2016), 089.
  * [56] S. M. Ryan et al. [Hadron Spectrum], JHEP 02 (2021), 214 [arXiv:2008.02656 [hep-lat]].
  * [57] G. Ray and C. McNeile, [arXiv:2110.14101 [hep-lat]].
  * [58] J. Soto and J. Tarrús Castellà, Phys. Rev. D 102 (2020) no.1, 014012.
  * [59] R. Oncala and J. Soto, Phys. Rev. D 96 (2017) no.1, 014004.
  * [60] J. Tarrús Castellà, EPJ Web Conf. 202 (2019), 01005.
  * [61] M. Sadl and S. Prelovsek, [arXiv:2109.08560 [hep-lat]].
  * [62] S. Prelovsek, H. Bahtiyar and J. Petkovic, Phys. Lett. B 805 (2020), 135467
  * [63] P. Bicudo, K. Cichy, A. Peters and M. Wagner, Phys. Rev. D 93 (2016) no.3, 034501
  * [64] G. S. Bali et al. [SESAM], Phys. Rev. D 71 (2005), 114513
  * [65] J. Bulava, B. Hörz, F. Knechtli, V. Koch, G. Moir, C. Morningstar and M. Peardon, Phys. Lett. B 793 (2019), 493-498
  * [66] R. Bruschini and P. González, [arXiv:2111.07653 [hep-ph]].
  * [67] V. Leino, N. Brambilla, J. Mayer-Steudte, A. Vairo "The static force from generalized Wilson loops using gradient flow", TUM-EFT 157/21