SU(N ) gauge theories in 2+1 dimensions: glueball spectra and k-string tensions

We calculate the low-lying glueball spectrum and various string tensions in SU(N ) lattice gauge theories in 2 + 1 dimensions, and extrapolate the results to the continuum limit. We do so for for the range N ∈ [2, 16] so as to control the N -dependence with a useful precision. We observe a number of striking near-degeneracies in the various JPC sectors of the glueball spectrum, in particular between C = + and C = − states. We calculate the string tensions of flux tubes in a number of representations, and provide evidence that the leading correction to the N -dependence of the k-string tensions is ∝ 1/N rather than ∝ 1/N2, and that the dominant binding of k fundamental flux tubes into a k-string is via pairwise interactions. We comment on the possible implications of our results for the dynamics of these gauge theories.


Introduction
In this paper on SU(N ) gauge theories in 2 + 1 dimensions we provide a calculation of the light glueball spectrum and of the confining string tension for fluxes in various representations. This work improves significantly upon earlier work [1,2,5]; in particular, we cover a wider range of N , from N = 2 to N = 16, and we calculate the continuum limit of more excited states, by using a much larger basis of operators. And we do all this with greater accuracy than before. Since one of the motivations of such calculations is to look for interesting regularities in the spectrum, one needs accuracy and many states. Another motivation is to provide a testbed for analytic methods or models, e.g. closed flux tube [6][7][8][9] or constituent gluon [10] models, and this makes similar demands.

JHEP02(2017)015
In section 2 we begin with a brief description of the lattice setup and outline how we calculate the masses of glueballs as well as the energies of flux tubes from which we extract string tensions. We also list the main systematic errors one needs to be aware of. We then move on, in section 3, to the results of our string tension calculations and, in section 4, to our results for the glueball spectrum. Our calculations point to a number of interesting regularities which we discuss in the concluding section of the paper.
Our results for the glueball spectrum in this paper are largely consistent with the earlier work in [1,2,5] but are much more precise and extensive. Some of our preliminary results appeared in [11] where it was shown how a comparison of SO(N ) and SU(N ) gauge theories can provide an explanation for the very weak N -dependence of a number of glueball mass ratios. While the spectrum in this paper supersedes that quoted in [11], none of the conclusions of that paper are affected by the changes.

Calculating energies on a lattice
Here we outline how our lattice calculations are performed. The methods we use are standard and so we shall be brief.

Lattice setup
Our lattice field variables are SU(N ) matrices, U l , residing on the links l of the periodic L 2 s L t lattice, with lattice spacing a. The Euclidean path integral is where DU is the Haar measure over all the lattice link matrices and we use the standard plaquette action, Here U p is the ordered product of link matrices around the plaquette p. We write β = 2N/ag 2 , where g 2 has dimensions of mass and this becomes the continuum coupling when a → 0. Monte Carlo calculations are performed using a standard Cabibbo-Marinari heat bath plus over-relaxation algorithm. (See e.g. [1] for more details.)

Calculating glueball masses
The quantum numbers of our particle ('glueball') states are as follows [1][2][3][4] (We ignore momentum since we always take p = 0). In continuum D = 2 + 1 the rotation group is Abelian and the spins are J = 0, ±1, ±2, . . .. (More exotic values -anyons -are possible in two spatial dimensions, but we shall assume that they do not occur in SU(N ) gauge theories.) Parity, P , does not commute with rotations: J P → −J. So particle ('glueball') states can be labelled by parity P = ± and |J| as well as charge conjugation C. (For notational simplicity we shall use the label J P C rather than the more correct |J| P C in this paper.) Now, for J = 0 the state |J will have an orthogonal partner | − J that is degenerate with it. So the states ∝ {|J ± | − J |}, which are degenerate and non-null, JHEP02(2017)015 have parities P = ± respectively. That is to say we have the well known parity doubling for states with J = 0 in 2 space dimensions For J = 0 the argument breaks down because the P = − state can clearly be null. We also recall [1,2]. that on a square lattice, where the full rotation group is broken down to the group of rotations by π/2, states with J P = 1 + , 1 − (and odd spins in general) will still be degenerate, while states with J P = 2 + , 2 − (and even non-zero spins in general) need not be -although any non-degeneracy is due to lattice spacing correction and so will usually be very small. Note also that the finite volume periodic boundary conditions may also break the 2 P =± degeneracy, even in the continuum limit, but not that of odd spins. The fact that on a square lattice we have exact rotational invariance only under rotations of π/2 implies that states that tend to continuum states with J = 0 or J = 4 or J = 8,. . . will all fall into the same representation of the lattice rotation group, and similarly for J = 2, 6, . . . and for J = 1, 3, . . .. For simplicity, we shall label states belonging to these representations as J = 0, J = 2 and J = 1. What are the real spins of these states has been investigated in [2][3][4] and we shall refer to the detailed implications of that work for our results later on. However the reader should be aware that these were exploratory calculations, performed with smaller operator bases than those used in the present paper, and so there is some uncertainty in matching those results to the spectrum we obtain here. For now we merely remark that in some channels the actual ground states are not those with the lowest spin J.
Ground state masses M are calculated from the asymptotic time dependence of correlators, i.e. φ † (t)φ(0) = n | vac|φ † |n | 2 e −Ent t→∞ ∝ e −M t , (2.4) where M is the mass of the lightest state with the quantum numbers of the operator φ. (If φ has the quantum numbers of the vacuum then we use vacuum-subtracted operator φ − φ , so as to remove the contribution of the vacuum in eq. (2.4).) The operator φ will be the product of SU(N ) link matrices around some contractible closed path, with the trace then taken. We will use zero momentum operators, p = 0, so that there is no momentum integral on the right side of eq. (2.4). To calculate the excited states E n in eq. (2.4), one starts with a number of operators, φ i , and calculates all their (cross-)correlators φ † i (t)φ j (0) . One uses this as the basis for a systematic variational calculation in e −Ht 0 where H is the Hamiltonian (corresponding to our lattice transfer matrix) and t 0 is some convenient distance. Typically we choose t 0 = a. The procedure produces a set of orthonormal operators ψ i that are our best variational estimates for the true eigenstates. We then improve the energy estimate by using eq. (2.4) with the correlators C i (t) = ψ † i (t)ψ i (0) . With these operators we can hope to have good overlaps onto the desired states, so that one can evaluate their masses at values of t where the signal has not yet disappeared into the statistical noise. To achieve this in practice, we use blocked and smeared operators (for details see e.g. [1,12]) and a large number of such operators -in our case O(100) for each set of quantum numbers. Some examples of the basic loops we employ are displayed in table 1.

JHEP02(2017)015
An important aside. The accuracy of such a calculation 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 a simple exponential. 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.
We extract the energy of the state by going to vale of t = an t that is large enough that C i (t = an t ) ∝ exp{−aE i n t } within our errors. For illustrative purposes it is convenient to define the effective energy .
(2.5) Ift = añ t is the lowest value of t for which E eff,i (n t ) is independent of t ≥t within our errors, then we can say that a single exponential dominates and E eff,i (t) provides an estimate of the energy E i of the state. (In practice we extract our E i from a fit over a range of t.) That is to say, we search for a 'plateau' in the values of E eff (n t ) against n t . As a concrete example we show in figure 1 the values of E eff,i (n t ) for our six lightest 0 ++ glueballs in a calculation in SU(6) on a 60 3 lattice at β = 206.84, which corresponds to our smallest value of a. As we can see, the overlaps are not far from ∼ 100% so that it is easy to identify a plateau in E eff even for the highly excited states. In figure 2 we show the ground states for various quantum numbers, as obtained in SU(8) at a comparably small value of a.
Here the overlaps of some higher excited states are not quite so good, but an E eff plateau can still be plausibly identified. All this to indicate that while the mass estimates given in this paper should be largely reliable, this is less certain for the most massive states.

Calculating flux tube energies
We calculate the string tension by calculating the energy of a flux tube of length l that closes on itself by winding once around a spatial torus of size l. We use exactly the same technique as for the glueballs, except that the operator φ is now the product of SU(N ) 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 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, t) = (an x , an y , an t ), and we sum over n y to produce an operator with zero transverse momentum, p ⊥ = p y = 0. (We only consider flux tubes with zero transverse momentum, since we do not expect to learn anything qualitatively new from p ⊥ = 0.) In addition to this simple straight-line operator we construct other winding operators with a variety of kinks and loops extending from the original straight line and we also use smeared and blocked operators just as for the glueballs. Using this large basis JHEP02(2017)015 of operators we apply a variational procedure as for the glueballs, and this allows us to calculate not only the ground state energy but also excitation energies of the flux tube. We can label the flux tubes states by their transverse parities and their longitudinal momenta, and by their longitudinal parities if the momentum is zero. We can do all this for flux tubes carrying flux in any representation R by taking the trace in eq. (2.6) in that representation, i.e. by using Tr R . All this is described in detail in for example [13,14]. Since we will focus in this paper on the string tension rather than on the full flux tube spectrum [14] we only need to calculate the ground state energy, E gs (l), of the flux tube. We extract the string tension σ from E gs (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 [14,15] for reasons that have now become well understood [16][17][18][19][20] (see also [21]). We calculate in this way not only the fundamental string tension, σ f , but also string tensions, σ R , for the flux in a number of other representations R as discussed in section 3.
The reliability of such calculations depends on how well we can calculate E gs (l). As R becomes larger, so does the corresponding ground state energy E gs,R (l) which makes the calculation less reliable. To give some indication of the uncertainties involved we show in figure 3 the effective energies, as defined in eq. (2.5), for the representations of interest in this paper. We do this for the case of SU(8) at our smallest lattice spacing. We also show the estimates of E gs,R (l) that we obtain from these effective energies. We see that the identification of the effective energy 'plateau' appears to be unambiguous in all cases.

Ambiguities and systematic errors
A calculation of the glueball mass spectrum is subject to a number of systematic errors that one needs to control. (For flux tubes this is also the case, but because we are only interested in their ground states here, they are less important.) We discuss in this section how one can control some of these systematic errors. Our approach is to then to try to choose lattice sizes and operator bases such that the systematic errors are within our statistical errors. The errors we then quote on our masses and energies are always just the statistical errors.

Missing states
Our variational bases are, of course, incomplete and as we calculate ever higher excited states we will, at some point, begin to miss states. As it happens we have two sets of operators that were independently produced and which differ substantially, and comparing the spectra obtained using these bases provides some check on whether we are in fact missing any states in the range of interest to us here. One set of operators was used for our N = 3, 6 calculations, while the other was used for N = 2, 4, 8, 12, 16. However we have also performed some extra calculations to compare the results obtained using the two bases. One comparison was in SU(4) at β = 51 on a 40 2 48 lattice. Another was in SU(6), comparing the spectrum on a 54 2 60 lattice at β = 206.0 with that from the second basis on JHEP02(2017)015 a 60 3 lattice at β = 206.84. This latter comparison has a much smaller lattice spacing and hence a more reliable estimate of the masses of the heavier states, although the differing lattice sizes and values of β make it less than ideal. (But since to leading order aM ∝ 1/β, the expected shift in the masses should be a negligible ≈ 0.4%.) We compare the spectra obtained with the two bases of operators in tables 2, 3. We include for each set of quantum numbers as many states as we will eventually include in our continuum extrapolations. We see that the spectra are broadly consistent within errors. (This provides us with some confidence that we are not missing any states in the range of energies being considered. (The consistency of the spectra obtained with the two bases is of course a prerequisite for the large N extrapolations we perform later in this paper.)

Multi-glueball and unstable states
The full glueball spectrum will contain multi-glueball states. Such states are represented by multi-trace operators and so as N increases the overlap of such a state onto our basis of single trace operators will vanish by standard large-N counting. However at modest values of N such states may be present. Here two J P C = 0 ++ glueballs with angular momentum J can produce |J| ±+ states (with P = − for J = 0). Similarly |J| ±− states can arise from a 0 ++ and 0 −− glueball with angular momentum J. Some of our more massive states certainly fall into this 'dangerous' mass range. However the overlap suppression tells us that such a state will manifest itself through an effective energy that is much larger at small t and then decreases, asymptoting to 2m 0 ++ or m 0 ++ + m 0 −− at large t. Since we order our states by the value of E eff (n t ) at t = a, it is unlikely that any such states would appear in our energy range, and indeed we see no obvious sign of any such states.
A closely related issue concerns genuine single glueball states that are heavy enough to decay into lighter glueballs. Again the decay width will vanish as N → ∞ and at moderate N we expect the state to resemble a narrow resonance. So the effective energy should have something like an approximate plateau at small t, but at large t it will drift down to the threshold value of its decay products. However because the error/signal ratio grows exponentially with t, this large-t behaviour is unlikely to be visible and the state will appear just like any stable state. Of course things may be different at small N . To see what happens at smaller N we turn to our SU(2) calculation, where the decay widths should be largest, and to our smallest lattice spacing, which corresponds to β = 30.0, performed on a 120 2 90 lattice, We show in figure 4 the effective energies of the lightest few 2 ++ states. The horizontal line corresponds to the threshold energy of the potential decay products, i.e. aE th = 2am 0 ++ ≃ 0.43. We see that while several of the states are heavier than this theshold energy, there is no sign, within our errors, of a drift in the effective mass towards this threshold energy as t increases. There is of course a decrease at small t, as the higher excited state contributions to the correlator die away, but that does not persist to larger t, as one would expect it to do if the decay width was substantial.
All this encourages us to conclude that any systemtic errors due to multiglueball states or the instability of the states we consider are unlikely to be substantial in the calculations of this paper.

Finite V corrections
The leading finite volume correction, at large volumes, to the mass of a glueball comes from the emission of the lightest glueball, of mass m, which then winds around the spatial torus of length l, before being reabsorbed by the glueball. So these contributions are typically ∝ exp(−ml) and they are suppressed at large N because the triple glueball coupling is g 2 GGG ∼ 1/N . Now for our calculations in this paper ml varies between ∼ 25 at small N and ∼ 12 at our largest values of N , so we can assume that these leading finite volume effects are completely insignificant.
In practice, as is well known, the important finite volume corrections are quite different and arise from the presence of finite volume states composed of, for example, a winding flux tube together with a conjugate winding flux tube. These states are described by doubletrace operators so their overlap onto our basis of single-trace glueball operators will vanish at large N but is known to be substantial at small values of N . Now, when both flux tubes are in their ground states, the energy of such a 'winding' state is E T ∼ 2σ f l, and this can contribute to both 0 ++ and 2 ++ . For J = 1 or P = − or C = −, one needs to have one or both of the flux tubes in a suitable excited state and this has a much larger energy. Our strategy to control these finite volume corrections is therefore the very simple one of making l large enough at small N that all the 0 ++ and 2 ++ glueball states that we consider have energies no greater than E T . We then rely on the large-N suppression of the overlaps to allow ourselves smaller values of l at larger N . (A decreasing overlap means higher effective masses at small t, pushing the state out of our range of masses, which is determined by the effective mass at t = a.) Our specific choices are listed in table 4. Of course one needs to ask whether the large-N suppression is actually sufficient to eliminate such extra finite volume states at our 'large' values of N . This can only be answered convincingly by explicit finite volume studies. So in table 5 we compare the spectrum obtained at β = 306.25 on a 44 2 48 lattice (our standard size) with the spectrum on a much larger 60 2 48 lattice. Any finite volume state should be apparent from a large upward shift in its mass as we go to the larger lattice. As we see in table 5 the spectra are consistent, with no sign of any finite volume states. The smallest volumes we use are in SU(16) so we perform another comparison in that case. We list in table 6 the spectrum on a 22 2 30 lattice (our standard volume) and that on a larger 26 2 30 lattice, both at β = 800. We see that these two spectra are entirely consistent at the level of 2 standard deviations. All this strongly suggests that these finite volume corrections are insignificant in the calculations presented in this paper.
Returning briefly to the N -dependence of the finite volume corrections, it is interesting to contrast what we find for SU (8) in table 5 with what one finds in SU(2) when one performs calculations on lattice volumes that are similar to the SU(8) ones (in units of the string tension). The results are displayed in table 7 and we observe very significant finite volume corrections. Which of these are due to mixing with the extra torelon states (the lightest of which will be ∼ 2aE f , up to a shift due to interaction) and which are caused by the emission and reabsorption of glueballs that wind around the torus, is something we cannot address here. (We have also performed a calculation on an 104 2 80 lattice to check that the 80 2 64 lattice is large enough for the states of interest in this paper.) Nonetheless this displays that the overall finite volume corrections decrease with increasing N , as expected.

Nearly degenerate states
One of the interesting features of our results is that there are a number of nearly degenerate states. When these have the same quantum numbers and therefore arise from a single variational calculation (as described above) the variational procedure can induce an extra splitting, which is driven by statistical fluctuations and is therefore on the order of the errors. Moreover this may also lead to the two states being ordered differently in different bins, leading possibly to biases in the estimated errors and eventually to unsatisfactory fits when performing continuum or large N extrapolations. Any ambiguities here are small as long as the errors are small, but some of the most interesting near-degeneracies occur in the J = 1 sector where all the states are massive and the statistical errors are substantial.

Fundamental string tension and N -dependence
The ground state of the fundamental flux tube of length l that winds once around our spatial torus is usually the lightest of all the states (on the volumes we typically use) and hence it is the state whose energy, E gs,f (l), we can calculate most accurately and most reliably. We extract the string tension σ f from the energy using the 'Nambu-Goto' formula in eq. (2.7).
We list in tables 8-14 the resulting values of a √ σ f for our various lattice calculations. To obtain values in the continuum limit we need to express √ σ f in units of a quantity with dimensions of mass, and an obvious choice is the coupling g 2 which in 2 + 1 dimensions has dimensions [m 1 ]. The lattice coupling scheme that we choose to use is the 'mean-field improved' coupling g 2 I defined by [22] Since different choices of lattice coupling differ at O(g 2 ), the leading correction to the continuum limit will be O(a) rather than O(a 2 ). We therefore extrapolate our string tensions to the continuum limit, at each value of N , using We have used the 't Hooft coupling, g 2 N , since that is what needs to keep fixed for a smooth N → ∞ limit. We provide some examples of such continuum extrapolations in figure 5. As is apparent from the figure, we can get good fits with just the leading O(a) correction. We list in table 19 the continuum limit for each N , together with the coefficients of the linear correction, and the goodness of fit as measured by the total χ 2 and the number of degrees of freedom, n dof , of the fit.
We can now extrapolate our results to N = ∞. (For N = 2 and N = 4 we use the values on lattices of medium size, where the flux tube energies are small enough that one expects to avoid the systematic error associated with the large energies one obtains on large JHEP02(2017)015 lattices.) We expect the leading correction to be O(1/N 2 ) [23] and indeed we find that we get a marginally acceptable fit to all our calculated values with just this correction, which is displayed in figure 6. In order to see how robust this result is to the inclusion of higher order corrections in 1/N 2 , or to dropping the lowest N data point, we perform the corresponding fits, giving In the rest of the paper we shall use the result obtained in eq. (3.4).

k-string tensions and N -dependence
We calculate string tensions in other representations in the same way as for the fundamental. A flux tube carrying a flux that transforms under the Z N centre like the product of k fundamental fluxes is called a k-string. For N ≥ 2k the lightest flux tube in each k-sector is stable against decay as one can infer for k = 2, 3, 4 from the energies of the various ground states listed in tables 15-18 and from the continuum string tension ratios listed in table 20. So there is no ambiguity in extracting the flux tube energy, and hence the string tension using eq. (2.7). The ground states of these flux tubes fall into the totally antisymmetric representation to an excellent approximation [13,15] and so we label them as 2A, 3A, 4A respectively. One might ask how this feature can survive screening by gluons, but it turns out that for long flux tubes the effects of screening are extremely weak [15]. The flux tubes in the other representations listed in table 20 are in principle not stable; for example σ 2S > 2σ f and σ 3M > σ 2A + σ f (albeit not by much in this case). Since the least ambiguous string tensions are those that are associated with stable flux tubes, we will focus here on the 2A, 3A, 4A string tensions. We obtain the continuum limit of √ σ k / √ σ f from our calculated values using the where we include only an O(a 2 ) correction term since that suffices to obtain an acceptable fit in all cases. In figure 7 we show the fits for k = 2, 3, 4 in the case of SU (8). We see that the a-dependence is so small as to be consistent with zero. Our continuum extrapolations are listed in table 20. We begin by fitting σ 2A /σ f for N ≥ 4. In addition to the values in table 20 we also add the constraint lim N →∞ σ 2A = 2σ f from the theoretical expectation that lim N →∞ σ k = kσ f . We attempt separate fits: one that is in powers of 1/N , and one in powers of 1/N 2 . These

JHEP02(2017)015
give the best fits Since n df = 3, the first fit is entirely acceptable, but the second is not. The second fit is not much improved if we restrict ourselves to N ≥ 6. It is clear that our calculations imply that the leading correction to σ 2A /σ f is ∝ 1/N rather than ∝ 1/N 2 . We repeat the exercise for σ 3A /σ f , this time adding the constraint lim N →∞ σ 3A = 3σ f . Fitting to N ≥ 6 we obtain Again this clearly points to the leading correction being ∝ 1/N and not ∝ 1/N 2 . Finally while the best fit with N → N 2 has an unacceptable χ 2 /n df = 12.1.
In the above we have restricted our fits to N ≥ 2k since the center symmetry allows a k-string to mix with a (N − k)-string, and the latter will have a lower string tension if N < 2k, and so it will then provide the lightest flux tube in the k-sector. So, for example, if we extrapolate eq. (3.7) to SU(3), we expect to find σ 2A su3 = σ f . Interestingly enough, substituting N = 3 in eq. (3.7) does indeed give a value very close to unity. We can extend the argument to SU(2), where k = 2 can mix with the k = 0 vacuum, and indeed the O(1/N ) fit in eq. (3.7) gives us σ 2A ∼ σ k=2 (where we estimate the value of σ k=2 in SU(5) using eq. (3.7)). Of course this rough agreement is only indicative: as we go to smaller N we must expect that the omitted higher order terms in 1/N will become significant. Indeed a constructive way to look at this is to use the constraints from N < 2k to fix, with no extra work, the next higher order terms in the expansion in powers of N , which our fitting for N ≥ 2k is not sensitive to.
Another interesting feature becomes apparent when we look at the coefficient c 1 (k) of the leading 1/N correction in the above equations. We observe that i.e. c 1 (k) ∝ k(k − 1)/2, which is just the number of pairwise interactions amongst the k fundamental strings which are bound into the k-string. A similar rough proportionality appears to hold when we look at the coefficient c 2 (k) of the second, 1/N 2 , correction term.

JHEP02(2017)015
This appears to suggest that the interactions binding the k fundamental strings into the k-string are predominantly pairwise. It is also interesting to compare our results to the Karabali-Kim-Nair (KKN) analysis [24,25] which provides a prediction not only for the k and N dependence of σ k , but also predicts its absolute magnitude: Although we know from earlier studies [26] that eq. (3.13) does not fit the calculated values of σ f /g 2 N perfectly, it does come within ∼ 2% for N ≥ 3, and within 3% even if we include SU (2). The prefactor in eq. (3.13) is simply the quadratic Casimir of the represention which, for reasons obvious from the above discussion, we take to be the totally antisymmetric one.
From the values of √ σ f /g 2 N in table 19 and table 20 we form the ratio √ σ k /g 2 N which we plot in figure 8. We observe very close agreement between the string tensions predicted by eq. (3.13) and our calculated values. Indeed the agreement is clearly better for σ k≥2 than for σ f . Of course, since lim N →∞ σ k = kσ f , any disagreement will be the same for σ k and for σ f in the large N limit. Nonetheless it is clear that the simple formula in eq. (3.13) provides a remarkably good approximation to all our k-string tensions for all values of N .

Glueball spectra
We calculate glueball masses as described in section 2.2. A particular concern is the finite volume states which might appear in our spectrum of excited glueball states, as discussed in section 2.4. The lightest such state is composed of a ground state flux loop and its conjugate, so it will have a mass m T ∼ 2σ f l, and where it might appear is in the 0 ++ and 2 ++ spectra. Since our glueball operators are single trace while these finite volume states are double trace, large-N counting tells us that such states will decouple from our correlators at larger N . So our strategy is to have a large enough volume l 2 at small N so that m T is heavier than the states we calculate, and then to relax the volume to a smaller value at larger N where we expect such states to become invisible. A summary of the 'standard' volumes we have chosen to use, expressed in the relevant physical units, is presented in table 4. We check these choices with some explicit comparisons between smaller and larger volumes. For example in SU(8) at β = 306.25 we compare the spectrum on our 'standard' 44 2 48 volume to the one on a larger 60 2 48 volume. The resulting glueball masses are given in table 5, where we include all the excitations for which we will attempt to obtain continuum limits. We see that the spectrum on the l = 44a volume matches that on the l = 60a volume within say 2σ, suggesting that for this moderately large value of N the choice of l √ σ f ∼ 3.8 is sufficient to exclude finite volume corrections (of whatever source) at the level of our statistical errors. In table 6 we perform a similar comparison in SU(16) at β = 800 where our 'standard' volume is l = 22a, corresponding to to l √ σ f ∼ 3.1.
Again we see that the glueball spectra agree, within say 2σ, indicating that our reduction in l √ σ f with increasing N (based on the idea of the large-N suppression of finite volume

JHEP02(2017)015
corrections) is indeed appropriate. The specific lattices and β-values used are listed in tables 8-14, along with the values of the string tension and the mass gap. For N ∈ [2,6] there are two sets of volumes listed: the larger are used for glueball calculations while the smaller are used for string tension calculations. The reason for doing this is that the energy of the flux loop grows with l, so we can perform a (statistically) more precise calculation with smaller l. This is only possible, of course, because we have a very precise theoretical control of finite volume corrections in this case. Having obtained our 'infinite' volume glueball masses in this way, for various lattice spacings, we extrapolate the results to the continuum limit. Mostly we do so for the dimensionless ratio am/a √ σ f = m/ √ σ f since in general a √ σ f is our most accurately calculated physical quantity on the lattice. We perform the continuum extrapolation in the standard way and we find in nearly all cases that the first O(a 2 ) correction suffices for a good fit. The results of these continuum extrapolations are listed in tables 21-27, together with the χ 2 and number of degrees of freedom n dof of the fit. In table 25 we also show the values one finds for the coefficient c 1 in eq. (4.1); they are very similar for other values of N . We see that the lattice corrections are very modest even at our coarsest lattice spacing, where a 2 σ f ∼ 0.1. Note that at the coarser values of a, the masses of higher excited states can be too large for them to be identified, and then the value of n dof is smaller than for the ground states.
In addition to these continuum fits we also calculate in one or two cases the continuum limit of the ratio m/g 2 N just as we did for the string tension in eq. (3.2). This might seem an attractive way to calculate continuum physics because the error on the denominator g 2 I N comes from the average plaquette and is negligible. In practice, however, this advantage is more than outweighed by the fact that the leading correction is O(a) rather than O(a 2 ), so we generally focus on the extrapolations using eq. (4.1). We will also look at some of the masses in units of the mass gap, which is of interest for the reasons discussed in [11].
Finally we extrapolate our continuum mass ratios to N = ∞ using  state and the lowest few 0 −− masses can also be calculated accurately. The lightest 0 −+ and 0 +− are quite heavy and the mass estimates for them are correspondingly less reliable.
In figure 9 we plot the lattice values of m/ √ σ f , calculated in SU (6), against the value of a 2 σ f for the lightest six 0 ++ states, the lightest three 0 −− states, and for the 0 −+ ground state. We also show the continuum extrapolations, and note that in all cases a leading O(a 2 ) correction suffices to give an acceptable fit. We observe some striking regularities in figure 9. Firstly, each of the three 0 −− states is nearly (but not exactly) degenerate with an excited 0 ++ state: to be specific, the first, second and fifth excited 0 ++ states. Secondly, the ground state 0 −+ is degenerate, within errors, with the fourth excited 0 ++ . In fact earlier analyses [2] have revealed that the lightest '0 −+ ' state is in fact a 4 −+ state so, by parity doubling, there should be a degenerate 4 ++ state which will appear in what we label as the set of '0 ++ ' states. In other words, we believe that what we have called the fourth 0 ++ excitation is in fact a 4 ++ state. Moreover this 4 ++ state is nearly degenerate with the fifth 0 ++ excitation (and also with the second excited 0 −− ). We also note that the ground state 0 +− appears to be nearly degenerate with the first excited 0 −+ , and that in those cases where we have calculations (i.e. for N = 4, 12, 16) the first excited 0 +− appears to be nearly degenerate with the second excited 0 −+ . (With the caveat that for these very massive states the errors are large.) As we can see in figure 9 the gaps between the states that are 'nearly degenerate' are much smaller than the typical gaps between other states. That is to say, these regularities appear to be real rather than statistical coincidences.
While we have chosen to use SU (6)  As an example we plot in figure 11 our values in SU (12) of the lightest few 2 ++ and 2 −− glueball masses as a function of a 2 σ f together with the corresponding linear continuum extrapolations. Just as for J = 0 we observe a striking pattern of near-degeneracies: the first, second and fourth 2 ++ excited states appear to be nearly degenerate with the lightest three 2 −− states. We also see from tables 21-27 that this near-degeneracy holds for all N except in SU(3) and SU(8) for those excited states where the parity doubling is poor, and, in any case, in those cases the near-degeneracy holds between the 2 −+ and 2 ±− .

JHEP02(2017)015
The large-N extrapolations of our J = 2 spectra are listed in table 29. In the 2 ++ sector the second and fourth excited states have such poor fits that we do not attempt to provide any estimate based on a large-N extrapolation. Instead we provide an estimate based on the average of the SU(12) and SU (16) values (with errors enhanced to encompass the two values.) Since the fits to the supposedly degenerate 2 −+ parity partners are mostly acceptable, it is these that we plot in figure 12 against 1/N 2 , together with the lightest 2 −− states. We clearly see the near-degeneracies discussed above. We also see that the large-N limits of the second and third excited 2 −+ states are very similar. All this is very similar to what we observed in the J = 0 sector, displayed in figure 10.
One might ask if there is a J = 6 state located in the J = 2 sector, just like the J = 4 state in the J = 0 sector. This has been addressed in [3] for SU (2) where one finds (see tables 2 and 3 therein) that m J=6 gs ∼ 1.5m J=2 gs . This suggests that the lightest five J = 2 states which we consider in this paper (see our table 21) are indeed J = 2 and not J = 6 because their masses are within ∼ 1.4 of the J = 2 ground state mass. For larger N we do not have detailed guidance, but we will assume that the lightest five J = 2 states are also really J = 2, just as in SU (2), even though they are now up to ∼ 1.55 times heavier than the ground state. Needless to say this identification should not be taken as being definitive.

|J | = 1, 3, . . . glueballs
For J = 1 (and indeed any odd J) we expect to have exact parity doubling, not only in the continuum limit but even on a square lattice in a finite volume as long as the latter respects π/2 rotational invariance. Our calculated masses are consistent with parity doubling, as we see for example in tables 2, 3 and tables 5, 6.
We list in tables 21-27 the values that we obtain for the J = 1 glueball masses expressed in units of the string tension, after an extrapolation to the continuum limit. For simplicity we refer to the |J| = odd states as J = 1 both here and in the tables, but will shortly return to the question of their true spin.) A striking feature of the J = 1 spectrum is that for N ≥ 3 the lightest two 1 ++ states are nearly degenerate, within our errors, and that the next two states are also nearly degenerate. (By parity doubling the same is true for the 1 −+ states.) The gap between these two pairs of states is much larger than any gap within each pair, making it plausible that the near-degeneracy is not just a statistical accident but has a dynamical origin. Moreover the fact that all these states are very massive, with correspondingly large statistical errors, creates some technical problems for the variational procedure, where the ordering of the states may differ in different (jack-knife) data bins so that the observed small splittings between the states in each pair may be enhanced, or even entirely driven, by the statistical fluctuations. This may be the reason for the fact that many of the continuum fits are quite poor, and that we occasionally see discrepancies between the P = ± J = 1 masses.
Another feature of the J = 1 spectrum, that differs from the J = 0, 2 spectra, is that the lightest state is the C = − ground state. Moreover the first excited C = − state appears to be roughly degenerate with the pair of nearly degenerate C = + ground states and, albeit now within larger errors, the third excited C = − state is nearly degenerate

JHEP02(2017)015
with the next pair of nearly degenerate C = + states. All this is very reminiscent of what we saw for J = 0, 2, but with C = − and C = + interchanged.
To illustrate these features of the J = 1 ±+ and J = 1 ±− spectra, we plot in figure 13 the masses of the J = 1 states against a 2 σ f for our SU(8) calculation. We have averaged the P = ± masses since they should be degenerate, and this should help to suppress fluctuations. The fact that the lightest state is C = − is unambiguous. We also see good evidence that the lightest two C = + states are (nearly) degenerate, and that the first excited C = − state is (nearly) degenerate with them. The next two C = + states have similar masses to each other and may be degenerate although the errors on these massive states are so large that this is something of a conjecture. Equally, the next two C = − states have similar masses to each other and also, interestingly, to the second pair of C = + states.
These features are characteristic of all our SU(N ≥ 3) spectra and hence also of the N → ∞ extrapolations listed in table 30. These extrapolations have reasonably small errors, and the (very near) degeneracy of the lightest two C = + states is convincing as is the approximate degeneracy of the next two excited C = + states. It is also clear that the C = − ground state is the lightest J = 1 state and that it is not degenerate with any other state. However the first excited C = − state appears to be (nearly) degenerate with the lightest pair of C = + states. The second C = − excitation does not appear to be very close to any other states, but the third C = − excitation appears to be nearly degenerate with the second pair of C = + states. One can, as for SU (8), reduce the fluctuations a little by averaging the P = ± values. We show the result of doing so in figure 14 where we plot these averaged values for all our SU(N ) continuum limits against 1/N 2 . Apart from SU (2) where there are no C = − states, and where the lightest four C = + states definitely do not form nearly degenerate pairs, all the other SU(N ) spectra display the features emphasised above.
Given these striking regularities, it is interesting to ask whether there is any evidence that some of these states might be |J| = 3, 5, . . . rather than J = 1. Now the analysis for SU(2) summarised in table 2 of [3] and the extension to SU(3) and SU(5) in section 6.4 of [4] does in fact provide evidence that the lightest two states in the '1 ++ ' sector are J = 3 and not J = 1, while in the '1 −− ' sector the ground state is J = 1 but the first excited state is J = 3. If we take these assignments as correct then this will alter our conclusions as follows: the lightest J = 1 P =± state has C = − and is lighter than the lightest J = 3 state and is much lighter than the lightest C = + J = 1 state. The lightest J = 3 P =± states consist of three nearly degenerate states: two with C = + and one with C = −. For the higher excited states we have no evidence concerning their spin, other than that it is odd.

Discussion
A striking feature of the fundamental string tension expressed in units of the coupling, is how small are the lattice corrections, as we can see in figure 5. Even at the coarsest lattice spacings, where a √ σ f ∼ 0.3, the correction is less than 10%. Indeed a simple O(ag 2 I N ) correction term is all that is needed despite the precision of our lattice calculations of a 2 σ f and the substantial range of a being fitted. Moreover the same is true of the mass gap, as we see in table 19. The lattice corrections are even smaller when we consider the masses JHEP02(2017)015 of the lightest glueballs expressed in units of the string tension, m G / √ σ f , as we see in figures 9, 11 and 13. The same is true for the tensions of the stable higher representation flux tubes, as we see in figure 7 in the case of SU (8). This remarkably precocious scaling may have to do with the fact that the theory is super-renormalisable: the dimensionless running coupling decreases linearly as the distance scale is reduced,g 2 (l) = lg 2 , so that corrections that are higher powers of g 2 are necessarily accompanied by higher powers in a. (In contrast to 4 dimensions where the leading O(a 2 ) lattice correction is a power series in g 2 .) A striking feature of our (continuum extrapolated) values of m G / √ σ f , m G /g 2 N and √ σ f /g 2 N , is how weak is their variation with N , over the whole range of N ≥ 2, as we see in figures 10, 12 and 14. A possible explanation for this has been given in our earlier paper [11] where we discuss the constraints on the N -dependence that arise when we consider SO(N ) as well as SU(N ) gauge theories.
In any case all this suggests that there may be some underlying simplicity in the dynamics of D = 2 + 1 SU(N ) gauge theories, and this motivates a close examination of our results for unexpected regularities. We have indeed found evidence for a number of these (some already remarked upon in earlier less precise calculations) which we shall now summarise, beginning with the string tensions.
In addition to the stable fundamental flux tube, it is known that the ground states of k-strings (i.e. flux tubes which carry the flux of local static sources consisting of k fundamental charges) are also stable, and that these flux tubes are in the totally anti-symmetric representation (when not very short). We confirm all this with our more accurate calculations. Furthermore, our large range of N allows us to make an unambiguous statement that the leading corrections to the N = ∞ limit of σ k /σ f are O(1/N ) rather than O(1/N 2 ). Moreover the O(1/N ) binding energy is consistent with being produced by a pairwise interaction between the k fundamental strings that are bound into a k-string. And we confirm previous observations that to a very good approximation the string tensions are proportional to the quadratic Casimir of the representation. It is intriguing that the absolute value of √ σ k /g 2 N is very close to that predicted in [24,25] as we see in figure 8.
Turning now to our results for glueballs, here the striking feature is a number of unexpected near-degeneracies amongst glueballs with different J P C quantum numbers. (This is in addition to the expected parity doubling for J = 0 which we also observe.) For J = 0 we saw that the first, second and fifth excited 0 ++ glueballs are nearly degenerate with the lightest three 0 −− glueball states. The fourth excited '0 ++ ' state is presumably the 4 ++ ground state given its near degeneracy with the '0 −+ ' ground state which is believed to be, in reality, the 4 −+ ground state [2]. The only 0 ++ states that are not nearly degenerate with another J = 0 state (amongst the lightest 6 states) are the ground state and the third excited state. We note that the mass of the latter is very nearly twice the mass of the former. (Although this may well be accidental.) In the J = 2 glueball sector we observe a nearly identical pattern of degeneracies. The first, second and fourth excited 2 ++ glueballs are nearly degenerate with the lightest three 2 −− glueball states. Only the ground and third excited of these lightest 2 ++ states are not nearly degenerate with some other J = 2 state.

JHEP02(2017)015
In the J = 1 sector the lightest glueball is a 1 −− . There are two very nearly degenerate '1 ++ ' ground states which are nearly degenerate with the first excited '1 −− ' state. But there is in fact good evidence [2,4] that these three nearly degenerate states are J = 3 rather than J = 1. It also appears that the next two pairs of excited states in the 1 ++ and 1 −− sectors are nearly degenerate and nearly degenerate with each other, although their large masses mean that this observation is more speculative.
We finish by recalling that in the simplest closed flux tube model of glueballs [6][7][8][9], the C = + and C = − glueball states in two spatial dimensions are degenerate. We also note that while we have a number of near-degeneracies, they do not appear to be exact. That is to say, even if they point to some kind of relatively simple dynamics, this is likely to be the property of a field theory from which the SU(N ) gauge theory is -at the very least -a small perturbation. Helping to identify such a neighbouring field theory which one may hope to be analytically tractable, is a motivation for the present study.

Acknowledgments
AA has been partially supported by an internal program of the University of Cyprus under the name of BARYONS. In addition, AA acknowledges the hospitality of the Cyprus Institute where part of this work was carried out. MT acknowledges partial support under STFC grant ST/L000474/1. The numerical computations were carried out on the computing cluster in Oxford Theoretical Physics.