Glueball Spins in $ D=3$ Yang-Mills

We determine spins of more than 100 low lying glueball states in $D=2+1$ dimensional $SU(4)$ gluodynamics by a lattice calculation. We go up to $J=8$ in the spin value. We compare the resulting spectrum with predictions of the Axionic String Ansatz (ASA). We find a perfect match for 39 lightest states, corresponding to the first four string levels. In particular, this resolves tensions between the ASA predictions and earlier spin determinations. The observed spins of heavier glueballs are also in a good agreement with the ASA. We did not identify any sharp tension between lattice data and the ASA, but more work is needed to fully test the ASA predictions for the spins of 64 states at the fifth string level.


Introduction
SU (N c ) Yang-Mills (YM) theory turns into a free string theory in the 't Hooft planar limit N c → ∞ [1]. Understanding these confining strings remains an outstanding challenge. Lattice Monte-Carlo simulations deliver precious experimental guidance for this enterprise. Lattice observables sensitive to the stringy origin of YM glueballs belong to two large distinct categories. The most direct way to probe the worldsheet dynamics of confining strings on the lattice is to measure the spectrum of tolerons-long flux tubes wound around one of the spatial directions [2][3][4][5][6]. Applying the effective string theory methods, such as large radius expansion [7][8][9][10][11][12], Thermodynamic Bethe Ansatz (TBA) [13,14] or flux tube S-matrix bootstrap [15] one can then extract properties of the worldsheet action and S-matrix from the spectral data. These studies resulted in the Axionic String Ansatz (ASA) for the worldsheet theory in both D = 3 and D = 4 gluodynamics [16,17].
The second set of lattice data comes from glueball spectroscopy [18][19][20][21][22]. Glueball spectra have an advantage of representing traditional physical observables, directly measurable in the lab. However, relating these data to the worldsheet properties is hard. In particular, effective string theory is only useful in the vicinity of the leading Regge trajectory [23]. Hopefully, one can do better with the ASA, which postulates that the worldsheet theory can be treated as a deformation of an integrable model (see [24,25] for a further discussion of the underlying physics). To predict the glueball spectrum one needs then to solve for the short string spectrum in the integrable approximation, and to develop perturbative methods accounting for violations of integrability.
A very first step towards implementation of this program has been made in [17] for D = 3 YM. The ASA string is particularly simple at D = 3, where it belongs to the same equivalence class as the minimal Nambu-Goto string with no additional massless or massive worldsheet degrees of freedom present. In the integrable approximation the ASA worldsheet S-matrix is the same as the worldsheet S-matrix of critical strings [26]. Building on this similarity, [17] came up with a prediction of glueball quantum numbers. At D = 3 these are glueball spin J, spatial parity P and charge conjugation parity C. J > 0 glueball states come out as doublets of states carrying opposite spatial parities P , so that P carries non-trivial information only for J = 0 glueballs. We will refer to the resulting spectrum as the ASA spectrum, but it is important to keep in mind that some additional assumptions have been made in [17], which may be not a part of the ASA.
In the present paper we report results of a dedicated lattice simulation whose goal is to test the predictions of [17]. The major difficulty for such a measurement is that the physical SO(2) rotation symmetry is broken to its Z 4 subgroup on a cubic lattice. As a result, the conventional glueball mass measurements, such as the ones presented in [22], deliver only Z 4 transformation properties of the glueball states, rather than the physical glueball spins.
A technique to measure the actual physical spins J was described in [18][19][20], where it was applied to several low lying glueball states. Here we build on this proposal and perform continuous spin measurements for all accessible glueball states (in practice, this is more than 100 states). Let us stress that our main goal is to provide spin determinations rather than to perform a high precision mass measurement. As a result, the mass values which came out from our simulation are less precise than those reported in [22]. Also we restrict ourselves to the SU (4) gauge group. This value of N c should be large enough for our purposes, given that the typical mass shifts between N c = 4 and N c = ∞ measured in [22] are all within 10%.
The rest of the paper is organized as follows. In Section 2 we present a short summary of the ASA predictions for the glueball spectrum. We also compare them to the previously existing lattice results. In Section 3 we describe in detail our methods. The results are presented in Section 4, and also in Tables 2a- 18 and Figures 2-11 which can be found at the end of the paper. We compare the measured glueball spectra with the ASA predictions and find good agreement. Future directions are discussed in the concluding Section 5.

ASA spectrum and earlier lattice data
Let us review the ASA predictions for the spectrum of D = 3 glueballs. It is instructive to start with a brief reminder of what the spectrum of closed critical bosonic strings looks like (see, e.g., [27]). The Hilbert space H c of closed critical strings is a direct sum of subspaces labeled by a non-negative integer N , For critical strings the masses of states at the same level N are all equal to each other. This mass degeneracy follows from exact integrability of the worldsheet theory. Furthermore, each subspace H N is a tensor product of spaces describing left-and right-moving worldsheet perturbations, The level number N counts the total left-moving momentum along the string. The right-moving momentum necessarily takes the same value as a consequence of the worldsheet reparametrization invariance. This is known as a level matching condition. In order to define charge conjugation one needs H L (N ) = H R (N ), then A geometric meaning of (3) is that charge conjugation acts by flipping the string orientation. Note that the tensor square structure (2) follows from the absence of massive perturbations on the critical string worldsheet. If present, these would allow for zero momentum excitations, modifying the tensor square structure.
Finally, for critical strings one finds that the spectrum of closed strings is related to the spectrum of open strings via where H o (N ) is the Hilbert space of open string states at level N . According to the ASA confining strings at D = 3 do not carry any massive excitations. Unlike in the critical case, the worldsheet theory is no longer integrable [14,28], so that the exact mass degeneracies at each level do not hold. However, as a consequence of approximate integrability, one still expects to find approximate level degeneracies. This leads to the first (weak) ASA prediction-glueballs are expected to come in approximately degenerate groups (levels), where each group can be decomposed as a tensor square (2). As discussed in [17], this "weak" prediction is already quite restrictive. In particular, it implies that multiplicities of each level are equal to exact squares. Also, odd spin glueballs at each level appear in pairs of opposite charge parity and equal spin.
A stronger set of predictions follow from an assumption that the relation (4) continues to hold for confining strings. The structure of the open string Hilbert space H o can then be described using the following semiclassical considerations. One starts with decomposing H o as a sum of subspaces H j of a definite angular momentum j, Note that the sum in (5) runs over all integer values of j. Open string states of spin J > 0 are doublets of states with j = ±J. At large J the minimum energy state in the corresponding H j is described by a classical rotating rod solution. These states give rise to the leading Regge trajectory and satisfy N = J. Low lying excitations in H j can be constructed by quantizing small perturbations around the rotating rod. This results in a prediction for quantum numbers and multiplicities of states close to the leading Regge trajectory. The strongest set of the ASA predictions follow from an additional assumption that this perturbative analysis applies at all j (including j = 0 !) and for arbitrarily high excited states in each H j . We will see that this somewhat ad hoc assumption works surprisingly well. Its possible justification may be that it is often the case that semiclassics produces exact results in integrable models. Confining strings are not integrable, however, according to the ASA they can be described as a deformation of an integrable theory. Discrete data, such as state's multiplicities and quantum numbers, cannot change under a continuous deformation.
Adopting this assumption allows one to fully predict the multiplicities of states in each H j . Namely, the multiplicity of states at the level N in H j is given by the N -th coefficient in the Taylor expansion of the following generating function where The open string Hilbert space does not have direct physical meaning for pure gluodynamics 1 . To obtain glueball quantum numbers one makes use of the tensor square decomposition (2). The resulting predictions of the glueball quantum numbers for the first five levels (i.e., for N = 0, . . . , 4) are presented in Table 1. We see that the "strong" ASA completely fixes quantum numbers of the first 39 glueball states (up to N = 3).
The ASA also fixes the quantum numbers of 64 glueball states at the N = 4 level, up to a spatial parity assignment for two (out of ten) J = 0 glueballs. The latter ambiguity arises because the arguments above do not fix parities of j = 0 open string states. This leaves parities of some J = 0 glueballs undetermined, if there is more than one j = 0 state at the corresponding open string level. N = 4 is the first level when this happens. Independently of how this ambiguity is resolved, one of these states will be an exotic CP = −1 state, i.e. either 0 +− or 0 −+ . This makes the N = 4 level especially interesting-this is the first place where such exotic state appears in the ASA spectrum.
To conclude this section and to set the stage for our further study let us review the status of the ASA spectrum as compared to the lattice data prior to the current work. A high precision mass determination for a large number of glueball states (and for up to N c = 16) has been recently performed in [22]. However, SO(2) spin values were not measured in this simulation. Still, some limited information about glueball spins has been available. First, a number of low lying Z 4 singlet states observed in [22] do not have nearly degenerate partners of opposite parity. This implies that these are true J = 0 states. In addition, the lightest J = 4 doublet has been identified in the earlier dedicated study presented in [18,20]. Furthermore, combining the results of [22] with the earlier ones presented in [19,20] it was argued in [22] that the preferred J C assignments for the four lightest odd spin glueballs are 1 − and 3 + , 3 + , 3 − .
These results are summarized in Figure 1. One immediately identifies in this figure wellpronounced N = 0, 1, 2 levels perfectly matching the ASA predictions. This agreement is nontrivial, but should not be considered as a definitive confirmation of the ASA logic, given that the results presented in Figure 1 were available prior to [17] and served as a guide for the analysis presented in [17]. Hence, the crucial test of the ASA should come from spin determinations for heavier glueballs with masses in the N = 3, 4 regions. The goal of the present paper is to perform such a test. In particular, if the preliminary spin assignments for the four lightest odd J doublets in Figure 1 were correct, they would sharply contradict even the weak ASA. Addressing this tension is one of the primary motivations for the current study.

Lattice simulations
We generate configurations of SU (4) matrices representing lattice link variables, corresponding to exponentiated gauge fields, in accordance with the standard Euclidean path integral where DU is the Haar measure over the lattice link matrices, a the lattice spacing and S is the Wilson plaquette action Here U p is ordered product of link matrices around a plaquette at location p.
We use a 80 × 80 × 96 lattice with β = 98.25; 18 million lattice configurations were generated using a standard Cabibbo-Marinari heat bath plus over-relaxation algortihm, with operators being evaluated every 100 configurations.

Spin ambiguity and parity degeneracy on a cubic lattice
As is well known, on a cubic D = 3 lattice spatial Euclidean SO(2) symmetry is broken to its Z 4 subgroup [20,29]. Using only rotations of loop operators restricted to the lattice we cannot distinguish the spins of eigenstates which have J = 0 mod 4, nor among those that have odds spins J = 1 mod 2, nor among those with J = 2 mod 4. As discussed in detail in Section 3.3 we get around this problem by implementing up to 16 approximate rotations of loop operators, so that we can distinguish states of specific spins up to J = 8.
In three spatial dimensions P parity can be defined as a reflection w.r.t. the origin. This transformation commutes with spatial rotations, so that the extended orthogonal group is simply a direct product O(3) = SO(3) × Z 2 , where the Z 2 factor corresponds to P . On the other hand, in two spatial dimensions a reflection w.r.t. the origin is equivalent to a conventional rotation. P parity has then to be defined as a reflection w.r.t. one of the axis, for example the y-axis. As a result, P does not commute with continuous rotations and the extended orthogonal group is a semidirect (rather than a direct) product O(2) = SO(2) Z 2 . In other words, the parity and angular momentum operators do not commute. A parity transformation switches clockwise and counter-clockwise rotations, so that P |j = |−j . Now if |j is a mass eigenstate we can form two opposite parity eigenstates P |j, ± = |j, ± of the form assuming both |j, ± are non-zero. Since H and P commute, these would then be degenerate mass eigenstates of opposite parity. For j = 0 the states |−j and |j are orthogonal so that j, ±|j, ± = 1 .
Hence, |j, ± are indeed both non-zero at j = 0. This continuum argument gets weaker on the lattice. It still applies for J = odd states because Z 4 rotations distinguish j = 1 from j = −1 states, so that there is an exact parity doubling for odd spin states. However, 2 = −2 mod 4, so the argument does not apply to even spin states, and the parity doubling is broken by the lattice for these states. Of course, it does get restored in the continuum limit β → ∞ and a → 0.

Approximately rotated operators
In order to get around the breaking of the full rotational symmetry on a cubic lattice, we make use of approximately rotated (AR) operators as introduced in [18]. Figure 2 shows six basic operators that we use, while Figure 3 shows an example of an approximate rotation for one of these. An overall physical scale of these operators is chosen to be similar to the size (after blocking) of the operators used in [22] which had large overlaps with light glueball mass eigenstates.
In order to construct AR operators we make use of additional gauge field matrices corresponding to diagonal links, as illustrated in Figure 3. For instance, the diagonal link matrix in the (1, 1) direction is built by forming U xy (n x , n y , n t ) = 1 2 (U x (n x , n y , n t ) U y (n x + 1, n y , n t ) + U y (n x , n y , n t ) U x (n x , n y + 1, n t )) and by projecting the result back into the space of SU (N c ) matrices exactly as it is done in the smeared link construction. Here U x(y) (n x , n y , n t ) is the gauge field matrix on the link, which starts at lattice site (n x , n y , n t ) and points in the spatial direction determined by the subscript. The definition of a diagonal link matrix in the (1,-1) direction is completely analogous. Note that similarly to the smeared links, we don't use the diagonal links in the definition of the action, but only to construct a basis of AR operators. Using these links we calculated sixteen rotated versions of each basic operator at angular increments of π/8. This allows construction of operators with (approximate) spin values J from 0 up to 8. Clearly, the spin fidelity of rotated Wilson loops is better for shapes that are larger in lattice units. Keeping the physical size fixed, this requires smaller lattice spacing a. An important consequence is that it is not feasible to construct small operators and then scale them up using blocking. We do use smeared operators. However, the inability to implement blocking significantly restricts the quality of the overlap of the AR operator basis with the light glueball mass eigenstates.
Additionally, without blocking the calculation of the loop trace for operators of this size in lattice units adds up to 100 multiplications of SU (N c ) matrices at each lattice point. In [22] glueball masses were extrapolated to N c = ∞ values using masses obtained for N c ∈ [2,16]. However, the authors of [22] made use of blocking and also reduced the lattice size depending on N c to as small values as 26 2 × 30 for N c = 16. In the current work, the spin fidelity of approximate lattice rotations would be unacceptable using a 26 2 × 30 lattice size. In addition, each matrix multiplication requires O N 3 c floating point operations. As a result, our mass determinations are limited to the SU (4) gauge group.
We further enlarge the AR operator basis by attaching pairs of basic operators at the origin as shown in Figure 4. Since we already calculated the product of link matrices around each individual basic loop we can combine them at the cost of just one more matrix multiplication per lattice site for each double loop operator. For any pair of distinct operators we can construct 16 distinct shapes corresponding to different angles between the individual loops. For each resulting double loop we then construct 16 rigid rotations of that double loop.
This procedure allows us to generate a large number of double loops at relatively small computational costs. Given constraints on memory and processing time we restrict ourselves to using double loop operators obtained by combining two of the loops in Figure 2a or two of the loops in Figure 2b. This gives us 2 × 96 = 192 double loop operators before implementing rigid π/8 rotations.
Every operator is calculated after three and five smearings of lattice links. Since we don't perform any blocking, after a few smearings the values of the loop trace for a particular operator at different smearing levels at a given lattice site become highly correlated.
Finally, in order to partially compensate for distortions due to the approximate nature of rotations, we renormalize each rotated operator by replacing it with

Filtering of operators
In [18] the rotational quality of operators was assessed by examining the correlation matrix for rotations of a particular shape and checking that it is well approximated by a circulant matrix, which it becomes in the continuum theory. We take a different approach, dictated by a different use of AR operators in our analysis. Namely, we use them to construct operators with a definite value of J via the angular Fourier transform. These are used then to measure glueball masses in different J subsectors. Hence to test the rotational quality we look at correlations between different J components in the same Z 4 representation. For every Z 4 representation we reject operators for which the absolute value of the overlap between any pair of different J harmonics is greater than 0.02. After variational analysis this choice of threshold is found to result in mass eigenstates that in most cases do not have significant overlaps between states of different J P C .
Nevertheless, we noticed several exceptions from this rule. For example Table 3 shows overlaps between 0 ++ and 8 ++ states (for time separation 4 lattice units; see Equation 10, Section 4.2 for a precise definition of what we mean by overlap at some time separation), which are still statistically significant. At some point in our analysis this effect gave rise to a large systematic underestimation of 8 ++ masses. To avoid this in the variational estimates for the mass eigenstates for each J P C we project out components with smaller J values in the same Z 4 representation and with the same P C values. For example in the case of 8 ++ we project out the overlap with 0 ++ and 4 ++ .
As an illustration of the rejection procedure, Table 2a shows the overlaps between the odd J ++ components of the second operator in Figure 2a. Since at least one of these is greater than 0.02, this operator is excluded from the variational analyses in the 1 ++ , 3 ++ , 5 ++ and 7 ++ representations. On the hand, the overlaps shown in Table 2b for another operator are deemed satisfactory and the corresponding operator is included the 1 ++ , 3 ++ , 5 ++ , 7 ++ variational analyses. .

Mass determination
After constructing a set of M operators φ i with zero momentum and the same definite J P C , which are obtained via angular Fourier transforms, charge conjugation, reflections and spatial averaging of particular loop shapes, we define the cross-correlation matrix where the sum is over mass eigenstates with the same J P C , t is time separation and E n s are eigenstate energies. Given lattice measurements for C (t) for different time lags t we can then recover the mass eigenstates by solving the generalized eigenvalue problem [30] for n = 1, . . . M . If we write .
Of course for any J P C there is an infinite tower of mass eigenstates, while we have only a finite set of operators. As a consequence, each generalized eigenstate ψ n obtained from (9) has non-vanishing overlaps with more than one true mass eigenstate. However, as the time lag t increases, the higher mass components of a particular ψ n get exponentially suppressed. As a result at large t the effective mass E ef f,n (t) approaches the mass of the lightest component of ψ n (of course the ψ n are obtained by variational analysis at some short time separation and may include small components of lighter states, but these will become material in correlators at larger time separations than we effectively fit in this paper).
For each ψ n we fit the normalized correlators c n (t)/c n (0) with a sum of two exponentials and take the minimum of two masses min (m 1 , m 2 ) as an estimate for the mass of the corresponding eigenstate. The fit is performed by minimizing the inner product where Ω is the jackknife covariance matrix of the normalized correlators. As discussed later, a comparison of masses obtained with the AR operators with those obtained with the conventional operators shows that the conventional masses are systematically lower. This indicates that the AR operators have larger residual overlaps with heavy eigenstates. This is not surprising, given that no blocking is implemented for AR operators. In the future this problem can be addressed, for example, by fitting the normalized correlators by a sum of three or more exponentials. Unfortunately, with the currently available amount of data, this gives rise to very large statistical errors which make it impossible to extract the masses to a useful precision. In Tables 4, 5 and 6 we present the mass values in lattice units for 15 lightest states in each channel where we perform the analysis. In Figure 5 we show the same glueball spectrum in a Chew-Frautschi plot (with M 2 and J axes swapped as compared to how it is usually presented), where the squared mass has been expressed relative to the string tension σ f = l −2 s , where l s is the string width. σ f is determined in a conventional manner from the ground state temporal correlator and corresponding energy for fundamental flux tubes winding once round the lattice in the spatial direction [22], [29]) for the value of β = 98.25 used here and the square root is found to be 0.06603 (8). For J > 0 this Figure shows for each J C only the masses averaged over two P values 2 . Also, in this Figure we cut the spectrum at high masses at m = 1 in lattice units. The errors indicate 1σ statistical uncertainties only. As we discuss below, in most cases systematic errors are quite a bit larger than the statistical ones (typically, systematic errors result in masses presented in Figure 5 being heavier than the actual physical masses).
For now, let us neglect possible systematic errors and compare these results with the ASA predictions (see Table 1) for up to N = 3 string level, taking the mass values (or rather, the ordering of masses) presented in Figure 5 at face value. Most conveniently this can be done using Figure 6, where we colored the states according to the ASA level assignment up to N = 3, and cut the heavier states by keeping only one heavier N = 4 state in each J P C channel.
Remarkably, we observe a perfect match to the ASA predictions for the N = 3 states, in addition to the previously observed match for N = 0, 1, 2. In particular, the sharp conflict between earlier preliminary spin determinations for four low lying odd J states and the weak ASA predictions got resolved-our spin values for these states are in a complete agreement with the ASA. We see a well-pronounced group of 1 ± and 3 ± ground states separated by a large gap from heavier states in the respective channels. In Figure 7 we also presented the effective mass for each of these four states as a function of the time lag in lattice units together with the corresponding fits with two exponentials.
This achieves the minimal goal for the current study. We see that the ASA does provide a convincing description of glueball quantum numbers below the N = 4 level. Now, given the large amount of spin determinations for heavier states, it is natural to check the ASA predictions for the N = 4 states as well. Before doing this, let us take a closer look at the systematic errors in our data.

Spectrum with conventional operators and systematics
As a way to cross-check our results and to assess the involved systematics it is natural to compare the AR spectrum with the one based on a conventional set of operators. To achieve this we measure glueball masses from the same set of gauge field configurations using one of the operator bases used in [22]. Of course, spin identifications for the latter states are subject to the usual limitations of square lattices.
The resulting masses are presented in Tables 11 to 18, where we show first 20 states for each set of Z 4 and P, C quantum numbers. In addition we show the normalized overlaps between each corresponding generalized eigenvector and the generalized eigenvectors in different J channels obtained with the AR operator set. Here by a normalized overlap between operators ψ and φ at time separation t we understand As a consequence of this definition at late time one finds if the lowest energy eigenstates appearing in the decompositions of ψ|0 and ψ|0 in the energy eigenbasis are the same and non-degenerate. Overlaps shown in Tables 11 to 18 are defined by where the sum goes over 20 lightest states in each corresponding J channel 3 . For example, in Table 11 we show overlaps between each A 1 , C = +, conventional state with the 0 ++ , 4 ++ and 8 ++ AR states and finally the total overlap with all the latter. In all cases we show overlaps at two different times, t = 0 and t = 4 in lattice units. We see that for most states the total overlap does not add up to unity. Most likely the main cause of this is that both AR and conventional operators have a substantial admixture of heavier states, which still contribute at t = 4. In agreement with this interpretation total overlaps get considerably closer to unity at t = 4 as compared to t = 0. Note that the masses obtained with conventional operators are in a good agreement with the results of [22], while the masses for the same states obtained with the AR operators are quite a bit higher. This indicates that the AR operators have larger admixture of heavy states, as expected, because no blocking was applied for them.
For all states with mass less than 0.80 in lattice units we make a tentative assignment of spin based on overlaps at separation t = 4. In the Tables, the corresponding "dominant" overlap is denoted in colored bold or italic. In many cases this assignment is unambiguous, meaning that the corresponding state indeed contains a clear dominant component with some value of J. However, for a number of states this is not the case, and generalized eigenstates have comparable components with different values of J when decomposed in the AR basis. This happens when there is a group of states with different spins and close values of the mass. In these cases the generalized eigenstates obtained with conventional operators are actually admixtures of several states with different J. Even in these cases usually the overall count for the number of states is unambiguous, even though an assignment of a definite J value to a specific conventional mass eigenstate does not have well-defined meaning.
For example, a situation like this arises for the two lightest E + states. As follows from Table  15 one finds here two nearly degenerate states with a significant admixture of both J = 1 and J = 3 (even though in both cases one of the components is considerably larger than another).
Still, the conclusion that there is one true eigenstate with J = 1 and another with J = 3 at this mass is robust, in agreement with the ASA predictions and the analysis in the AR basis. On the other hand, as can be seen from Table 16, the mass gap between the two lightest E − states is larger and no ambiguity of this kind arises there.
At high spins the conventional basis starts missing states. Neither a clean J = 8 candidate in Tables 11, 12, 13, 14 nor a clean J = 7 candidate in Tables 15, 16 are found. This problem is present even at J = 6-what we marked as the lightest J = 6 + state in Table 17 is one of two nearly degenerate states with a significant J = 6 component (with another such state also being closeby), but for all of these states the J = 2 component dominates. This is not surprising given that conventional operators have only been rotated by multiples of π/2, so their overlaps with highly oscillating large J states are suppressed.
With all these caveats in mind, we have then counted a number of states in each spin channel from the lowest to the highest mass and for N = 0, 1, 2, 3, 4 identified the corresponding ASA level by color coding superimposed on bold font. States in italic have N > 4 according to the ASA. Figure 9 shows all these states in the Chew-Frautschi plot. We again find a well pronounced N = 3 level with the same spin content as observed with the AR basis (barring the above subtlety with identification of the J = 6 state).
As we already said, the AR masses are systematically heavier than the conventional ones. It is natural to use this as an estimate for systematics in the AR mass determinations. To illustrate this effect in Figure 10 we plotted the AR masses with additional dashed bars extending down to the corresponding point in Figure 9. As expected, systematic errors increase with mass. By the time one reaches the N = 4 level they become larger than the level separation.
Another characterization of systematic errors may be obtained by fitting the leading Regge trajectories for both sets of masses. These can be approximated by for the conventional ones. The greater slope in the AR case reflects the increase in systematic errors at larger masses. Note that at large J the trajectory should be well described by semiclassical folded rotating rod solution and the slope should approach 4π. Even with conventional masses the slope is larger than this value by a factor 1.19. Of course, one does not expect the exact agreement at lowest spins, but it is apparent from Figure 9 that the shape of the observed trajectory does not go in the right direction as the spin increases from J = 4 to J = 6. The most likely explanation for this is the presence of a substantial systematic error in the measured mass of the J = 6 state even with the conventional basis of operators. We show effective masses and sum of exponential fits for the leading Regge trajectory states obtained from the AR operators in Figure 8

Finite volume and discretisation errors
Our calculations have been performed on a single volume and at a single value of the lattice spacing. So an important question is whether the finite volume and lattice spacing corrections are small enough that they do not affect our conclusions. Our choice of volume was guided by the finite volume checks in earlier work [22] and our lattice spacing is smaller than the smallest lattice spacing used in that work. However the calculations in [22] involved a smaller number of states than in the present work and there was no attempt to identify the real continuum spins of the states involved. In particular one might expect that states of higher spin, which require higher angular resolution for their identification, might have larger O(a 2 ) discretisation corrections, and this is one reason that we performed our calculations at such a small value of a. In addition to this, states of higher spin should be larger (as will be some of the more massive states of lower spin) and this may lead to larger finite volume corrections. The only unambiguous way to deal with these issues is to repeat our calculations for a range of values of β ∝ 1/a, and also for a range of spatial volumes. Unfortunately this ideal approach would involve a much larger calculation which is beyond the scope of the present paper. Instead what we shall do here is to return to some of the calculations in [22] and focus on a selection of states that have particular relevance to our present work.
We focus on the SU (4) calculations at the three largest values of β in [22], i.e. β = 63, 74, 86. Since masses vary as am ∝ O(1/β) the O(a 2 ) lattice spacing corrections should decrease by a factor ∼ 0.5 when we go from β = 63 to β = 86, so any large lattice corrections will lead to visible changes in ratios of masses as we vary β over this range. The values of the mass gap and string tension are listed in Table 19. The larger spatial size, l = 70a, at β = 86 satisfies lM 0 ++ 21.8 as do the lattice sizes at β = 74 and β = 63. The lattice size in the present paper, l = 80a at β = 98.25, is roughly the same. In Table 19 we have also included calculations at β = 74 on a smaller lattice size so as to investigate finite volume corrections. All these calculations have been performed with the same "conventional" basis of operators used in the present paper. These operators project onto states in the irreducible representations of the square symmetry, so we do not know what their continuum spins are, but we can use what we have learned in Tables 11 to 18 to determine which of these states are interesting for our purposes.
To compare the mass of a state at different β (and volumes) we compare the ratio of this mass to that of the lightest scalar glueball (the mass gap). (As is evident from Table 19 the mass gap, expressed in units of the string tension, is independent of β within the small statistical errors.) The first states we choose to focus upon are the lightest 6 in the A 4 C = + representation. The predictions of the ASA in Table 1 are that these states include 4 J P C = 2 −+ states and 1 J P C = 6 −+ state from the N = 1, 2, 3 levels, and one further state from the N = 4 level. The calculated overlaps in Table 17 support this spin identification (with some possible mixing). We observe from Table 19 that for all of these states any finite volume corrections cannot be much larger than the small statistical errors. Any lattice spacing corrections are also invisible except possibly for the first excited state, but even here any shift is far too small to alter our conclusions about the ordering of the lightest levels.
The second set of states we consider are the lightest 2 in the A2 C = + representation, which we know from Table 12 to be states with J P C = 4 −+ . Here we see no significant finite volume corrections, and any finite lattice spacing corrections are small and within about two standard deviations.
Finally we consider the lightest 3 states in each of the E C = + and C = − representations. From Table 1 we expect the lightest two states to be J = 1 and J = 3 and that they will be close in mass since they both belong to the N = 3 level. The gap to the next state should be larger since it belongs to the N = 4 level. This expectation is supported by the overlap analysis in Tables 15, 16. We see from Table 19 that these states suffer no significant finite volume or lattice spacing corrections.
This comparison of a number of states selected from the N = 1, 2, 3, 4 levels suggests that our calculations and conclusions in this paper are unaffected by significant finite volume and lattice spacing corrections. This check includes a few states with J = 4, 6 but cannot, of course, reassure us about the states of higher spin considered in this paper. For example, J = 8 appears as too high an excitation in the A 1 or A 2 spectra to allow an unambiguous matching between spectra at various β or volumes using just our conventional operators. A larger, dedicated calculation is clearly desirable.

Lattice data vs N = 4 ASA
As should be clear from the results presented in Section 4.2 and summarized in Figure 10, at the current stage the measurements with the AR operator basis are not to be considered as high precision determinations of glueball masses. On the other hand, these measurements do provide a reliable determinations of glueball spins and of the relative order of glueball states, at least for N < 4 levels. Figures 9 and 10 confirm our earlier conclusion that the ASA provides the correct description of glueball quantum numbers at N = 0, . . . , 3 levels.
Let us now discuss what can be said about the N = 4 states. Of course, if one literally took the shifts in Figure 10 as an estimate of systematic errors, it would be very hard to distinguish between N = 4 and N = 5 states. However, these shifts move all the states downwards. Even though the shifts may vary depending on a state or a symmetry channel, still it looks plausible that the AR spectrum correctly captures relative ordering of glueball states in the range of masses corresponding to the N = 4 level and may be used there to assess the validity of the ASA. This is most conveniently done by using Figure 11, which is a direct analogue of Figure 6, but moving one level higher.
We feel that Figure 11 and Table 10 demonstrate an overall good agreement between the N = 4 ASA predictions and lattice data. At some spin values (such as J = 8, 7, 5, 3 and perhaps J = 6) the agreement is in fact very convincing. At the remaining spin values one finds some "contamination" of the N = 4 level by putative N = 5 states. It seems premature to conclude that these contaminations are indicative of a strong tension between the ASA and the lattice data before more precise mass determinations with the AR basis become available.
One should also keep in mind that the number of glueball states rapidly grows with N , while the level separation stays constant in m 2 units. For instance, at N = 5 level the ASA predicts (0 P 1 + 0 P 2 + 1 + 1 + 2 + 3 + 5) 2 = 144 (12) states (their quantum numbers can be easily worked out using the O(2) "multiplication table" presented in Eqs. (9)-(13) of [17]). As a result, one expects the levels to become broader and to start overlap at large N , making the ASA quantum numbers predictions not very useful for many channels. Note that the "clean" spin values in Figure 11 are also the ones with a smaller number of states, suggesting that we are starting to see the onset of this effect. This expected broadening of high levels does not contradict the approximate integrability of the worldsheet theory at high energies assumed in the ASA. Indeed, a typical large N glueball correspronds to a long string state with a broad spectrum of excitations, so that the spectrum of heavy glueballs does not provide a direct probe of the high energy worldsheet dynamics. The broadening does not make ASA completely useless at large N , because one can still follow its validity in exclusive less "crowded" channels. A very interesting example of such a situation happens at N = 4 level. This is the first time the ASA predicts an appearance of (one) exotic spin 0 state with CP = −1. One indeed finds an isolated 0 +− state in Figure 11. Its AR mass is somewhat heavier than other N = 4 masses, but not badly so. With the current quality of data we find that the very appearance of such isolated exotic state in the right ballpark of masses is encouraging.
Note that using the conventional operators in Figure 9 one also finds an isolated 0 +− state but it appears to be separated further from the rest of N = 4 states (as far as counting of "contaminating" N = 5 states goes). However, by inspecting Table 13 one finds that this state is a member of an extended group of A − 1 states with relatively large and overlapping statistical error bars, so it well may be that the 0 +− state is misidentified in Figure 9 and a lower member of that group should be used instead.

Conclusions and Future Directions
To summarize, the main conclusion from the presented results is that the ASA provides an accurate description of the glueball spectrum in D = 3 YM at least at the N = 0, . . . 3 string levels, which correspond to the first 39 glueball states. It also provides, at the very least, a good organizing principle to describe the glueball spectrum as a whole. There is a number of directions to improve and extend these results both on a lattice and on a theory side.
On a lattice side it will be important to reach the stage when the AR basis can be used for high precision mass spectroscopy. This requires a significant reduction in the systematic overestimate of the glueballs masses which we observe with the present AR basis. This goal may be achieved with future simulations with a larger basis of AR operators and with higher statistics. In addition, for a cleaner test of the ASA, it is important to extend the analysis to a larger number of colors N c , similarly to how this has been done in the conventional spectral measurements in [22]. When these goals are achieved, it will be possible to check that the shape of the leading Regge trajectory agrees with the effective string theory predictions at large spins (this may require extending our analysis to J > 8).
As we discussed, due to level broadening, the ASA stops being a useful description of the glueball spectrum at higher N in all channels. Still, it will be interesting to check its predictions in some clean exclusive channels at N > 4 with the improved precision of mass determination. A particularly interesting class of heavy states to study are the exotic CP = −1 J = 0 states. At N = 5 the ASA predicts another isolated exotic state (which can be either 0 +− or 0 −+ , while at N = 6 there should be six such states (they can be either 6 0 +− s or 4 0 −+ s and 2 0 +− s or evenly split between 0 +− s and 0 −+ s).
In addition to looking at the exclusive channels the inclusive counting of states is still interesting even if different levels cannot be clearly separated as N grows. Indeed, perhaps the most striking prediction of the ASA in D = 3 YM is the absence of any massive excitations on the string worldsheet. If present, massive modes would violate the tensor square structure (2). It is hard to test (2) directly, when the levels start to overlap. However, the presence of massive modes would also imply a much faster growth in the number of the glueball states, compared to the ASA, which can still be tested. On a theory side it is important to push the ASA spectrum beyond the ansatz for quantum numbers. In general, this requires a comprehensive understanding of the integrable approximation, which is lacking at the moment. However, certain progress can be achieved, especially at large J, using the effective string theory technology. For instance, as clear from the glueball spectra based both on the AR and conventional basis, the mass of glueballs at the same level tend to decrease as the spin gets smaller. This fits well with the presence of an additional attractive correction to the Nambu-Goto phase shift extracted from long fluxtube data [14,28]. Indeed, the smaller J glueball states are obtained by adding short wavelength perturbation to the classical rotating rod solution, describing the leading Regge trajectory. Additional attraction between these excitations makes their energy smaller compared to the Nambu-Goto prediction (the latter corresponds to an exactly degenerate level).
A further inspection of Figure 11 suggests additional patterns of glueball spectrum at large J. Namely, it appears that at each N the highest spin (2N ) + state is followed by a pair of (2N − 2) + and (2N − 2) − states, where the latter is considerably lighter than the former, and then, perhaps, by a pair of nearly degenerate (2N − 3) ± states (at N > 2). If correct, this pattern should be relatively straightforward to calculate within the effective string theory using the phase shift correction from [14,28] as an input.
Finally, it will be most important to extend these results to D = 4 YM. On a theory side this is more complicated due to the presence of a massive excitation on the worldsheet of D = 4 confining strings [2,13]. Hence, it is likely that the lattice input will play an important role in achieving this goal.