SU(N) gauge theories in 3+1 dimensions:
glueball spectrum, string tensions and topology

Andreas Athenodoroua,b and Michael Teperc,d

a Dipartimento di Fisica, Università di Pisa and INFN, Sezione di Pisa,

Largo Pontecorvo 3, 56127 Pisa, Italy.

b Computation-based Science and Technology Research Center, The Cyprus
Institute,

20 Kavafi Str., Nicosia 2121, Cyprus

cRudolf Peierls Centre for Theoretical Physics, University of Oxford,

Parks Road, Oxford OX1 3PU, UK

dAll Souls College, University of Oxford,

High Street, Oxford OX1 4AL, UK

Abstract

We calculate the low-lying glueball spectrum, several string tensions and some
properties of topology and the running coupling for $SU(N)$ lattice gauge
theories in $3+1$ dimensions. We do so for $2\leq N\leq 12$, using lattice
simulations with the Wilson plaquette action, and for glueball states in all
the representations of the cubic rotation group, for both values of parity and
charge conjugation. We extrapolate these results to the continuum limit of
each theory and then to $N=\infty$. For a number of these states we are able
to identify their continuum spins with very little ambiguity. We calculate the
fundamental string tension and $k=2$ string tension and investigate the $N$
dependence of the ratio. Using the string tension as the scale, we calculate
the running of a lattice coupling and confirm that $g^{2}(a)\propto 1/N$ for
constant physics as $N\to\infty$. We fit our calculated values of
$a\surd\sigma$ with the 3-loop $\beta$-function, and extract a value for
$\Lambda_{\overline{MS}}$, in units of the string tension, for all our values
of $N$, including $SU(3)$. We use these fits to provide analytic formulae for
estimating the string tension at a given lattice coupling. We calculate the
topological charge $Q$ for $N\leq 6$ where it fluctuates sufficiently for a
plausible estimate of the continuum topological susceptibility. We also
calculate the renormalisation of the lattice topological charge,
$Z_{Q}(\beta)$, for all our $SU(N)$ gauge theories, using a standard
definition of the charge, and we provide interpolating formulae, which may be
useful in estimating the renormalisation of the lattice $\theta$ parameter. We
provide quantitative results for how the topological charge ‘freezes’ with
decreasing lattice spacing and with increasing $N$. Although we are able to
show that within our typical errors our glueball and string tension results
are insensitive to the freezing of $Q$ at larger $N$ and $\beta$, we choose to
perform our calculations with a typical distribution of $Q$ imposed upon the
fields so as to further reduce any potential systematic errors.

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Contents

  1. 1 Introduction
  2. 2 Calculating on a lattice
    1. 2.1 lattice setup
    2. 2.2 energies and correlators
    3. 2.3 systematic errors
    4. 2.4 topological freezing
  3. 3 String tensions
    1. 3.1 finite volume corrections
    2. 3.2 fundamental ($k=1$) string tensions
    3. 3.3 $k=2$ string tensions
    4. 3.4 $N\to\infty$ extrapolations
  4. 4 Running coupling
    1. 4.1 scaling with N
    2. 4.2 perturbative running and ${\mathrm{\Lambda_{\overline{MS}}}}$
    3. 4.3 interpolating and extrapolating functions for $a(\beta)$
  5. 5 Glueball masses
    1. 5.1 quantum numbers
    2. 5.2 finite volume effects
    3. 5.3 lattice masses
    4. 5.4 strong-to-weak coupling transition
    5. 5.5 continuum masses
    6. 5.6 $N\to\infty$ extrapolation
    7. 5.7 continuum spins
    8. 5.8 scattering states
    9. 5.9 some comparisons
  6. 6 Topological fluctuations
    1. 6.1 topology and cooling
    2. 6.2 tunneling between topological sectors
    3. 6.3 topological susceptibility
    4. 6.4 ${Z_{Q}(\beta)}$ and lattice $\theta$ parameter
  7. 7 Conclusions
  8. A Lattice running couplings
  9. B Scattering states

##  1 Introduction

In this paper we calculate various physical properties of $SU(N)$ gauge
theories in $3+1$ dimensions. We do so by performing calculations in the
corresponding lattice gauge theories over a sufficient range of lattice
spacings, $a$, and with enough precision that we can obtain plausible
continuum extrapolations. We also wish to be able to extrapolate to the
theoretically interesting $N=\infty$ limit and to compare this to the
physically interesting $SU(3)$ theory. To do so we have performed our
calculations for $N=2,3,4,5,6,8,10,12$ gauge theories.

Our main aim is to provide a calculation of the low-lying ‘glueball’ mass
spectrum for all quantum numbers and all $N$. This means calculating the
lowest states in all the irreducible representions, $R$, of the rotation group
of a cubic lattice, and for both values of parity $P$ and charge conjugation
parity $C$. We also calculate the confining string tension, $\sigma$, so as to
provide a useful scale for the glueball masses. In addition we calculate a
number of other interesting quantities. However these are side-products of our
glueball calculations rather than being dedicated calculations of these
quantities. The first is the $k=2$ string tension, $\sigma_{k=2}$ and the
nature of its approach to $N=\infty$. The second is an estimate of the scale
parameter $\Lambda_{\overline{MS}}$ for all our $SU(N)$ gauge theories.
Thirdly, in order to monitor the topological ‘freezing’ as $a$ becomes small
at fixed $N$, or $N$ becomes large at fixed $a$, we calculate the interval (in
units of our Monte Carlo steps) between changes in the topological charge, as
well as the renormalisation of that charge, and also the topological
susceptibility.

Our calculations of the glueball spectrum are intended to make necessary
improvements to previous work. We recall that the pioneering glueball
calculations in [1, 2] were restricted to obtaining the continuum masses of
the lightest and first excited $J^{PC}=0^{++}$ glueballs and the lightest
$2^{++}$ glueball, in units of the string tension, with the continuum spin
assignments being based on plausible assumptions. This sufficed to provide the
first fully non-perturbative demonstration that the basic physical quantities
of the $SU(3)$ gauge theory are ‘close to’ those of $SU(\infty)$, as had long
been hoped, but it did not provide the kind of detailed $N=\infty$ glueball
spectrum that would be useful in testing theoretical approaches. Such a
detailed spectrum was subsequently provided in [3] where the lightest masses
in all the irreducible representations of the rotation group of the cubic
lattice were calculated for $N\in[2,8]$ and extrapolated to $N=\infty$. In
addition a serious effort was made in that work to identify those states that
might be multi-glueball ‘scattering’ states or finite volume artifacts (based
on conjugate pairs of winding flux tubes), rather than being the single
glueballs that one would obtain in an infinite volume. The important drawback
of this calculation is that it was made at a fixed value of the lattice
spacing $a$ (‘fixed’ in units of the deconfinement temperature), corresponding
to a value of $\beta\simeq 5.895$ in $SU(3)$. This corresponds to a coarse
lattice discretisation, typically between the two largest lattice spacings
used in our present work. This has two important adverse consequences. The
first is that it leaves uncertain the values of the masses in the desired
continuum limit. (Although the calculations in [3] and earlier $SU(3)$
calculations left room for being optimistic about the limited size of any
lattice corrections.) The second problem is that since the lattice masses $aM$
are large when $a$ is large, as in that paper, it makes the extraction of
heavier glueball masses much more ambiguous and the statistical errors much
larger than they would be at smaller values of $a$. A corollary is that it
renders continuum spin assignments more ambiguous, except for the very
lightest glueballs, since the spin assignment depends on observing near-
degeneracies amongst states in appropriate irreducible representations of the
lattice rotation group, and for this it is clearly essential to achieve an
adequate precision in the various mass estimates. In addition to this earlier
work, there recently appeared a potentially more relevant paper [4] when our
work was largely completed. This paper provides a pioneering calculation of
the masses of the lightest glueballs in all the irreducible representations,
of a number of continuum $Sp(2N)$ gauge theories with an extrapolation to
$N=\infty$. In principle this extrapolation is equally relevant to the
$N\to\infty$ limit of $SU(N)$ gauge theories, as pointed out in [4], since the
$Sp(2N)$ and $SU(N)$ gauge theories share a common (perturbative) planar
limit. In practice, however, the usefulness for $SU(N\to\infty)$ turns out to
be very limited. Firstly, $Sp(2N)$ has no $C=-$ sector, so there are no
predictions for half the $SU(N\to\infty)$ spectrum. Secondly the limited range
of $N$ in [4], coupled with the fact that in $Sp(2N)$ the leading correction
to $N=\infty$ is $O(1/N)$ renders the large-$N$ extrapolations less convincing
and much less precise when compared to $SU(N)$ where the leading correction is
$O(1/N^{2})$. The net result is that the errors on the $SU(N\to\infty)$
glueball masses obtained in [4], are larger by a factor of $\sim 5-20$
compared to the values obtained in this paper, and the mass estimates are
mainly for just the ground states in each channel. This means that the
evidence for the assignment of a continuum spin $J$ only appears plausible for
the $J^{P}=0^{+},0^{-},2^{+},2^{-}$ ground states. Of course, none of these
comments detract from the success of [4] in their primary aim of elucidating
the mass spectrum of $Sp(2N)$ gauge theories.

In this paper we provide a calculation of the masses of the ground states and
some excited states in all the irreducible representations $R$ of the rotation
group of our cubic lattice with an extrapolation of these masses (in units of
the string tension) to the continuum limit. To make this extrapolation more
reliable and more precise we extend the range of our calculations to much
smaller lattice spacings than earlier work. This enables us to extract the
masses of the heavier ground states and most of the first few excited states
much more reliably than earlier calculations. As an important by-product, all
this will allow us to make a significant number of continuum spin assignments
after the extrapolation to the continuum limit. We recall that the assignment
of a continuum spin typically depends on observing near-degeneracies amongst
both ground and excited states in the various irreducible representations of
the lattice rotation group, and this requires both precision in the mass
estimates and small enough lattice spacings for the masses of the heavier
excited states to be plausibly estimated. We also extend the range of our
calculations to larger values of $N$. This is not only to make the
extrapolation to $N=\infty$ somewhat more reliable and precise but also, and
perhaps more importantly, to help in excluding from our $N\to\infty$ glueball
spectrum any multiglueball scattering states and finite volume states, since
these states will decouple from our single trace operators as $N\uparrow$.

This kind of calculation is of course standard in the case of $SU(3)$, see for
example [5, 6, 7, 9, 8, 10, 11], and what we have attempted to do is to bring
the $SU(N)$ calculations towards a similar level of sophistication. The
computational cost of performing calculations at larger $N$ means that further
improvement to our work is still desirable. (It would also be very interesting
to see the existing calculations at very large $N$ using space-time reduction
[12] extended to calculations of the glueball spectrum, just as they have been
in $2+1$ dimensions [13].) Nonetheless we are able in this paper to provide
the first calculation of the masses of the ground states in all the $R^{PC}$
channels of the continuum $SU(N\to\infty)$ gauge theory, as well as some
excited states in most channels.

The plan of the paper is as follows. In the next section we introduce our
lattice setup and describe how we calculate energies from correlators. We
discuss some of the main systematic errors affecting these calculations and
how we deal with them, with a particular focus on the rapid loss of tunneling
between sectors of differing topology as $N\uparrow$. In Section 3 we describe
in detail our calculation of the confining string tension, $\sigma$, which we
will later use as the physical scale in which to express our glueball masses.
As a side product we also calculate the string tension, $\sigma_{k=2}$, of the
lightest flux tubes carrying $k=2$ units of fundamental flux. In Section 3.4
we study the approach of the ratio $\sigma_{k=2}/\sigma$ to the $N=\infty$
limit so as to address the old controversy concerning the power of $1/N$ of
the leading correction. In Section 4 we show how our precise calculation of
the fundamental string tension, $a^{2}\sigma$, as a function of the lattice
coupling enables us to confirm the expected scaling of $g^{2}$ with $N$, and
motivates (lattice improved) perturbative fits that allow us, in Section 4.2,
to estimate a value for $\Lambda_{\overline{MS}}$ as a function of $N$. We
also provide, in Section 4.3, some analytic interpolation/extrapolation
formulae for the variation of the string tension with the coupling, which may
be of use in other calculations. We then turn, in Section 5, to our main
calculation in this paper, which is that of the low-lying glueball spectrum.
We calculate the masses on the lattice, extrapolate to the continuum limit,
and then extrapolate to $N=\infty$. Although these states are classified
according to the representations of the rotation group of our cubic lattice,
we are able to identify the continuum spins in many cases, as described in
Section 5.7. In doing all this we need to address the problem of the extra
states that winding modes introduce into the glueball spectrum and also the
possible presence of multi-glueball ’scattering’ states. We complete Section 5
with a brief comparison of our results with those of some earlier
calculations. We then return, in Section 6, to some of the properties of the
topological fluctuations in our lattice fields. After illustrating in Section
6.1 how our cooling algorithm reproduces the topological charge of a lattice
field, we calculate in Section 6.2 the rate of topological freezing with
increasing $N$ and with decreasing $a(\beta)$ and compare our results to the
simplest theoretical expectations. In Section 6.3 we provide our results for
the topological susceptibility in the continuum limit of $SU(N\leq 6)$ gauge
theories and in the $N\to\infty$ limit. Then, in Section 6.4, we calculate the
multiplicative renormalisation of our lattice topological charge for each
$SU(N)$ group and provide interpolating formulae which may be useful in
calculations with a $\theta$ parameter in the lattice action. Section 7
summarises our main results.

Finally we remark that in parallel with the present calculations most of our
$SU(3)$ calculations, which are of particular physical interest, have recently
been published separately [11].

##  2 Calculating on a lattice

###  2.1 lattice setup

We work on hypercubic lattices of size $L_{s}^{3}L_{t}$ with lattice spacing
$a$ and with periodic boundary conditions on the fields. The Euclidean time
extent, $aL_{t}$, is always chosen large enough that we are in the confining
phase of the theory, at a temperature that is well below the deconfining phase
transition [14]. Our fields are $SU(N)$ matrices, $U_{l}$, assigned to the
links $l$ of the lattice. The Euclidean path integral is

$Z=\int{\cal{D}}U\exp\\{-\beta S[U]\\},$ (1)

where ${\cal{D}}U$ is the Haar measure and we use the standard plaquette
action,

$\beta
S=\beta\sum_{p}\left\\{1-\frac{1}{N}{\text{ReTr}}U_{p}\right\\}\quad;\quad\beta=\frac{2N}{g^{2}_{L}}.$
(2)

Here $U_{p}$ is the ordered product of link matrices around the plaquette $p$.
We write $\beta=2N/g^{2}_{L}$ since in this way we recover the usual continuum
action when we take the continuum limit of the lattice theory and replace
$g^{2}_{L}\to g^{2}$. The subscript $L$ reminds us that this coupling is
defined in a specific coupling scheme corresponding to the lattice and the
plaquette action. Since $g^{2}_{L}$ is the bare coupling corresponding to a
short distance cut-off $a$, it provides a definition of the running coupling
$g^{2}_{L}(a)=g^{2}_{L}$ on the length scale $a$. Since the theory is
asymptotically free $g^{2}_{L}(a)\to 0$ as $a\to 0$ and hence $\beta\to\infty$
as $a\to 0$ and so we can decrease the cut-off $a$ and so approach the desired
continuum limit of the theory by increasing $\beta$. Our calculations of this
lattice path integral are carried out via a standard Monte Carlo using a
mixture of Cabibbo-Marinari heat bath and over-relaxation sweeps through the
lattice. We typically perform $\sim 2\times 10^{6}$ sweeps at each value of
$\beta$ at each lattice size, and we typically calculate correlators every
$\sim 25$ sweeps and the topological charge every $50$ or $100$ sweeps.
Naturally we choose values of $\beta$ that place us on the weak coupling
branch of the lattice theory. We discuss the ‘bulk’ transition that separates
weak coupling from strong coupling, and becomes a first order transition for
$N\geq 5$ [15], in Section 5.4. There are of course other possible choices for
the action, and we briefly comment upon our choice in Section 5.9, where we
compare our results to those of other recent calculations.

We simulate $SU(N)$ lattice gauge theories for $N=2,3,4,5,6,8,10,12$ over a
range of values of $\beta$ so as to be able to plausibly obtain, by
extrapolation, the glueball spectra and string tensions of the corresponding
continuum gauge theories. A summary of the basic parameters of our
calculations is given in Tables 1-8. For each of our $SU(N)$ calculations we
show the values of $\beta$, the lattice sizes, the average plaquette, the
string tension, $a^{2}\sigma$, and the mass gap, $am_{G}$. In the header of
each table we show the approximate spatial size in units of the string
tension. Since finite volume corrections are expected to decrease with
increasing $N$ we decrease the lattice volume as we increase $N$. This
expectation needs to be confirmed by explicit calculations which we shall
provide in Section 5.2.

###  2.2 energies and correlators

We calculate energies from correlation functions in the standard way. Suppose
we wish to calculate glueball masses in some representation $R^{PC}$. If one
picks an operator $\phi(t)$ with quantum numbers $R^{PC}$ and momentum $p=0$
then the correlator will provide us with the energies $E_{i}$ of states with
those quantum numbers

$C(t)=\langle\phi^{\dagger}(t)\phi(0)\rangle=\sum_{n=0}|c_{n}|^{2}\exp\\{-E_{n}t\\}\stackrel{{\scriptstyle
t\to\infty}}{{\longrightarrow}}|c_{0}|^{2}\exp\\{-E_{0}t\\},$ (3)

where $E_{n}\leq E_{n+1}$ and $E_{0}$ is the lightest state with $R^{PC}$
quantum numbers – which will often be the lightest glueball in that sector.
(If $\phi$ has vacuum quantum numbers then one uses the vacuum subtracted
operator.) Since on the lattice time is measured in lattice units, $t=an_{t}$,
what we obtain is the value, $aE_{0}$, of the mass in lattice units. If we
calculate two masses, $aM$ and $a\mu$, in this way then the lattice spacing
drops out of the ratio and if we calculate the ratio for several values of
$a(\beta)$ we can extrapolate to the continuum limit in the standard way,
using

$\frac{aM(a)}{a\mu(a)}=\frac{M(a)}{\mu(a)}\simeq\frac{M(0)}{\mu(0)}+ca^{2}\mu^{2}$
(4)

once $a(\beta)$ is small enough. (This standard tree level extrapolation could
be improved with perturbative corrections [16] and with higher order power
corrections but our calculations are not so extensive and precise as to
motivate such modifications.)

The starting point for an operator $\phi$ is the trace of a closed loop on the
lattice. For an operator that projects onto glueballs the loop should be
contractible. A non-contractible loop that closes across the periodic spatial
boundary will project onto a confining flux tube that winds around that
spatial torus. For the contractible glueball loop one takes a suitable linear
combination of rotations of that loop for it to be in the desired
representation $R$ of the lattice rotation group and together with the parity
inverse this allows us to form an operator with parity $P$. The real and
imaginary parts of the original loop will correspond to $C=+$ and $C=-$
respectively. Summing this linear combination over all spatial sites at time
$t$ gives us the operator with $p=0$.

The statistical errors on $C(t)$ are roughly independent of $t$ while its
value decreases exponentially with $t$, so if we are to estimate $aE_{0}$ from
$C(t)$ in eqn(3) we need to be able to do this at small $t$. (The fluctuations
of $\phi^{\dagger}(t)\phi(0)$ involve the higher order correlator
$(\phi^{\dagger}(t)\phi(0))^{2}$ which typically has a vacuum channel and this
is independent of $t$.) For this to be possible we need the operator $\phi$ to
have a large overlap onto the state $|n=0\rangle$ that corresponds to $E_{0}$,
i.e. that $|c_{0}|^{2}/\sum_{n}|c_{n}|^{2}\sim 1$. One can achieve this using
iteratively ‘blocked’ link matrices and loops [17, 18] as described in detail
in, for example, [2]. To monitor the approach of $C(t)$ to the asymptotic
exponential decay in eqn(3) it is useful to define an effective energy

$\frac{C(an_{t})}{C(a(n_{t}-1))}=\exp\\{-aE_{eff}(n_{t})\\}.$ (5)

(In practice we use a cosh modification of this definition to take into
account the periodicity in the temporal direction.) If, within errors,
$aE_{eff}(n_{t})=\mathrm{const}$ for $n_{t}\geq\tilde{n}_{t}$ then we can use
$aE_{eff}(\tilde{n}_{t})$ as an estimate for $aE_{0}$, or we can do a simple
exponential fit to $C(t)$ for $n_{t}\geq\tilde{n}_{t}-1$ to estimate $aE_{0}$.

To calculate not just the ground state energy but also some usefully precise
excited state energies from $C(t)$ in eqn(3) would require a precision that is
at present unachievable. Instead the standard strategy is to use a variational
calculation. One chooses some basis of $n_{0}$ operators
$\\{\phi_{i}(t):i=1,..,n_{0}\\}$, calculates the cross-correlators
$C_{ij}(t)=\langle\phi^{\dagger}_{i}(t)\phi_{j}(0)\rangle$ and finds the
linear combination, $\Phi=\Phi_{0}$, that maximises
$C(t)=\langle\Phi^{\dagger}(t)\Phi(0)\rangle$ for some suitable small
$t=t_{0}$. This is then our best estimate of the wave-functional of the ground
state $|n=0\rangle$. We then extract $aE_{0}$ from the asymptotic exponential
decay of $C(t)=\langle\Phi^{\dagger}_{0}(t)\Phi_{0}(0)\rangle$. To calculate
the first excited energy $aE_{1}$ we consider the subspace of
$\\{\phi_{i}:i=1,..,n_{0}\\}$ that is orthogonal to $\Phi_{0}$ and repeat the
above steps. This gives our best estimate $\Phi_{1}$ for $|n=1\rangle$ and we
estimate $aE_{1}$ from the asymptotic exponential decay of the correlator of
$\Phi_{1}(t)$. Similarly by considering the subspace of the operator vector
space that is orthogonal to $\Phi_{0}$ and $\Phi_{1}$ we can obtain estimates
for the next excited state and so on. As with any variational calculation, the
accuracy of such estimates will depend on having a large enough basis of
operators.

###  2.3 systematic errors

In addition to the statistical errors in the calculation of energies and
masses, which can be estimated quite reliably, there are systematic errors
that are harder to control. We will now briefly point to some of these,
leaving a more detailed discussion till later on.

Since the error to signal ratio on $C(n_{t})$ increases at least as fast as
$\propto\exp\\{+aE_{0}n_{t}\\}$, the range of $n_{t}$ which is useful rapidly
decreases when the energy $aE_{0}$ increases, and an obvious systematic error
is that we begin our fit of the asymptotic exponential decay at too small a
value of $t=an_{t}$. Given the positivity properties of our correlators, this
means that as the true value of $E_{0}$ increases, so does our overestimate of
its value. For any particular state this problem will become less severe as we
decrease $a$, and this provides some check on this error.

The correlator of an excited state, using our best variational choice of
operator, may contain small contributions from lighter excitations since our
basis is far from complete. So at large enough $n_{t}$ these will dominate the
correlator and the effective mass will drop below that of the mass of the
state of interest. This does not appear to be a severe problem for us because
the pattern of our statistical errors, which grow rapidly with $n_{t}$,
prevents us from obtaining values at large $n_{t}$. In any case we attempt to
identify at most an intermediate effective mass plateau from which we extract
the desired mass.

If a state has too small an overlap onto our basis of operators then, given
the exponential growth with $t$ of the statistical error on the effective
energy, the state will be, at best, assigned a mass that is much too high and
so it will typically not appear in its correct order in the mass spectrum. If
so this means that there will be a missing state in our calculated glueball
spectrum and this will not be helpful if, for example, one wishes to check
theoretical models of the spectrum against the lattice spectrum. To reduce the
possibility of this occurring one needs to use a large basis of operators,
just as in any variational calculation.

The states we are primarily interested in are single particle glueball states
rather than multi-glueball scattering states. (In our finite spatial volume
multiglueball states are interacting at all times but for simplicity we refer
to them as scattering states.) We use single trace operators which should
ensure that at large $N$ our correlators receive no multi-glueball
contributions but this is not necessarily the case at smaller $N$. This is of
course related to the fact that at smaller $N$ heavier glueball states may
have substantial decay widths. Later in the paper we will provide an
exploratory analysis of this issue by combining single and double trace
operators.

The spectrum in a finite spatial volume will differ from the spectrum in an
infinite volume. One type of correction to a glueball propagator comes, for
example, from a virtual glueball loop where one line of the two in the loop
encircles the spatial torus. Such corrections are exponentially small in,
typically, $M_{G}L$ where $M_{G}$ is the mass gap (the lightest scalar
glueball) and $L$ is the spatial size [19]. These corrections will be very
small in our calculations. For the intermediate volumes on which we work a
much more relevant finite volume correction comes from a confining flux tube
that winds once around one of the spatial tori. Such a flux tube is often
referred to as a torelon. By itself a torelon has zero overlap onto local
particle states. However a state composed of two torelons, where one is
conjugated, will in general have a nonzero overlap. Such a ‘ditorelon’ will
not only shift the energies of the usual particle glueballs but will introduce
extra states into the glueball spectrum. These extra states can be identified
through their volume dependence and in some other ways as we shall see later
on.

In addition to the usual gradual loss of ergodicity as the lattice spacing
decreases, the tunnelling between sectors of differing topological charge $Q$
suffers a much more rapid supression not only as $a\to 0$ at fixed $N$ but
also as $N$ increases at fixed $a$. This affects a significant part of our
calculations and we discuss next how we deal with it in practice.

###  2.4 topological freezing

We begin by briefly summarising the reason that one loses ergodicity in the
topology of lattice gauge fields much faster than in other typical quantities
as $a\to 0$ at fixed $N$ and as $N\to\infty$ at fixed $a$, when using a local
Monte Carlo algorithm,. We then describe where this impacts on our particular
calculations and why this inpact should be minor on theoretical grounds. We
provide explicit calculations to confirm this expectation and, finally, we
describe the additional procedure that we follow to minimise any remaining
systematic errors. It is important to note that all these arguments are for
pure gauge theories and would have to be reconsidered in the presence of light
fermions as in $QCD$.

In a periodic four-volume a continuum gauge field of topological charge $Q$
cannot be smoothly deformed into a gauge field of charge $Q^{\prime}\neq Q$.
So in a sequence of lattice fields that have been generated using a Monte
Carlo algorithm that is local any change in the value of $Q$ should be
increasingly suppressed as we approach the continuum limit. This statement can
be made quantitative as we do below. But before that we note that there is a
separate suppression that occurs at any fixed small value of $a$ when we
increase $N$. To smoothly deform a lattice field from, say, $Q=1$ to $Q=0$ in
a finite volume with periodic gauge potentials, the core of an original
extended instanton must gradually shrink until it disappears within a
hypercube, leaving a simple gauge singularity centered on that hypercube.
Before it does so it will be a small instanton (assuming $a$ is small) which
we can describe using standard semiclassical methods. So its effective action
will be $S_{I}(\rho)\sim 8\pi^{2}/g^{2}(\rho)$ where $\rho$ is the size of the
core and provides the relevant scale, and the probability of it existing in a
field will be
$\propto\exp\\{-S_{I}(\rho)\\}\sim\exp\\{-8\pi^{2}/g^{2}(\rho)\\}$ per unit
volume up to a $1/\rho^{4}$ volume factor and some less important factors [20,
21]. Now at large $N$ the ’t Hooft coupling $\lambda=g^{2}N$ needs to be kept
constant for a smooth limit [22, 23]. In terms of the running coupling this
means $\lambda(l)=g^{2}(l)N$ is independent of $N$ up to corrections that are
powers of $1/N$, where $l$ is measured in some physical units, e.g. the mass
gap. So the weight of a small instanton vanishes as
$\propto\exp\\{-8\pi^{2}N/\lambda(\rho)\\}$, i.e. it vanishes exponentially
with $N$ [25, 26]. Since the process of going ‘smoothly’ from a field with
$Q=1$ to a field with $Q=0$ necessarily involves at an intermediate stage of
the process the presence of fields containing such small instantons, we
immediately infer that the change in topology must be at least exponentially
suppressed with increasing $N$ once $a$ is small. Such a strong suppression
has indeed been observed in earlier lattice calculations [1, 27]. Returning to
the $a\to 0$ limit at fixed $N$ we note that one can use the one-loop
expression for $g^{2}(\rho)$ and this leads to a suppression of the tunnelling
between topological sectors that is the appropriate power of $\rho$ and hence
of $a$, once $\rho\sim O(a)$.

In our calculations we normally start a Monte Carlo sequence with the trivial
$Q=0$ gauge field $U_{l}=I,\forall l$, and ‘thermalise’ with some tens of
thousands of sweeps to reach the ‘equilibriated’ lattice gauge fields which we
can begin to use for our calculations. As described above, for larger $N$ and
smaller $a$ the topological tunnelling can become sufficiently rare that the
fields remain at $Q=0$ even after our attempted equilibriation. In practice
this proves not to be an issue for any of our $SU(2)$, $SU(3)$, $SU(4)$ and
$SU(5)$ calculations. For $SU(6)$ it becomes a problem for the smallest values
of $a$ and for $SU(N\geq 8)$ it is a problem for all but our largest values of
$a$. That is to say, within our statistics and our range of $a$ topological
freezing is largely an issue to do with $N\uparrow$ rather than with
$a\downarrow$. This distinction is important since the currently accepted
theoretical picture [28] tells us that we are sitting in a vacuum with
periodicity in the $\theta$ parameter of $2\pi N$ rather than $2\pi$. (It is
other interlaced $\theta$-vacua that provide the expected $2\pi$ periodicity.)
This implies that higher moments in $Q$ disappear as inverse powers of $N$,
and so does the correlation between the topological charge and gluonic
operators. Lattice calculations [29] have confirmed this expectation. This
tells us that at large enough $N$ the loss of ergodicity should not affect
physical quantities such as glueball masses.

Of course we do not know in advance what value of $N$ is ‘large enough’ and so
it is useful to perform an explicit check of how such freezing affects the
physics we are interested in. For this we choose to work in $SU(8)$ at
$\beta=45.50$ where changes in topology are extremely rare, but which is close
to the value $N=6$ where we begin to encounter a serious loss of topological
ergodicity in our calculations. We begin with a comparison of two sets of
about $2\times 10^{6}$ fields on a $14^{3}20$ lattice, with glueball and flux
tube ‘measurements’ taken every 25’th update. One set of fields has $Q=0$
throughout while the other has what one expects to be the correct distribution
of $Q$ imposed through about 40 different starting configurations as described
above. We calculate glueball and flux tube correlators in exactly the same way
in both ensembles of fields and perform the usual variational calculation for
each representation $R^{PC}$. Note that here we include all the states that we
obtain and make no attempt to identify any of the finite volume ditorelon
states alluded to in Section 2.3. In most, but not all, cases the effective
energy plateau in the effective energies, $aE_{eff}(t)$, defined in eqn(5)
begin at $t=2a$. So we show in Table 9 the values of $aE_{eff}(t=2a)$ that we
obtain in each of the two sets of fields, for the ground states and some of
the lowest lying excited states. We see that the two sets of values are
remarkably consistent within errors, with only a small handful of heavier
states differing by slightly more than $2\sigma$. In addition to the glueball
effective energies we also show those for the lightest two states of a flux
tube that winds around the spatial torus. We do so for the flux in the
fundamental representation ($k=1$) and for the flux in the
$(k=2)=(k=1)\otimes(k=1)$ representation. Here too we see consistency within
statistical errors. In addition to this study on the $14^{3}20$ lattice we
perform a similar one on a smaller $12^{3}20$ lattices at the same $\beta$.
This is the lattice size, in physical units, that is used in our calculations
later on in this paper, and we recall that in general the local effects of
constraining the topological charge should increase with decreasing volume
(see e.g. [30]). Again we have consistency within errors for almost all the
states. The conclusion of these explicit calculations is that within our
typical statistical errors the physical quantities we calculate in this paper
are not affected by constraining the total topological charge to $Q=0$.

Although both the theory and our above explicit checks suggest that the loss
of topological ergodicity in our calculations will have very little effect on
our measured glueball spectrum, we implement the following additional
procedure that is designed to reduce any remaining biases. Suppose we are
confronted with topological freezing on a lattice of size $L_{s}^{3}L_{t}$ at
some value of $\beta$ in some $SU(N)$ theory. Now instead of starting with a
$Q=0$ field with $U_{l}=I,\forall l$, we produce a set of sequences in
parallel which start with fields that have various values of $Q$. We choose
the ensemble of these starting fields to roughly reproduce the expected
distribution. What is ‘expected’ we infer from what we observe at the lower
values of $N$ where there is no topological freezing and also from the values
obtained at the same value of $N$ but at the larger values of $a(\beta)$ (if
any) where the freezing has not yet set in. This is a plausible strategy since
we find that in those cases where the distribution of $Q$ can be determined
there is little variation in the distribution of $Q$ with either $N$ or $a$ on
an equal physical volume. To produce such a starting field with some $Q$ we
generate a sequence of $SU(3)$ fields on the same size lattice as our $SU(N)$
one and we pick out a field with the desired charge $Q$. We ‘cool’ this field
to produce a field with minimal fluctuations apart from the net topological
charge. (See Section 6 for details.) We then embed this lattice field of
$3\times 3$ complex matrices into a suitable corner of our set of $N\times N$
unit matrices, so as to produce a set of $SU(N)$ matrices with charge $Q$. If
we have topological freezing then this starting distribution of $Q$ will be
maintained in our set of Monte Carlo sequences. This method has some
disadvantages but also some advantages. One advantage is that the lack of
topological tunneling is a better representation of the continuum theory: we
are sampling disconnected sectors. After all, the field configurations near
the point of topological tunnelling have no continuum correspondence and are a
lattice artifact. The major disadvantage is that we necessarily have a rather
limited set of parallel sequences (typically 20 or 40) and hence values of
$Q$, and so only the most probable values of $Q$ are properly represented.
Lattice fields with values of $Q$ that have a low probability simply do not
appear. That is to say we are missing the tail of the $Q$ distribution, and so
any physics that depends on higher moments of $Q$ will be something we cannot
address here. However theory tells us that the correlation between such higher
moments of $Q$ and our gluonic operators will be suppressed as a
correpondingly higher power of $N$, so any bias should be a small fraction of
what, as we have seen above, will be at most a very small shift to the
glueball masses.

In conclusion, our explicit calculations, the theoretical arguments, and our
additional procedure of imposing upon our fields something quite close to the
expected distribution of $Q$ makes us confident that within our typical
statistical errors any systematic error from the freezing of topology is
negligible in our glueball calculations.

##  3 String tensions

Let us label the sites of the hypercubic lattice by the integers
$n_{x},n_{y},n_{z},n_{t}$. To project onto a fundamental ($k=1$) flux tube the
simplest operator we can use is $\mathrm{Tr}(l_{f})$ where $l_{f}$ is the
product of link matrices along a minimal path encircling the spatial torus.
This is just the usual spatial Polyakov loop. For example, if the spatial
lattice size is $L$ then the Polyakov loop in the $x$ direction is

$l_{f}(n_{y},n_{z},n_{t})=\prod_{n_{x}=1}^{n_{x}=L}U_{x}(n_{x},n_{y},n_{z},n_{t}).$
(6)

We also use such Polyakov loops composed of blocked links and we translate the
operator along the $x$ direction. Clearly the parities of such states are
positive. This loop is translationally invariant along the flux tube direction
and so the momentum along the flux tube is zero; $p_{\parallel}=0$. It is
rotationally invariant around its axis, so $J=0$. Because of the usual centre
symmetry (see below) the $C=+$ and $C=-$ flux tubes are degenerate. We then
sum $l_{f}(n_{y},n_{z},n_{t})$ over $n_{y}$ and $n_{z}$ so that the transverse
momentum $p_{\perp}=0$. Thus the lightest state onto which this final operator
projects should be the ground state of the winding flux tube that carries
fundamental flux. To project onto a $k=2$ flux tube we do the same using
linear combinations of the operators $\mathrm{Tr}(l_{f}l_{f})$ and
$\mathrm{Tr}(l_{f})\mathrm{Tr}(l_{f})$. For $k=3$ we use linear combinations
of $\mathrm{Tr}(l_{f}l_{f}l_{f})$, $\mathrm{Tr}(l_{f}l_{f})\mathrm{Tr}(l_{f})$
and $\mathrm{Tr}(l_{f})\mathrm{Tr}(l_{f})\mathrm{Tr}(l_{f})$, and so on for
higher $k$. The defining property of a $k$-string that winds around the
spatial torus in say the $x$ direction is that under a transformation of the
link matrices
$U_{x}(n_{x},n_{y},n_{z},n_{t})\to\exp\\{i2\pi/N\\}U_{x}(n_{x},n_{y},n_{z},n_{t})$
for all $n_{y},n_{z},n_{t}$ and at some fixed value of $n_{x}$, it will
acquire a factor $\exp\\{i2k\pi/N\\}$. Since this transformation is a symmetry
of the theory (the action and measure are unchanged) and at low temperatures
the symmetry is not spontaneously broken, any two operators for which this
factor of $\exp\\{i2k\pi/N\\}$ differ will have zero mutual overlap. So we
have separate sectors of flux tubes labelled by $k$ with the fundamental $k=1$
strings existing for $N\geq 2$, the $k=2$ strings for $N\geq 4$, the $k=3$
strings for $N\geq 6$, and so on. We note that our basis of operators for a
$k$-string could easily be extended by multiplying the operators by factors
such as $\mathrm{Tr}(l_{f})\mathrm{Tr}(l_{f}^{\dagger})$, but while this might
be useful to improve the overlap onto massive excitations it would be very
surprising if it improved the overlap onto the ground state which is what we
are interested in here.

###  3.1 finite volume corrections

As described above, we calculate the ground state energy of a flux tube that
winds once around a spatial torus. If the lattice has size $l^{3}l_{t}$, this
flux tube has length $l$ and we generically denote the energy by $E_{k}(l)$,
where $k=1$ corresponds to a flux tube carrying fundamental flux and $k$ to a
$k$-string (carrying $k$ units of fundamental flux). The string tension
$\sigma_{k}$ is the energy per unit length of a very long flux tube

$\sigma_{k}=\lim_{l\to\infty}\frac{E_{k}(l)}{l}.$ (7)

In practice our tori are finite and so we expect that there will be
corrections to the leading linear dependence of $E_{k}(l)$. In our
calculations we will normally estimate the asymptotic string tension
$\sigma_{k}$ using the ‘Nambu-Goto’ formula

$E_{k}(l)=\sigma_{k}l\left(1-\frac{2\pi}{3\sigma_{k}l^{2}}\right)^{\frac{1}{2}}.$
(8)

When expanded in powers of $1/\sigma_{k}l^{2}$ this formula generates all the
known universal terms [31, 32, 33, 34] and, at least for fundamental flux,
proves to accurately describe numerical calculations for the range of $l$
relevant to our calculations [35]. However, strictly speaking, these results
hold for flux tubes that are effectively in an infinite transverse volume,
while in our case, in order to maintain the rotational symmetry needed for our
glueball calculations, the transverse size is $l\times l$ and so decreases
when we decrease $l$. We therefore need to perform some explicit checks, i.e.
to check whether eqn(8) encodes all the finite-$l$ corrections that are
visible within our typical statistical errors.

We see from eqn(8) that the relevant scale is $l$ in units of the string
tension, i.e. $l\surd\sigma$. These scales are listed in the headers of Tables
1-8 and they vary from $l\surd\sigma\sim 4$ for lower $N$ to $l\surd\sigma\sim
2.6$ for our highest values of $N$. As far as glueball masses are concerned
the theoretical justification for this decrease in the volume is that we
expect finite volume corrections to decrease as $N$ grows and the practical
reason is that the expense of the $SU(N)$ matrix product calculations grows as
$\propto N^{3}$. In order to check whether there are significant corrections
to eqn(8) at the values of $l$ we use, we calculate $\sigma$ on larger
(sometimes smaller) lattices at selected values of $\beta$. These calculations
are summarised in Table 10. For $N=8,10,12$ the scale $l\surd\sigma\sim 2.6$
corresponds to $l=12$ in the Table, and we see that the value of
$a\surd\sigma$ is the same for the longer $l=14$ flux tube, suggesting that
there are no significant finite $l$ corrections to eqn(8) for these values of
$N$. The same conclusion holds for the other groups in Table 10. All this
provides evidence that the fundamental string tensions quoted in Tables 1-8 do
not suffer significant finite volume errors, at least within our typical
statistical errors.

The energies of $k\geq 2$ flux tubes are larger and therefore so are the
statistical and systematic errors in estimating them. Since the energies grow
with $k$, we have only attempted a finite volume study for $k=2$ flux tubes.
This is presented in Table 11. The message here is more nuanced than for
$k=1$: while any finite $l$ corrections to eqn(8) appear to be insignificant
when $l\surd\sigma\gtrsim 3$, and hence for our $N\leq 6$ string tension
calculations, there is strong evidence that there is a correction of $\simeq
1.9(5)\%$ to the values obtained with $l\surd\sigma\sim 2.6$, which are the
values we emply for $N\geq 8$. This is something we will need to consider when
we estimate the errors on our $k=2$ string tensions. Presumably the $k=3$ and
$k=4$ string tensions show corrections at least as large, but we do not have
any finite $l$ study of these.

###  3.2 fundamental ($k=1$) string tensions

We calculate the ground state energies of the fundamental flux tube for the
lattices, couplings and gauge groups listed in Tables 1-8. We then use eqn(8)
to extract our estimates of the infinite volume string tensions. The energies
are obtained by identifying effective energy plateaux as described in Section
2.2. As an example we show in Fig.1 the effective energies that lead to the
$SU(8)$ string tension estimates in Table 6. As one increases $\beta$ and
hence decreases $a$, the energy in lattice units decreases so that we can
obtain precise effective energies over a larger range of $n_{t}$ which in turn
makes it easier to identify a plateau and so estimate the corresponding value
of $E_{eff}(n_{t}\to\infty)$. Although we attempt to show in Fig.1 the error
bands on our estimates of the corresponding flux tube energies, these are
almost invisible because the errors are too small on the scale used in the
plot. We therefore replot, in Fig.2, the effective energies for $\beta=47.75$,
which is the calculation that is the closest to the desired continuum limit,
with a sufficient rescaling to expose the errors on the effective masses. On
the plot we show the best estimate for the energy obtained from a fit to the
correlation function, together with the energies corresponding to $\pm 1$
standard deviations.

For the corresponding $SU(3)$ plots we refer the reader to Fig.1 of [11]. The
fact that the $SU(3)$ lattices are (much) larger than the $SU(8)$ ones at
similar values of $a\surd\sigma$ means that the flux tubes are more massive,
so the correlators decrease faster with $n_{t}$ and the extraction of
$E_{eff}(n_{t}\to\infty)$ is less compelling – despite the fact that the
calculations extend closer to the continuum limit in the case of $SU(3)$.
These $SU(3)$ plots are typical of our $N\leq 4$ calculations, where the
lattices are chosen to be particularly large, while the $SU(8)$ plot is
typical of $N\geq 8$ where the spatial volumes used are smaller.

In the case of $SU(2)$ and $SU(3)$ we have also calculated the string tension
on smaller volumes which are still large enough for the flux tube calculations
even if too small for reliable glueball calculations. This enables us to
extend the calculations to smaller values of $a\surd\sigma$ at modest
computational cost. The purpose of these calculations is to feed into our
analysis of the running coupling later on in this paper. We list the results
of the calculations in Tables 12,13.

###  3.3 $k=2$ string tensions

One expects that as $N\to\infty$ the $k=2$ flux tube will become two non-
interacting $k=1$ flux tubes and we will have $\sigma_{k=2}\to 2\sigma_{k=1}$.
However, as earlier calculations have shown [36, 35, 27], at lower values of
$N$ one finds that $\sigma_{k=2}$ is substantially less than $2\sigma_{k=1}$,
so that one can think of it as behaving like a bound state of two fundamental
flux tubes. This suggests that at lower values of $N$ the $k=2$ flux tube can
be treated as a single ‘string’ with finite volume corrections well described
by eqn(8). As we increase $N$ the $k=2$ flux tube should increasingly look
like two fundamental flux tubes that are loosely bound with the binding energy
vanishing as $N\to\infty$. Here one might expect each of these two flux tubes
to have finite volume corrections given by eqn(8) so that the overall finite
volume behaviour of the energy of a $k=2$ flux tube becomes significantly
different:

$E_{k=2}(l)\stackrel{{\scriptstyle\mathrm{small}\,N}}{{=}}\sigma_{k=2}l\left(1-\frac{2\pi}{3\sigma_{k=2}l^{2}}\right)^{\frac{1}{2}}\stackrel{{\scriptstyle
N\to\infty}}{{\longrightarrow}}\sigma_{k=2}l\left(1-\frac{4\pi}{3\sigma_{k=2}l^{2}}\right)^{\frac{1}{2}}$
(9)

with $\sigma_{k=2}=2\sigma_{k=1}$. One can expect a smooth transition between
these two behaviours which clearly creates ambiguities in extracting the
string tension $a^{2}\sigma_{k=2}$ from the flux tube energy $aE_{k=2}(l)$.

Assuming for now that we can treat the $k=2$ flux tube as a single string with
finite volume corrections as given by eqn(8) we calculate $\sigma_{k=2}$ on
our various ensembles of lattice fields. Using the previously calculated
values of $\sigma_{f}\equiv\sigma_{k=1}$ we form the dimensionless ratio
$\sigma_{k=2}/\sigma_{f}$ which we then extrapolate to the continuum limit
using eqn(4). We show the extrapolations in Fig. 3. For $N\leq 6$ these appear
to be under good control, but for larger $N$ that is less clear. One would
expect the slope of the extrapolation to vary smoothly with $N$, and that is
certainly what one sees for $N=4,5,6$. But then there is a violent break in
the behaviour between $N=6$ and $N=8$ and wild oscillations when comparing the
$N=8,10,12$ slopes. It may be that this is due to our use of smaller spatial
volumes for $N\geq 8$. In any case, we calculate the continuum limits from
these fits and the results for our various gauge groups are presented in the
second column of Table 14. As remarked earlier, there is good evidence from
Table 11 that for the smaller lattices that we have used for $N\geq 8$ we
should apply an additional finite volume correction of $\simeq 1.9(5)\%$.
Doing this leads to the values in the third column of Table 14, which are the
values we consider to be more reliable. There is a caveat here: at our largest
values of $N$ the values of $\sigma_{k=2}/\sigma_{f}$ are quite close to the
asymptotic value of two, and so we might expect that we are in the range of
$N$ where the stronger finite volume behaviour displayed in eqn(9) is setting
in. To settle this question would require a dedicated calculation that is
beyond the scope of this paper.

We have also obtained continuum extrapolations of the $k=3$ and $k=4$ string
tensions. These have been obtained from the corresponding flux tube energies
using eqn(8) and this undoubtedly underestimates the finite volume
corrections. Indeed as $N\to\infty$ we expect $\sigma_{k}\to k\sigma_{f}$ as
the $k$-string becomes $k$ non-interacting fundamental strings, and then we
expect a version of eqn(9),

$E_{k}(l)\stackrel{{\scriptstyle
N\to\infty}}{{\longrightarrow}}k\sigma_{f}l\left(1+\frac{2\pi}{3\sigma_{f}l^{2}}\right)^{\frac{1}{2}}=\sigma_{k}l\left(1+\frac{2k\pi}{3\sigma_{k}l^{2}}\right)^{\frac{1}{2}}.$
(10)

Given this uncertainty and the increasing ambiguity with increasing $k$ of
extracting an effective energy plateau (because of their increasing energies),
we do not discuss these values any further here.

###  3.4 $N\to\infty$ extrapolations

An interesting question about the ratio $\sigma_{k=2}/\sigma_{f}$ concerns its
dependence on $N$. On general grounds we expect the leading correction to the
asymptotic value to be $O(1/N^{2})$. However there are old ideas under the
label of ‘Casimir Scaling’, that the (string) tension of a flux tube should be
proportional to the smallest quadratic Casimir of the representations that
contribute to that flux, as though the flux tube joining two sources behaved,
in this respect, just like one gluon exchange. For the $k$ flux tube the
relevant representation is the totally antisymmetric one and this predicts

$\frac{\sigma_{k}}{\sigma_{f}}\stackrel{{\scriptstyle
CS}}{{=}}\frac{k(N-k)}{N-1}\stackrel{{\scriptstyle
N\to\infty}}{{=}}k-\frac{k(k-1)}{N}-\frac{k(k-1)}{N^{2}}+O\left(\frac{1}{N^{3}}\right)$
(11)

so that the leading correction is $O(1/N)$ rather than $O(1/N^{2})$. Previous
lattice calculations have often favoured Casimir scaling as a good
approximation [39] so it is interesting to test this idea against our $k=2$
values listed in Table 14. In Fig.4 we plot our continuum values of
$\sigma_{k=2}/\sigma_{f}$ against $1/N^{2}$. In addition we know that
$\sigma_{k=2}/\sigma_{f}=2$ at $N=\infty$ so we impose this as a constraint in
our fits. In Fig.4 we show our best fit in powers of $1/N^{2}$

$\frac{\sigma_{k=2}}{\sigma_{f}}=2.0-\frac{1.28(19)}{N}-\frac{4.78(90)}{N^{2}}$
(12)

with a $\chi^{2}$ per degree of freedom of $\sim 0.5$ and an alternative best
fit in powers of $1/N^{2}$

$\frac{\sigma_{k=2}}{\sigma_{f}}=2.0-\frac{14.43(60)}{N^{2}}+\frac{73.8(12.1)}{N^{4}}$
(13)

with a $\chi^{2}$ per degree of freedom of $\sim 2.2$. The former fit is
clearly better, but the latter cannot be entirely excluded. On the other hand
while the fit in powers of $1/N$ is very good, the coefficients are very
different from those obtained if we expand the Casimir scaling prediction, as
in eqn(11). One may speculate that for the range of $N$ where the $k=2$ flux
tube is a bound state the $N$ dependence of the string tension is best
described in powers of $1/N$ while once $N$ is large enough that it has become
a weakly interacting pair of fundamental flux tubes, the dependance will be
best described in powers of $1/N^{2}$. At which $N$ this transition occurs
will depend on the length $l$ of the flux tube. If we use the formula in
eqn(8) for the fundamental and $k=2$ flux tubes and if we assume that
$\sigma_{k=2}=2\sigma-O(1/N)$, then one finds that the lightest state will be
the weakly interacting pair of fundamental flux tubes for $l\leq l_{c}$ where
$l_{c}\surd\sigma\propto N^{\frac{1}{2}}$. (This essentially reflects a
competition between the leading linear terms and the $O(1/l)$ Luscher
corrections.) That is to say, for any fixed length $l$, the (very) asymptotic
leading correction is $O(1/N^{2})$ as expected by the usual large-$N$
counting. However one has to be careful in which order one takes the
$l\to\infty$ and $N\to\infty$ limits.

##  4 Running coupling

A question of theoretical interest is whether our calculations show that we
should keep the ’t Hooft coupling $g^{2}N$ fixed as $N\to\infty$ in order to
have a smooth large-$N$ limit. A question of phenomenological interest is
whether we can estimate the scale $\Lambda$ of our $SU(N)$ gauge theories,
from the calculated running of the gauge coupling, and in particular whether
we can do so for $SU(3)$. Finally something that is often useful in lattice
calculations is to have an interpolation function for $a(\beta)$. These are
the three issues we address in this section. A summary of the technical
background to the perturbative calculations has been placed in Appendix A.

###  4.1 scaling with N

Since the lattice coupling $g^{2}_{L}$ defined in eqn(2) provides a definition
of the running coupling on the scale $a$, and our above calculation of the
(fundamental) string tension at various $\beta$ enables us to express $a$ in
physical units, i.e. as $a\surd\sigma$, we can use these calculations to
address some questions about the properties of the running coupling in $SU(N)$
gauge theories. Before doing so we recall that this lattice coupling,
corresponding to the particular coupling scheme defined by the lattice and the
plaquette action, is well-known to be a ‘poor’ definition of a running
coupling in the sense that higher order corrections will typically be very
large. This is indicated, for example, by the relationship between the scale
parameters $\Lambda_{L}$ and $\Lambda_{\overline{MS}}$ in this scheme and the
standard ${\overline{MS}}$ scheme [37, 38]:

$\frac{\Lambda_{\overline{MS}}}{\Lambda_{L}}=38.853\exp\left\\{-\frac{3\pi^{2}}{11N^{2}}\right\\}.$
(14)

This is a long-standing issue that has led to the formulation of a number of
improved couplings. (For a review see [40].) Here we shall use the ‘mean-
field’ improved coupling of Parisi [41],

$\frac{1}{g^{2}_{I}}=\frac{1}{g^{2}_{L}}\langle\frac{1}{N}\mathrm{Tr}U_{p}\rangle,$
(15)

which has a nice physical motivation as the effective coupling experienced by
a background field (in a simple approximation). Denoting the corresponding
scale by $\Lambda_{I}$, one finds

$\frac{\Lambda_{\overline{MS}}}{\Lambda_{I}}=\frac{\Lambda_{L}}{\Lambda_{I}}\frac{\Lambda_{\overline{MS}}}{\Lambda_{L}}=\exp\left\\{-\frac{w_{1}}{2b_{0}}\right\\}\times
38.853\exp\left\\{-\frac{3\pi^{2}}{11N^{2}}\right\\}\simeq 2.633$ (16)

using eqn(14) and eqns(58,63), with the value of $\Lambda_{L}/\Lambda_{I}$
being obtained using eqn(53) for each $\Lambda$, with $g^{2}_{I}$ and
$g^{2}_{L}$ related by eqns(62,63), and then taking $g^{2}\to 0$. This already
suggests that $g^{2}_{I}$ has the potential to be more-or-less as good as
$g^{2}_{\overline{MS}}$. One can also show [42] that this coupling tracks
quite accurately the Schrodinger-functional coupling [43, 44] over a very wide
range of scales. (See [42] for a detailed discussion.) So from now on, in this
section, we shall use $g^{2}_{I}(a)$ as our lattice running coupling.

The main question we address concerns the $N$ dependence of $g^{2}_{I}(a)$.
The usual large $N$ counting tells us that we expect to approach constant
physics as $N\to\infty$ if we keep $g^{2}N$ fixed. For the running coupling
this means that if we plot $g^{2}_{I}(a)N$ against the calculated values of
$a\surd\sigma$ in our $SU(N)$ gauge theories, they should approach a common
envelope as $N\to\infty$. The interesting question is whether the approach is
slow or fast, indicating that our values of $N$ are ‘far from’ or ‘near to’
$N=\infty$ respectively. In Fig. 5 we plot our calculated values of the
running ’t Hooft coupling $g^{2}_{I}(a)N$ against $a\surd\sigma$ for all our
values of $N$. The results are quite remarkable: even $SU(2)$ is very close to
$SU(\infty)$ in this respect, and it is only when the lattice spacing becomes
large that appreciable differences appear. It is also interesting to see what
happens if one uses the poor $g^{2}_{L}(a)$ coupling instead of
$g^{2}_{I}(a)$, and one finds that while the convergence to a large $N$ limit
is still evident, the corrections at lower $N$ are substantial indicating
that, not surprisingly, higher order non-planar contributions are important
using this coupling scheme. Finally we comment that these conclusions are not
unexpected or novel: similar analyses have appeared in for example [1, 42],
albeit usually with less accuracy and over more limited ranges of $N$ and
$a\surd\sigma$.

###  4.2 perturbative running and ${\mathrm{\Lambda_{\overline{MS}}}}$

The second question we address is whether the running coupling dependence
displayed in Fig. 5 can be described by the usual perturbative
$\beta$-function for at least some $N$ once $a$ is small. If so we can extract
a value of $\Lambda_{I}$ for each such $N$ and a corresponding value of
$\Lambda_{\overline{MS}}$ using eqn(16). Our analysis will broadly follow that
of [42] and we refer the reader to that paper for background and context; here
we merely outline the calculation, with more details in Appendix A. To fix our
notation we begin with the standard $\beta$-function for the lattice bare
coupling $g_{I}^{2}$:

$\beta(g_{I})=-\frac{\partial g^{2}_{I}}{\partial\log
a^{2}}=-b_{0}g^{4}_{I}-b_{1}g^{6}_{I}-b^{I}_{2}g^{8}_{I}+...+O(a^{2}),$ (17)

where the scheme independent coefficients $b_{0},b_{1}$ and the scheme
dependent $b^{I}_{2}$ are given in eqn(58) and eqn(64) of Appendix A. As shown
in Appendix A this motivates the following 3-loop expression for $a$:

$a\sqrt{\sigma}(a)\stackrel{{\scriptstyle
3loop}}{{=}}\frac{\sqrt{\sigma}(0)}{\Lambda_{I}}\left(1+c_{\sigma}a^{2}\sigma\right)\left(b_{0}g_{I}^{2}(a)\right)^{-\frac{b_{1}}{2b^{2}_{0}}}e^{-\frac{1}{2b_{0}g_{I}^{2}(a)}}e^{-\frac{1}{2}\int^{g_{I}^{2}(a)}_{0}dg^{2}\left(\frac{b_{0}b^{I}_{2}-b^{2}_{1}-b_{1}b^{I}_{2}g^{2}}{b^{3}_{0}+b^{2}_{0}b_{1}g^{2}+b^{2}_{0}b^{I}_{2}g^{4}}\right)}.$
(18)

Here the first factor on the right after the coefficient is a non-perturbative
tree-level lattice spacing correction. The motivation is that we could just as
well use some physical mass, $\mu(a)$, as a scale in place of $\surd\sigma(a)$
and since $\mu(a)/\surd\sigma(a)=\mu(a=0)/\surd\sigma(a=0)(1+O(a^{2}))$ it
must be the case that we have, in general, a factor $(1+O(a^{2}))$ multiplying
the perturbative expression [42, 45]. Here we choose to use $\sigma$ as our
scale for $a^{2}$: using some other $\mu^{2}$ would only change the $O(a^{4})$
term and we will assume that $a$ is small enough that we can neglect any
$O(a^{4})$ terms. We evaluate numerically the integral in the exponential for
any given value of $g_{I}(a)$: for such a smooth integrand any simple
technique will be able to give accurate results. We now describe the result of
fitting this function to our values of $a\surd\sigma$, as listed in Tables
12,13 for $SU(2)$ and $SU(3)$ and in Tables 3–8 for $SU(4)$ to $SU(12)$. In
performing these fits one can either use the measured values of $a^{2}\sigma$
in the $(1+c_{\sigma}a^{2}\sigma)$ lattice correction factor, or the value
calculated from the formula itself. In the latter case eqn(18) becomes a
quadratic equation for $a\surd\sigma$. Not surprisingly, for the range of
$g_{I}^{2}(a)$ where eqn(18) provides a good fit to our calculated string
tensions, the difference between the two methods is insignificant, and we
choose here to use the measured value in the correction term.

The results of our fits using eqn(18) are shown in Table 15. For each value of
$N$ we fit our lattice values of $a\surd\sigma$ discarding the largest values
until we obtain a reasonably acceptable fit. The range of lattice spacings for
the fits is listed for each $N$, together with the $\chi^{2}$ per degree of
freedom, $\chi^{2}/n_{df}$. In the few cases labelled by $\ast$ in Table 15
the fits are very poor with $\chi^{2}/n_{df}\geq 3$, and so we quadruple the
stated error in those cases in the hope that this encodes the increased
uncertainty. Referring to Tables 1–8 we see that the number of values used
within a given fit varies from 6 in $SU(3)$ to 4 in $SU(10)$ and $SU(12)$.
Since we have 2 parameters in our fit, 4 points is a very small number for a
fit. Moreover, for $SU(10)$ and $SU(12)$ the lattice spacings do not extend to
values as small as for lower $N$. Hence one should treat the resulting fits
with extra caution.

In Table 15 we also list the fitted values of the perturbative scale
parameter, $\Lambda_{I}$, in units of the continuum string tension. We then
convert this to the corresponding $\Lambda_{\overline{MS}}$ scale in the
widely used $\overline{MS}$ coupling scheme, using the $N$-independent
relation $\Lambda_{\overline{MS}}/\Lambda_{I}\simeq 2.633$ from eqn(16) and
the results are listed in the last column of Table 15, again in units of the
string tension. The quoted errors on our values of $\Lambda_{I}$ are very
small, but they are purely statistical and their smallness reflects the very
small errors on the string tensions that are fitted. The systematic errors
involved in, for example, the truncation of the beta-function may be much
larger. To get some measure of this error we also perform an extra fit using
the exact 2-loop running,

$a\sqrt{\sigma}(a)\stackrel{{\scriptstyle
2loop}}{{=}}c_{I}\left(1+c_{\sigma,I}a^{2}\sigma\right)e^{-\frac{1}{2b_{0}g_{I}^{2}}}\left(\frac{b_{1}}{b_{0}^{2}}+\frac{1}{b_{0}g_{I}^{2}}\right)^{\frac{b_{1}}{2b_{0}^{2}}}.$
(19)

The resulting 2-loop values of $\Lambda_{I}$ are also listed in Table 15. We
see that $\Lambda^{2loop}_{I}$ is typically about $10\%$ smaller than the
$\Lambda^{3loop}_{I}$. We have decided to use half of this difference as a
measure of the error associated with dropping 4-loop and higher coefficients
from the 3-loop calculation. This is added as a second error, within square
brackets, to the resulting 3-loop values for $\Lambda_{\overline{MS}}$ that
are listed in the last right-hand column of Table 15. Unlike the statistical
errors, this error will affect the results for all values of $N$ in the same
direction.

In Fig. 6 we display the values that we obtain from our 3-loop fits for
$\Lambda_{\overline{MS}}/\surd\sigma$ as a function of $N$. (We use only the
statistical errors in these fits.) It is clear that the variation with $N$ is
very weak. If we fit all our values for $N\leq 12$ we find the fit,

$\frac{\Lambda_{\overline{MS}}}{\surd\sigma}=0.5055(7)[250]+\frac{0.306(12)}{N^{2}},\quad\chi^{2}/n_{df}=2.70,\qquad
N\in[2,12],$ (20)

which is displayed in Fig. 6. Although the $\chi^{2}/n_{df}$ is not good, the
calculated values appear to be scattered around the fit in a random pattern.
Other possible fits include

$\displaystyle\frac{\Lambda_{\overline{MS}}}{\surd\sigma}$ $\displaystyle=$
$\displaystyle
0.5067(11)[250]+\frac{0.258(12)}{N^{2}},\quad\chi^{2}/n_{df}=1.65,\qquad
N\in[4,12],$ $\displaystyle\frac{\Lambda_{\overline{MS}}}{\surd\sigma}$
$\displaystyle=$ $\displaystyle
0.5060(8)[250]+\frac{0.303(13)}{N^{2}},\quad\chi^{2}/n_{df}=3.81,\qquad
N\in[2,8],$ $\displaystyle\frac{\Lambda_{\overline{MS}}}{\surd\sigma}$
$\displaystyle=$ $\displaystyle
0.5093(16)[250]+\frac{0.206(38)}{N^{2}},\quad\chi^{2}/n_{df}=0.44,\qquad
N\in[4,8].$ (21)

In all the above fits we have included within square brackets an estimate of
$\sim 5\%$ for the systematic error associated with the truncation of the
perturbative expansion. We note that this systematic error is much larger than
the differences between the various fits above, and is the dominant source of
uncertainty.

To finish we turn briefly to the case of $SU(3)$ where the value of
$\Lambda_{\overline{MS}}$ has phenomenological interest. Here we can attempt
to transform our value into $\mathrm{MeV}$ units as follows. We begin with its
value in units of the string tension, as listed in Table 15, and we then
transform this into a value in terms of the Sommer length scale $r_{0}$ [46]
using a recent calculation [11] that gives $r_{0}\surd\sigma=1.160(6)$. A
recent review [47] of calculations in lattice QCD with light quarks concludes
that $r_{0}=0.472(5)\mathrm{fm}=1/418(5)\mathrm{MeV}^{-1}$. The usual
expectation is that the value of $r_{0}$ is not very sensitive to the
inclusion of light quarks, so we use this $\mathrm{MeV}$ value in the pure
gauge theory. (This is of course the arguable step.) Doing so we arrive at

$\left.\frac{\Lambda_{\overline{MS}}}{\surd\sigma}\right|_{SU3}=0.5424(13)[185]\Longrightarrow
r_{0}\Lambda_{\overline{MS}}=0.629(4)[22]\Longrightarrow\Lambda_{\overline{MS}}\stackrel{{\scriptstyle
SU3}}{{=}}263(4)[9]\mathrm{MeV},$ (22)

where the first error is statistical and the second is systematic. This value
is consistent with the values recently obtained in the dedicated calculations
[48, 49] that use very different methods, and our value has a similar
accuracy. This adds confidence that our calculations of
$\Lambda_{\overline{MS}}$ for the other $SU(N)$ groups are also reliable.

###  4.3 interpolating and extrapolating functions for $a(\beta)$

It can often be useful to know the value of $a(\beta)$, in physical units, at
some value of $\beta$. This can be provided, for example, by the value of
$a\surd\sigma$. However calculations of $a\surd\sigma$ are obtained at a
number of discrete values of $\beta$ within some finite range
$\beta_{0}\leq\beta\leq\beta_{1}$ and the $\beta$ value of interest may lie
outside this range, or may be within this range but not at one of the discrete
values where $a\surd\sigma$ has been calculated. In the latter case one needs
to find an interpolating function that will work in the range
$\beta\in[\beta_{0},\beta_{1}]$ and since $a\surd\sigma$ typically varies
smoothly in this range this is easy to do: a few sensibly chosen terms from
almost any complete set of functions will work adequately. However if the
$\beta$ value of interest lies outside the range $[\beta_{0},\beta_{1}]$ and,
in particular, if it is at some weaker coupling, then extrapolating such a
‘random’ interpolating fit will invariably work badly unless the original
interpolating function has a form motivated by weak coupling perturbation
theory, in which case it should (in principle) be reasonably accurate. In the
section above we have used precisely such functions. We will here present a
fit in an explicit form so that the reader can readily employ it.

Our fit will be in terms of the mean-field improved coupling $g^{2}_{I}$
defined in eqn(15). This requires that one know the value of the average
plaquette; this is always calculated in a Monte Carlo but one needs to record
the average, which is normally done so as to provide a first check on the
calculation. (And if not, one can obtain a value with adequate accuracy very
quickly using small lattices and modest statistics.) The interpolating
function we will use is a variation on the ones used above and is as follows:

$a\sqrt{\sigma}(a)=c_{0}\left(1+c_{\sigma}a^{2}\sigma(a)\right)F_{2l+}(g_{I}),$
(23)

where

$F_{2l+}(g_{I})=e^{-\frac{1}{2b_{0}g_{I}^{2}}}\left(\frac{b_{1}}{b_{0}^{2}}+\frac{1}{b_{0}g_{I}^{2}}\right)^{\frac{b_{1}}{2b_{0}^{2}}}e^{-\frac{b^{I}_{2}}{2b_{0}^{2}}g_{I}^{2}}.$
(24)

The first two factors of $F_{2l+}(g_{I})$ constitute the exact dependence when
$\beta(g)$ is truncated to the first two terms. The last factor on the right
is the extra dependence if one keeps in the exponent the leading
$O(b^{I}_{2})$ term to $O(g^{2}_{I})$. That is to say, it is ‘more’ than 2
loops, but ‘less’ than 3 loops. Hence the subscript $2l+$ on $F(g)$. This is
of course an arbitrary truncation with no guarantee that it does indeed do
better than the 2-loop one; however we use it because, as we shall see below,
it turns out to be very close to the 3-loop result, but without the need to
perform any numerical integrations. Note that here we choose to use in the
$\propto a^{2}\sigma$ correction term the value given by fitting the formula
to our data since this is the mode in which it needs to be used in an
extrapolation.

To make use of eqn(24) we need to fit the constants $c_{0}$ and $c_{\sigma}$
to our calculated values of $a\surd\sigma$, at each $N$. Once the constants
$c_{0}$ and $c_{\sigma}$ have been fitted we can solve eqn(23), which is
quadratic in $a\surd\sigma$, at any value of $\beta_{I}$:

$a\sqrt{\sigma}(a)=\frac{1}{2c_{0}c_{\sigma}F_{2l+}(\beta_{I})}\left(1-\left[1-4c^{2}_{0}c_{\sigma}F_{2l+}(\beta_{I})^{2}\right]\right).$
(25)

In eqn(24) the values of $b_{0},b_{1},b^{I}_{2}$ need to be specified. The
values of $b_{0}$ and $b_{1}$ are universal and are given in eqn(58). The
value of $b^{I}_{2}$ is as given in eqn(64)

$b^{I}_{2}=b^{L}_{2}+w_{2}b_{0}-w_{1}b_{1},$ (26)

where we use eqns(59,61) in eqn(60) to obtain the explicit expression for
$b^{L}_{2}$ and then inserting that together with the functions in eqn(63) and
eqn(58) into eqn(26) we obtain the explicit expression for $b^{I}_{2}$ for any
$N$.

It only remains now to give our fitted values of $c_{0}$ and $c_{\sigma}$ for
each of the $SU(N)$ groups for which we have calculated the string tension.
This we do in Table 16. There we show for each $N$ our best fits to these two
parameters as well as the fitted range and the $\chi^{2}$ per degree of
freedom of the fit. In performing these fits we systematically drop the
largest values of $a\surd\sigma$ from the fit until we obtain an acceptable
$\chi^{2}$. This is appropriate since our interpolating function is based on
weak-coupling perturbation theory. Note that being based on weak coupling the
resulting function is not designed to work for values of $a$ or $g^{2}$
greater than the range within which it provides a good fit.

As an aside it needs to be emphasised that all the above fits are only
relevant if one uses the Wilson plaquette action. As a second aside, we note
that if one wants $a(\beta)$ in terms of some physical scale $\mu$ instead of
$\surd\sigma$ then one can straighforwardly modify the above formula for
$a\surd\sigma$ to one for $a\mu$ by using the relation
$\mu/\surd\sigma=d_{0}+d_{1}a^{2}\sigma$ if it accurately holds in the range
of couplings of interest.

##  5 Glueball masses

###  5.1 quantum numbers

The glueballs are colour singlets and so our glueball operator is obtained by
taking the ordered product of $SU(N)$ link matrices around a contractible loop
and then taking the trace. To retain the exact positivity of the correlators
we use loops that contain only spatial links. The real part of the trace
projects on $C=+$ and the imaginary part on $C=-$. We sum all spatial
translations of the loop so as to obtain an operator with momentum $p=0$. We
take all rotations of the loop and construct the linear combinations that
transform according to the irreducible representations, $R$, of the rotational
symmetry group of our cubic spatial lattice. We always choose to use a cubic
lattice volume that respects these symmetries. For each loop we also construct
its parity inverse so that taking linear combinations we can construct
operators of both parities, $P=\pm$. The correlators of such operators will
project onto glueballs with $p=0$ and the $R^{PC}$ quantum numbers of the
operators concerned.

The irreducible representations $R$ of our subgroup of the full rotation group
are usually labelled as $A_{1},A_{2},E,T_{1},T_{2}$. The $A_{1}$ is a singlet
and rotationally symmetric, so it will contain the $J=0$ state in the
continuum limit. The $A_{2}$ is also a singlet, while the $E$ is a doublet and
$T_{1}$ and $T_{2}$ are both triplets. In Section 5.7 we will outline the
detailed relationship between these representations and the continuum spin
$J$. Since, for example, the three states transforming as the triplet of
$T_{2}$ are degenerate on the lattice, we average their values and treat them
as one state in our tables of glueball masses and we do the same with the
$T_{1}$ triplets and the $E$ doublets. (Just as we would treat the 5 states of
a continuum $J=2$ glueball as one entry.)

###  5.2 finite volume effects

For reasons of computational economy we wish to calculate on lattice sizes
that are small but, at the same time, large enough that any finite volume
corrections remain smaller than our typical statistical errors. Since the
computional cost of calculating in $SU(N)$ gauge theories grows roughly
$\propto N^{3}$ (the multiplication of two $N\times N$ matrices) and since
finite volume corrections are expected to decrease as powers of $1/N$, we
reduce the size in physical units of our lattices as we increase $N$, as shown
in Tables 1-8.

There are two important types of finite volume corrections. The first can be
thought of as arising when the propagating glueball emits a virtual glueball
which propagates around the spatial torus. The resulting shift in the mass of
the propagating glueball decreases exponentially in $m_{G}l$ where $m_{G}$ is
the mass gap and $l$ is the length of the spatial torus [19]. As we see from
Tables 1-8 the value of $am_{G}\times l/a$ is quite large in all of our
calculations, so we can expect this correction to be small.

The second type of finite volume correction consists of states composed of
multiple flux tubes winding around a spatial torus in a (centre) singlet
state. The lightest of these will be a state composed of one winding flux tube
together with a conjugate winding flux tube, which we refer to as a
‘ditorelon’. (A single winding flux tube is usually referred to as a
‘torelon’.) Since it can have a non-zero overlap onto the contractible loops
that we use as our glueball operators, it can appear as a state in our
calculated glueball spectrum. Neglecting interactions, the lightest ditorelon
will consist of each flux tube in its ground state with zero momentum and will
have an energy, $E_{d}$, that is twice that of the flux tube ground state,
$E_{d}=2E_{f}$. Interactions will shift the energy but this shift should be
small on our volumes so we shall use $E_{d}\simeq 2E_{f}$ as a rough guide in
searching for these states. This ground state ditorelon has simple rotational
properties and only contribute to the $A_{1}^{++}$ and $E^{++}$
representations. If we allow one or both of the component flux tubes to be
excited and/or to have non-zero equal and opposite transverse momenta we can
populate other representations and produce towers of states. However these
excited ditorelon states will be considerably heavier on the lattice volumes
we employ and so we will not consider them any further in this paper, although
they certainly warrant further study.

The first of the above corrections leads to small shifts in the masses of the
glueballs. The second leads to extra states in the glueball spectrum. The
signature of such an extra ditorelon state is that its mass grows roughly
linearly with the lattice size: $aE_{D}\simeq 2aE_{f}\simeq 2a^{2}\sigma_{f}L$
where $L$ is the relevant spatial size in lattice units and $\sigma_{f}$ is
the (fundamental) string tension. So to test for finite volume effects we
perform calculations at the same value of $\beta$ on different lattice sizes
and compare the glueball spectra. To identify any ditorelon states we look for
extra states in the $A_{1}^{++}$ and $E^{++}$ spectra whose masses increase
roughly linearly with the volume. Since the mass shift associated with the
first kind of correction decreases exponentially with the lattice size any
shift on a significantly larger volume should be much smaller than on the
smaller volume; so to check that it is negligible compared to our statistical
errors we simply compare the masses of the states that are not ditorelons on
the different volumes. As we see from Tables 1-8 the lattice sizes we use fall
into three groups: $l\surd\sigma\sim 4.0$ for $SU(2),SU(3),SU(4)$;
$l\surd\sigma\sim 3.1$ for $SU(5),SU(6)$; $l\surd\sigma\sim 2.6$ for
$SU(8),SU(10),SU(12)$. We will exhibit a finite volume analysis for a
representative of each of these three groups.

We begin with $SU(2)$. In Table 17 we list the low-lying glueball spectra that
we obtain on 3 lattice sizes, together with the energies of the ground and
first excited states of the winding flux tubes. In physical units the
$L=12,14,20$ spatial lattice sizes correspond to $l\surd\sigma\simeq
2.9,3.4,4.8$ respectively and we recall that the typical lattice size we use
in $SU(2)$ is $l\surd\sigma\simeq 4$. For the glueball states we list the
effective energy at $t=2a$ or, where applicable, the energy of the fit to the
effective energy plateau when that begins at $t=2a$. (These two measures
differ very little in practice.) These are obtained from the correlators of
our variationally selected best operators. For most of our lighter states the
value of $aE_{eff}(t=2a)$ is very close to our best estimate of the mass. For
these finite volume comparisons we prefer this measure to the mass itself
because it serves to minimise the statistical errors and makes any finite
volume corrections more visible. Putting aside the $R^{P}=A_{1}^{+}$ and the
$E^{+}$ spectra for the moment, we see that most of the glueball energies on
the $L=14$ and $L=20$ lattices are consistent within one standard deviation
and all within two standard deviations. The energies from the $L=12$ lattice
are broadly consistent but there are now some examples, such as the
$A_{2}^{+}$ and $E^{-}$ ground states, where there appear to be significant
differences. We now return to the $A_{1}^{+}$ and the $E^{+}$ spectra listed
in Table 17 to see if there is any evidence of the extra ditorelon states. We
do indeed see these on the $L=12$ and $L=14$ lattices in both the $A_{1}^{+}$
and the $E^{+}$ spectra: these states are displayed in the Table with no
corresponding entries at other lattice sizes. Their effective energies
increase with $L$ and are just a little heavier than twice the flux tube
energy. For the $L=20$ lattice any ditorelon state would be much heavier than
any of the states shown: at such energies the spectrum is denser, the errors
are larger, and so identifying an ‘extra’ state becomes ambiguous. The
remaining states in the $A_{1}^{+}$ and the $E^{+}$ spectra are broadly
consistent across the lattice sizes, except for the lightest $E^{+}$ on $L=12$
and the first excited $A_{1}^{+}$ on $L=14$, both of which are quite close to
their respective ditorelons and possibly mixing with them. Also the 4th state
in the $L=12$ $A_{1}^{+}$ spectrum shows a shift. Since the typical lattice
size we use in $SU(2)$ is $l\surd\sigma\sim 4$ (except for the relatively
unimportant calculations at the largest values of $a(\beta)$) we can estimate
the ditorelon states to have energies $E_{D}\sim 2\sigma l\sim
2(l\surd\sigma)\surd\sigma\sim 8\surd\sigma$ and we shall therefore only
perform continuum extrapolations of $A_{1}^{+}$ and $E^{+}$ states that are
lighter than $\sim 8\surd\sigma$, although it is still the case that the
heaviest states in these channels may be perturbed by ditorelon contributions.
As for the other channels, since the size $l\surd\sigma\sim 4$ falls between
our $L=14$ and $L=20$ lattices at $\beta=2.427$, we can conclude from the
above that there should be no finite volume corrections that would be visible
outside our statistical errors. We have performed similar finite volume checks
in $SU(3)$ in our earlier paper [11] and we refer the reader to that paper for
details, in particular for a check of the $C=-$ states. In the case of $SU(4)$
we have performed a finite volume analysis comparing the spectra on $18^{3}20$
and $22^{4}$ lattices at $\beta=11.02$, which has helped us identify the
positions of the ditorelon states in the spectra at other values of $\beta$,
and to remove these states from our spectra. With all these checks we have
some confidence that our $SU(2)$, $SU(3)$ and $SU(4)$ continuum spectra will
not be afflicted by significant finite volume corrections.

Our $SU(5)$ and $SU(6)$ calculations are on significantly smaller volumes,
typically with $l\surd\sigma\sim 3.1$, and therefore deserve a separate finite
volume study to the one above. To do so we compare the low-lying glueball
spectra that we obtain in $SU(5)$ on $14^{3}20$ and $18^{3}18$ lattices at
$\beta=17.46$, together with the energies of the ground and first excited
states of the fundamental and $k=2$ winding flux tubes. In physical units
these $L=14$ and $L=18$ spatial lattice sizes correspond to
$l\surd\sigma\simeq 3.06$ and $3.93$ respectively and we recall that the
typical lattice size we use in $SU(5)$ and $SU(6)$ is $l\surd\sigma\simeq 3.1$
i.e. close to that of the smaller of our two lattices. As in $SU(2)$ we
readily identify the extra ditorelon states in the $A_{1}^{++}$ and $E^{++}$
representations, which have energies very close to twice that of the winding
fundamental flux tube. (The ditorelon on the $L=18$ lattice will lie outside
our energy range.) They appear to be well separated from neighbouring glueball
states and so we can readily identify and exclude them from our $SU(5)$ and
$SU(6)$ glueball spectra and exclude them from our continuum limits. The other
energies on our two lattice sizes are mostly within errors of each other, with
a few out by up to two standard deviations, and two states a little more than
that. Given the large number of states being compared, such rare discrepancies
are inevitable and we can claim that we see no significant finite volume
corrections to any of our listed states.

We now turn to our largest $N$ calculations where we use the even smaller
spatial lattice size $l\surd\sigma\simeq 2.6$. Our finite volume study is in
$SU(12)$ on $12^{3}20$ and $14^{3}20$ lattices at $\beta=103.03$, where the
$12^{3}20$ lattice represents our typical physical size, and the energies are
listed in Table 18. In contrast to the previous tables, we list under the
$A_{1}^{++}$ and $E^{++}$ representations our best estimate of the true
glueball spectrum with the ditorelons removed. We will return to the
identification of the latter shortly. We observe that most of the energies are
the same within errors and all are within two standard deviations. That is to
say, any finite volume shifts in the energies are within our statistical
errors at these values of $N$, despite the fact that the $L=12$ lattice volume
is quite small.

The situation with respect to the ditorelons is more complex. The first
complication is that on the $12^{3}20$ lattice the $A_{1}^{++}$ ditorelon is
nearly degenerate with the first excited $A_{1}^{++}$ glueball, so that these
states may well mix even if the overlaps are small due to the large-$N$
suppression. Moreover the same occurs in the $E^{++}$ representation where the
ditorelon is nearly degenerate with the ground state $E^{++}$ glueball. The
second complication arises from the fact that the overlap of the ditorelon
double trace operators onto the single trace operators that we used to
calculate the glueball spectrum is suppressed since $N$ is large. This should
be an advantage and indeed if the overlap is small enough then the ditorelon
will not appear in the spectrum obtained from the single trace operators so
that there is no issue. However it may well appear as a minor component of a
state that at first sight is quite massive but then its $aE_{eff}(t)$ drops
towards the ditorelon energy as $t$ increases and the more massive components
die away. This can lead to ambiguities and indeed does on both our lattice
sizes.

Since all our $SU(12)$ glueball calculations, and also those in $SU(10)$ and
$SU(8)$, are on lattices of roughly the same physical size as our $12^{3}20$
one at $\beta=103.03$, we need to address these problems with the ditorelons.
We do so as follows. First we find that if in constructing our glueball
operators we use blocked links whose extent is smaller than the lattice size
then the overlap of ditorelons onto these operators is small enough that we do
not see any ditorelon state in the $A_{1}^{++}$ spectrum but we do see it
embedded as a small component in more massive states in the $E^{++}$ spectrum.
Since a link at blocking level $bl$ joins sites that are $2^{bl-1}a$ apart,
this means keeping to blocking levels $bl=1-4$ on our $l=12a$ and $l=14a$
lattices. If we now do something that may appear less reasonable and include
$bl=5$ blocked links, which join lattice sites $16a$ apart, so that the
operators formed out of these links wrap multiply around the spatial torus in
all spatial directions, we find that the ditorelon overlaps are much larger
and the ditorelon states now appear quite clearly in the resulting spectrum.
The results of these various calculations are shown in Table 19. As usual the
energies shown are $aE_{eff}(t=2a)$ and so are often slightly higher than the
value of the effective energy plateau. Consider first the $A_{1}^{++}$ spectra
on the $12^{3}20$ lattice. We see that the second and third states in the
$bl=1-5$ spectrum are nearly degenerate, with masses close to what one might
expect for the ditorelon, i.e. $2aE_{f}\sim 0.96$, but only one of them
appears in the $bl=1-4$ spectrum. Thus we infer that one of those two (or a
mixture of the two) is a ditorelon. On the $14^{3}20$ lattice we also see an
extra state in the $bl=1-5$ spectrum, as compared to the $bl=1-4$ spectrum,
but now the state is heavier, as we would expect for a ditorelon since
$2aE_{f}\sim 1.2$ on this lattice. Turning to the $E^{++}$ representation on
the $14^{3}20$ lattice we find an extra state using the $bl=1-5$ basis, as
compared to using the $bl=1-4$ basis, with roughly the expected mass of a
ditorelon. On the $12^{3}20$ lattice things are a bit different: the second
state in the $bl=1-5$ spectrum has no obvious partner in the $bl=1-4$
spectrum, but there is a state in the latter spectrum where $aE_{ff}(t)$
decreases rapidly with increasing $t$ to a comparable value. This is the
second entry in the table and the mass from the ‘plateau’, indicated in square
brackets, is similar to the energy of the second state in the $bl=1-5$
spectrum. We interpret this as follows: the $E^{++}$ ditorelon has a
substantial overlap onto the $bl=1-5$ basis and only a small overlap onto the
$bl=1-4$ basis, but large enough to appear within what appears to be a higher
excited state. This mass is somewhat smaller than that of the state we
identified as a ditorelon on the $14^{3}20$ lattice which is what one expects.
We note that unlike the states we identify as ditorelons, most of the other
states have nearly the same energies on the two lattice sizes, as one expects
for states that are not ditorelons. Recalling that the lattice sizes we use in
our $N\geq 8$ calculations are of the same physical size as our $12^{3}20$
lattice at $\beta=103.03$, we can use the above results in all those cases to
identify ditorelon states and remove them from our glueball spectra. The same
type of technique is useful for $N=5,6$.

A final comment is that the rotational symmetries of the ditorelon states –
and indeed multitorelon states – remain those of our finite volume, i.e.
$\pi/2$ rotations, even in the continuum limit. That is to say, even in that
limit they will fall into the representations of the cubic subgroup of the
continuum rotation group. The genuine glueballs, on the other hand, will fall
into representations of the full continuum rotation group in the continuum
limit once the volume is large enough, up to corrections that are
exponentially small in the spatial size, and this difference can also be
useful in distinguishing ditorelons from genuine glueballs.

###  5.3 lattice masses

As described earlier, we calculate glueball masses from the correlators of
suitable $p=0$ operators. These operators are chosen to have the desired
$R^{PC}$ quantum numbers, where $R\in\\{A_{1},A_{2},E,T_{1},T_{2}\\}$ labels
the irreducible representation of the cubic subgroup of the rotation group and
$P=\pm$ and $C=\pm$ label parity and charge conjugation. We typically start
with a set of 12 different closed loops on the lattice. (For $N=2,3$ we used
27, but we then observed that we can get almost equally good spectra by using
a suitable subset of 12 loops.) We calculate all 24 rotations of the loop and
construct linear combinations of the traces that tranform as $R$. We do so
separately using the real and imaginary parts of the traces, which gives us
operators with $C=\pm$ respectively. We also calculate the parity inverses of
each of these twelve loops, and of their rotations, and by adding and
subtracting appropriate operators from these two sets we form operators for
each $R$ with $P=\pm$. To ensure that we have non-trivial operators for each
set of quantum numbers it is useful to include some loops that have no
symmetry under rotations and parity inversion. With our particular choice of
12 (or 27) closed loops we are able to construct the number of independent
operators shown in Table 20. This is the number of operators at each blocking
level and we typically use 4 or 5 blocking levels in our calculations. This
makes for quite a large basis of operators for all $R^{PC}$ quantum numbers,
and indeed a very large number for most of them. However this large number is
slightly deceptive as at any given $\beta$ only two or perhaps three blocking
levels make an important contribution to the low-lying spectrum. In addition
the states in the $E$ representation are doubly degenerate and in the $T_{1}$
and $T_{2}$ are triply degenerate, so one should divide the corresponding
numbers in Table 20 by 2 or 3 in order to estimate the number of different
energy levels accessible to the basis.

For completeness we list our 12 basic loops in Table 21. The links are
labelled by $1,2,3$ to indicate their spatial directions, with negative signs
for backward going links, and the path ordering is from left to right. So, for
example, the plaquette in the $2-3$ plane would appear as $\\{2,3,-2,-3\\}$.
Our list includes the plaquette, three 6 link loops, four 8 link loops and
four 10 link loops. We also show for each loop which representations it
contributes to once we include all the rotations and the inverses, doing so
for the real ($C=+$) and imaginary ($C=-$) parts of the traces separately.
This basis is used for $N\geq 4$. For $SU(2)$ and $SU(3)$ we use a larger
basis of 27 loops that includes the 12 listed here. For some finite size
studies we have used a reduced set of 8 loops.

We use these bases of operators in our variational calculations of the
glueball spectra. In each $R^{PC}$ sector we calculate a number of the lowest
masses from the correlation functions of the operators that our variational
procedure selects as being the best operators for those states within our
basis. This relies on identifying an effective mass plateau in the correlator,
as described in Section 2.2, and performing an exponential fit to the data
points on that plateau. So an important question is: how reliably can we
identify such a plateau? We will illustrate this with our calculations in
$SU(8)$ on the $20^{3}30$ lattice at $\beta=47.75$. This is at our smallest
value of $a(\beta)$ which, being the closest to the continuum limit, is one of
the most important values in our subsequent extrapolations to the continuum
theory. Also it provides the best resolution in $t$ of our correlation
functions. We begin with the effective masses, $aE_{eff}(t=an_{t})$, of the
three lightest $A_{1}^{++}$ states and of the two lightest $A_{1}^{-+}$
states, as shown in Fig.7. As we will see below, the lightest two $A_{1}^{++}$
and $A_{1}^{-+}$ states become the lightest two $J^{PC}=0^{++}$ and $0^{-+}$
glueballs in the continuum limit. The third $A_{1}^{++}$ state is most likely
part of the nonet of spin states making up the $4^{++}$ glueball ground state
(see below), although if not then it would probably be the third scalar
glueball. It is clear from Fig.7 that we have plausible effective mass
plateaux for all the states, with the weakest case being the excited
$A_{1}^{-+}$ state. The solid lines show our best mass estimates as obtained
from fits to the correlation functions, together with their error bands. The
error band on the lightest $A_{1}^{++}$ state, which is the mass gap of the
gauge theory, is invisible on this plot, so given the importance of this state
we replot it in Fig.8 with an axis rescaling that exposes the errors. Here the
solid line is the best mass estimateand the dahsed lines bound the $\pm 1$
stadard deviation error band. Fig.7 illustrates the obvious fact that as the
masses get larger, the error to signal ratio grows at any fixed value of
$n_{t}$, so that at large enough masses the effective energies becomes too
imprecise to unambiguously indicate the values of $n_{t}$ where the effective
mass plateau begins. In this case we can turn to our calculations of the
excited $A_{1}^{-+}$ state at the smallest lattice spacings in $SU(2)$ and
$SU(3)$, where there is much better evidence for the effective mass plateau
beginning at $aM_{eff}(n_{t}=3a)$, and use that in estimating that the
effective mass plateau in this $SU(8)$ calculations begins at $n_{t}=3a$. This
is the type of argument we use for a number of the heavier states at larger
values of $N$ where our calculations do not extend to very small values of
$a(\beta)$.

As an aside, we also show in this plot the effective mass plot for the
ditorelon. It shows every sign of plateauing to a value not far from $\sim
2aE_{f}\simeq 0.59$ as one would expect. It is clearly important to make sure
that one excludes such a state from the scalar glueball spectrum, as we have
done here, since it is located near the first excited glueball state and would
create a false level ordering if included.

We turn next to states which become $J=2$ glueballs in the continuum limit. As
we shall see below, the five components of a $J=2$ state are spread over the
two components of an $E$ energy level and the three of a $T_{2}$. As always we
average the masses of the ‘degenerate’ doublets of $R=E$ and the triplets of
$R=T_{1}$ and $T_{2}$ to provide single effective masses in each case. What we
show in Fig.9 are the resulting effective mass plots for a number of such
pairs of $E$ and $T_{2}$ states. The equality of the $E^{++}$ and $T_{2}^{++}$
masses for the lightest two pairs is convincing, as it is for the lightest
$E^{-+}$ and $T_{2}^{-+}$. It is also plausible for the lightest
$E^{--},T_{2}^{--}$ pair, and the lightest $E^{+-},T_{2}^{+-}$ pair, and there
is some convergence for the second $E^{-+},T_{2}^{-+}$ pair, even if these
more massive states do not show unambiguous plateaux. Again we show the
ditorelon, this time in the $E^{++}$ representation. In this spectrum it is
nearly degenerate with the $E^{++}$ ground state. The fact that it does not
appear in a nearly degenerate $T_{2}^{++}$ state confirms that it is a finite
volume state, reflecting the limited rotational symmetries of the spatial
volume.

If we examine the effective mass plots for the $T_{1}^{PC}$ states that we
shall later argue approach $J=1$ in the continuum limit, we see that while the
lightest $T_{1}^{+-}$ glueball has a well-defined effective mass plateau, this
is less evident for the more massive states. Finally we remark that if we plot
the effective masses of our heaviest states, such as those of the lightest
five $A_{1}^{--}$ states, the evidence for the effective mass plateaux is not
strong, although one can speculate on the presence of some the plateaux. This
will clearly result in substantial systematic errors on the corresponding mass
estimates, which we cannot readily quantify.

One lesson of the effective mass plots is that the heavier the state the less
reliable will be our error estimate. However these plots also tell us at what
$t=an_{t}$ the effective mass plateau begins for those states where this can
be identified. Since the iterative blocking means that when we vary $\beta$
our variationally selected operators are of roughly constant physical size and
shape, we can assume that the overlap of a given state onto our basis will be
roughly independent of $\beta$ and hence that the effective mass plateau will
begin at a value of $t=a(\beta)n_{t}$ that is roughly constant in physical
units, i.e. at smaller $n_{t}$ as $a(\beta)$ grows. So at larger values of
$a(\beta)$ where $aE_{eff}(t)$ may be too large for us to identify a plateau,
we can nonetheless use the value of $aE_{eff}(t_{0})$, where $t_{0}$ is the
value where our above calculations tell us that a plateau begins, as an
estimate of the mass. This is something we do in our calculations, where
appropriate.

We turn now to the glueball spectra that we obtain by the methods described
above. We obtain the spectra for the lightest glueballs in each $R^{PC}$
sector for each of our gauge groups. As an example we display in 22 our
results for $SU(8)$. Here we have removed the finite volume ditorelon states
in the $A_{1}^{++}$ and $E^{++}$ sectors whenever they are present, as
discussed in Section 5.2, so that what we present in the Table is our best
estimate of the infinite volume glueball spectrum. A similar Table for $SU(3)$
has been published separately in [11]. All this assumes of course that none of
these states is a multiglueball state. Since we use single trace operators,
their overlap onto multiglueball states should decrease with increasing $N$
and in any case they would be quite heavy. Some explicit calculations that
provide evidence that such mutliglueball states do not appear in our spectra
are described in Section 5.8.

The reader will note that in some cases what we list as a higher excited state
has a lower mass than what is listed as a lower excited state. This typically
involves states that are nearly degenerate. Our procedure is to order the
states according to the values of their effective masses at $t=a$ which is the
value of $t$ at which our variational calculation typically operates. This
nearly always corresponds to final mass estimates (from the plateaux) that are
in the same order, but very occasionally this not the case. In the very few
cases where the difference is well outside the statistical errors we invert
the ordering, but otherwise we do not. In addition it may be that the
continuum extrapolation leads to a level inversion; again we only take that on
board if the difference is well outside the errors. We have a good number of
states that are nearly degenerate with a neighbouring state and in these cases
the variational procedure itself is likely to mix them, and if that is done
differently at different values of $\beta$ it can lead to poor extrapolations
to the continuum limit.

Finally we remark that the ensembles of lattice fields used at each $N$ and
$\beta$ have a reasonable distribution of topological charge, either because
the tunnelling between topological charge sectors is sufficiently frequent, or
because we have imposed such a distribution on the initial lattice fields of
the collection of sequences that make up the total ensemble, as discussed in
Section 2.4.

###  5.4 strong-to-weak coupling transition

Before moving on to the continuum extrapolation of the lattice masses that we
have calculated above, we briefly remark on the ‘bulk’ transition that
interpolates in $\beta\propto 1/g^{2}$ between the strong and weak coupling
regimes. From eqns(1,2) we see that the naive expectation is that at strong
coupling, $\beta\to 0$, the natural expansion of a lattice quantity is in
positive powers of $\beta$ while at weak coupling, $\beta\to\infty$, the
natural expansion is in positive powers of $1/\beta$. On the weak coupling
side asymptotic freedom promotes the dependence of physical quantities to
exponentials, $\propto\exp\\{-c\beta\\}$, and it is important to make sure
that our calculations are indeed in the weak coupling regime.

As remarked in Section 2.1, for $SU(N\leq 4)$ the transition is known to be a
smooth crossover, while for $SU(5)$ it is weakly first order, and for
$SU(N\geq 6)$ it is strongly first order. For $SU(N\leq 4)$ the location of
the crossover coincides with a dip in the mass gap which then drives a peak in
the specific heat. (The specific heat is proportional to the sum of the
plaquette-plaquette correlator and so peaks where the correlation length has a
peak, i.e. where the mass gap has a dip.) For $SU(2)$ this peak is around
$\beta\sim 2.15$ as we see, for example, in Fig.8 of [50], while for $SU(3)$
the peak is around $\beta\sim 5.4$ as we see, for example, in Fig.4a of [51].
For $SU(4)$ the dip in the mass gap can be seen in Fig.1 of [1] and is located
around $\beta\sim 10.45$. For larger $N$ the transition has been shown in [14]
to be first order, and strongly first order for $N\geq 6$. In that case we
have a strong hysteresis, i.e. if we lower $\beta$ slowly from large values
then the transition occurs at $\beta=\beta_{b}^{\downarrow}$, while if we
increase $\beta$ slowly from small values then the transition occurs at
$\beta=\beta_{b}^{\uparrow}$, with $\beta_{b}^{\uparrow}$ significantly larger
than $\beta_{b}^{\downarrow}$. Values for these transitions can be found
listed in Table 16 of [14]. As described in that paper, calculations performed
in the weak coupling false vacuum just above $\beta=\beta_{b}^{\downarrow}$,
and well below $\beta=\beta_{b}^{\uparrow}$, show no sign of being affected by
the simultaneous presence of the deeper true vacuum elsewhere in field space.
This fact has been exploited in previous large $N$ calculations. For example
all but one of the $\beta$ values in the $SU(8)$ calculations in [2] lie
within this hysteresis window, as does the $SU(8)$ calculation in [3]. As is
clear from Tables 1-8 all our calculations have been performed on the weak
coupling side of the cross-over or transition. In addition we note that none
of our $SU(5)$ or $SU(6)$ calculations fall within the hysteresis window, only
one of our $SU(8)$ calculations lies within this window, while some of the
$SU(10)$ and about half of the $SU(12)$ calculations lie within the window.

As remarked above, the main effect of the crossover at smaller $N$ is an
anomalous dip in the value of the mass gap, i.e. that of the lightest scalar
glueball $am_{0^{++}}$. This dip is caused by a nearby critical point which
lies at the end of a first-order transition line in an extended fundamental
and adjoint coupling plane, at which critical point the scalar glueball mass
vanishes. For a detailed calculation in $SU(3)$ see for example [52]. (The
fundamental coupling is the one that which appears in our path integral.) As
$N$ increases the first-order line intersects the fundamental axis and
provides the first-order bulk transition separating weak and strong coupling.
However until the critical point has moved far away from the fundamental axis
there may still be a dip in the mass gap near the first-order transition.
Indeed one can see some plausible evidence for this occuring in the behaviour
of the $SU(6)$ and $SU(8)$ mass gaps listed in [2]. Since there is no
theoretical reason to expect that this dip can be encoded in the weak coupling
expansion of the lattice action in powers of $a^{2}$, we minimise the risk to
the continuum extrapolation of the scalar glueball mass by simply excluding
from the continuum fit the value obtained at the largest value of $\beta$ when
including this value would require adding an additional $a^{4}$ correction to
the fit. As it turns out, we need to do this for all our $SU(N)$ groups except
for $SU(2)$.

###  5.5 continuum masses

For each of our $SU(N)$ lattice gauge theories we now have the low-lying
glueball spectra for a range of values of $a(\beta)$. These are all given in
lattice units as $aM$, and to transform that to physical units we can take the
ratio to the string tension, $a\surd\sigma$, that we have simultaneously
calculated. We can then extrapolate this ratio to the continuum limit using
the standard Symanzik effective action analysis [53] that tells us that for
our lattice action the leading correction at tree-level is $O(a^{2})$:

$\frac{aM(a)}{a\surd\sigma(a)}=\frac{M(a)}{\surd\sigma(a)}=\frac{M(0)}{\surd\sigma(0)}+a^{2}\sigma(a)+O(a^{4}).$
(27)

Here we have used the calculated string tension, $a^{2}\sigma(a)$, as the
$O(a^{2})$ correction. Clearly we could use any other calculated energy, and
this would differ at $O(a^{4})$ in eqn(27). It is convenient to use
$a^{2}\sigma(a)$ since its calculated value has very small errors.

The results of these continuum extrapolations are listed in Tables 23-30 for
the gauge groups ranging from $SU2)$ to $SU(12)$. (As always, each $E$
doublet, and each $T_{2}$ or $T_{1}$ triplet appears as a single state in our
tables and discussions.) A few entries are accompanied by a star denoting a
poor fit, $2.5<\chi^{2}/n_{df}<3.5$, or a double star denoting a very poor
fit, $\chi^{2}/n_{df}\geq 3.5$. For $SU(10)$ and $SU(12)$ we have only 5
values of $\beta$ and since the masses at the coarsest value of $\beta$ often
have to be discarded from the fit (not surprisingly) and since we are fitting
2 parameters, we often have only 2 degrees of freedom. For lower $N$ we
usually have 3 and sometimes 4 degrees of freedom, which gives much more
confidence in the extrapolations. All of which is to say that while the fits
for larger $N$ are certainly not trivial, it would be good to be able to do
better in future calculations.

We now illustrate the quality of these linear continuum extrapolations for our
most interesting glueball states. We do so for $SU(4)$. In Fig.10 we show our
extrapolations of the lightest two $A_{1}^{++}$, $E^{++}$ and $T_{2}^{++}$
states. These states are of particular importance because, as we shall see
latter on, they correspond to the lightest two $J^{PC}=0^{++}$ and $2^{++}$
states. We see that all the linear fits are convincing, even if in some cases
we have to exclude from the fit the value at the largest $a(\beta)$. In Fig.11
we show the corresponding plot for $P=-$ which, as we shall see latter on,
correspond to the lightest two $J^{PC}=0^{-+}$ and $2^{-+}$ states. The
lightest states have very plausible continuum extrapolations, although the
excited states, which are heavier than those for $P=+$, begin to show a large
scatter indicating a poor fit. In Fig.12 we show the extrapolations of various
$T_{1}^{PC}$ states that we shall later argue correspond to $J=1$, and again
we see fits that appear convincing for the lighter states and quite plausible
for the heavier states. Even for the heaviest states to which we can plausibly
assign a continuum spin, such as the states that pair up to give the
$J^{PC}=2^{--}$ and $2^{+-}$ ground states and the states that provide the
seven components of the $3^{+-}$ ground state, where the larger errors lead to
a greater scatter of points, continuum fits linear in $a^{2}\sigma$ are still
plausible.

###  5.6 $N\to\infty$ extrapolation

Amongst the various $SU(N)$ calculations, the one that is most interesting
from a phenomenological point of view is the $SU(3)$ one, and that is why we
devoted a separate paper to that case [11]. From a theoretical point of view
however the most interesting glueball spectra are those of the
$SU(N\to\infty)$ theory since the theoretical simplifications in that limit
make it the most likely case to be accessible to analytic solution, whether
complete or partial.

To obtain the $N=\infty$ spectrum from our results so far, we use the fact
that in the pure gauge theory the leading correction is $O(1/N^{2})$. So we
can extrapolate the continuum mass ratios in Tables 23-30 using

$\left.\frac{M_{i}}{\surd\sigma}\right|_{N}=\left.\frac{M_{i}}{\surd\sigma}\right|_{\infty}+\frac{c_{i}}{N^{2}}+O\left(\frac{1}{N^{4}}\right).$
(28)

The results of these extrapolations are presented in Table 31. In this table
the stars point to poor fits exactly as described earlier for the continuum
fits. The very poor fit for the first excited $A_{1}^{++}$ is due to a large
mismatch between the $SU(5)$ and $SU(6)$ mass estimates which may well be due
to an inadequate treatment of the ditorelon influence. (We can obtain a good
fit to $N\geq 6$ and this gives a mass $\sim 5.85(9)$ which is the same within
errors.) Most of the fits are to $N\geq 2$ or $N\geq 3$ but some are over a
more restricted range of $N$ and this is indicated by a dagger.These states
are the $A_{2}^{++}$ ground state, fitted to $N\geq 4$, the $T_{2}^{-+}$
second excited state, also fitted to $N\geq 4$, the $T_{2}^{--}$ ground state,
again fitted to $N\geq 4$, and finally the $A_{2}^{+-}$ ground state which is
more arguable since it was fitted to $N\leq 8$. From the practical point of
view the most important extrapolations are for those states to which we are
able to assign a continuum spin (in the next section). We therefore show these
extrapolations in Figs.13,14,15 for states with $J=0,2,1$ respectively.

###  5.7 continuum spins

So far we have used the representations of the rotational symmetry of our
cubic spatial lattice to label our glueball states. However as we approach the
continuum limit these states will approach the continuum glueball states and
these belong to representations of the continuum rotational symmetry, i.e.
they fall into degenerate multiplets of $2J+1$ states where $J$ is the spin.
In determining the continuum limit of the low-lying glueball spectrum, it is
clearly more useful to be able to assign the states to a given spin $J$,
rather than to the representations of the cubic subgroup which have a much
less fine ‘resolution’ since all of $J=0,1,2,...,\infty$ are mapped to just 5
cubic representations. How the $2J+1$ states for a given $J$ are distributed
amongst the representations of the cubic symmetry subgroup is given, for the
relevant low values of $J$, in Table 32. So, for example, we see that the
seven states corresponding to a $J=3$ glueball will be distributed over a
singlet $A_{2}$, a degenerate triplet $T_{1}$ and a degenerate triplet
$T_{2}$, so seven states in total. These $A_{2}$, $T_{1}$ and $T_{2}$ states
will be split by $O(a^{2})$ lattice spacing corrections, generated by
irrelevant operators. So once $a$ is small enough these states will be nearly
degenerate and one can use this near-degeneracy to identify the continuum
spin.

This strategy works for the lightest states but becomes rapidly unrealistic
for heavier states. The latter will have larger statistical errors and will
fall amongst other states that become more densely packed as the energy
increases, so that identifying apparent near-degeneracies between states in
different representations becomes highly ambiguous. There exist more sytematic
ways of assigning continuum spin to lattice states, such as [8, 9, 54], but
these are beyond the scope of this work. Accordingly we shall limit ourselves
to the lightest states where any ambiguity in identifying near-degeneracies is
either small or non-existent.

The states whose spin $J$ we feel confident in identifying are listed in Table
33 and the resulting continuum masses are listed in Table 34 and Table 35. We
will now briefly illustrate the argument for the assignments, taking $SU(8)$
as a typical example. Consider the masses in the $SU(8)$ column of Table 35.
We have obtained these by interpreting the masses listed in Table 28 as
follows. The ground state of $R^{PC}=A_{1}^{++}$ is much lighter than any
other mass and so must be the singlet $J^{PC}=0^{++}$ ground state. The first
excited $A_{1}^{++}$ state certainly has no ‘nearly-degenerate’ partners
amongst all the other $P,C=+,+$ states and so it must be the first excited
$0^{++}$ glueball. The ground states of the doubly degenerate $E^{++}$ and the
triply degenerate $T_{2}^{++}$ states are nearly degenerate and there are no
other states with similar masses, so together they must provide the five
components of the $J^{PC}=2^{++}$ ground state. Similarly the first excited
$E^{++}$ doublet and $T_{2}^{++}$ triplet have no nearly degenerate companion
states elsewhere, and so they provide the first excited $J^{PC}=2^{++}$ state.
The argument for the ground and first excited $0^{-+}$ states is equally
straightforward, as it is for the $E^{-+}$ and the $T_{2}^{-+}$ ground states
forming the $2^{-+}$ ground state. The first excited $E^{-+}$ and the
$T_{2}^{-+}$ states are consistent with forming the first excited $2^{-+}$
state, but here there is some ambiguity: it is also possible that these states
pair with the second excited $A_{1}^{-+}$ state and the $T_{1}^{-+}$ ground
state to make up the $4^{-+}$ ground state. This ambiguity arises because the
errors are sufficiently large that both possibilities can be entertained. Here
we can resolve the ambiguity by observing that for neighbouring values of $N$,
$SU(6)$ and $SU(10)$, the ambiguity is not present and the first excited
$E^{-+}$ and the $T_{2}^{-+}$ states form the first excited $2^{-+}$ state,
strongly suggesting that this is also the case for $SU(8)$. This choice is not
as solid as the earlier choices and fortunately it is not an argument we need
to make very often in coming to our choices of $J$. With this choice we can
plausibly assign the lightest $T_{1}^{-+}$ state to the $1^{-+}$ ground state.
(The second excited $A_{1}^{-+}$ state has a similar mass but there are no
further partners allowing us to assign it elsewhere according to Table 32.)
Similar arguments to the above lead to the $J^{+-}$ and $J^{--}$ choices in
Table 33 and the corresponding masses in Table 35. Determining the spin $J$ in
the above way, clearly requires us not to miss any intermediate states in the
mass range of interest. For this one needs a good enough overlap onto all the
low-lying states and for that one needs a large basis of operators, which has
been our goal in these calculations.

Using the arguments sketched above we infer the $J^{PC}$ glueball masses
listed in Tables 34,35, from the various masses listed in Tables 23-30. In
addition by inspecting the remaining states in Tables 23-30 we can infer some
lower bounds on the ground states in some other channels, and these are listed
in Table 36. These estimates are fairly rough and at that level apply to all
our $SU(N)$ theories.

The most interesting glueball spectra are those for $SU(3)$ and for
$SU(\infty)$: the former for its potential phenomenological relevance and the
latter for its potential theoretical accessibility. The $SU(3)$ theory is
sufficiently close to the real world of $QCD$ that it makes sense to attempt
to present the masses in physical units, using scales that one believes to be
relatively insensitive to the presence of quarks. We refer the reader to [11]
for a more detailed discussion. Here we simply recall from our earlier
discussion that the Sommer scale $r_{0}$, which is defined by the heavy quark
potential at intermediate distances, is believed to vary weakly with quark
masses, with a physical value of $r_{0}\simeq 0.472(5)\mathrm{fm}$ [47]. One
can extract from the published data the value $r_{0}\surd\sigma=1.160(6)$ [11]
This allows us to translate masses in units of $\surd\sigma$ to units of
$r_{0}$ and finally to units of $\mathrm{GeV}$. Doing so we obtain the masses
listed in Table 37. This Table is similar to that in [11] except that we have
included lower bounds for most of the masses where we previously lacked
entries.

We now turn to our main focus: the $N\to\infty$ limit. We extrapolate the
masses in Tables 34,35 with $O(1/N^{2})$ corrections, as in eqn(28). The
extrapolations for most of these states are displayed in Figs.13,14,15. This
leads to the masses (in units of the string tension) listed in Table 38. In
making these extrapolations we have sometimes excluded the mass for the
smallest value of $N$ ($N=2$ for $C=+$ and $N=3$ for $C=-$), in order to
achieve a significantly better fit. This we have done for the $2^{++}$,
$1^{-+}$ and $1^{--}$ ground states. This presumably reflects the need for at
least a further $1/N^{4}$ correction, but our number of data points is too
small for this. In some cases we obtain better fits by excluding the $N=12$
values, and although this could be argued for on the basis that our range of
$a$ for the continuum extrapolation is relatively limited in that case, the
danger of ‘cherry-picking’ has led us to avoid doing so.

The main features of these continuum spectra are that the lightest glueball is
the $J^{PC}=0^{++}$ scalar, with the $J^{PC}=2^{++}$ tensor about $50\%$
heavier, and with the $J^{PC}=0^{-+}$ pseudoscalar just slightly above that.
The next state is the $J^{PC}=1^{+-}$ vector, which is special in that it is
the only relatively light $C=-$ glueball, with most $C=-$ states being very
heavy. Just a little heavier than this vector is the first excited $0^{++}$
and then the $2^{-+}$ ground state. The $N$ dependence of nearly all these
states is weak and readily absorbed into a $O(1/N^{2})$ correction down to at
least $SU(3)$; all this confirming that the $SU(3)$ gauge theory is indeed
‘close to’ $SU(\infty)$.

###  5.8 scattering states

A given operator will project onto all states with the quantum numbers of that
operator. In particular our single trace glueball operators will project onto
multi-glueball states in addition to the glueballs that we are interested in.
This means that some of the states in the ‘glueball’ spectra that we have
calculated may in fact be multi-glueball states. In the ideal case of a very
large spatial volume these will be scattering states where at sufficiently
large times the states are far apart

There is some ambiguity here of course. Sufficiently heavy glueballs will be
unstable, and will decay into multiglueball states but to the extent that
their decay width is not too large, we still regard them as glueball states.
Indeed as $N$ increases the decay widths and the overlap of multiglueball
states onto our single trace basis of glueball operators should decrease. So
although we expect this to be an issue primarily at smaller $N$, it is clearly
important to address. We do so here in some detail because it has been
addressed only briefly and occasionally in past glueball calculations.

A full analysis of scattering states and glueball decays would involve
different techniques to those used in this paper and so are outside the scope
of this work. Here we will, instead, perform some exploratory calculations to
obtain an indication of the impact of multiglueball states on our calculated
glueball spectra. We limit ourselves to considering states consisting of two
glueballs since these are the lightest ones and our calculated glueball
spectra do not extend to masses that are much higher. If there is weak mixing
between double trace and single trace operators then we can introduce two
glueball scattering states by using double trace operators, such as

$\phi_{ab}(t)=(\phi_{a}(t)-\langle\phi_{a}\rangle)(\phi_{b}(t)-\langle\phi_{b}\rangle)-\langle(\phi_{a}-\langle\phi_{a}\rangle)(\phi_{b}-\langle\phi_{b}\rangle)\rangle.$
(29)

Here $\phi_{a}$ and $\phi_{b}$ are single trace operators. The vacuum
subtractions ensure that $\phi_{ab}$ does not have a projection onto the
vacuum or a ‘trivial’ projection onto single glueball states through terms
such as $\langle\phi_{a}\rangle\phi_{b}(t)$. One can calculate the mass
spectrum obtained when one adds such double trace operators to the basis of
single trace operators and compare it to that obtained using just the single
trace operators. If the resulting spectra are the same then this strongly
suggests that the spectrum obtained using single trace operators already
includes two glueball states. If, on the other hand, the spectrum obtained
using the expanded basis produces extra states that can be plausibly
interpreted as two glueball scattering states, then this suggests that the
states obtained with our original basis of single trace operators does not
contain such states and that categorising them as single glueball states is
probably correct.

The potential number of operators $\phi_{ab}$ is clearly too huge to be
practical here and so we severely limit the number for our study as follows.
Firstly, we take both $\phi_{a}$ and $\phi_{b}$ to have zero momentum, so both
the total and relative momenta are zero. Of course this does not mean that the
relative momentum of the two glueballs has to be zero, since the relative
momentum is not a conserved quantum number in an interacting system, but one
naively expects that the main overlap will be onto zero relative momentum.
Secondly we only keep the 2 or 3 blocking levels that are most important at
the $\beta$ we calculate. (We systematically avoid using the largest blocking
levels so as to exclude unwanted ditorelon states.) Thirdly, we take the same
blocking levels for $\phi_{a}$ and $\phi_{b}$. For $\phi_{a}$ we use only the
rotationally invariant sum of (blocked) plaquettes. That is to say it is in
the $A_{1}^{++}$ representation. So $\phi_{ab}$ will be in the same
representation as $\phi_{b}$. For $\phi_{b}$ we take 3 different loops chosen
so that we can have projections onto all representations including $P=\pm$ and
$C=\pm$. The lightest energy of the corresponding asymptotic two glueball
state would be the sum of the lightest $A_{1}^{++}$ mass plus the lightest
mass in the representation of $\phi_{b}$ if we were in an infinite spatial
volume. Since our spatial volume is finite the glueballs are interacting at
all times and there will be a shift in the total energy. However, to the
extent that our spatial volume is not very small, this shift should be small
and we will use the naive sum as the quantity against which to compare our
supposed scattering states.

We perform calculations in $SU(3)$, which is the most physically interesting
case amongst our lower $N$ calculations, and in $SU(8)$ which is
representative of our large $N$ calculations. In $SU(3)$ we work at
$\beta=6.235$ on the same $26^{4}$ lattice used in our above glueball
calculations. For completeness we have also carried out some calculations on
$18^{3}26$ and $34^{3}26$ lattices at the same value of $\beta$. In $SU(8)$ we
work at $\beta=46.70$ on a $16^{3}24$ lattice which, again, is the same as
that used in our glueball calculations. The spatial volumes at larger $N$ are
smaller than those at smaller $N$, taking advantage of the expected
suppression of finite volume corrections with increasing $N$, and this is the
reason for this additional calculation in $SU(8)$.

We calculate correlators of the double trace operators with each other as well
as with the single trace operators, since all of these are needed in the
variational calculation using the basis that combines single and double trace
operators. Because we need high statistics for glueball calculations, we
calculate the correlators at the same time as we generate the lattice fields.
At that stage we do not know the values of the vacuum subtractions in eqn(29)
and so we need to calculate a number of other correlators so that at the later
analysis stage we can reconstruct the desired subtracted correlators. The
relevant expressions are given in Appendix B.

We begin with our $SU(3)$ calculation. In Fig. 16 we show the effective masses
of the lightest few states in the $A_{1}^{++}$ representations obtained using
the basis of single trace operators (open circles) and the same basis extended
with our double trace operators (filled points). We see quite clearly that the
two sets of states match each other well except for one state, whose effective
masses are represented by a $\blacklozenge$, and which clearly has no partner
amongst the states from the single trace basis. This state appears to
asymptote, at larger $t$, to an energy that is close to twice the lightest
$A_{1}^{++}$ mass, i.e. that of the lightest free two glueball state, and so
we infer that this is indeed the lightest two glueball state. Apart from this
state, the lightest 4 states in the two bases match each other quite
precisely. That is to say, we have good evidence that the $A_{1}^{++}$
continuum masses listed in Table 24 are indeed those of single glueball
states.

It is interesting to ask what kind of $A_{1}^{++}$ energy spectrum one obtains
if one uses a basis consisting only of the double trace operators. The result
is shown in Fig. 17. We see that the state that is lightest at small $t$ is
consistent with being the lightest 2 glueball state, and that our variational
wave-functional has a good overlap onto that state. What is equally
interesting is the state that is the first excited state at small $t$ and
whose effective mass drops monotonically below that of our two glueball state
as $t$ increases, consistent with heading towards the mass of the lightest
glueball. This presumably represents the overlap of our double trace basis
onto the lightest glueball. A rough estimate suggests that this variationally
selected operator has an overlap (squared) of $\sim 10\%$ which again points
to little mixing of scattering states into our low-lying glueball spectrum.

Clearly we would like to repeat our above analysis of the $A_{1}^{++}$
spectrum for all the other representations. However once the lightest state in
some $R^{PC}$ representation is much heavier than in the above case the energy
of the predicted scattering state lies in a dense part of the spectrum, all
the energies have much larger statistical errors and we cannot follow the
effective energies beyond the smallest values of $t$. That is to say, the
analysis becomes hopelessly ambiguous. We will therefore restrict ourselves to
a few cases where the energies are still manageably small. That is to say, we
will look at the $T_{2}^{++}$ which contains 3 components of the lightest
$J=2^{++}$ state, the $A_{1}^{-+}$ where the lightest state is the interesting