The glueball spectrum of SU(3) gauge theory in 3+1 dimension

We calculate the low-lying glueball spectrum of the SU(3) lattice gauge theory in 3+1 dimensions for the range of beta up to beta=6.50 using the standard plaquette action. We do so for states in all the representations R of the cubic rotation group, and for both values of parity P and charge conjugation C. We extrapolate these results to the continuum limit of the theory using the confining string tension as our energy scale. We also present our results in units of the r0 scale and, from that, in terms of physical `GeV' units. For a number of these states we are able to identify their continuum spins J with very little ambiguity. We also calculate the topological charge Q of the lattice gauge fields so as to show that we have sufficient ergodicity throughout our range of beta, and we calculate the multiplicative renormalisation of Q as a function of beta. We also obtain the continuum limit of the SU(3) topological susceptibility.


Introduction
The spectrum of the SU(3) gauge theory in 3+1 dimensions is of interest for both phenomenological and theoretical reasons. On the phenomenological side it may inform us about the approximate location of glueballs in QCD if the mixing with qq states is as modest as suggested by the OZI rule. Here it is the lighter glueballs, with masses below say 3 GeV, that are of greatest interest given the experimental constraints in observing and identifying glueballs. The masses of such glueballs, some 4 or 5 in total, have been reasonably well determined for a long time: see for example [1]. Subsequent calculations, for example [3,4,5,6], have remained consistent with these earlier estimates, albeit with much reduced errors, indicating that this part of the spectrum is well understood. Slight shifts in the mass estimates are mostly due to differences in importing the physical 'GeV' scale into the unphysical pure gauge theory. Since the two theories are not the same this reflects an intrinsic ambiguity.
On the theoretical side the spectrum represents the obvious challenge for attempts to obtain an analytic control of the gauge theory. (This is also of indirect phenomenological importance since much of the difficulty in obtaining analytic control of the long distance physics of QCD arises from the non-perturbative physics of the underlying gauge theory.) To use the glueball spectrum as a test of theoretical approaches it is important to have a substantial number of states for a variety of spins J, parities P and charge conjugations C, with excited as well as ground states. Many of these states will be quite heavy and calculating their masses on the lattice is more challenging. Moreover here we need accurate mass estimates, since the errors cannot be hidden within the uncertainty associated with importing a physical 'GeV' scale. Such calculations are still open to considerable improvement and the purpose of the present paper is to take some further steps in that direction.
The first step is to try to improve the reliability of our calculation of lattice masses. Since the usual way to calculate masses on the lattice is from the time decay of correlators of lattice operators, we use a large basis of operators so as to decrease the probability that we miss intermediate states in the spectrum. To improve our extrapolation of these masses to the continuum limit we perform calculations over a wide range of lattice spacings, a, down to a √ σ ∼ 0.10 where σ is the confining string tension (which we also calculate). We also perform a finite volume comparison at one of our smaller lattice spacings so as to identify any finite volume corrections to our calculations. The calculations at the smallest values of a are of course the most important, and here one needs to be concerned about the ergodicity of the Monte Carlo and in particular about the possible freezing of the fields into sectors of unchanging topological charge Q. To investigate this we calculate the decorrelation of Q for all our values of a. In the process we need to demonstrate that our method of calculating Q is reliable and this also allows us to provide a precise estimate of the continuum topological susceptibility χ t . The glueball states that we calculate will belong to the representations R of the subgroup of rotations that leave the cubic lattice invariant. That is to say, they are labelled by R P C rather than the continuum J P C . Since one knows how the 2J + 1 dimensional vector space associated with a glueball of spin J is partitioned amongst the various cubic representations R, and since these 2J + 1 states should be nearly degenerate at small a, one can attempt to identify the value of J through the observation of the appropriate near-degeneracies. This is limited by the accuracy of our calculation of enough excited states as well as the ground states but we are still able to identify J for a substantial number of states. It is however worth emphasising that if one wants to compare the predicted spectrum of a theoretical approach or model to our lattice data, it may make sense to convert the predictions from those for J P C to those for R P C . This is simple to do using the conversion that we will present later in this paper.
Our paper is as follows. In the next section we give a brief overview of our lattice setup and how we calculate energies. We then move on to present our calculations of the confining string tension σ which we shall use as our scale for the calculated glueball masses in our extrapolation of the latter to the continuum limit. We then turn to the calculation of the glueball masses m G on our cubic lattices. Here we discuss the representations of the cubic group and present the relation between these and the continuum rotation group. We also provide a detailed discussion and check of finite volume effects. In the following section we present our extrapolations to the continuum limit of the dimensionless ratios m G / √ σ and the conversion of the masses into 'GeV' units. Where possible we also identify the continuum spin J of the states. We follow this section with a comparison of our results to those of some of the most comprehensive earlier spectrum calculations. We then return to the issue of the freezing of the topological charge at the smallest lattice spacings. We calculate Q on large ensembles of lattice fields, to check for ergodicity, and as a byproduct obtain a value for the topological susceptibility in the continuum limit. We then turn to pinpointing the ways in which our calculations need to be improved. We summarise our findings in the concluding section.

Calculating on a lattice
Our lattice setup is entirely standard. We work on hypercubic lattices of size L 3 s L t with lattice spacing a and with periodic boundary conditions on the fields. Our fields are SU(3) matrices, U l , assigned to the links l of the lattice, The Euclidean path integral is where DU is the Haar measure and we use the standard plaquette action, Here U p is the ordered product of link matrices around the plaquette p. We write β = 6/g 2 , where g 2 is the bare coupling and this provides a definition of the running coupling on the length scale a. Since the theory is asymptotically free g 2 → 0 and hence β → ∞ as a → 0. Monte Carlo calculations are performed using a mixture of standard Cabibbo-Marinari heat bath and over-relaxation sweeps through the lattice in the ratio 1:4. We typically perform 2 × 10 6 sweeps at each value of β at each lattice size, and we typically calculate correlators every 25 sweeps. We calculate masses M and energies E from the exponential decay of appropriate correlation functions. For example, the energy E gs of the lightest state with the quantum numbers of the operator φ is given by These masses and energies are obtained in lattice units, i.e. as aE gs , and are relative to the vacuum once L t is large enough, so that the vacuum energy is zero. If φ has the quantum numbers of the vacuum then the lightest state of interest is the first excited state above the vacuum. In this case it can be convenient to use the vacuum-subtracted operator φ − φ , which will remove the contribution of the vacuum in eqn (3), so that the lightest non-trivial state appears as the leading term in the expansion of states.
The accuracy of such a calculation of aE gs is constrained by the fact that the statistical errors are roughly independent of t (for pure gauge theories) while the desired 'signal' is decreasing exponentially with t. So if the overlap of the desired state onto our basis is small, the relevant correlator will disappear into the statistical 'noise' before we get to large enough t for the correlator to be dominated by the simple exponential from which we can extract the ground state energy. Similarly the correlator of a more massive state will disappear into the 'noise' at smaller t = an t and this may make ambiguous the judgement of whether it is dominated by a single exponential. All this may provide an important source of systematic error.
The glueball masses calculated using eqn (3) are in lattice units, i.e. aM, and to remove this factor of a we can take a ratio of masses or energies. For this purpose we simultaneously calculate the confining string tension, a 2 σ, and form the dimensionless ratio aM/a √ σ = M/ √ σ from which the lattice scale has disappeared. The leading lattice corrections to such a ratio are O(a 2 ) so we extrapolate to the continuum limit using We will usually only need the leading correction. In reality the coefficients c i are power series in g 2 but this dependence on a is too weak (logarithmic) to be visible in our data, so we follow common practice in treating the c i as constants.
The values of β we use are listed in Table 1. The lattice sizes next to these β values are the ones we use for our string tension and glueball mass calculations. The lattices listed on the right are used for calculating topology and also, in some cases, to check for finite volume corrections. We also show the average plaquette at each value of β as well as the string tension (see Section 3) and the mass gap (see Section 4).

String tension
We calculate the string tension by calculating the ground state energy E(l) of a flux tube of length l that closes on itself by winding once around a spatial torus of size l. We use eqn (3) where the operator φ is the product of SU(3) link matrices taken around a non-contractible closed path that winds once around the spatial torus. The simplest such operator is the spatial Polyakov loop φ(n t ) = l p (n t ) = ny,nz Tr Lx nx=1 U x (n x , n y , n z , n t ) (5) where we take the product of the link matrices in the x-direction around the x-torus of length l = aL x . Here (x, y, z, t) = (an x , an y , an z , an t ), and the sum over n y and n z produces an operator and hence a state with zero transverse momentum, p y = p z = 0. Since φ is clearly invariant under translations in x, the state also possesses p x = 0. We do the same for loops over the y and z spatial tori and average over all three. We do the same replacing the link matrices by blocked link matrices, which is the essential step in obtaining a useful overlap onto the ground state. This provides a vector space of operators on which we perform a variational calculation to extract the operator that is our best approximation to the ground state wavefunctional. From the correlator of this operator we extract a value for E(l).
We extract the string tension σ from E(l) using the 'Nambu-Goto' formula which arises from the light-cone quantisation of the bosonic string and is known to provide an excellent approximation to the lattice calculations [7] for reasons that have now become well understood. (See for example [8,9].) To extract the ground state energy of the flux tube using eqn(3) we need to be confident that we have gone to large enough t = n t that the correlator is dominated by a single exponential within errors, and moreover that these errors are small enough for this observation to be significant. A useful way to test for a single exponential is to extract an effective energy If a single exponential dominates C(n t ) for n t ≥ñ t , then E ef f (n t ) should become independent of n t for n t ≥ñ t , and in that case E ef f (ñ t ) provides an estimate of the energy E(l) . (In practice we extract our E from a fit over a range of n t .) That is to say, we search for a 'plateau' in the values of E ef f (n t ) against n t . An important point is that as long as our operators are entirely spacelike, the reflection positivity of our action ensures that E ef f (n t ) ≥ E(l) for any n t , up to statistical errors, and indeed that the value of E ef f (n t ) must be monotonically decreasing with increasing n t . In Fig. 1 we plot the values of aE ef f (t) for each of our β values together with the ±1σ error band for our best estimate of the energy E(l). At our smallest value of β, a and hence aE are large and the exponential disappears into the statistical errors before we see an unambiguous effective energy plateau -hence the very wide error band around our estimate of E(l). At other values of β there is little ambiguity and the error bands are correspondingly narrow. We also note that the asymptotic exponential typically begins at t = a at smaller β and at t = 2a at larger β. Since we expect that our space of multiply blocked operators will have an overlap onto the ground state that is weakly dependent on a(β), this is as expected. It is also the reason that for our lowest value of β, where no effective energy plateau apears, we are confident in choosing a value of E that encompasses, within errors, the values of aE ef f (t) for t = a and t = 2a.

Glueball spectrum
We calculate glueball masses from Euclidean correlators as in eqn (3). Since only states |n such that vac|φ † |n = 0 will contribute to the correlator, we can restrict the contributing states to have particular quantum numbers by choosing operators φ with those quantum numbers. Since the error to signal ratio of such a correlator increases exponentially in t = an t , we can only hope to identify an effective mass plateau using eqn (7) if the overlap of the state onto the operator φ is large enough. So in order to avoid missing intermediate states in the spectrum one needs a sufficiently large basis of operators.
In this paper we employ a basis of 27 different operators. The operators are contractible closed loops of links, and we take their traces. We include all loops of length 4 and length 6 as well as fifteen of length 8 and eight of length 10, with all their rotations and parity inverses. In addition to constructing the loops out of elementary links we also construct them out of link matrices that have been 'blocked' [10] using the algorithm described for example in [11]. The length of a blocked link is 2 n b a, where n b is the blocking level, and we use all blocked link matrices up to a maximum level that depends on the lattice size and is 6 for our largest 38 4 lattice which corresponds to blocked links of length 32a. Our loops are of length 4, 6, 8, 10 in units of the blocked links being employed in their construction.
As we increase β, and hence decrease a(β), we increase the spatial lattice size L s so that the lattice volume is approximately constant in physical units. We choose for our physical units our calculated values of the string tension, σ, so that we keep a √ σL s roughly constant. Since our blocked operators cover all length scales on such a lattice, we expect that their overlap onto the desired glueball states will only change weakly with β and this indeed appears to be the case. We calculate all the cross-correlators, , of all the operators φ i , φ j with given quantum numbers (see below) and this means that we can calculate the correlator of any linear combination in the vector space generated by the {φ i } and so we can perform a variational calculation to find the linear combination Φ 0 that minimises the energy (and hence maximises exp{−Et}). We then look for an effective mass plateau in the correlator of Φ 0 to provide our best estimate of the ground state energy. We now take the subspace of {φ i } that is orthogonal to Φ 0 and repeat the variational procedure to find the linear combination Φ 1 that minimises the energy in this subspace and hence provides, from its correlator, our best estimate of the mass of the first excited state. And so on. Note that one needs to ensure that the φ i one keeps are linearly independent so as to avoid singularities in the standard variational procedure.
The number of our operators φ i in a given quantum number sector can be very large (up to nearly 800) and we sometimes do not use all of our operators where doing so would make the analysis very unwieldy. In a few quantum number sectors the number is not large, but in any case the minimum is 30 operators, which should be large enough for the lightest couple of states, which is all we need in those cases.
In this section we start by describing the rotational quantum numbers of our lattice states and how these relate to the continuum spins. We then address the important issue of corrections due to the finite volume of our lattice. Finally we provide our results for the glueball spectrum with some demonstrations of how reliable are our mass estimates.

Quantum numbers
In the continuum theory the relevant quantum numbers of glueballs are their spin J, parity P and their sign under charge conjugation, C. (They also have momentum but since we are interested in the masses we choose to calculate states with zero momentum.) On a cubic lattice the rotational symmetry is restricted to the octahedral subgroup of the full rotation group and there are only five irreducible representations which we shall label generically as R. These five representations are often labelled as A 1 , A 2 , E, T 1 , T 2 . The dimensions of these are 1, 1, 2, 3, 3 respectively. The states that we calculate will belong to one of these representations and we will generically label them by R P C or more specifically as A P C 1 etc. As a(β) → 0 the glueball size diverges in lattice units and we recover the full rotational invariance where the states fall into the 2J + 1 multiplets labelled by J. Thus the lattice states can also be labelled by the continuum spin J to which they tend continuously as a(β) → 0. (Ignoring ambiguities due to level crossing, which should be appropriate once we are close enough to the continuum limit.) Of course at non-zero a the energy degeneracy will be broken except for that part of it that is enforced by the dimensionality of the lattice representations.
The way that the states labelled by the spin J of the full rotation group are partitioned amongst the representations of the subgroup of the symmetry rotations of the cube is wellknown (see for example Section 9.4 and Table 9.4 of [12]) and we reproduce the result for the lowest spins in Table 2. One can readily check that for each J the total number of states in the corresponding cubic representations, taking into account their degeneracies, is indeed 2J + 1. Where these states fall into different cubic representations they will not be degenerate, but one expects that the breaking of the degeneracy will be small once a(β) is small. Thus one can try to identify the spin J of a lattice state by seeing if there are a set of near-degenerate states in the corresponding cubic represntations listed in Table 2. In practice this works well for the lightest states, but breaks down for heavier states where the spectrum of states becomes sufficiently dense that an apparent near-degeneracy (within errors) is increasingly likely to be accidental. This method has been used to good effect, within its limitations, in for example [5], and it is the method we will employ in the present paper. There also exists a more powerful method which consists in ascertaining the approximate rotational properties of the states using an appropriate variation of a Fourier transform. This approach has been used successfully in [6,15] and more recently in [13] but requires a dedicated calculation of a kind that is beyong the scope of the present paper.

Finite volume corrections
The finite volume dependence of the glueball spectrum has two main components. The first is due to self-energy corrections that are modified by the periodic boundary conditions [14]. For example when a glueball emits a virtual glueball, one of the loop propagators can wind around the spatial torus before being reabsorbed. These corrections are typically exponentially suppressed in m G l s where m G is the mass gap and l s is the spatial periodicity, and since in our typical calculations m G l s ∼ 10 the resulting mass shift will be invisible within our statistical errors.
The second finite volume influence on the glueball spectrum consists of extra states whose energy → ∞ as l s → ∞, so that they disappear in the thermodynamic limit, but can have the same quantum numbers as glueballs and so can mix with the latter and be mistaken for them in a finite volume. These states are composed of flux tubes that wind around a spatial torus, and are generically called torelons. A single flux tube has non-trivial quantum numbers under a certain centre group transformation under which glueballs tranform trivially, and hence has zero overlap with contractible glueball operators and so is irrelevant here. (Although it is of course relevant to calculating the string tension as in Section 3.) The simplest non-trivial such finite volume state will consist of a torelon with a conjugate torelon closed around the same spatial torus. If the ground state energy of the torelon is E T , then the minimum energy of such a state will be E TT ∼ 2E T ∼ 2σl s and it will in principle appear as an extra state in our glueball correlator. Such a 'ditorelon' state may be a bound state or a scattering state. There will be a set of scattering states where each of the two torelons has an equal and opposite momentum (since our glueball operators have zero total momentum). One can estimate their energies using the allowed lattice momenta, e.g. ap x = 2πn/L x , although the actual momenta may have substantial shifts from these values due to the interactions. (It is only the total momentum of the state that has these quantised values.) Since the ditorelon states correspond to double trace operators while our glueball operators are single trace, we expect the overlap of the ditorelons onto our glueball operator basis to have some 'large-N' suppression which may mean that some of these finite volume states may become invisible in our actual calculations. However it has long been known that the lightest of these ditorelon states will certainly appear as the lightest glueball states in the A ++ 1 and E ++ representations once the volume is small enough. (See e.g. [16].) Since this simplest ditorelon operator is real, invariant under reflections and also invariant under rotations of π, it will only contribute to the A ++ 1 and E ++ representations. Nonetheless this poses a potential problem since the glueball states that are lightest and hence of greatest phenomenological interest are the J P C = 0 ++ , 2 ++ ground states and these fall into the A ++ 1 and E ++ representations respectively. So it will be important to make sure that the spatial volume is large enough for the ditorelon not to influence these mass estimates. For the J P C = 2 ++ state 2 of the 2J + 1 = 5 components are in the E ++ and 3 are in the T ++ 2 and since the ditorelon does not contribute to the latter this provides a useful crosscheck. However the 0 ++ only appears in the A ++ 1 . Since the ditorelon mass is ∝ l s one can make the spatial volume large enough for the ditorelon not to influence these important mass estimates. Then by comparing calculations on different volumes we can identify the location of any such ditorelon states in the glueball spectra. If we keep l s roughly constant in units of σ as we vary β we can keep such a state in the same location in the ordered spectrum and so obtain reliable continuum limits of the other 'true' glueballs.
In addition to the ditorelon one can form a tritorelon composed of three winding flux tubes. In SU(3) this has trivial centre symmetry quantum numbers and so can mix with local glueball states. The operator has an imaginary part and so it can contribute to some C = − states. However its energy is ∼ 3σl s and so it will not influence the lightest states given that we will work on lattices where the ditorelons are much heavier than the lightest states. Moreover this is a triple trace operator whose overlap onto the single trace operators we use in our glueball calculations is likely to be small enough that these states make no impact on our spectrum.
We also note that the above ditorelon is the simplest of a larger family of ditorelon states. A winding flux tube has a spectrum of excitations [7] and one can have a ditorelon consisting of an excited torelon as well as a ground state (conjugate) torelon, or indeed both torelons being excited. Such states may contribute to various R P C representations but will be heavier and will be less easy to identify since their volume dependence is much weaker at intermediate volumes. It is however not clear that they will have a large enough overlap onto our glueball operator basis to be relevant and we do not attempt to identify them explicitly in the present work.
To establish some explicit control over these finite volume ambiguities we have performed calculations on three different volumes at β = 6.235. This corresponds to one of our smaller values of a(β) and that enables us to obtain reasonably reliable mass estimates for our heavier states as well as our lightest states. In addition to the 26 3 26 lattice, which corresponds to the standard physical volume that we will use in our glueball calculations below, we also perform calculations on a smaller 18 3 26 lattice and a larger 34 3 26 lattice. In this comparison we will normally extract the masses from the same range of n t in the correlators on these three lattices so as to minimise enhancing the statistical fluctuations betweeen the calculations. This means that some of these masses may differ slightly from those presented in the following section.
We begin with the P, C = +, + sector of states, our mass estimates being listed in Table 3. The states are ordered by the effective mass aM ef f (t) of the corresponding correlator at t = a, which will normally, but not always, correspond to being ordered in the value of the mass taken from the effective mass 'plateau'. To make the comparison clearer we have left gaps where a state on a smaller lattice has no partners on the larger lattices. At the bottom of the table is listed the energy of a pair of non-interacting torelons on each lattice; presumably the mass of a ditorelon would be close to this. In the A ++ 1 sector of the 18 3 26 lattice we clearly see a state at aM ≃ 0.70 which does not appear in the spectra of the larger lattices. This has a mass close to twice that of the torelon and is undoubtedly a ditorelon. Similarly it is clear that the spectrum on the 26 3 26 lattice has an extra state: between the very well defined first excited state at aM ≃ 0.86 and the well defined states at aM ≃ 1.23 and aM ≃ 1.33 there are 4 states on the 26 3 26 lattice but only 3 states on the 34 3 26 lattice. Which of these states is the extra one is ambiguous since all 4 states have quite similar masses. We note that the third of these four states on the 26 3 26 lattice has a lower mass, aM ≃ 1.054, than the preceeding two; this is because it has a smaller overlap so that there is a larger gap between the value of the mass and the value of aM ef f (t = a). We have chosen this to be the extra state, mainly because we expect the ditorelon to have a reduced overlap onto the single trace operators of our basis. However from a purely matching point of view, we could equally well have chosen the next higher state at aM ≃ 1. 19. In both cases this extra state has a mass close to twice that of that lattice's torelon, as expected. In the E representation on the 18 3 26 lattice it is the ground state with aM ≃ 0.66 that is definitely an extra state and hence the candidate for a ditorelon, and on the 26 3 26 lattice it is the excited state with aM ≃ 1.206 that appears to be the extra state and hence a ditorelon. If we now compare the other states from the 34 3 26 and 26 3 26 lattices we see that they agree well with each other: there are no differences that are greater than 2 standard deviations, except 2.6σ for the A 2 ground state. More surprisingly, the same is true for most (but not all) states on the small 18 3 26 lattice.
We perform a similar comparison for the other values of P and C. In Table 4 we do so for the P, C = −, + sector of states. We see that all the masses on the 26 3 26 and 34 3 26 lattices are within 2σ of each other, with most differences much less than that. So that it appears that the 26 3 26 lattice is large enough within our statistics. This is not the case for the 18 3 26 lattice which shows large differences for a number of states. The conclusion is the same for the P, C = +, − and P, C = −, − sectors of states, as we seefrom Table 5 and Table 6 respectively.
While we can provide a plausible identification of the ditorelon contribution to the spectrum, as we have done above, it would be useful to complement that with an alternative calculation. One way to do this is to include operators that are explicitly ditorelon. We shall now provide a first limited step in that direction at β = 6.235 where we have performed our above finite volume comparison. We take operators that are a product of a zero momentum Polyakov loop times its conjugate. (Note that in so far as the torelons interact with each other this mometum is not conserved so this operator should couple to ditorelons with non-zero relative momenta.) We include only the three largest blocking levels since, as we have checked, it is these that provide essentially the whole of the flux tube ground state wave function at this value of β. These ditorelon operators should contribute primarily to the A1 ++ and E ++ representations. Taking into account the fact that the two torelon operators in the ditorelon can have different blocking levels this means that we have 6 distinct double trace operators in the A1 ++ representation and 12 in the E ++ (since it has dimension two). Note that to contribute efficiently to other representations we would need to include Polyakov loop operators with non-trivial quantum numbers, employing suitable transverse deformations, and with non-zero momenta (equal and opposite so as to give zero total momentum). This would also give a more complete basis for the A1 ++ and E ++ representations, but would go well beyond the scope of the present work. In the same exploratory spirit of this calculation we include only 12 different loops in our glueball basis, rather than the 27 used in most of our other calculations, and we have lower statistics (by a factor of two) and we restrict ourselves to the three largest blocking levels. That is to say, we have 36 and 108 single trace glueball operators for A1 ++ and E ++ respectively. We perform the calculations on our smaller 18 3 26 lattice as well as on our standard 26 3 26 lattice. We begin by showing in Fig. 2 the effective mass plots for the lightest three A1 ++ and E ++ states that we obtain on the 18 3 26 lattice using only our ditorelon operators. We compare these to the lightest three states obtained with two free torelons on the same lattice. (We also indicate the lightest free tritorelon energy, although its relevance is not clear given that it should have a reduced overlap onto our ditorelon basis.) In interpreting this plot one must bear in mind that these ditorelon operators should have a non-trivial projection onto glueball states, and where these are lighter they will drive the behaviour of the effective masses at larger t. We begin by noting that the E ++ ground state has the same energy as the free ditorelon suggesting a negligible interaction energy. And it is the lightest E ++ state overall. By contrast the lightest A1 ++ ditorelon shows a nice plateau at intermediate t but one which is about 15% below the free ditorelon energy, indicating a significant interaction energy in this representation: presumably a genuine bound state. At larger t the effective energy appears to drift lower, presumably towards the mass of the lightest A1 ++ glueball. Finally we remark that while the excited ditorelon states begin roughly within striking range of the corresponding free ditorelon energies, they appear to be drifting rapidly lower, presumably driven by their overlaps onto the lighter glueball states, and there is no evidence that they harbour any bound states. We turn next to Fig. 3 where we show effective mass plots for our lightest four A1 ++ states on the 18 3 26 lattice. We show the results using our basis of 36 single trace operators and also when this basis is augmented by the 6 ditorelon operators. We see that the ground state and the third excited state are unaffected by the ditorelon operators; presumably because they are far enough from the energy of the free ditorelon. In contrast, the overlaps of the first and second excited states onto the operator basis are both significantly increased when the ditorelon operators are included. This is consistent with our earlier conclusion that the first excited state is the extra ditorelon state when compared to the spectra on larger volumes, but is not so specific. In Fig. 4 we repeat the plot for the four lightest E ++ states, and the picture is very similar. The overlap of the two lightest states is increased by the inclusion of the ditorelon operators, but the higher states, which are much further in energy from the free ditorelon energy, are unaffected. Again while this is consistent with our earlier conclusion that the lightest E ++ state is in fact a ditorelon, it is more ambiguous. The ambiguity is due to the fact that the ditorelon state, if any, already appears in the spectrum obtained with our single trace glueball operators. That is to say there is a strong mixing between the single and double trace operators, at least on this smaller volume. The corresponding spectra obtained on the larger 26 3 26 lattice are shown in Fig. 5 and Fig. 6. The picture here is quite different in two main ways. The first and trivial difference is that the individual flux tubes are longer and more massive, so that the ditorelon is much more massive and appears amongst the higher excited states. This should, in principle, make its identification more difficult. The second and more interesting difference is that the inclusion of the ditorelon operators appears to lead to an extra state in both the A1 ++ and E ++ representations, i.e. the second excited state when using the ditorelon extended operator basis, and these states are close in energy to that of a free ditorelon. There thus appears to be little ambiguity in identifying these as ditorelon states. Since we use the same single trace basis on the 18 3 26 and 26 3 26 lattices the fact that these states appear with that basis on the former lattice but not on the latter suggests that the overlap of the ditorelon onto the conventional glueball basis decreases as its length increases. This means that the value of E ef f (t = a) will increase relative to the ditorelon energy, so that it will appear higher in the spectrum and with larger errors. If we make the single trace operator basis much larger, then the overlap should increase, perhaps sufficiently for the ditorelon states to reappear as low-lying states in that spectrum. This is indeed what appears to happen: our finite volume comparison between the 26 3 26 and 34 3 26 lattices used a much larger basis and this appears to resurrect the ditorelon states in the conventional glueball spectrum. We can conclude that even though these ditorelon calculations are intended to be exploratory, they already reveal some usefulness in complementing our earlier comparison of spectra across different volumes.
We conclude that the important finite volume effects on a lattice that has the physical spatial volume of a 26 3 lattice at β = 6.235 are the extra ditorelon states in the A ++ 1 and E ++ representation. As we see in Table 3 this is probably the 5'th A ++ 1 state and the 5'th E ++ state. We will work on lattices at other values of β which are chosen so as to be roughly the same physical size (in units of the string tension) and so provided there is no strong lattice correction to the level ordering, we can assume that the location of the ditorelons in the spectrum of these extra states is the same. When we come to perform continuum extrapolations of our lattice masses we shall limit ourselves to the lightest 4 A ++ 1 states and the lightest 4 E ++ states, which means that if our identification of the ditorelon states in Table 3 is correct, then we should not have any contamination from these ditorelon states in our final results. This is not guaranteed of course: there are several states close together in the relevant energy range which may enhance mixing and changes in level ordering on slightly different volumes. Nonetheless, by the same token this means that our energy estimates of the genuine glueball states should be quite accurate.

Glueball masses
We calculate a number of the lightest glueball masses for each set of R P C quantum numbers. Since our mass estimates becomes less reliable as the the mass becomes larger, our first task is to provide some evidence for their reliability. This we do by providing effective mass plots for all the states to which we will eventually assign a cotinuum spin J, i.e those states that are of greatest phenomenological interest. We do so for our β = 6.50 calculation on a 38 4 lattice, because it has the smallest lattice spacing and so is the closest to the continuum limit and has the finest resolution of the glueball correlators.
We begin, in Fig. 7, with the lightest and first excited J P C = 0 ++ , 2 ++ glueballs. The 0 ++ ground state is the lightest of all the glueballs and constitutes the theory's mass gap. We show also, as horizontal lines, our best estimates of the various masses. We see that we have well defined effective mass plateaux in all cases and that these begin from either t = a or from t = 2a. In the case of the 0 ++ ground state we see a noticeable increase in the aM ef f (t) at larger values of t, but since we know that the exact value of aM ef f (t) must be monotonically decreasing as t increases this rise must be a statistical fluctuation. In general the effecive masses do show fluctuations up or down at large t which is presumably due to the incompleteness of our straightforward error analyses. In addition for excited states the incompleteness of our variational basis means that our best operator for an excited state may contain components of lower mass states. When these components are small, which we try to achieve with our operator basis being large, one will see the desired effective mass plateau over an intermediate t range followed by a slow drift at larger t to a secondary plateau (and so on). So sometimes a downward drift at large t will be real rather than a fluctuation. In any case our interest is in the initial effective mass plateau. Returning to Fig. 7 what we see is that in the cases shown the mass estimates are robust. It is also clear that the T 2 and E states that make up the five components of the J = 2 state are degenerate within errors as we would expect once a(β) is small enough and the spatial volume is large enough.
Note that since the states of T 1 and T 2 come in degenerate triples and those of E in doublets, we shall average the corresponding correlators to produce single values for the masses. That is to say, whenever in this paper we refer to a state of one these representations we are actually referring to an average of the 3 or 2 states that would be exactly degenerate with infinite statistics. This is what we have done in Fig. 7 for the E and T 2 masses. So the former state actually represents two states and the latter three giving a total of five as needed for J = 2.
In Fig. 8 we display the effective masses of some P, C = −, + states. The lightest and first excited 0 −+ have well defined plateaux starting at around t = 2a. As do the ground and first excited 2 −+ , with a near-degeneracy between the E and T 2 components. The lightest three T −+ 1 triplets are heavier and have worse overlaps: for two of the states there is evidence for a plateau starting near t = 2a but less evidence for the third. We assign all three to J = 1 (see Section 5.3) but with modest confidence for some. In Fig. 9 we display the effective masses of some P, C = +, − states. The lightest and first excited 1 +− have well defined plateaux starting at around t = 2a. The 2 +− are reasonably convincing. The near degeneracies needed to identify the J = 3 state are convincing, and the effective mass plateau is quite plausible. Finally in Fig. 10 we display the 1 −− and 2 −− effective masses and it is clear that here we can only make an estimate by assuming that, as in the previous states, the plateau starts from t = 2a.
The message here is that our lighter states are robust but that as the mass becomes higher the extraction of a mass will typically become more ambiguous. In addition we have seen that for most states the effective mass plateau begins at t = 2a, or even t = a, in this calculation at β = 6.50. Since our operator basis is designed so that the correlators scale approximately with the physical mass scale, this implies that for β ≤ 6.0 the effective mass plateau will usually begin at t = a. Since at our smaller β values the lattice spacing a and hence the lattice mass aM are considerably larger, an effective mass plateau is harder to establish except for the lightest states, and in these cases the fact that we can assume the plateau typically begins from t = a becomes useful.
We list our resulting mass estimates in Tables 7 to 13.

Continuum limit
We take the masses in Tables 7 to 13 and express them as M/ √ σ in units of the string tension that we have already calculated. It is this dimensionless ratio that we extrapolate to the continuum limit in Section 5.1 using eqn (4). However while expressing the continuum masses in terms of the string tension is theoretically elegant it is not ideal if we are thinking of phenomenological applications where we want to express masses in GeV units. For this purpose one typically expresses masses in units of the r 0 scale [17] and hence into 'GeV' units, as discussed in Section 5.2. Of course this still involves some ambiguity because the SU(3) gauge theory and QCD are two different theories and there can be no single rescaling that will match their physics. In addition to expressing the masses in physical units we will obviously want to identify the continuum spins of as many states as possible, which we do in Section 5.3.

Continuum mass ratios
We use eqn(4) to extrapolate our mass ratios, M/ √ σ, to the continuum limit. We choose to use just the leading O(a 2 σ) correction and where necessary we drop masses at the largest values of β in order to obtain acceptable fits. We list in Table 14 the resulting continuum extrapolations for states with various quantum numbers. Some of the entries in Table 14 are starred, indicating particularly poor fits where the χ 2 per degree of freedom is greater than 2.5. How good the fits are is cleary an important issue. This involves not only the χ 2 per degree of freedom but also the range fitted. So to provide evidence for the reliability of these continuum extrapolations we shall show explicitly a number of them below. Since the most important states are those whose continuum spin we are able to identify, these are the ones we shall show.
Note that by deliberately limiting ourselves in Table 14 to the lightest four A ++ 1 and E ++ states we expect to exclude from our spectra the finite-volume ditorelon states that were discussed in Section 4.2.
We begin by showing in Fig. 11 two extrapolations of the very lightest glueball mass, the A ++ 1 ground state which has J = 0. One extrapolation keeps just the leading O(a 2 σ) correction, while the other includes an additional O(a 4 σ 2 ) correction. As is well known this state has a deep dip around β ∼ 5.5 where one has a crossover between the regimes of weak and strong bare coupling. This dip is a reason for attempting to improve the linear fit with a higher order correction. However, since our calculations go to quite large β the continuum limits obtained with these two extrapolations turn out to be identical within errors: (21) and 3.391(23) respectively. We shall employ the first value, primarily because it makes no difference, but also because the strong to weak coupling transition involves a change in the functional form of the expansion in β, and there is no guarantee that simply including a higher order term derived in weak-coupling is the right way to parametrise it. In Fig. 12 we display the continuum extrapolations of the lightest four A ++ 1 states. The lightest two have J = 0; and the very lightest has already been shown in Fig. 11. The third state probably belongs to J = 4 (see below) and we have no J assignment for the fourth state.
As we see, the continuum limits appear to be quite unambiguous in all cases.
In Fig. 13 we move on to the lightest and first excited E ++ and T ++ 2 states. Recall that since the states in the E come in degenerate pairs (up to fluctuations), we average the correlators to produce a single mass. Similarly for the triply degenerate T 2 (and T 1 ). So, for example, the lightest T 2 mass corresponds to the lightest three states and the first excited T 2 mass corresponds to the next three states. This labelling of the states is used throughout the paper. Returning to Fig. 13 we see that for these states, which we shall see below correspond to J = 2, there is remarkably little lattice dependence over our whole range of β, leading to continuum extrapolations that are unambiguous.
Some heavier states with P = + and C = + are displayed in Fig. 14 and Fig. 15. Here one sometimes sees a scatter of points indicating worse fits, and this is presumably due to the increasing uncertainty in extracting the masses of heavier states. Nonetheless the extrapolations are plausible within the errors quoted in Table 14.
Moving now to states with P = − and C = +. The lightest two A −+ 1 states plotted in Fig. 16 show small O(a 2 ) corrections and are simple to extrapolate to the continuum. The lightest two E −+ and T −+ 2 states, plotted in Fig. 17, also show small corrections and hence easy extrapolations. In Fig. 18 we plot the lightest three T −+

Physical units
We have used the string tension as the energy scale for our continuum extrapolation of the glueball masses both because it is a quantity that is of fundamental dynamical significance in the theory and because one can calculate it very accurately in our lattice calculation, as described in Section 3. However for phenomenological applications one would like to present the spectrum in 'GeV' units. Strictly speaking this is of course not possible: the pure gauge theory and QCD, with its various quark masses at their physical values, are two different theories. Nonetheless phenomena such as the OZI rule and lattice calculations showing that SU (3) is 'close to' SU(∞), both for the pure gauge theory [18] and also when quarks are added [19], motivates attempting a translation into GeV units, something which has traditionally been done in past lattice calculations.
One possibility is to use the QCD value for the string tension σ. Since the QCD string will contain virtual qq pairs -which eventually cause it to break -this is not unambiguous. For a rough estimate one might look to the long-distance linearly rising piece of the potential describing heavy quarkonia, but since the dynamics of these heavy states is primarily determined by the interactions at short distances this is not very constraining. The traditional approach is to take the slope of the observed Regge trajectories and to use a simple asymptotic rotating string model, and this gives something like √ σ ∼ 440 MeV . How good this interpretation is for the relatively low J states that in practice determine these Regge trajectories is of course uncertain. An alternative scale that has been widely used is to determine the distance r at which the slope of the static (heavy) quark potential, the force F (r), has some value. For example r 2 F (r)| r=r 0 = 1.65 determines the scale r 0 [17] and there are arguments that this should be relatively insensitive to the inclusion of quarks. An analysis of various lattice QCD calculations with at least 3 light quarks leads to the estimate [20] We have not attempted to calculate r 0 /a in our calculations but there are calculations covering our range of β in [21,22] which provide interpolating formulae for r 0 /a| β . We have used these formulae to interpolate to our values of β and have extrapolated the resulting values of r 0 √ σ to the continuum limit. We find good fits which agree within errors and give us using the value for r 0 in eqn (8). Using this scale we can convert the continuum mass ratios in Table 14 to GeV units and this results in the masses listed in Table 15.

Continuum glueball spins
In the continuum theory a glueball transforms according to an irreducible representation of the (proper) rotation group SO(3) and so will possess some continuum spin J, and it will be part of a set of 2J + 1 degenerate states. All such states will in principle be included in the set of states we have obtained by extrapolating our lattice calculations, albeit with the obvious caveats arising from the limitations of our basis, our statistics etc. However our calculated states are all labelled by the irreducible representations A 1 , A 2 , E, T 1 , T 2 of the octahedral subgroup O of SO (3), so we need to map back from these to the multiplets corresponding to a given J. The mapping can be readily derived (see Section 9.4 of [12]), and is summarised for the lowest values of J in Table 2. (There is a simple relation allowing one to calculate the mapping for any J once one has the mappings for J ≤ 5, as explained in Section 9.4 of [12]). Taking into account the triple degeneracy of the T 1 and T 2 representations, and the double degeneracy of E, we can see that the sum of states is indeed 2J + 1 in each case, as it should be. We have, optimistically included the mapping up to J = 8 although, as we will see, we are not able to calculate masses beyond J = 4. In a calculation with infinite statistics identifying states of spin J from our continuum calculations would be unambiguous: we would simply search for 2J + 1 degenerate states located in the towers of states in each of the appropriate cubic representations listed in Table 2. In practice, with finite statistical and systematic errors and with the energy gaps between neighbouring states in each tower decreasing with increasing energy, this strategy has severe limitations. An alternative approach would be to exploit the fact that even on a cubic lattice the full SO(3) rotational invariance becomes a better approximation on physical length scales as a → 0. Once one calculates the variationally selected glueball wave-functional at each β one can analyse its approximate rotational properties so as to determine the probable continuum spin and see if it approaches the corresponding continuum angular wave-function as a(β) decreases. A pioneering calculation of this type [6,15] in D = 3 + 1 gauge theories appears promising. Moreover in D = 2 + 1, where the rotation group is Abelian so that there is no J-dependent degeneracy to exploit, this technique has proved quite powerful in recent work [13] that builds upon the earlier exploratory studies [6].
In the present paper we attempt to identify states of various spin J using the decomposition in Table 2 applied to our results in Table 14. The arguments will be as follows, leading to the identifications presented in Table 16. We will also find it useful to refer to the spectra obtained at various values of a, on the assumption that lattice corrections are small, as appears to be the case for most states.
We begin with the P, C = +, + sector of states. First we note that the lightest two masses in R P C = A ++ 1 are very much lighter (as always, compared to errors) than any states in other R ++ . So these are, unambiguously, the ground and first excited J P C = 0 ++ states. The lightest two E ++ states (each doubly degenerate) are nearly degenerate with the lightest two T ++ 2 states (each triply degenerate) but not with any other states (by a large margin). Moreover this near-degeneracy is the case at all values of β as we see in Fig. 13. So these are unambiguously, the ground and first excited J P C = 2 ++ states. One might wonder if the lightest T ++ 1 state in Table 14 is J = 1? We note however that it is essentially degenerate with the first excited T ++ 1 and also with the lightest A ++ 2 and with either the second ot third excited T ++ 1 state, which are nearly degenerate. In addition the second excited A ++ 1 and E ++ have similar masses. This suggests that we have a 3 ++ ground state and a 4 ++ ground state as in Table 16, which happen to have very similar masses. If we look at the states as a function of β, with the 3 ++ in Fig. 14 and the 4 ++ in Fig. 15, then the assignments look plausible, particularly in the former case. However it is clear that these choices are no more than plausible, because all the states shown in Fig. 14 and Fig. 15 are very close in mass. An alternative interpretation could be that the lightest T ++ 1 is J = 1 and that the remaining states listed above make up the J = 6 ++ ground state. However this would mean that the J = 6 ground state is lighter than the J = 4 ground state, which appears less likely to us. Hence our choice. However this discussion illustrates the limitations of our approach.
We now turn to the P, C = −, + sector of states. The lightest two A −+ 1 states have no partners in other representations and so are unambiguously the ground and first excited 0 −+ states. The lightest two E −+ states are nearly degenerate with the lightest two T −+ 2 states but not with any other states, so it is unambiguous that they are the ground and first excited 2 −+ states. Again this is reinforced by noting the near-degeneracy at all values of β as shown in Fig. 13. The lightest three T −+ 1 states (each triply degenerate) have no obvious partners in other representations, so we have labelled all of them as being 1 −+ states. We do so with some unease: why the apparent degeneracy of these three states? (See also Fig. 18.) Perhaps only one is truly 1 −+ and the other two are nearly degenerate because they appear together as components of a higher spin state? Unfortunately we have no plausible partners in other representations to realise this suggestion.
In the P, C = +, − sector of states the lightest glueball is the T +− 1 and it has no partners elsewhere: hence it is unambigously the 1 +− ground state. The second excited T +− 1 has no partners elsewhere so we are confident that it is the first excited 1 +− glueball. Both their continuum limits are robust as we see from Fig. 20. The A +− 2 ground state is nearly degenerate with the T +− 2 ground state and with the first excited T +− 1 , with no nearby states elsewhere. And, as we see in Fig. 21, this is so at all β. Hence we are confident that this is the 3 +− ground state. There appears to be little ambiguity in identifying the 3 +− ground state with the A +− 2 and T +− 2 ground states and the first excited T +− 1 , since these states are not very heavy and it is clear that there are no partner states elsewhere. We further note that the ground state E +− and the first excited T +− 2 are nearly degenerate and so we assign them to the 2 +− ground state. This is also the case throughout our range of β as we see in Fig. 19. However the mass and errors are large, and there are other states with masses about 2σ away, so the assignment has some uncertainty. However looking at the corresponding spectra at β = 6.50, 6.338, 6.235 does increases the evidence for no partner states elsewhere. Finally we can speculate that the A +− 1 ground state partners the first excited E +− , the second excited T +− 2 and the third excited T +− 1 to form the 4 +− ground state. This claim is slightly weakened by what we observe in the spectra at smaller β, as shown in Fig. 22 and is somewhat flimsy.
Finally, in the P, C = −, − sector the states are all heavy making it difficult to calculate much of the spectrum with sufficient accuracy to be able to search for near-degeneracies. The lightest three states are the E −− , T −− 1 and T −− 2 ground states and there are no potential partners to these. Given the possibilies listed in Table 2, there is no ambiguity in assigning the E −− and T −− 2 ground states to give the 2 −− ground state, and the T −− 1 ground state to be the 1 −− ground state.
Having made these identifications we average the masses of the states in Table 14 that correspond to a state of given J P C , taking some account of any splittings between the masses in the different cubic representations. This leads to the masses listed in Table 17, for a number of continuum states with various values of J P C . Where we have some hesitation about the spin assignment, as discussed above, we qualify the value with a single star. Two stars indicate a speculation.
Our predictions for the glueball masses are given in terms of the string tension. If one wishes one may express the masses in units of the mass gap, the mass of the lightest 0 ++ glueball, simply by multiplying the various mass ratios by the inverse of m 0 ++ / √ σ. Since many of these states may be of phenomenological interest it is useful to provide an estimate in physical GeV units and this we do in Table 18 using the value of √ σ estimated in Section 5.2.

Comparisons
It is clearly useful to compare our calculations to earlier ones. Since the number of these is too large to attempt a comprehensive comparison, we shall choose a few representative calculations (within which one may find references to other calculations).
One of the earlier realistic attempts to calculate the ground states of all the representations of the cubic group was in [1], and this was extended to a smaller value of a(β) in [3]. The values for the masses of the states of greatest phenomenological interest, i.e. those that have masses < 3GeV, are consistent with ours (within two standard deviations) albeit with statistical errors at least ten times larger than ours. The point here is that the approximate masses of these states, and their level ordering, have been established for a very long time.
When one reduces the statistical errors, one also needs to reduce the systematic errors so that they do not become the dominant error. An important example is the need to improve the extrapolation to the continuum limit: in [1]. the mass ratios were consistent within errors with being independent of β, so that a constant extrapolation was justified, but with smaller errors a variation with β becomes visible and one needs to include at least an O(a 2 σ) correction. One also needs to go to smaller values of a(β) to better control the continuum extrapolation, as for example in [3].
Another systematic error arises in the calculation of heavier glueballs. One needs the mass in lattice units, aM, to be not too large if one is to observe some kind of effective mass plateau in the correlation function. One way to achieve this is to go to small enough a(β). (At the same time reaping benefits in the continuum extrapolation.) An alternative is to note that it is the temporal lattice spacing a t that appears in the correlator, so if one works on lattices where the spatial, a s , and temporal lattice spacings differ with a t ≪ a s one can achieve modest values for a t M even if M is large, while working on L 3 s L t lattices with L s ≪ L t which are of a manageable size. This technique proved invaluable in some of the earliest calculations of the lightest scalar and tensor glueballs [23] in providing evidence for effective mass plateaux. It has been used subsequently to good effect in a comprehensive calculation of the mass spectrum that includes many heavy glueball states [4,5]. One weakness with such calculations is that they typically involve larger spatial lattice spacings a s than those on symmetric lattices, and hence risk larger sytematic errors in the continuum extrapolation where the leading correction is O(a 2 s ). So, for example, the smallest value of a s in [4,5] is a s M g ≃ 3a t M g ≃ 0.84, where M g is the mass gap, whereas the calculations in the present paper go down to a s M g ≃ 0.35. A second issue is that it is difficult to estimate the corrections to the inputted tree level value of a s /a t , and the introduction of a physical scale such as r 0 or √ σ will depend on this ratio.
While the resulting uncertainty may be only at the (few) percent level, this becomes significant in precise calculations. A further important systematic error arises once one is more ambitious and wishes to calculate some significant number of excited states in each R P C sector. The problem arises from the finite size of one's operator basis: one will completely miss states which do not have a substantial overlap onto the basis. Such missing states can lead to problems in a comparison with theoretical predictions and in the identification of the near degenerate multiplets associated with a given continuum spin J. Here the best one can do is to make the operator basis as large as practical, as we have done in the present work.
If one is interested in the spin J that a lattice state extrapolates to, the identification of the 2J + 1 nearly degenerate states in the appropriate cubic representations becomes rapidly problematic as the mass increases. Typically the density of states increases with increasing mass, and with finite errors that also increase with the mass, any apparent degeneracies are likely to be accidental. An alternative approach to determining J is to use sets of loop operators that are approximate rotations of each other by angles that are some fraction of π/2. Since the important operators have an extension that is constant in physical units, the approximate rotation invariance becomes better as a(β) deacreases, and one can resolve higher spins. This method has been used to good effect in [6] where the spins of an extensive set of states in various R P C representations were determined.
We now turn to a comparison with the specific results of [5,6]. In our calculations we have a large basis of 27 loops, with each typically constructed out of 5 or 6 different levels of blocked links, which is a considerably larger basis of operators than that used in the earlier calculations we have discussed above. This has enabled us to calculate more excited states within each quantum number sector. Having these extra states also helps in the identification of their continuum J using the 2J + 1 degeneracy partitioned amongst the lattice states as in Table 2. In practice we are able to so identify more states than in [5]. However apart from that we do broadly agree: for the 12 states listed in Table XXI of [5] we agree with the spin assignment and mass (within errors) of 10 and only have serious doubts about 2 of the heaviest states, the J P C = 0 +− and the 3 −− . The former state is based on the very heavy A1 +− state which is quite possibly nearly degenerate with excited E +− , T 1 +− and T 2 +− states making it 4 +− rather than 0 +− . The claimed 3 −− is based on the A2 −− ground state and we do not see an excited T 1 +− with a mass that is plausibly degenerate. But only better calculations will be able to tell. As for the lightest phenomologically interesting states, our masses agree remarkably well except for a significant discrepancy for the 0 ++ ground state. (Which is slightly increased by the fact that the scale used in [5] is r −1 0 = 410MeV while we use a later determination, r −1 0 = 418MeV .) Since this state has large lattice spacing corrections, we put this down to our better control of the O(a 2 s ) corrections in taking the continuum limit. The very different method used in [6] for determining the spin J of lattice states is more powerful and allows more states to be categorised than we are able to do. This despite our much higher statistics and larger range of β which lead to much smaller statistical errors as well as a better control of the continuum limit. Within these large errors we are in agreement for most spins and masses, as we see by comparing our Table 14 to Tables 7.10,7.11 in [6]. There are some differences: in the P C = ++ sector our candidate for the J = 4 ground state seems too light; in the P C = −+ sector two of our three nearly degenerate T 1 states might be part of the J = 5 ground state according to [6]. We note that with this method one could identify states with spin up to J = 6 in the now rather old calculations of [6], something we could not hope to do despite our statistical errors being a factor of 3 or 4 smaller and our range of a(β) being significantly larger.

Topological fluctuations
An interesting property of Euclidean non-Abelian gauge fields is that they possess non-trivial topological properties, characterised by a topological charge Q which is integer-valued in a space-time volume with periodic boundary conditions. This charge can be expressed as the integral over Euclidean space-time of a topological charge density, Q(x), where Recalling that the plaquette matrix U µν (x) = 1 + a 2 F µν (x) + .... on smooth fields, one can write a lattice topological charge density Q L (x) as While this definition is adequate for sufficiently smooth fields, it lacks the explicit reflection properties of the continuum operator in eqn (10), since all the plaquettes U µν (x) are defined as forward going in terms of our coordinate basis. So it needs to be extended to include terms with backward going plaquettes as well, if it is to be used on the 'rough' Monte Carlo generated fields, and so we employ the version of this operator that is antisymmetrised with respect to forward and backward directions [24]. The fluctuations of realistic lattice gauge fields affect the observed value of Q L (x) in two important ways. Firstly we note that the fluctuations of Q L (x) are related to expectation value of the composite operator Q 2 L (x) whose operator product expansion contains the unit operator [24]. So these fluctuations are powerlike in β while the average of Q L (x) is O(a 4 ) and hence exponentially suppressed in β. Thus as β increases the fluctuations around Q L (x) and Q L = x Q L (x) diverge compared to the physically interesting mean values. In addition the lattice composite operator also receives a multiplicative renormalisation Z(β) such that Z(β → ∞) = 1 but which at accessible values of β strongly suppresses the charge [25].
In practice all this means that one cannot extract the topological charge of a typical lattice gauge field by directly calculating Q L = x Q L (x) on that gauge field. However we note that the fluctuations obscuring the value of Q are ultraviolet, while the physically relevant topological charge is on physical length scales. Thus if we perform a very limited local smoothening of the fields to suppress the ultraviolet fluctuations, this should not affect physics on long distance scales, and the value of Q L = x Q L (x) calculated on these smoothened fields should provide a reliable estimate of Q. Moreover, recalling that the total topological charge of a gauge field is unchanged under smooth deformations, we can expect that even under a moderately large amount of continued smoothening the value of Q L will not change, even though Q L (x) itself does gradually change. One convenient way to smoothen the gauge fields is to locally minimise the action. Such a 'cooling' of the original 'hot' lattice gauge field [26] involves sweeping through the lattice one link at a time, precisely like the Monte Carlo except that one chooses the link matrix that minimises the total action of the plaquettes containing that link matrix. This is a standard technique that one can find described in more detail in, for example, [27]. An alternative and attractive smoothing method with purturbatively proven renormalisation properties is the gradient flow [28,29,30] which has been shown to be numerically equivalent to cooling [31,32,33]. Cooling appears to perform nearly two orders of magnitude faster than the gradient flow and since we are aiming for large statistics we adopted cooling. After the first couple of cooling sweeps the fields are already quite smooth, as we shall see below. Since we are minimising the action and since in the continuum the minimum action field with a given Q is a multi-instanton field, we expect that under systematic cooling the lattice field will be driven to become some multi-instanton field, which one can see by calculating the distribution Q L (x) on such a field. Of course, because of the discretisation of space-time the topological properties of lattice fields are only approximately like those of a continuum field. One can deform a large instanton by gradually shrinking its non-trivial core and on a lattice this core can shrink to within a hypercube. At this point what was an instanton has been transformed into a gauge singularity and the value of Q will now differ from its original value by one unit. (Equally, one can gradually grow an instanton out of a hypercube.) This process can occur during cooling but it can equally occur during the course of our Monte Carlo. In the latter case, it is these changes in Q that allow us to sample all possible values of Q and hence maintain the ergodicity in Q of our Markov process. As β ↑ the distance between physical and ultraviolet scales grows and these changes in Q become increasingly suppressed -in fact more strongly than the usual critical slowing down (see below).
One of our reasons for studying the topology of our SU(3) gauge fields is to establish that in our range of β we have not in fact lost ergodicity in the topology Q of the fields. In the process we shall provide some insights into how reliable is the cooling method for establishing the topological charge of the lattice gauge fields. As a side-product we will also provide an accurate calculation of the topological susceptibility, which has been much studied because of its appearance in sum rules for the lightest pseudoscalars in QCD. For a more detailed discussion of these issues we refer to, for example, [27].

Topology on the lattice
Since we calculate the topological charge Q by cooling the lattice gauge fields, it is relevant to ask how (un)ambiguous is the identification of Q on those cooled fields. (Note that from now on we will use the labels Q and Q L interchangeably unless specified otherwise.) The answer will depend on the value of β and since the intrinsic ambiguity in defining topology is a lattice artifact we might expect that any ambiguity is greatest at the smallest values of β where the lattice spacing is largest. We therefore look at our smallest value of β, β = 5.6924, an intermediate value, β = 6.0625, and our largest value, β = 6.50. We note that the final values of Q that we shall quote are those that are obtained after 20 cooling sweeps.
We begin with β = 5.6924, which we expect to be our 'worst' case. In Fig. 25 we show the distributions of Q obtained after 4, 8 and 20 cooling sweeps on a total of 96000 8 3 16 lattice gauge fields. We see that after 20 cooling sweeps any ambiguity in identifying Q -arising from the overlap of the tails between the peaks -is at most at the percent level. The peaks are shifted away from integer values due to the fact that the typical instanton size, ρ/a in units of the lattice spacing, is not large at this β. (One can calibrate the shifts using classical instanton solutions discretised on a lattice.) However, while the value of Q is well defined after 20 cooling sweeps, it is clearly quite ambiguous after 8 cooling sweeps, and totally invisible after 4 cooling sweeps.
We turn now to β = 6.0625 where we calculate Q on a total of 40000 14 3 20 lattice gauge fields. In Fig. 26 we show the distributions of Q obtained after 4 and 20 cooling sweeps. It is clear that after 20 cooling sweeps any ambiguity in identifying the value of Q is totally negligible. More interestingly, the value after only 4 cooling sweeps is reasonably well defined with only a minor ambiguity. We also note that the peaks have moved closer to integer values, as should be the case since we expect the instanton sizes to be larger in lattice units at this larger β.
Finally we show in Fig. 27 the distributions of Q obtained after 4 and 20 cooling sweeps on 44000 26 3 38 lattice gauge fields generated at β = 6.50. Now any ambiguity is completely negligible even after only 4 cooling sweeps.
We see from the above that after 20 cooling sweeps the value of Q is very well defined over our whole range of β. Now, although it is plausible that the value of Q after 20 cools reflects the value of Q of the original 'hot' field configuration, it is worth asking how well that expectation holds up. The first obvious comment is that the smaller the number of cooling sweeps, the more convincing is the calculated value of Q. This is because the cooling is a local process to the same extent as the Monte Carlo heat bath, so a few cooling sweeps will deform the original fields only on correspondingly small distance scales. However one can ask if the fields with a certain value of Q after a very few cooling sweeps have the same value of Q after 20 cooling sweeps. To address this kind of question better we have taken a sample of 2750 fields generated at β = 6.50 and measured the value of Q for 1,2,3,4,5,10,15,....,50 cooling sweeps. In Fig. 28 we plot the values of Q after 2 and after 20 cooling sweeps. We see that even after only 2 cooling sweeps the values of Q are peaked around near-integer values with little overlap so that we can assign each value of Q to an integer, Q I , with a high level of confidence. And since 2 cooling sweeps can only deform the fields over the shortest distance scales we can be confident that this topological charge provides a very accurate representation of the topological charge of the original lattice gauge fields. The interesting question now is whether the fields with a given Q I after 2 cooling sweeps maintain that same value when we continue cooling the lattice gauge fields to, say, 20 cooling sweeps. The answer is yes to a very good approximation. As an example we show in Fig. 29 the values of Q obtained after 2 cooling sweeps which lie between the minima that define the peak near Q = 1, as in Fig. 28, and which we therefore label by Q I = 1. We follow these same fields to 20 cooling sweeps and plot their values of Q in the figure. They nearly all have Q ≃ 1. In fact of the 665 lattice fields that we assigned to Q I = 1 at 2 cooling sweeps, 658 had Q ≃ 1 after 20 cooling sweeps with 4 at Q ≃ 0 and 3 at Q ≃ 2. Moreover, we observe at intermediate cooling sweeps that the 7 fields that end up with Q I = 1 migrate gradually from Q I = 1 to Q I = 1 through intermediate values of Q in a way that is consistent with a narrow (anti)instanton shrinking out of the lattice with increased cooling.
We have given the fields with Q I = 1 as an example. In Table 19 we show what happens to those fields possessing other values of Q I at 2 cooling sweeps when one continues to cool these fields to a total of 50 cooling sweeps. We indicate the range of Q that we assign to each integer Q I . That such a assignment should possess little ambiguity should be evident from the histogram of Q after 2 cooling sweeps in Fig. 28. We observe that the value of Q I is almost completely unchanged between 2 and 50 cooling sweeps in all cases.
One might ask why we start with fields after 2 cooling sweeps rather than after just 1 cooling sweep. The answer is indicated in Fig. 29 where we take the fields that we assign to Q I = 1 after 2 cooling sweeps and we display the values of Q for those same fields after only one cooling sweep. While we do indeed observe a distribution of values that has a mean not too far from what one might expect, the width of the distribution is so large that it will clearly overlap heavily with fields that belong to neighbouring values of Q I (and which are not shown). So after one cooling sweep we cannot make a plausible estimate of Q I unless we run more highly cooled fields backwards as we have done here. This is even more so if we consider the fields prior to any cooling at all.
What this study demonstrates is that calculating Q after 20 cooling sweeps provides a reliable method for calculating the true integer-valued topological charge of our Monte Carlo generated lattice fields -at least once β is not too small.
The shift of the value of Q as calculated after 20 cooling sweeps from the integer value that it would have in the continuum limit is primarily due to the finite size in lattice units of the typical instanton in the cooled fields. That is to say, it is a lattice artifact. So if we assign the corresponding integer Q I to a given field -an assignment which we have seen is essentially unambiguous over our range of β -we can think of it as partially compensating for the lattice artifacts. Of course since different lattice artifacts may mutually cancel, there is no guarantee that this will hasten the overall approach to the continuum limit. Nonetheless this integer-valued topological charge Q I has at least as much validity as the measured charge Q and we shall present results for both.
The possibility of assigning an integer Q I to a lattice field depends on there being a clear multi-peak structure in the histogram of values of Q with a very small overlap between neighboring peaks. As we see from Fig. 26 and Fig. 27 this is the case not only after 20 cooling sweeps but even after only 4 cooling sweeps for these values of β. So we also calculate an integer Q I for all cooling sweeps where a reliable assignment is possible. (Variations in how one does this will make a negligible difference as long as one more-or-less assigns fields between neighbouring minima in the histogram of Q to the same value of Q I .) Examples of such an assignment for β = 6.0625 and β = 6.50 are shown in Table 20 where we compare the values of Q 2 and Q 2 I against the number of cooling sweeps. As we see, a property of Q 2 I is that it is essentially independent of the number of cooling sweeps, unlike the value of Q 2 . (Note that the statistical errors shown are almost completely correlated when cooling, which is why we display them separately at the bottom of the Table.) This is a consequence of the fact that, as we have seen above, to a very good approximation the topological charge does not change with cooling. Beyond a very few cooling sweeps the topological charge only changes through an instanton shrinking out of the lattice (or the unlikely reverse), a process which becomes rapidly much rarer as β increases for any given number of cooling sweeps. Note that the fact that Q 2 I does not vary with the number of cooling sweeps means that we can evaluate it at some fixed number of cooling sweeps for any of our values of β. We shall choose to use 20 cooling sweeps unless otherwise stated.
As an aside, it is interesting to see how the average value of the lattice topological charge Q hot on the lattice fields prior to cooling, is related to the true (cooled) Q I . We show the result of the comparison in Fig. 30. We see that the relation is a simple multiplicative constant to a good approximation, i.e. Q hot = Z(β)Q I . Here we have indicated that we expect this constant to vary with the coupling β. Such a linear dependence is expected if one performs a perturbative calculation of Z(β) which at 1 loop leads to Z 1−loop (β) = 1 − 5.451/β in SU(3) [25]. Note that the correction is O(1/β) in our range of β which means that we expect Q hot ≪ Q I and this is certainly what we observe in Fig. 30. We have repeated this analysis for our other values of β and we list the values obtained for Z(β) in Table 21, together with the 1 loop perturbative predictions. We note that the latter are, given the large size of the leading correction, surprisingly close to our values.
As we remarked earlier, the fluctuations of the topological charge in the pure gauge theory are of particular interest since they are related, through the U A (1) anomaly to the masses of the lightest pseudoscalars, up to 1/N colour corrections [34,35]. In Table 22 we list our values for both Q 2 and Q 2 I at our values of β. The interesting quantities are these values divided by the lattice volumes, i.e. the topological susceptibility a 4 χ t = Q 2 /L 3 s L t , where L 3 s L t is the lattice volume, and the susceptibility expressed in physical units, e.g.
σ where a 2 σ is the string tension at that β. These values are also shown in Table 22. We extrapolate these values to the continuum limit with a leading O(a 2 σ) correction obtaining the continuum values given in the Table. The best fits are displayed in Fig. 31. Not surprisingly the fits deviate from our measurements at the two coarsest lattice spacings. At the largest β value we are beginning to lose ergodicity in topology (see below) so that our errors may be underestimates, but the overall fits are still statistically acceptable. In Table 22 we show the result of fitting to two different ranges of β. The results are consistent within errors as are the continuum limits using Q and Q I . If we take the most concervative value, based on the integer valued charge in the range 5.99 ≤ β ≤ 6.50, we obtain the continuum limit To convert the scale from √ σ to r 0 one can use the conversion factor given in Table 17 giving We note that this agrees, within errors, with the value r 0 χ 1 4 t = 0.4928(62) obtained in [36] where the calculation of Q is performed by the very different method of calculating the number of exact fermionic zero-modes using the Neuberger overlap-Dirac operator [37,38]. One can now use eqn (8) to convert to 'MeV' units. If we do so we obtain which agrees with the value χ 1 4 t = 208(6)MeV obtained in [39] using the gradient (Wilson) flow technique. (Note that this value has been rescaled from that in the paper using the more recent value of r 0 given in eqn (8). Note also that true systematic error in obtaining these physical 'MeV' units is no doubt underestimated.) Our above calculation assumes that finite volume corrections to the topological susceptibility are negligible when working on the lattice volumes shown in Table 22. To check that this is so we have calculated the susceptibility on different lattice sizes at β = 5.8941, β = 5.99 and β = 6.235. When expressed in units of the string tension the lattice volumes in Table 22 are spanned by the volumes 10 3 16 and 12 3 Table 23 and we see that there are no significant finite volume corrections within our small errors.

Critical slowing down
Since our Monte Carlo exploration of the phase space of all fields proceeds through local steps, by changing one link matrix at a time, the value of Q only changes when the core of a topological charge shrinks down to the scale of a lattice hypercube where it becomes nothing more than a gauge singularity. (Or the reverse process.) Once the lattice spacing a is small enough, a topological charge whose core is a few lattice spacings across will be embedded in a background gauge field that is relatively smooth, and one can calculate its weight in the path integral using standard semiclassical methods. To be specific the weight of an instanton of size ρ is where we omit factors varying weakly with ρ. We note the scale-invariant integration measure, a factor to account for the fact that a ball of volume ρ 4 can be placed in 1/ρ 4 different ways in a unit volume; and a factor arising from the classical instanton action, S I = 8π 2 /g 2 (ρ), where g 2 (ρ) is the running coupling on the scale ρ. When we insert the asymptotically free form of the coupling, we obtain where ξ is a physical length scale of the theory, such as r 0 , or √ σ, or Λ MOM . This formula should be reliable as long as ρ ≪ ξ. On the lattice it will also break down once ρ ∼ a in a way that depends on the particular lattice action, but it should be accurate once ρ f ew × a. In order for Q to change on the lattice, an instanton will have to pass through sizes ρ ∼ f ew × a and as we see from eqn (17) such fields are suppressed by a factor ∼ (a/ξ) 6 . That is to say, our sequence of Monte Carlo generated lattice fields will be locked into given values of Q once we are at large enough β and we will have lost ergodicity in Q. Whether any such loss of ergodicity in Q would have a significant impact on the physics that we are interested in is not clear. (In QCD the situation is quite different due to the associated zero-modes of the Dirac operator.) Nonetheless it is worth monitoring this potential loss of ergodicity, as we shall now do. We focus on two quantities. The first is to calculate directly the number of Monte Carlo sweeps, ξ Q , that it takes to 'decorrelate' the value of Q, as defined by where is and is + ξ Q label the number of MC sweeps and is is averaged over. This is of course not a physical quantity: it depends on our particular implementation of the Monte Carlo algorithm. However by comparing the value of ξ Q (β) to the total number of sweeps at a given β we can judge how ergodic is our calculation. The values of ξ Q (β) that we obtain in our calculations are plotted in Fig. 32. (We restrict our calculation to β ≥ 5.99 since we only calculate Q every 25 sweeps and so our determination of ξ Q can only be reliable once ξ Q ≫ 25.) Given that our calculations are typically based on ∼ 10 6 sweeps it is clear that we do not have a serious breakdown of ergodicity at any β although it is also clear that at β = 6.50 we are at the borderline. It is interesting to note that the dependence of ξ Q (β) on a(β) is well given by ξ Q ∝ a −6 , as shown in Fig. 32.
The second quantity we calculate is the number of sweeps, τ |∆Q|=1 , that is takes to change Q by one unit. We list our calculated values in Table 24. Clearly the probability of producing or annihilating an instanton as described above will be proportional to the space-time volume, so it is useful to give its value for a standard physical volume which we choose to be V = (3/ √ σ) 4 , which is very close to the volumes we actually used for calculating Q. The resulting values, labelledτ |∆Q|=1 , are also listed in Table 32 and, not surprisingly,they are close to the values of τ |∆Q|=1 . Now, since these changes in Q are randomly ±1 it takes roughly Q 2 such changes to disorder a charge of Q. On our lattices we have Q 2 ∼ 3 so the number of sweeps to disorder the lattice charge is ∼ Q 2 τ |∆Q|=1 ∼ 3τ |∆Q|=1 . This comes to about ∼ 4 × 10 3 for β = 6.50, which is not very far from our other measure, displayed in Fig. 32. We note that since Q 2 ∝ V and since τ |∆Q|=1 ∝ 1/V , we expect this measure of decorrelation to be approximately volume independent. In any case we confirm that our calculations in this paper are reasonably ergodic in Q.

Possible improvements
The calculations in this paper cover a large enough range of a(β) for the continuum extrapolations to be under reasonable control. And our large basis of loop operators means that it is unlikely that we are missing any glueball states in the part of the mass spectrum that we have presented. In addition we have checked that even at our largest values of β we maintain reasonable ergodicity as far as the topology of the gauge fields is concerned. In this section we briefly reflect upon the main ways in which our calculations can be improved.
While our calculation of the lightest states is adequate, the calculation of heavier states is less so. Here using lattices with an asymmetric discretisation, a t ≪ a s can be very useful. However since the glueball masses are calculated as a t M while the physical scale, whether provided by σ or r 0 , also involves the value of a s , we need to know the value of the ratio ξ = a s /a t to form dimensionless ratios such as M 0 ++ / √ σ or M 0 ++ r 0 . To overcome any issues with the shift in the value of the asymmetry ξ = a s /a t from its inputted tree-level value, one can perform two calculations in parallel. The first is with a s = a t and is much like the one presented in this paper but one that can be done with a much smaller basis of loop operators. It would be designed to calculate the lightest 0 ++ , 2 ++ and 0 −+ states as well as the string tension. (And perhaps r 0 and the spectrum of string excitations.) This would give us accurate values of the continuum limit of ratios such as M 0 ++ / √ σ and M 2 ++ / √ σ. Our second calculation would be with some well chosen value of ξ = a s /a t > 1 with a large basis of loop operators and would be designed to provide accurate values of the heavy as well as the light glueball masses. This latter calculation would give the continuum limit of, for example, M G /M 0 ++ for the glueball G. Then one can, for example, multiply it by the value of M 0 ++ / √ σ obtained in the first calculation to obtain M G / √ σ.
The main reason to calculate the heavier states is to provide a target for theoretical and model calculations. However with only 5 rotational representations on the lattice, each will typically end up having a densely packed set of excited states making the challenge ambiguous. So it is (almost) essential to be able to identify the spin, J, towards which each lattice state tends as a → 0, and also its 2J + 1 (near) degenerate partners. Doing so by simply identifying nearly degenerate states partitioned appropriately amongst the 5 rotational representations of the lattice is not an option once the states become dense enough. This means incorporating something like the approach in [6] into the glueball calculations. Doing so could in addition have the benefit of identifying the lightest high spin glueballs and hence the leading glueball Regge trajectory and so test the old idea that this trajectory is the Pomeron.
Amongst the more massive states will be multiglueball scattering states and we need to know if any of our supposed single glueball states are in fact of this type. (This is closely related to the issue of the decay widths of heavier glueballs.) While one expects there to be some suppression of the overlap between our single trace operators and the multi-trace operators that are the natural wave-functionals for multiglueball states, this is something that needs to be examined explicitly. A way to do this is to include some appropriate multitrace operators in our basis. For example we can take our variationally selected lightest 0 ++ glueball operator, give it a non-zero momentum, then take a product of two such operators with equal and opposite momenta to give a zero-momentum state. Do so for various momenta and possibly different glueballs. States in the resulting spectrum that have a small overlap onto such multi-glueball states can be identified with more confidence as being single glueballs. For a first step in this direction see for example Section 7.2.3 of [6] where small overlaps are found indicating that the multi-glueball states are likely to be invisible to current glueball calculations. However this needs to be studied more systematically and should form part of an improved calculation.
The above discussion also extends to the question of identifying ditorelon states. For the lightest such states their large variation with the lattice volume makes them easy to detect. However when such states are heavier and are located in the dense part of the spectrum, identifying finite volume shifts becomes ambiguous and a direct study of overlaps, as for multiglueball states, becomes useful and should be one of the incorporated improvements.
While do not claim the above improvements to be exhaustive, they would represent a major qualitative improvement to this and earlier glueball spectrum calculations.

Discussion
In this paper we provide a new calculation of the glueball spectrum of SU(3) gauge theories which is designed to improve upon some limitations of earlier calculations. We have a very large basis of glueball operators which gives us confidence that we are not missing any states in the important low-lying part of the spectrum. We perform high statistics calculations out to lattice spacings as small as a(β = 6.50) ≃ 0.042f m, which both improves our control of our extrapolations to the continuum limit and gives us a finer resolution of the glueball correlators, which is particularly useful for the heavier glueballs. We perform a high statistics calculation to check that our spatial volume is large enough and identify the most important multi-torelon states so that we can exclude them from the spectrum. We perform a high statistics study of the topological fluctuations of the gauge fields, with particular care in establishing the reliability of our measure of topology, so as to establish the ergodicity of our Monte Carlo sequences throughout our calculations. Given our cubic lattice, our glueball states fall into representations of the octahedral subgroup of the full rotation group, but using the observed near-degeneracies of various states we are able to assign a continuum spin to more states than previous calculations. At the same time, such an extensive calculation highlights important areas where improvements are needed: a parallel calculation with a lattice asymmetry allowing a much finer discretisation of the correlators [23,5] would be useful for the massive states; to systematically establish the continuum spins of the glueball states one needs to incorporate methods such as those used in [6]; to be sure one that one can distinguish multiglueball scattering states in the spectrum one needs to include the relevant operators into the calculations; and the same applies in controlling the contribution of the finite volume multi-torelon states.
Our results for the glueball spectrum of the continuum theory are listed in Tables 14, 15  and 18. Table 14 lists the masses in units of the confining string tension, and labels them by their parity P , charge conjugation C, and the irreducible representation of the rotation group of our cubic lattice. Table 15 translates these masses into 'GeV' units as explained in Section 5.2. In Table 18 we list those glueballs whose continuum spins J we can identify.
If we assume that the phenomenologically relevant states are those with masses 3GeV then our calculations in Table 18 confirm what has been established for a long time: these are the 0 ++ , 2 ++ , 0 −+ , 2 −+ , 1 +− ground states and the first excited 0 ++ . Of course our calculations are much more accurate than those early calculations, and improve significantly upon the more recent ones as well, but at this stage the greatest phenomenological uncertainty lies in the translation of the masses into GeV units and that is something that a calculation in the pure gauge theory cannot deal with by itself; here we need glueball calculations in full QCD to be further pursued.
We have emphasised that an important role for the glueball spectrum lies in providing a benchmark challenge for theoretical and model approaches that attempt to understand the non-perturbative physics of the SU(3) gauge theory and hence of QCD, usually in some approximation. For this purpose our excited states in Table 18 are useful, as indeed is the more extensive set of states listed in Table 14 since a calculation in the continuum theory can always project the states into the representations of the octahedral subgroup of the full rotation group.