Glueball Spectrum with four light dynamical fermions

A. Athenodorou(a,b), J. Finkenrath(c,a), A. Lantos(d), M. Teper(e)

(a) Computation-based Science and Technology Research Center, The Cyprus
Institute, 20 Kavafi Str., Nicosia 2121, Cyprus
(b) Dipartimento di Fisica, Università di Pisa and INFN, Sezione di Pisa,
Largo Pontecorvo 3, 56127 Pisa, Italy
(c) Bergische Universität Wuppertal, Gaußstraße 20, Wuppertal, Germany.
(d) Cyprus University of Technology, Archiepiskopou Kyprianou 30, Limassol
3036, Cyprus
(e) All Souls College and Department of Physics, University of Oxford, Oxford
OX1 2JD, United Kingdom

We perform a calculation of the glueball spectrum for $N_{f}=4$ degenerate
dynamical fermions with masses corresponding to light pions. We do so by
making use of ensembles produced within the framework of maximally twisted
fermions by the Extended Twisted Mass Collaboration (ETMC). We obtain masses
of states that fall into the irreducible representations of the octahedral
group of rotations in combination with the quantum numbers of charge
conjugation $C$ and parity $P$; the above quantum numbers result in 20
distinct irreducible representations. We implement the Generalized Eigenvalue
Problem (GEVP) using a basis that consists only of gluonic operators. The
purpose of this work is to investigate the effect of light dynamical quarks on
the glueball spectrum and how this compares to the statistically more accurate
spectrum of $SU(3)$ pure gauge theory. Given that glueball states may have
broad widths and thus need to be disentangled from all the relevant mixings,
we use large ensembles of the order of ${\sim{~{}\cal O}}(20{\rm K})$
configurations. Despite the large ensembles, the statistical uncertainties
allow us to extract the masses for only a few irreducible representations;
namely $A_{1}^{++}$, $A_{1}^{-+}$, $E^{++}$ as well as $T_{2}^{++}$. The
results for the scalar $A_{1}^{++}$ representation show that an additional
state appears as the lightest state in the scalar $A_{1}^{++}$ channel of the
glueball spectrum, while the next two excited states are consistent with the
lightest two states of the pure gauge theory. To further elucidate the nature
of this additional state we perform a calculation using $N_{f}=2+1+1$
configurations and this demonstrates that it possesses a large quark content.
Finally, the ground states of the $E^{++}$ and $T_{2}^{++}$ tensor channels
and of the $A_{1}^{-+}$ pseudoscalar channel show, at most, minor effects due
to the inclusion of dynamical quarks.

##  1 Introduction

The extraction of the spectrum of glueballs in full QCD at physical quark
masses and in the continuum limit is an open question which requires a
coordinated effort from the Lattice community. A number of experiments such as
PANDA [1] and BESIII [2] are currently looking to identify such states, and so
the extraction of the associated spectrum from first principles as well as the
understanding of how light dynamical quarks affect the spectrum of glueballs
is a timely question that merits careful investigation. Recent reviews on
searches for glueballs can be found in the Lattice 2022 plenary presentation
by Davide Vadacchino [3, 4] as well as in the review in Ref. [5] by Eberhard
Klempt. Glueball states have large fluctuations and possibly broad resonance
widths and so the precise determination of their spectrum requires high
statistics, at least of the order of ${{\cal O}}(10-100{\rm K})$
configurations [6]. Hadron physics performed using current Lattice methodology
and machinery, on the other hand, makes use of sets of typically a few hundred
configurations. Nevertheless, for a number of measurements on the lattice such
as the extraction of renormalisation factors, large sets of configurations are
produced and can thus be used for exploratory investigations of the glueball
spectrum. This work results from the availability of accessing such large sets
of configurations.

In this work we invested effort in the investigation of the effect of light
dynamical quarks on the spectrum of glueballs. The spectrum of QCD in the
presence of light quarks is of interest for both phenomenological and
theoretical reasons. On the phenomenological side it would provide important
information regarding the mixing of ${q{\bar{q}}}$ states with glueball
states. This is expected to be modest by the OZI rule [7], which could mean
that, in practice, states with quark content could actually be undetectable in
the glueball spectrum, as calculated using purely gluonic correlators, or they
might appear with extremely small overlaps. For this purpose we use
configurations produced with $N_{f}=4$ light twisted mass quarks, since we
expect that such a set up will enhance the impact of light dynamical quarks on
the spectrum of glueballs.

We extract the glueball spectrum and then compare it with the one calculated
using pure gauge $SU(3)$ configurations [8]. Using decoupling arguments, we
expect that for massive dynamical quarks the glueball spectrum becomes similar
to that of the pure gauge theory. The important question which arises here is
what happens if one includes dynamical fermions with low masses close to
physical values.

Our investigation within $N_{f}=4$ QCD with four light degenerate fermions
reveals the presence in the scalar channel ($R^{PC}=A_{1}^{++}$) of an extra
state which is also the lightest state, in addition to excited states that
match the lightest and first excited states in the pure gauge theory. It is
possible that this additional state involves some mixing with $q{\bar{q}}$
states. To understand what is occurring in the low-lying spectrum of the
scalar channel and whether this additional state does indeed include a large
$q{\bar{q}}$ component we turn to some available $N_{f}=2+1+1$ ensembles.
Within this setup we extract the ground state energy for the scalar channel
and find that, in contrast to the first excited state, this indeed exhibits a
strong dependence on the pion mass.

Overall, our main findings can be summarized in the following points.

  1. 1.

In the scalar channel we obtain an additional state when we introduce light
dynamical quarks; this is the lightest state with the next two states having
masses consistent with the lightest two glueballs in the pure gauge theory.

  2. 2.

In the additional calculations that we perform using the $N_{f}=2+1+1$ gauge
ensembles, we observe that the additional ground state in the scalar channel
depends strongly on the pion mass, unlike the first excited state.

  3. 3.

The tensor glueball mass appears to be insensitive to the presence of light
dynamical quarks.

  4. 4.

The pseudoscalar glueball mass is affected only slightly by the inclusion of
light dynamical quarks and is close to the mass of the tensor glueball.

  5. 5.

The string tension is suppressed by the inclusion of four light dynamical
quarks by 50-70%.

This article is organized as follows: In Section 2 we provide a detailed
presentation of the lattice set up which has been used for the production of
the configurations with $N_{f}=4$ and $N_{f}=2+1+1$ twisted mass fermions as
well as those with the $SU(3)$ pure gauge action. Then, in Section 3 we
explain how one can extract the spectrum of glueballs in Lattice QCD by making
use of the Generalized Eigenvalue Problem. In the same manner, in Section 4 we
describe how to evaluate the masses of torelons. Following that, in Section 5,
we describe the calculation of the topological charge which is used as a
measure of the ergodicity of the system. Furthermore we explain how we
evaluate the energy scale $t_{0}$ using the gradient flow. Subsequently in
Section 6, we move to the presentation of the results by focussing on the
scalar glueballs in the $R^{PC}=A_{1}^{++}$ channel, the lightest tensor
glueball whose components are split between the $R^{PC}=E^{++}$ and the
$T_{2}^{++}$ channels, as well as the pseudoscalar glueball obtained in the
$R^{PC}=A_{1}^{-+}$ channel. Then, in Section 7 we focus our discussion on the
ground state obtained in the scalar channel for $N_{f}=2+1+1$, where it
appears to display a dominant $q{\bar{q}}$ component. Subsequently, in Section
8, we present our results for the string tension obtained from the energies of
the torelons. We then compare our glueball results with other older as well as
recent results from the literature, before presenting our conclusions in
Section 10.

##  2 Lattice Set-up

For the evaluation of the glueball spectrum we use gauge ensembles of clover
improved twisted mass fermions produced with $4$ degenerate light flavours
($N_{f}=4$) at two different lattice spacings as well as two ensembles with
$2$ degenerate light flavours and a strange and a charm quark ($N_{f}=2+1+1$).
All the configurations with dynamical quarks have been generated within the
context of the Extended Twisted Mass collaboration. They were also used in a
parallel work for extracting renormalisation constants [9].

For the gluonic part of the action we use the Iwasaki improved gauge action
[10, 11, 12], which is given by the expression

$S_{G}=\frac{\beta}{3}\sum_{x}\left(c_{0}\sum_{\begin{subarray}{c}\mu,\nu=1\\\
\mu<\nu\end{subarray}}^{4}\left[1-\mathrm{Re}\,{\rm
Tr}\left(U_{x,\mu\nu}^{1\times
1}\right)\right]+c_{1}\sum_{\begin{subarray}{c}\mu,\nu=1\\\
\mu\neq\nu\end{subarray}}^{4}\left[1-\mathrm{Re}\,{\rm
Tr}\left(U_{x,\mu\nu}^{1\times 2}\right)\right]\right)\,,$ (1)

where $\beta=6/g^{2}$, $U^{1\times 1}$ is a plaquette and $U^{1\times 2}$ is a
rectangular Wilson loop. The Symanzik coefficients are $c_{0}=3.648$ and
$c_{1}=(1-c_{0})/8$. The fermionic action is given [13, 14] by:

$S_{F}^{l}=a^{4}\sum_{x}\bar{\chi}^{(l)}(x)\left(D_{W}[U]+\frac{i}{4}c_{SW}\sigma^{\mu\nu}\mathcal{F}^{\mu\nu}[U]+m_{0,l}+i\mu_{l}\gamma_{5}\tau^{3}\right)\chi^{(l)}(x)\,.$
(2)

In the equation above, $\chi^{(l)}$ is the field representing the quark
doublets, expressed in the twisted basis, $m_{0,l}$ and $\mu_{l}$ are
respectively the untwisted and twisted mass parameters, $\tau^{3}$ is the
third Pauli matrix acting in flavor space and $D_{W}$ is the massless Wilson-
Dirac operator. The clover term $\propto\sigma^{\mu\nu}\mathcal{F}^{\mu\nu}$
is included in the action to suppress cut-off effects, reducing the difference
between the mass of the charged and neutral pions [15].

Moving on to the $N_{f}=2+1+1$ case, the additional strange and charm quarks
are included as a non-degenerate twisted doublet $\chi^{(h)}=(s,c)^{t}$, with
the action [16]

$S_{F}^{h}=a^{4}\sum_{x}\bar{\chi}^{(h)}(x)\left(D_{W}[U]+\frac{i}{4}c_{SW}\sigma^{\mu\nu}\mathcal{F}^{\mu\nu}[U]+m_{0,h}-\mu_{\delta}\tau^{1}+i\mu_{\sigma}\gamma_{5}\tau^{3}\right)\chi^{(h)}(x)\,,$
(3)

where $m_{0,h}$ is the bare untwisted quark mass for the heavy doublet,
$\mu_{\delta}$ the bare twisted mass along the $\tau^{1}$ direction and
$\mu_{\sigma}$ the mass splitting in the $\tau^{3}$ direction.

The partial conserved axial current (PCAC) mass is tuned to zero in order to
achieve maximal twist. This ensures automatic $\mathcal{O}(a)$ improvement for
the expectation values of the observables of interest [17]. The simulation
parameters of the gauge ensembles with $N_{f}=4$ as well as $N_{f}=2+1+1$
quarks are given in Tables 1 and 2 respectively.

| $\beta$ | $\quad c_{SW}$ | $\mu_{l}$ | $L$ | $am_{PS}$ | $t_{0}/a^{2}$
---|---|---|---|---|---|---
cB4.06.16 | $1.778$ | $1.69$ | $0.006$ | 16 | $0.2652(53)$ | $4.947(62)$
cB4.06.24 | $1.778$ | $1.69$ | $0.006$ | 24 | $0.1580(\phantom{0}8)$ | $4.667(17)$
cC4.05.24 | $1.836$ | $1.6452$ | $0.005$ | 24 | $0.1546(20)$ | $6.422(48)$
Table 1:  Simulation parameters of the $N_{f}=4$ gauge ensembles [9] used in this work. | $\beta$ | $\quad c_{SW}$ | $\mu_{l}$ | $L$ | $am_{PS}$ | $t_{0}/a^{2}$
---|---|---|---|---|---|---
cA211.53.24 | $1.726$ | $1.74$ | $0.0053$ | 24 | $0.1661(\phantom{0}4)$ | $2.342(\phantom{0}6)$
cA211.25.32 | $1.726$ | $1.74$ | $0.0025$ | 32 | $0.1253(\phantom{0}1)$ | $2.392(\phantom{0}4)$
Table 2:  Simulation parameters of the $N_{f}=2+1+1$ gauge ensembles [15, 18]
used in this work.

For the purpose of comparison with the $SU(3)$ pure gauge theory, we also
simulate that theory using the standard Wilson action which is given by Eq. 1
with Symanzik coefficients $c_{0}=1$ and $c_{1}=0$. The simulation algorithm
combines standard heat-bath and over-relaxation steps in the ratio 1:4; these
are implemented by updating $SU(2)$ subgroups using the Cabibbo-Marinari
algorithm [19]. The parameters of the pure gauge runs are shown in Table 3.

$\beta$ | $L$ | $\sqrt{\sigma}$ | $am_{G}$ | $t_{0}/a^{2}$
---|---|---|---|---
$6.222$ | $30$ | $0.1533(6)$ | $0.499(6)$ | $6.422$
$6.135$ | $26$ | $0.1750(9)$ | $0.578(4)$ | $4.947$
$6.117$ | $26$ | $0.1781(12)$ | $0.585(8)$ | $4.667$
Table 3: Simulation parameters for the $SU(3)$ pure gauge ensembles used in
this work. The last column represents the values of $t_{0}/a^{2}$ for
$N_{f}=4$ simulations corresponding to the given values of $\beta$ via a cubic
spline interpolation of data from Refs. [20, 21, 22, 23] and presented in
Figure 1. Figure 1: The cubic spline interpolation between $t_{0}/a^{2}$ and
$\beta$ extracted from results taken from Refs [20, 21, 22, 23].

The idea behind choosing the particular parameters is based on matching the
gradient flow times $t_{0}/a^{2}$ between the $N_{f}=4$ QCD and pure gauge
runs. Since $N_{f}=4$ QCD configurations had already been produced, we
calculate the corresponding values for $t_{0}/a^{2}$ and subsequently match
them to the right value of $\beta$ for the pure gauge theory using cubic
spline interpolations based on data from Refs [20, 21, 22, 23]. The curve
obtained from the cubic spline interpolation is given in Figure 1. The
motivation for this procedure is discussed in more depth in Section 6.

##  3 Calculation of Glueball Masses

Glueballs are colour singlet states and their masses can be determined by
employing the standard decomposition of a Euclidean correlator associated with
an operator denoted as $\phi(t)$. This correlator has an overlap with the
physical states expressed in the form of the eigenstates of the system’s
Hamiltonian, represented as $H$:

$\displaystyle\langle\phi^{\dagger}(t=an_{t})\phi(0)\rangle$ $\displaystyle=$
$\displaystyle\langle\phi^{\dagger}e^{-Han_{t}}\phi\rangle=\sum_{i}|c_{i}|^{2}e^{-aE_{i}n_{t}}$
(4) $\displaystyle\stackrel{{\scriptstyle t\to\infty}}{{=}}$
$\displaystyle|c_{0}|^{2}e^{-aE_{0}n_{t}}\,,$

where the energy levels are arranged in ascending order as $E_{i+1}\geq
E_{i}$, where $E_{0}$ represents the ground state energy. Only states with
non-zero overlaps, denoted as $c_{i}=\langle{\rm
vac}|\phi^{\dagger}|i\rangle\neq 0$, contribute to the summation above. Thus,
it is necessary to ensure that the quantum properties of the operator $\phi$
match those of the states we are interested in. In our study, our primary
focus is on glueballs and torelons, so we must embed the appropriate quantum
numbers within the operator $\phi$ to project onto these specific states.
These quantum numbers include the angular momentum-spin $J$, Parity-$P$, and
Charge conjugation-$C$, which define the properties of glueball states.

The determination of the ground state hinges on two key factors: the quality
of the overlap with this state and the rapidity of exponential decay with
respect to the time variable $t$, as stipulated by Eq.(4). Maximizing the
overlap involves constructing operators that effectively encapsulate the
essential properties of the target state. Achieving this entails projecting
onto the appropriate quantum characteristics and physical length scales
associated with the relevant state. To expedite the onset of exponential decay
with respect to $t$, it becomes imperative to minimize the influence of higher
excited states. To address this, we employ the Generalized Eigenvalue Problem
(GEVP)[24, 25], which is applied to a set of operators constructed using
various lattice loops at different blocking levels [26, 27]. This approach
mitigates contamination from excited states, enhancing the operators’
alignment with the physical length scales.

The glueballs are color singlets, and thus, an operator designed to represent
a glueball state is created by taking the trace of an ordered product of
$SU(3)$ link matrices along a contractible loop. To ensure that the
correlators maintain exact positivity, we restrict the use of loops to those
containing only spatial links. The real part of the trace corresponds to
Charge conjugation $C=+$, while the imaginary part corresponds to $C=-$. To
generate an operator with zero momentum ($p=0$), we sum over all spatial
translations of the loop. By considering all possible rotations of the loop,
we construct linear combinations that conform to the irreducible
representations $R$ of the rotational symmetry group associated with our cubic
spatial lattice. We consistently opt for a cubic spatial lattice volume
($L_{x}=L_{y}=L_{z}$) that respects these symmetries. For each loop, we also
create its parity-inverse counterpart, allowing us to form linear combinations
and construct operators that encompass both parities, denoted as $P=\pm$. The
correlators of these operators will project onto glueballs with $p=0$ and
possess the quantum numbers $R^{PC}$ as defined by the relevant operators. In
total, there are 12 paths used in constructing the glueball operators, which
are illustrated in Figure 2.

The irreducible representations denoted as $R$ within our subgroup of the
complete rotation group are conventionally labeled as $A_{1},A_{2},E,T_{1},$
and $T_{2}$. $A_{1}$ representation is a singlet and exhibits rotational
symmetry, thus, containing the $J=0$ state in the continuum limit. Similarly,
the $A_{2}$ representation is also a singlet, while the $E$ representation
forms a doublet, and both $T_{1}$ and $T_{2}$ are triplets. For instance, when
considering the three states that transform under the triplet of $T_{2}$, they
are degenerate on the lattice, so we average their values and treat them as a
single state in our calculations of glueball masses. The same approach is
applied to the $E$ doublets to maintain consistency in our estimates.

Glueball energy values are extracted through the utilization of correlation
matrices denoted as $C_{ij}=\langle\phi_{i}^{\dagger}(t)\phi_{j}(0)\rangle$,
where $i,j=1...N_{\rm op}$ and $N_{\rm op}$ representing the number of
operators. In the scalar channel $A_{1}^{++}$, there exists a non-zero
projection onto the vacuum. In such cases, it can be convenient to employ the
vacuum-subtracted operator $\phi_{i}-\langle\phi_{i}\rangle$. This subtraction
effectively eliminates the vacuum contribution from Equation 4, ensuring that
the lightest non-trivial state appearing in the summation becomes the dominant
term in the state expansion.

The representations of the rotational symmetry group described above are
specific to our cubic lattice formulation. As we approach the continuum limit,
these states will converge towards the continuum glueball states, which are
categorized into representations of the continuum rotational symmetry group.
In simpler terms, they fall into degenerate multiplets of states with $2J+1$
states each. When determining the continuum limit of the low-lying glueball
spectrum, it becomes more meaningful to assign the states to a particular spin
value $J$ rather than to the representations of the cubic subgroup, which
provide less detailed information, as they map all spins from
$J=1,2,3,\dots,\infty$ to just five cubic representations. The distribution of
$2J+1$ states for a given $J$ among the representations of the cubic symmetry
subgroup is outlined in Table 4 for relevant low values of $J.$ For example,
the five states corresponding to a $J=2$ glueball will be divided between a
degenerate doublet $E$ and a degenerate triplet $T_{2}$, resulting in a total
of five states. These $E$ and $T_{2}$ states will be affected by lattice
spacing corrections of the order $O(a^{2})$, but as the lattice spacing $a$
becomes sufficiently small, these states will become nearly degenerate. This
near-degeneracy can then be utilized to identify the continuum spin of the
states.

$J$ | $A_{1}$ | $A_{2}$ | $E$ | $T_{1}$ | $T_{2}$
---|---|---|---|---|---
0 | 1 | 0 | 0 | 0 | 0
1 | 0 | 0 | 0 | 1 | 0
2 | 0 | 0 | 1 | 0 | 1
3 | 0 | 1 | 0 | 1 | 1
4 | 1 | 0 | 1 | 1 | 1
5 | 0 | 0 | 1 | 2 | 1
6 | 1 | 1 | 1 | 1 | 2
7 | 0 | 1 | 1 | 2 | 2
8 | 1 | 0 | 2 | 2 | 2
Table 4: Subduced representations of continuum spin $J\downarrow R$ of the
octahedral group of rotations up to $J=8$, illustrating the spin content of
the representations $R$ in terms of the continuum $J$. Figure 2: All the
different closed loops used for the construction of the glueball operators.

##  4 Calculation of torelon masses

In the same manner as in the glueball mass calculation, we can extract the
mass of the closed spatial Polyakov loop, also called a torelon, which winds
around a spatial compactified direction. The quantum numbers which
characterize the torelon states are the angular momentum-spin $J$, Transverse
Parity-$P_{\perp}$ as well as the Longitudinal Parity-$P_{||}$. The ground
state of the torelon which provides the easiest way to extract the string
tension is described by $J^{P_{\perp},P_{||}}$=$0^{++}$. For the purposes of
this work we will focus on the $0^{++}$channel, leaving the extraction of the
torelon mass for several configurations of the quantum numbers
$J^{P_{\perp},P_{||}}$ a subject of future investigations.

We obtain the mass of the torelon by calculating the ground state energy
$m_{T}(L)$ of a flux-tube of length $L$ that closes on itself by winding once
around the spatial compactified torus. We use Eq. (4) where the operator
$\phi$ is the product of $SU(3)$ link matrices taken around a non-contractible
closed path that winds once around the spatial torus. The simplest such
operator is the elementary Polyakov loop:

$\phi(n_{t})=P(n_{t})=\sum_{n_{y},n_{z}}\mathrm{Tr}\left\\{\prod^{L_{x}}_{n_{x}=1}U_{x}(n_{x},n_{y},n_{z},n_{t})\right\\}\,,$
(5)

and, thus, has $0^{++}$ quantum numbers.

The above formula denotes the path ordered product of link matrices in the
$x$-direction winding once around the $x$-torus. Then we sum over translations
along the $x$-torus and a time slice so that we project onto zero longitudinal
as well as transverse momentum respectively i.e.
$(p_{x},p_{y},p_{z})=(0,0,0)$. Confining torelons, at least in pure gauge
theories, behave, to an adequate extent, as bosonic strings [28]. So far no
such observation has been reported in full QCD calculations for the torelon.
One would expect that the closed flux-tube behaves as a confining bosonic
string, the spectrum of which can be approximated to an adequate accuracy by
the Nambu-Goto string [29]. The ground state energy of the bosonic string is
given by the expression

$\displaystyle E(l)=\sigma l\sqrt{1-\frac{2\pi}{3\sigma l^{2}}}\,,$ (6)

where $\sigma$ the string tension and $l=aL_{x}$ the length of the torelon. By
imposing non-zero longitudinal momentum, the energy of the torelon is elevated
by a momentum factor as well as by a Kaluza-Klein mode due to the difference
between the total contribution of the left and right moving phonons, while by
imposing transverse momentum the energy level of the torelon is altered
according to the associated lattice dispersion behaviour. Hence, for the
purposes of this work we focus on zero longitudinal and transverse momenta.

As explained above, the operator of Eq. 5 is invariant under rotations about
its main axis and so has angular momentum $J=0$. It is also clearly invariant
under both transverse and longitudinal parity transformations. Therefore, this
operator is ideal for projecting onto the torelon absolute ground state
$0^{++}$. Once more we employ smearing and blocking techniques in order to
enhance the projection onto the physical states. In $SU(N)$ pure gauge
theories, torelon operators are not invariant under center symmetry
transformations and, thus, their vacuum expectation value is zero as long as
this is not broken. However, in full Lattice QCD, the action is not invariant
under center symmetry transformation and thus the vacuum expectation value
should be subtracted from the correlator.

In the presence of a quark-dynamical background, flux-tubes should decay to
pions as long as this is energetically favorable. That is to say, if the flux-
tube ground energy is equal or larger than twice the mass of the pion, it
should decay to two pions leading to a constant ground energy while the
linearly rising contribution ”jumps up” to the next excitation level. Of
course such a phenomenon has never been observed in the closed flux-tube
channel and would be extremely interesting to test this scenario.

##  5 Topological charge and scale setting

In the continuum theory, the topological charge is calculated as the integral
across the four-dimensional volume of the topological charge density,
specifically:

${\cal Q}=\frac{1}{32\pi^{2}}\int d^{4}x\>\epsilon_{\mu\nu\rho\sigma}{\rm
Tr}\left[F_{\mu\nu}(x)F_{\rho\sigma}(x)\right]\,.$ (7)

The discrete equivalent of this quantity can be derived by substituting the
gluonic field tensor with a lattice operator, ensuring that it converges
correctly to the continuum limit. This selection is not exclusive, and
operators with reduced discretization effects can be obtained by employing
${\cal O}(a)$-improved definitions of $F_{\mu\nu}$. Specifically, the ${\cal
Q}$ definition utilized in our study is the symmetric or ’clover’ definition,
initially introduced in Ref. [30]. It takes the following form:

$\mathcal{Q}_{L}=\frac{1}{32\pi^{2}}\sum_{x}\epsilon_{\mu\nu\rho\sigma}\text{Tr}\left[C_{\mu\nu}(x)C_{\rho\sigma}(x)\right]\,,$
(8)

and uses a discretization of the gauge strength tensor in terms of a ”clover
leaf” path $C_{\mu\nu}$, made by the sum of the plaquettes $P_{\mu\nu}(x)$
centered in $x$ and with all the possible orientations in the $\mu\nu$-plane,
i.e.

$C_{\mu\nu}(x)=\frac{1}{4}\text{Im}\left[P_{\mu\nu}(x)+P_{\nu,-\mu}(x)+P_{-\mu,-\nu}(x)+P_{-\nu,\mu}(x)\right]\,.$
(9)

This operator retains its even parity under transformations and displays
discretization effects of order $\mathcal{O}(a^{2})$. To smooth the
ultraviolet fluctuations of the gauge field and extract the topological
charge, we employ the gradient flow method as described in [20]. We utilize
the standard Wilson action for the smoothing action in the flow equation, with
an elementary integration step of $\epsilon=0.01$. The computation of the
topological charge is performed on the smoothed fields at intervals of
$\Delta\tau_{\rm flow}=0.1$. Selecting an appropriate flow time is crucial. It
should be large enough to cancel out discretization effects but should also
maintain the gauge field’s topological properties. Therefore, we investigate
the dependence of our final results on $\tau_{\rm flow}$ in search of a
plateau region. According to Ref. [20], we anticipate that this plateau
emerges when $a\sqrt{8\tau_{\rm flow}}\sim\mathcal{O}(0.1\text{fm})$.

In Figure 3, we illustrate the evolution of the topological charge and its
distribution across the three $N_{f}=4$ ensembles. The plots conspicuously do
not exhibit signs of topological freezing, strongly suggesting the ergodicity
of the simulated Markov-Chain. The topological charge distribution appears to
have a Gaussian-like shape, characterized by a relatively narrow spread for
the smaller physical volumes of the cB4.06.16 and cC4.05.24 ensembles, while
displaying a much broader distribution for the larger physical volume of
cB4.06.24. The presence of a topological charge distribution that closely
resembles a Gaussian distribution implies that our computations have
effectively probed the topological sectors of the theory. Similar behavior is
observed across all other ensembles, both in the context of pure gauge theory
and $N_{f}=2+1+1$.

The gradient flow technique also allows for the establishment of a precise
physical scale parameter denoted as $t_{0}$. The concept of this flow
observable was initially introduced in [31, 32]. The determination of the
gradient flow time $t_{0}$ follows this prescription: First, we establish

$\displaystyle F(t)=t^{2}\langle E(t)\rangle\,\ {\rm where}\
E(t)=\frac{1}{4}B^{2}_{\mu\nu}(t)\,,$ (10)

where the field strength $B_{\mu\nu}$ is obtained by evolving $F_{\mu\nu}$
along the direction of flow time. We establish the scale $t_{0}(c)$ by
determining the value of $t$ for which $F(t)|_{t=t_{0}(c)}=c$. The choice of
$c$ is critical and should satisfy the condition $a\ll\sqrt{8t_{0}}\ll L$,
where selecting a small $c$ can result in larger lattice artifacts, while
opting for a larger $c$ generally leads to increased autocorrelations, as
noted in [33]. In our case, we opt for $c=0.3$, a commonly used value in
lattice QCD calculations.

Figure 4 displays both the evolution of the gradient flow time $t_{0}/a^{2}$
and its corresponding distribution across the three $N_{f}=4$ ensembles. The
history plots reveal substantial fluctuations in $t_{0}$, while the
distribution notably approximates a Gaussian shape to a considerable extent.

Figure 3: The history of the topological charge as well as its distribution
for the three $N_{f}=4$ ensembles at $t=t_{0}$. Left column shows the history
while right column corresponds to the histogram.

Figure 4: The history of the gradient flow observable $t_{0}/a^{2}$ for the
three $N_{f}=4$ ensembles used for the purposes of this work. On the left we
demonstrate the actual measurements while on the right the corresponding
histogram.

##  6 Results

###  6.1 The spectrum of Glueballs in $N_{f}=4$

We have computed the low-lying spectrum for specific irreducible
representations, namely, $A_{1}^{++}$, $E^{++}$, $T_{2}^{++}$, and
$A_{1}^{-+}$. These representations correspond to the scalar ($A_{1}^{++}$),
tensor ($E^{++}+T_{2}^{++}$), and pseudoscalar ($A_{1}^{-+}$) channels. In
Figure 5, we present the effective mass plateaus for two of these irreducible
representations, namely, $A_{1}^{++}$ and $T_{2}^{++}$, specifically for the
ensemble labeled as cB4.06.24, which serves as our reference ensemble. One
notable observation from the calculation is that the effective mass plateaus
set relatively late in time compared to what one would typically observe in
pure gauge theory. This, results in an overlap of approximately 30 % to 50 %,
which is significantly lower compared to the rapid convergence seen in $SU(3)$
pure gauge theory, where overlaps often reach values of around 90 % to 100 %.
This behavior may be explained by the richness of the Hilbert space in the
presence of dynamical quarks, leading to mixings among our operators,
including states involving mesons. In an attempt to improve the overlaps with
the extracted states, we expanded the variational basis of operators by
introducing additional paths, but this led to only marginal improvements.

Figure 5:  Left panel: The effective masses corresponding to the ground and
first excited states within the irreducible representation $A_{1}^{++}$ for
the ensemble cB4.06.24. Right panel: The effective mass corresponding to the
ground state of the irreducible representation $T_{2}^{++}$ for the ensemble
cB4.06.24; states in $T_{2}$ come in triplets, thus, the three states that
appear as ground, first excited and second excited states within the
variational analysis.

In Figure 6, we display the results for glueball masses in $N_{f}=4$ for three
different scenarios: ${\beta=1.778,L=16}$, ${\beta=1.778,L=24}$, and
${\beta=1.836,L=24}$, presented in the left, middle, and right panels,
respectively. Additionally, we provide a comparison with the corresponding
values extracted from the pure gauge theory. These plotted data are also
tabulated in Tables 5, 6, 7 and 8 for the representations $A_{1}^{++}$,
$E^{++}$, $T_{2}^{++}$ and $A_{1}^{-+}$ respectively.

The middle panel of Figure 6, or Figure 7, serves as our reference plot
because we know that finite volume effects do not significantly impact the
pion mass for this specific ensemble. Nevertheless, our findings indicate that
glueball masses appear to be relatively insensitive to finite volume effects.
This observation becomes evident when comparing the middle panel, where the
spatial lattice size is $L=24$, with the left panel, where $L=16$, while
keeping the value of $\beta$ constant. However, it’s worth noting that
reducing the volume can introduce other finite volume effects, such as non-
glueball states known as di-torelons.

Indeed, in the left panel of Figure 6, it is apparent that for the
representation $E^{++}$, the ground state seems to be a di-torelon with a mass
slightly greater than twice that of the torelon, as corroborated by the data
in Tables 6 and 9. It’s important to note that this is a somewhat heuristic
way of identifying the nature of this state as a di-torelon. A more rigorous
approach to confirming whether this state is indeed a di-torelon involves the
following procedure. Initially, we expand our variational basis to include
operators deliberately designed to maximize their overlap with di-torelon
states. These operators are constructed by taking the product of a spatial
Polyakov loop with its charge conjugate, ensuring that the two lines have
opposite directions. We then form linear combinations of these spatial double
lines to encode the correct quantum numbers for the $E^{++}$ representation.
At first, we conduct the analysis using the complete basis of operators, which
includes the di-torelon operators. Surprisingly, the resulting spectrum shows
no discrepancies when compared to the spectrum obtained using the basis that
excludes the di-torelon operators. Subsequently, we carry out a similar
analysis using only the di-torelon operators. In this case, the spectrum
obtained confirms the existence of the suspected di-torelon states, but
intriguingly, the lightest state, typically identified as a glueball, does not
appear. As anticipated, as we increase the lattice’s spatial length from
$L=16$ to $L=24$, the di-torelon state vanishes.

Furthermore, we conducted an examination to determine whether decreasing the
lattice spacing would have an impact on the spectrum of $N_{f}=4$. This
investigation demonstrated that, particularly for the lower-energy spectrum,
the spectrum depicted in the middle plot of Figure 6 accurately represents the
continuum physics. This assertion becomes evident when comparing the spectrum
in the middle panel to that in the right panel of Figure 6, where we vary the
value of $\beta$.

$M\sqrt{t_{0}}$

Figure 6: The spectrum of glueballs for the representations $A_{1}^{++}$,
$E^{++}$, $T_{2}^{++}$, $A_{1}^{-+}$ for the ensembles cB4.06.16(left panel),
cB4.06.24(middle panel) and cC4.05.24(right panel). In purple we denote the
states of $N_{f}=4$ QCD while in green the states of $SU(3)$ pure gauge.

Let’s now focus to the comparison between $N_{f}=4$ QCD and pure $SU(3)$ gauge
theory. To achieve this, we begin by determining the value of $t_{0}/a^{2}$
for the $N_{f}=4$ configurations. Following this, we proceed to extract the
glueball spectrum for the pure gauge theory at a specific $\beta$ value that
corresponds to the same $t_{0}/a^{2}$. You can find the values of
$t_{0}/a^{2}$ in Table 2, while the corresponding $\beta$ values have been
determined through interpolations of data from references such as Refs [20,
21, 22, 23], and they are presented in Table 3 and graphically represented in
Figure 1.

In Figure 6, the glueball masses for pure $SU(3)$ gauge theory are depicted in
green, while those for $N_{f}=4$ are represented in purple. Notably, a
distinguishing feature is the presence of an additional state within the
spectrum of states for the $A_{1}^{++}$ channel. Specifically, the ground
state of the $SU(3)$ gauge theory in the $A_{1}^{++}$ channel seems to differ
from that of $N_{f}=4$ and instead aligns more closely with the first excited
state in the $N_{f}=4$, $A_{1}^{++}$ channel. This observation immediately
raises the hypothesis that the ground state in $N_{f}=4$ incorporates some
quark content, while the first excited state corresponds to a glueball state.
However, it’s important to note that this hypothesis requires further
investigation.

Shifting our attention to the ground states of $E^{++}$ and $T_{2}^{++}$,
which correspond to the ground state of $J^{PC}=2^{++}$, it’s noteworthy that
there is a consistent agreement between $N_{f}=4$ and pure gauge theory. This
suggests that conducting precise investigations of the tensor glueball in the
pure gauge context can serve as a reliable benchmark for comparisons with
experimental data. Lastly, similar trends have emerged for the pseudoscalar
ground state $A_{1}^{-+}$. In particular, it appears to exhibit a slightly
lower mass when dynamical quarks are included in the simulation, with a
deviation of one sigma. Additionally, the ground state of $2^{++}$ in the
$N_{f}=4$ context closely resembles that of $A_{1}^{-+}$, mirroring the
situation observed in pure gauge theory.

$M\sqrt{t_{0}}$$m_{\pi}\sqrt{t_{0}}$$2\times m_{\pi}\sqrt{t_{0}}$$4\times
m_{\pi}\sqrt{t_{0}}$

Figure 7: The spectrum of glueballs for the representations $A_{1}^{++}$,
$E^{++}$, $T_{2}^{++}$, $A_{1}^{-+}$. In purple we denote the states of
$N_{f}=4$ QCD for the ensemble cB4.06.24 while in green the states of $SU(3)$
pure gauge for $\beta=6.117$. The dashed lines correspond to 1, 2 and 4 times
the pion mass from bottom to top respectively.

##  7 On the content of the $A_{1}^{++}$ ground state

One of the most intriguing aspects of the obtained spectra is the appearance
of an additional state in the scalar channel, as indicated in Figure 6. To
gain a deeper understanding of the nature of this extra state within the
$A_{1}^{++}$ irreducible channel, we conducted a similar investigation for
$N_{f}=2+1+1$ configurations. Specifically, we extracted the low-lying
spectrum for $A_{1}^{++}$ under the more physically relevant conditions of
$N_{f}=2+1+1$ twisted mass fermions, considering two distinct values of pion
masses, approximately $260$ and $350$ MeV. Once again, our analysis involved
ensembles composed of a substantial number of configurations, although not as
numerous as those for $N_{f}=4$. For this investigation, we deemed the
ensembles cA211.53.24 (comprising 5000 configurations) and cA211.25.32
(consisting of 2500 configurations) to be suitable choices. You can find
detailed information about these two ensembles in Table 2.

Figure 8: The ground state and the first excited state of the scalar channel
$A_{1}^{++}$ for the ensembles cA211.53.24 and cA211.25.32. The blue band
represents the scalar glueball mass of the $SU(3)$ pure gauge theory at
$\beta=5.917$ $(t_{0}/a^{2}\sim 2.38)$ extracted for the purposes of this
work. The black dashed lines correspond to twice the pion mass.

Figure 8 illustrates the outcomes for both the ground and the first excited
states in the $A_{1}^{++}$ scalar channel for the $N_{f}=2+1+1$
configurations. In this visualization, the blue band represents the scalar
glueball mass in pure $SU(3)$ gauge theory at $\beta=5.917$
($t_{0}/a^{2}\approx 2.38$), which has been determined specifically for this
study. Meanwhile, the dashed black lines indicate approximately twice the mass
of the pion for each distinct ensemble.

The message conveyed by this plot is quite evident: the behavior of the ground
state strongly correlates with the pion mass, while the first excited state
remains relatively consistent. This strongly suggests that the ground state
contains quark content, possibly signifying a decay process involving a
glueball decaying into either two or four pions. Notably, the mass of the
ground state closely aligns with twice the mass of the pion (as indicated by
the dashed lines). In contrast, the first excited state, which appears to
remain stable, is indicative of representing the glueball mass. It’s important
to acknowledge that the coarseness of the lattice introduces substantial
influences from lattice artifacts, leading to a considerable deviation from
the pure gauge theory value.

##  8 The Torelon Ground State and the string tension

We have extracted the ground state of the spectrum of the torelon for the
three $N_{f}=4$ ensembles. The torelons are short enough, enabling us to
extract adequate results with good statistics. The effective mass plots
resulting from the torelon correlators are shown in Figure 9. As in the case
of glueballs, it is stiking that the mass plateau sets late in time, resulting
to a low overlap of $20\%$; nevertheless the plateaus are quite convincing,
leading to a precise estimation of the torelon masses. The closed flux-tube
masses are presented in the second column of Tab. 9. Assuming that the torelon
behaves approximately as a bosonic string with energy described by Eq. 6, we
extracted the associated string tension which is presented in the third column
of Tab. 9. The string tension for $L=16$ and $L=24$ and $\beta=1.778$ is
compatible within statistics suggesting that Eq. 6 provides a good
approximation for torelon energies for the intermediate values of lengths
between $L=16$ and $L=24$. In Fig. 10 we provide the prediction of Eq. 6 with
string tension fixed to the value extracted for $L=24$; both points lay on the
Nambu-Goto prediction.

Figure 9: The effective mass of the torelon ground state for the ensembles
cB4.06.16(left panel), cB4.06.24(middle panel) and cC4.05.24(right panel).

A careful look at the masses of the torelon as well as the pion masses shown
on the sixth column of Tab. 1, indicates that for cB4.06.24, the torelon
ground state is larger than twice the pion mass, and in accordance with Nambu-
Goto. Furthemore, the first excited state of torelon for the same ensemble is
much larger $\sim 0.9508(870)$. This disfavors the scenario of having a
complete string breaking of the flux-tube to two pions. Of course possible
string breaking phenomena could have an effect on the dynamics of the flux-
tube and drive the string tension to lower values. As a matter of fact, if we
compare the string tension extracted for $N_{f}=4$ QCD compared to the string
tension extracted in pure gauge theory for the same values of $t_{0}/a^{2}$
(Tab. 3) we observe that the dynamical fermions induce a suppression in
$a^{2}\sigma$ of 50-70 %. At this point we come across the subtlety on what we
should be using as a valid scale in order to give physical units in the
quantities being measured. The two options we can consider are the Wilson flow
time and the string tension; this question requires dedicated investigation.
However, for the purposes of this work it looks more suitable to use the
Wilson flow time which appears to result to gluonic quantities such as
glueball masses affected only negligibly due to the inclusion of dynamical
quarks.

$l$$aE$

Figure 10: The energy of the torelon ground state for $l=16a$ and $l=24a$ for
$N_{f}=4$ and $\beta=1.778$ as well as the Nambu-Goto formula of Eq. 6 with
string tension taken from the estimation extracted for $L=24a$.

Focusing now on the ensemble cB4.06.16, we observe that the value of the
torelon mass is $am_{T}\sim 0.22$ while the potential di-torelon state is
$am_{DT}\sim 0.51$ (Tab. 6). This is in accordance with the scenario of the
ground state of $E^{++}$ representation corresponding to the di-torelon, given
that it is slightly larger than twice the mass of the torelon.

##  9 Comparison with other results

In this section we provide a comparison of our results for the scalar, tensor
as well as pseudoscalar glueball masses with other results from the
literature. In this comparison we do not consider the $0^{++}$ additional
lightest state which appears to have quark content.

We begin with early studies on the unquenched glueball spectrum. In 2001 the
authors of Ref. [34] used for first time, in addition to pure gluonic
operators, also $q{\bar{q}}$ operators which were included in the variational
analysis and the mass of the flavor-singlet state was measured in a system
with $N_{f}=2$ flavors of clover improved Wilson fermions. The results
obtained out of this work demonstrated a suppression on the flavor-singlet
scalar mass compared to our work.

$J^{PC}$MeV

Figure 11: Comparison of our results for the ground states of scalar
($0^{++}$), tensor ($2^{++}$) and pseudoscalar ($0^{-+}$) denoted by the
magenta bands with other unquenched results [35, 6, 36, 37] as well as the
quenched results from [8].

In Ref. [6], the authors studied a system with $N_{f}=2+1$ flavors for ASQTAD
improved staggered fermions, at mass of $m_{\pi}=360$ MeV. They used the
variational technique in which in addition to the standard gluonic operators
they also included operators resembling glueball scattering states. There is
agreement with our results for the scalar and tensor glueball masses within
the bounds of statistical accuracy. The authors obtained no results for the
mass of the pseudoscalar glueball. These results are shown in yellow color in
Figure 11.

Moving on to Ref. [37], the work was carried out using $N_{f}=2$ flavors with
the masses of clover-improved Wilson fermions corresponding to pion masses of
$m_{\pi}=650$ MeV and $m_{\pi}=960$ MeV. Here, only gluonic operators built
from Wilson loops were used in the calculation of the corresponding
correlation functions. While the masses of the tensor and pseudoscalar
glueballs are in agreement with our findings, there is disagreement with
regards to the mass of the scalar glueball; this is true for both fermion
masses the authors used. More specifically, the author’s findings on the
scalar glueball mass show a suppression by $\sim 20\%$ compared to our
results. The above results for $m_{\pi}=650$ MeV and $m_{\pi}=960$ MeV are
denoted in black and purple respectively in Figure 11.

Furthermore, comparing our results to those of [38], we find an agreement on
the scalar and tensor glueball masses as well as an agreement within one sigma
for the pseudoscalar one. In this work, the authors used gluonic operators on
two $N_{f}=2+1$ RBC/UKQCD gauge ensembles at fermion masses corresponding to
physical pion masses of $m_{\pi}=140$ MeV. To enhance the statistics the
authors used cluster decomposition which breaks translation invariance
spontenaously. The above results are presented in yellow in Figure 11.

We now proceed to compare our results with data obtained on the quenched
glueball spectrum. First, we consider the work in [39], in which the authors
used improved lattice gauge actions on anisotropic lattices. Their findings
for the scalar, tensor and pseudoscalar glueball masses are all in agreement
with ours.

We now take a look at the results obtained in yet another work on the quenched
glueball spectrum [8]. As evident in Figure 11, while the mass of the tensor
glueball is in full agreement with our findings, the mass of the scalar as
well as pseudoscalar agree with our data, at the level of one sigma. The
results from Ref. [8] are presented in brown colour in Figure 11.

##  10 Conclusions

In this paper we have obtained the spectrum of glueball masses in the scalar,
tensor and pseudo-scalar channels, $J=0^{++}$, $2^{++}$ and $0^{-+}$
respectively, using lattice configurations produced with light dynamical
quarks. First, we obtained the spectrum of glueballs on configurations
produced with $N_{f}=4$ degenerate light quarks with masses corresponding to a
pion mass of $m_{\pi}\sim 250$ MeV. Furthermore, we investigated possible
finite volume as well as discretization effects. We adopted as a reference
investigation the spectrum calculated within the ensemble cB4.06.24 for which
the pion mass exhibits negligible finite volume effects. We compared our
results to the data obtained using pure gauge configurations by keeping
$t_{0}/a^{2}$ fixed, assuming that lattice artifacts are insignificant
compared to the large statistical uncertainties accompanying these results.
Strikingly, the $A_{1}^{++}$ spectrum appears to contain an additional
(lightest) state, whose mass is between 2 and 4 pion masses, suggesting that
it could be a multipion state; this is one of the most significant results of
this work. Given that the first excitation level of $A_{1}^{++}$, for
$N_{f}=4$, has a mass close to the ground state in pure gauge, one may
hypothesise that the ground state for $A_{1}^{++}$ in $N_{f}=4$ corresponds to
a decay state of the glueball. By contrast, the glueball $2^{++}$ ground state
in $N_{f}=4$ is in good agreement with the $2^{++}$ ground state in the
$SU(3)$ gauge theory and, in addition, the $0^{-+}$ ground state in $N_{f}=4$
is very close to the $0^{-+}$ ground state of the $SU(3)$ theory.

We attempted to shed more light on the additional state appearing in the
$A_{1}^{++}$ channel by extracting the spectrum of the $A_{1}^{++}$
irreducible representation with $N_{f}=2+1+1$ Twisted Mass fermions for a
fixed lattice spacing of $\sim 0.09$ fm and two different pion masses
revealing that the mass of the $A_{1}^{++}$ ground state depends strongly on
the pion mass, unlike the first excited state. This can be interpreted as a
confirmation that the $A_{1}^{++}$ additional state does indeed have a large
”quark content”. This suggests that in future work multimeson operators should
be included in the variational basis in order to identify mixings. Taken
altogether our findings, in combination with data obtained from other
projects, suggest that the glueball spectrum depends very little on the mass
of dynamical quarks.

Finally, we extracted the torelon mass, and from that the associated string
tension. The main outcome of this analysis is that (1) the energy of the
torelon increases with its length and is consistent with the Nambu-Goto
bosonic string expression (10) and that (2) the string tension is suppressed
by 50-70 % due to the inclusion of the $N_{f}=4$ dynamical quarks. Given that
string breaking can take place due to the inclusion of dynamical quarks, this
was perhaps unexpected.

Ensemble | ground state |  $1^{\rm st}$ excited |  $2^{\rm nd}$ excited
---|---|---|---
cB4.06.16 | 0.4257(297) | 0.6303(364) | 0.8462(525)
cB4.06.24 | 0.4449(264) | 0.6445(422) | 0.9916(743)
cC4.05.24 | 0.3666(213) | 0.5131(524) | 0.7210(766)
Table 5: The energies in lattice units $(aE)$ for the ground, first and second excited states for the $A_{1}^{++}$ irreducible representation and the three $N_{f}=4$ ensembles. Ensemble | ground state |  $1^{\rm st}$ excited
---|---|---
cB4.06.16 | 0.5146(118) | 0.8722(667)
cB4.06.24 | - | 0.8666(622)
cC4.05.24 | - | 0.7078(449)
Table 6: The energies in lattice units $(aE)$ for the ground state which corresponds to the di-torelon as well as first excited state, for the $E^{++}$ irreducible representation and the three $N_{f}=4$ ensembles. Ensemble | ground state
---|---
cB4.06.16 | 0.8680(627)
cB4.06.24 | 0.8699(740)
cC4.05.24 | 0.6890(410)
Table 7: The energies in lattice units $(aE)$ for the $T_{2}^{++}$ irreducible representation and the three $N_{f}=4$ ensembles. Ensemble | ground state
---|---
cB4.06.16 | 0.9276(435)
cB4.06.24 | 0.8252(764)
cC4.05.24 | 0.7182(498)
Table 8: The energies in lattice units $(aE)$ for the $A_{1}^{-+}$ irreducible representation and the three $N_{f}=4$ ensembles. Ensemble | ground state | $a^{2}\sigma$ | $a^{2}\sigma_{SU(3)}$
---|---|---|---
cB4.06.16 | 0.2167(100) | 0.01824(60) | 0.03063(32)
cB4.06.24 | 0.3655(241) | 0.01715(100) | 0.03172(43)
cC4.05.24 | 0.2247(183) | 0.011365(735) | 0.02350(18)
Table 9: The energies in lattice units $(aE)$ for the ground state of the
torelon for the three $N_{f}=4$ ensembles.

##  11 Acknowledgements

We express our gratitude to K. Hadjiyannakou and G. Spanoudes for providing
access to the $N_{f}=4$ and $N_{f}=2+1+1$ configurations. AA is indebted to
David Gross, Davide Vadacchino and Biagio Lucini for critical discussions over
the spectrum of glueballs. MT acknowledges support in part by the Simons
Collaboration on Confinement and QCD Strings. AA was supported by the Horizon
2020 European research infrastructures programme ”NI4OS-Europe” with grant
agreement no. 857645. Results were obtained using Cyclone High Performance
Computer at The Cyprus Institute, under access with id p014, the Oxford
Theoretical Physics cluster, the resources of computer cluster Isabella based
in SRCE - University of Zagreb University Computing Centre, as well as using
computing time granted by the John von Neumann Institute for Computing (NIC)
on the supercomputer JUWELS Cluster [40] at the Jülich Supercomputing Centre
(JSC). JF received financial support by the German Research Foundation (DFG)
research unit FOR5269 ”Future methods for studying confined gluons in QCD”, by
the H2020 project PRACE 6-IP (GA No. 82376) and by the EuroCC project (GA No.
951732) funded by the Deputy Ministry of Research, Innovation and Digital
Policy and the Cyprus Research and Innovation Foundation and the European
High-Performance Computing Joint Undertaking (JU) under grant agreement No
951732.

## References

  * [1] Denis Parganlija “Mesons, PANDA and the scalar glueball” In _J. Phys. Conf. Ser._ 503, 2014, pp. 012010 DOI: 10.1088/1742-6596/503/1/012010
  * [2] D. M. Asner “Physics at BES-III” In _Int. J. Mod. Phys. A_ 24, 2009, pp. S1–794 arXiv:0809.1869 [hep-ex]
  * [3] Davide Vadacchino “A review on Glueball hunting” Zenodo, 2022 DOI: 10.5281/zenodo.7338133
  * [4] Davide Vadacchino “A review on Glueball hunting” In _39th International Symposium on Lattice Field Theory_ , 2023 arXiv:2305.04869 [hep-lat]
  * [5] Eberhard Klempt “Glueballs, a fulfilled promise of QCD?”, 2022 arXiv:2211.12901 [hep-ph]
  * [6] E. Gregory et al. “Towards the glueball spectrum from unquenched lattice QCD” In _JHEP_ 10, 2012, pp. 170 DOI: 10.1007/JHEP10(2012)170
  * [7] S. Okubo “$\phi$-meson and unitary symmetry model” In _Physics Letters_ 5.2, 1963, pp. 165–168 DOI: https://doi.org/10.1016/S0375-9601(63)92548-9
  * [8] Andreas Athenodorou and Michael Teper “The glueball spectrum of SU(3) gauge theory in 3 + 1 dimensions” In _JHEP_ 11, 2020, pp. 172 DOI: 10.1007/JHEP11(2020)172
  * [9] C. Alexandrou et al. “Complete flavor decomposition of the spin and momentum fraction of the proton using lattice QCD simulations at physical pion mass” In _Phys. Rev. D_ 101.9, 2020, pp. 094513 DOI: 10.1103/PhysRevD.101.094513
  * [10] Y. Iwasaki and T. Yoshie “Renormalization Group Improved Action for SU(3) Lattice Gauge Theory and the String Tension” In _Phys. Lett. B_ 143, 1984, pp. 449–452 DOI: 10.1016/0370-2693(84)91500-4
  * [11] Y. Iwasaki and T. Yoshie “Improved Lattice Action for the SU(2) Gauge Theory and Monte Carlo Calculations” [Erratum: Phys.Lett.B 132, 467 (1983)] In _Phys. Lett. B_ 130, 1983, pp. 77 DOI: 10.1016/0370-2693(83)91067-5
  * [12] P. Weisz “Continuum Limit Improved Lattice Action for Pure Yang-Mills Theory. 1.” In _Nucl. Phys. B_ 212, 1983, pp. 1–17 DOI: 10.1016/0550-3213(83)90595-3
  * [13] Roberto Frezzotti, Pietro Antonio Grassi, Stefan Sint and Peter Weisz “Lattice QCD with a chirally twisted mass term” In _JHEP_ 08, 2001, pp. 058 arXiv:hep-lat/0101001
  * [14] R. Frezzotti and G.C. Rossi “Chirally improving Wilson fermions. 1. O(a) improvement” In _JHEP_ 08, 2004, pp. 007 DOI: 10.1088/1126-6708/2004/08/007
  * [15] Constantia Alexandrou “Simulating twisted mass fermions at physical light, strange and charm quark masses” In _Phys. Rev. D_ 98.5, 2018, pp. 054518 DOI: 10.1103/PhysRevD.98.054518
  * [16] R. Frezzotti and G.C. Rossi “Twisted mass lattice QCD with mass nondegenerate quarks” In _Nucl. Phys. B Proc. Suppl._ 128, 2004, pp. 193–202 DOI: 10.1016/S0920-5632(03)02477-0
  * [17] R. Frezzotti, G. Martinelli, M. Papinutto and G.C. Rossi “Reducing cutoff effects in maximally twisted lattice QCD close to the chiral limit” In _JHEP_ 04, 2006, pp. 038 DOI: 10.1088/1126-6708/2006/04/038
  * [18] C. Alexandrou et al. “Proton and neutron electromagnetic form factors from lattice QCD” In _Phys. Rev. D_ 100.1, 2019, pp. 014509 DOI: 10.1103/PhysRevD.100.014509
  * [19] N. Cabibbo and E. Marinari “A New Method for Updating SU(N) Matrices in Computer Simulations of Gauge Theories” In _Phys. Lett. B_ 119, 1982, pp. 387–390 DOI: 10.1016/0370-2693(82)90696-7
  * [20] Martin Lüscher “Properties and uses of the Wilson flow in lattice QCD” [Erratum: JHEP 03, 092 (2014)] In _JHEP_ 08, 2010, pp. 071 DOI: 10.1007/JHEP08(2010)071
  * [21] A. Francis et al. “Critical point and scale setting in SU(3) plasma: An update” In _Phys. Rev. D_ 91.9, 2015, pp. 096002 DOI: 10.1103/PhysRevD.91.096002
  * [22] Martin Luscher and Stefan Schaefer “Lattice QCD without topology barriers” In _JHEP_ 07, 2011, pp. 036 DOI: 10.1007/JHEP07(2011)036
  * [23] Marco Cè, Cristian Consonni, Georg P. Engel and Leonardo Giusti “Non-Gaussianities in the topological charge distribution of the SU(3) Yang–Mills theory” In _Phys. Rev. D_ 92.7, 2015, pp. 074502 DOI: 10.1103/PhysRevD.92.074502
  * [24] M. Luscher and P. Weisz “Definition and General Properties of the Transfer Matrix in Continuum Limit Improved Lattice Gauge Theories” In _Nucl. Phys. B_ 240, 1984, pp. 349–361 DOI: 10.1016/0550-3213(84)90270-0
  * [25] Martin Luscher and Ulli Wolff “How to Calculate the Elastic Scattering Matrix in Two-dimensional Quantum Field Theories by Numerical Simulation” In _Nucl. Phys. B_ 339, 1990, pp. 222–252 DOI: 10.1016/0550-3213(90)90540-T
  * [26] Biagio Lucini, Michael Teper and Urs Wenger “Glueballs and k-strings in SU(N) gauge theories: Calculations with improved operators” In _JHEP_ 06, 2004, pp. 012 DOI: 10.1088/1126-6708/2004/06/012
  * [27] M. Teper “An Improved Method for Lattice Glueball Calculations” In _Phys. Lett. B_ 183, 1987, pp. 345 DOI: 10.1016/0370-2693(87)90976-2
  * [28] Andreas Athenodorou, Barak Bringoltz and Michael Teper “Closed flux tubes and their string description in D=3+1 SU(N) gauge theories” In _JHEP_ 02, 2011, pp. 030 DOI: 10.1007/JHEP02(2011)030
  * [29] P. Goddard, J. Goldstone, C. Rebbi and C.B. Thorn “Quantum dynamics of a massless relativistic string” In _Nuclear Physics B_ 56.1, 1973, pp. 109–135 DOI: https://doi.org/10.1016/0550-3213(73)90223-X
  * [30] P. Di Vecchia, K. Fabricius, G.C. Rossi and G. Veneziano “Preliminary evidence for UA(1) breaking in QCD from lattice calculations” In _Nuclear Physics B_ 192.2, 1981, pp. 392–408 DOI: https://doi.org/10.1016/0550-3213(81)90432-6
  * [31] Martin Luscher “Trivializing maps, the Wilson flow and the HMC algorithm” In _Commun. Math. Phys._ 293, 2010, pp. 899–919 DOI: 10.1007/s00220-009-0953-7
  * [32] Szabolcs Borsanyi “High-precision scale setting in lattice QCD” In _JHEP_ 09, 2012, pp. 010 DOI: 10.1007/JHEP09(2012)010
  * [33] Georg Bergner et al. “Influence of topology on the scale setting” In _Eur. Phys. J. Plus_ 130.11, 2015, pp. 229 DOI: 10.1140/epjp/i2015-15229-7
  * [34] Craig McNeile and Christopher Michael “Mixing of scalar glueballs and flavor singlet scalar mesons” In _Phys. Rev. D_ 63, 2001, pp. 114503 DOI: 10.1103/PhysRevD.63.114503
  * [35] Yidian Chen, Danning Li and Mei Huang “The dynamical holographic QCD method for hadron physics and QCD matter”, 2022 arXiv:2206.00917 [hep-ph]
  * [36] Gunnar S. Bali et al. “Static potentials and glueball masses from QCD simulations with Wilson sea quarks” In _Phys. Rev. D_ 62, 2000, pp. 054503 DOI: 10.1103/PhysRevD.62.054503
  * [37] Wei Sun et al. “Glueball spectrum from $N_{f}=2$ lattice QCD study on anisotropic lattices” In _Chin. Phys. C_ 42.9, 2018, pp. 093103 DOI: 10.1088/1674-1137/42/9/093103
  * [38] Feiyu Chen et al. “Glueballs at Physical Pion Mass”, 2021 arXiv:2111.11929 [hep-lat]
  * [39] Y. Chen “Glueball spectrum and matrix elements on anisotropic lattices” In _Phys. Rev. D_ 73, 2006, pp. 014516 DOI: 10.1103/PhysRevD.73.014516
  * [40] Jülich Supercomputing Centre “JUWELS: Modular Tier-0/1 Supercomputer at the Jülich Supercomputing Centre” In _Journal of large-scale research facilities_ 5.A135, 2019 DOI: 10.17815/jlsrf-5-171