The spectrum of 2+1 dimensional Yang-Mills theory on a twisted spatial torus

We compute and analyse the low-lying spectrum of 2+1 dimensional SU(N) Yang-Mills theory on a spatial torus of size l × l with twisted boundary conditions. This paper extends our previous work [1]. In that paper we studied the sector with non-vanishing electric flux and concluded that the energies only depend on the parameters through two combinations: x = λN l/(4π) (with λ the ’t Hooft coupling) and the twist angle θ˜\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \tilde{\theta} $$\end{document} defined in terms of the magnetic flux piercing the two-dimensional box. Here we made a more complete study and we are able to condense our results, obtained by non-perturbative lattice methods, into a simple expression which has important implications for the absence of tachyonic instabilities, volume independence and non-commutative field theory. Then we extend our study to the sector of vanishing electric flux. We conclude that the onset of the would-be large-volume glueball states occurs at an approximately fixed value of x, much before the stringy torelon states have become very massive.


JHEP07(2018)169
The present paper extends our previous study [1][2][3] on the behaviour of pure SU(N) Yang-Mills theory in 2+1 dimensions, where space is compactified as a 2-dimensional torus with 't Hooft twisted boundary conditions. We focus upon the dependence of the spectrum of the theory on the parameters that define it, namely the torus size l (restricted to be rectangular of equal length in both directions), the number of colours N , and the magnetic flux k (a modulo N integer, coprime to N ) introduced by the boundary conditions. There are certain results and observations that suggest that these parameters conspire to produce a simpler result. The dependence on the value of the coupling constant λ ('t Hooft coupling) is implicitly contained, since it is used as the unit of energy or inverse length. In our first paper on the subject [1], we studied the spectrum for small torus sizes using perturbation theory. Our study showed that the individual energies depend only on two combinations of parameters: the product of N times the torus lateral size and the following angular variable: θ = 2πk/N , wherek is the modular multiplicative inverse of k (kk = 1 mod N ). The second parameterθ/(2π) is a rational number, but continuity of the energies with respect to it takes place in perturbation theory. Our analysis then used the lattice regularized version of the model to study how the energies evolve with the torus size from the perturbative results to the expected confinement behaviour. The results supported that the same two combinations described the dependence of the energies for all torus sizes and that continuity inθ still applied. The dependence on the product lN , or rather the dimensionless ratio x = λlN/(4π), is connected to the phenomenon of volume independence at large N , since in that limit x goes to infinity for any value of l. The question has been studied in 4 dimensions going back to the old observation of Eguchi and Kawai [4] that Schwinger-Dyson equations are independent of the volume at large N . However, the proof assumes that center symmetry remains unbroken, a feature that turns out to be wrong in the weak coupling regime [5]. Several ideas have been presented over the years to solve this problem. A simple proposal, which is adopted in this paper, is to use 't Hooft twisted boundary conditions [6,7]. Other alternatives have been proposed, such as that of adding fermions in the adjoint representation [8,9] and other modifications of the reduction idea [10]. Recent tests give compelling evidence that indeed these methods do rescue the idea of volume independence at least in 4 dimensions [11,12]. The old studies only required that the flux k associated to the twisted boundary conditions should be irreducible. However, more recent studies [13][14][15][16] showed problems and signs of symmetry breaking, indicating that the choice of the magnetic flux k is crucial to avoid instabilities and to reduce finite N corrections. In ref. [17] two of the present authors proposed that k/N andk/N have to be kept finite (and beyond a certain threshold) to preserve the validity of volume independence at large N . Although there is a rationale behind the bound on k/N , there is no rigorous argument and no precise estimate of its value. However, the practical attitude has been to show that volume independence can be used to obtain precise results about the large N infinite volume theory. This is similar to the use of perturbation theory in theories where no rigorous proof of summability is available. Our study in the simpler 2+1 dimensional case

JHEP07(2018)169
is aimed precisely at clarifying some of these points. Not much work has been done in this system in this context apart from that of Narayanan and Neuberger [18,19].
There is another perspective from which our work has a long-standing interest: noncommutative field theories [20]. Indeed, twisted boundary conditions and non-commutative field theories are intimately connected. The first appearance in the literature of the action and Feynman rules of gauge and non-gauge non-commutative field theories was as a generalization of the volume reduced twisted models [21]. The twisted Eguchi-Kawai model was also used as lattice-regularized version of non-commutative Yang-Mills theory [22,23]. As a matter of fact, Morita equivalence implies that gauge theories on the torus with twisted boundary conditions are particular cases of non-commutative field theories on the non-commutative torus with special values of the non-commutativity parameter. Indeed, one in which a dimensionless combination of the torus size and the non-commutativity parameter takes rational values. Our observations of refs. [1,2] appear natural within this context, since lN is the size of the non-commutative torus andθ/(2π) the rational dimensionless non-commutativity parameter. Several studies showed that although the torus size tames down the IR/UV mixing [24], new instabilities, called tachyonic instabilities, could appear [25,26]. These are of the type associated with spontaneous breaking of centre symmetry and breakdown of volume independence. Our perturbative analysis of refs. [1,2] agrees with this conclusion, but also shows that a suitable choice of the flux k can help to avoid them. Furthermore, our work also serves to extend this analysis to larger torus sizes for which the perturbative calculation breaks down. Hence, our analysis is capable of addressing some of the questions raised long time ago within non-commutative field theories beyond the domain of perturbation theory. Here we signal out very specially the ideas and hypothesis formulated in refs. [27,28].
Having set the general context of our present study, we go into an overview of the specific goals of the present paper. The spectrum of states can be split into sectors corresponding to different values of 't Hooft electric flux, a two-dimensional vector of integers modulo N . Our previous papers only studied certain sectors with non-vanishing values of this electric flux. In this paper we will extend the range of our study to cover all electric flux sectors for a wider interval of values of N , torus size and flux parameters. This will allow us to achieve a semi-quantitative phenomenological description of the dependence of these energies on the parameters. In addition, our study will also include the sector with vanishing electric flux, the glueballs. As we will see, the possibility of obtaining large-volume results from our large-N analysis, as implied by the volume independence hypothesis, seems challenging. The reason being the presence of an ever increasing number of torelon pair states whose mass only grows with l, not lN . Our results provide a consistent solution to this puzzle.
The layout of the paper is as follows. In section 2 we present the necessary background to the problem, including previous results and a summary of what is known and expected for the system. In the next two sections we present our main results, dealing with nonzero (section 3) and zero electric flux states (section 4) respectively. To facilitate reading we have divided each section into three main parts: a preamble, a part in which the nonperturbative results are presented, and another one where the conclusions from these results

JHEP07(2018)169
are extracted. Within the same spirit, we have concentrated all technical aspects about the methodology used, which is by itself quite challenging, to an appendix. A list of tables provide the actual measured numbers for the sake of other researchers who might want to analyse the huge wealth of information that we have obtained. The paper closes with a brief conclusion section which gives a big overview and lists open problems and paths for improvement.
2 SU(N ) gauge theory on a spatial two-dimensional twisted box We will be considering a SU(N ) gauge theory defined on T 2 × R. For simplicity, the spatial 2-dimensional torus is considered to be symmetric and of period l. The gauge potential satisfies twisted boundary conditions, given by with SU(N ) matrices Γ i subject to the consistency condition: derived from imposing univaluedness of the gauge potential under displacement along the two cycles of the torus. We will be considering the case of the so-called irreducible twists, in which the magnetic flux integer k is taken to be coprime with N [29]. In that case, eq. (2.2) defines the matrices Γ i uniquely modulo global gauge transformations. Twisted boundary condition on a torus were introduced by 't Hooft [30,31] as a way to induce topological (chromo-) electric and magnetic fluxes in Yang-Mills theories. This is best understood in the Hamiltonian formalism in the A 0 = 0 gauge. In this set-up, time-independent large gauge transformations with non-trivial periodicity: States are thus classified by a 2-dimensional vector of integers e defined modulo N : the electric flux vector. In this way, the Hilbert space decomposes into Z 2 N disjoint sectors parameterized by the value of the electric flux. The vacuum and the glueballs live in the sector with zero electric flux, while Polyakov loop operators with non trivial winding acting on the vacuum generate non-zero electric flux states, the torelons [32]. As mentioned in the introduction, the purpose of this paper is to present the results of a non-perturbative analysis of the volume and N dependence of the spectrum in the different sectors, extending to the zero electric flux sector the results obtained in ref. [1] for non-zero fluxes. This will be done in sections 3 and 4 where we will analyse the spectrum obtained from a lattice Monte-Carlo simulation of the 2+1 dimensional system. Before doing that, it is JHEP07(2018)169 instructive to discuss what is the expected volume dependence based on what we know from perturbation theory and from the large volume, confinement regime.
Asymptotic freedom implies that perturbation theory is a good approximation for small torus sizes. The calculations with twisted boundary conditions are easily performed when using an appropriate basis of the SU(N ) Lie algebra [7]. In this basis the vector potential can be expanded in a modified Fourier expansion: where the momentum dependent matricesΓ( p ) satisfy Using this formula it is easy to verify that the twisted boundary conditions amount to the quantization of momenta p = (n 1 , n 2 ) p 0 , n i ∈ Z, where the quantum of momentum p 0 = 2π/(lN ). This corresponds to the standard minimum momentum for an effective box sizel ≡ lN . TheΓ( p ) can be written explicitly as follows: wherek is the modular multiplicative inverse of k: For irreducible twists, there are N 2 independentΓ matrices of this sort that can be chosen as those with n i taking values from 0 to N −1. The one corresponding to n = 0 is proportional to the identity, while the remaining N 2 −1 matrices are traceless and provide a basis for the (complexified) SU(N ) Lie algebra. The Fourier coefficientsÂ i (t, p ) are complex and satisfy a hermiticity condition similar to the standard one, which restricts the vector potentials to live in the standard real SU(N ) Lie algebra. The primed momentum sum in eq. (2.5) runs over all momenta that lead to tracelessΓ matrices, excluding those with non-zero trace corresponding to n = 0 (mod N ). Notice that this restriction implies in particular that zero-momentum is forbidden in the twisted box, and the minimum value of | p| = p 0 . Equipped with the previous formalism it is very easy to compute the spectrum to leading order of perturbation theory. At this order the system can be described as a gas of free massless gluons whose energy is just given by the modulus of its momentum E = | p|. The ground state or vacuum is the state with no gluons and has zero energy. As a result of what we discussed in the previous paragraph the spectrum has a gap corresponding to a single gluon of minimum momentum p 0 (it is 4-fold degenerate). If we write the energy in units of 't Hooft dimensionful coupling λ this gap becomes 1/2x, where

JHEP07(2018)169
is the quantity introduced in ref. [1], and which we argued is the relevant dynamical variable, setting the scale both in the perturbative and in the non-perturbative regimes. We will review our arguments below. But, how does this leading-order spectrum correspond with the electric flux sectors mentioned earlier? To see this let us consider a spatial Polyakov loop with non-trivial winding in the spatial torus. With twisted boundary conditions the Polyakov loop adopts the form: where γ is a closed curve on the 2-torus and ω is the corresponding winding number. The state obtained by acting with this operator over the vacuum has an electric flux given by ω modulo N . Expanding the ordered exponential of the Polyakov loop and using the Fourier expansion of the vector potential given earlier, we conclude that a gluon of momentum p = np 0 carries electric flux given by e = (n 2 , −n 1 )k mod N. (2.12) The fact thatk is defined modulo N can be expressed by saying that the quantityθ defined as:θ is an angle. This is the other main quantity introduced in ref. [1] to describe the spectra. At this stage it must be said that our two main quantities have a natural interpretation within the non-commutative field theory description. The effective sizel is just the size parameter of the non-commutative torus andθ/(2π) =k/N is the dimensionless noncommutativity parameter.
Inverting the previous formula eq. (2.12) one can determine the minimum momentum corresponding to each electric flux e as follows p c ( e) = (±n 1 , ±n 2 ) p 0 with 14) where we have followed the notation in [33], with ||s|| denoting the distance of the real number s to the nearest integer. Hence, the lowest perturbative energy, corresponding to a one gluon state carrying momenta p = np 0 with | n| = 1, belongs to the electric flux sectors e = (±k, 0) and (0, ±k) -see eq. (2.12). For multiple gluon states the momenta add up and hence, so does the corresponding electric flux. Notice that conservation of electric flux directly follows from conservation of momenta. It can happen that there are multigluon states degenerate with the single gluon state. For example, the state with two gluons of momenta (p 0 , 0) each, is degenerate with the state of one gluon of momenta (2p 0 , 0). This degeneracy only occurs for collinear gluons. Excited states in the zero electric flux sector can be obtained to leading order in perturbation theory as multigluon states. The first excited state (the mass gap) is given by a pair of gluons of opposite momenta equal to the minimum one | p| = p 0 . In λ units the energy is just 1/x. Generically we might call these states glueballs, since they have the JHEP07(2018)169 same quantum numbers as the corresponding states at large volumes. However, later we will reserve that name to the states that are present in the infinite volume theory, while these other will be referred as torelon pairs. In the following sections we will present the results of our study of the spectrum separating the cases of non-vanishing and vanishing electric flux. 3 The torelon spectrum (non-zero electric flux)

General considerations
In this section we will deal with the case of non-vanishing electric flux. The corresponding energy eigenstates were called torelons in ref. [32]. This was already studied in ref. [1], but we will present additional results which will allow us to draw certain conclusions from them.
We will label the torelon energies in units of the 't Hooft coupling λ for each momentum value p = p c ( e) = np 0 by the symbol E n . In ref. [1] we computed the next to leading order perturbative contribution to these energies coming from self-energy gluon diagrams. This combines with the leading order result into the following expression: with the part in 1/x representing the gluon self-energy given in terms of the function: with θ 3 the Jacobi theta function [34]: One can see that for each value of n the energy square depends only on the two combinations of the arguments x andθ. Notice that the argument of the function G is just the (rotated) electric flux over N : ijθ n j /(2π) = e i /N . There are some interesting observations following from eq. (3.1). The first is that the simple form appears when writing the energy square and not the energy. Since the first term is the momentum of the gluon square, the second can be interpreted as the mass square. However, the latter is actually negative and would eventually drive negative the energy square at some finite value of x. This would signal a phase transition, that for obvious reasons is called a tachyonic instability. When the problem was studied many years ago within the context of non-commutative field theories [25,26], it was correctly understood that the new phase would be one in which there is condensation of electric flux sectors into the vacuum and centre symmetry spontaneously broken.
The arguments in favour of a tachyonic instability are non-rigorous because they are based on a truncated perturbative expansion. However, the function G( z) has a pole when the argument takes integer values, and hence the self-energies for large N and small electric fluxes grow, making the transition point occur at small values of x when the effect of higher order terms might be considered negligible. Hence, one of the goals of our study is that of JHEP07(2018)169 addressing whether such a phase transition exists or not. The result has implications for the non-commutative field theory program.
With or without transition, when the torus size is large the system should return to the confinement phase in which these electric flux energies would grow linearly with the torus period and eventually decouple from the rest of the system. The behaviour is expected to be that following from an effective string description. The two leading terms at large x for the energy square should take the following form The second term is the linear energy growth term whose magnitude is determined by the string tension σ in units λ 2 . The function φ gives the k-string spectrum and it depends on the electric flux. The first term in the right hand side of eq. (3.4) is the leading correction, related to the sometimes called Lüscher term [35], where the function χ is expected to be of order 1 and to have a mild electric flux dependence. It is interesting to point out that, although valid in the opposite regime than the perturbative formula, it shares some properties with it. First of all, that the formula simplifies when expressed in terms of the energy square. Indeed, the truncated formula is actually exact for the Nambu-Goto string. The second is that, once more, all the dependence comes through the two combinations x andθ. Our conjecture in ref. [1] was that actually, at all values of the torus size, the energies only depend (continuously) on these two quantities. Thus, one should obtain similar energies with two very different values of l and N provided the product remains constant and the parameterθ varies only slightly. This statement is stronger than the ordinary volume independence since it is valid also at finite N . However, it relies upon continuity inθ.
In the next subsection we will present the results of our non-perturbative study aiming at testing our conjecture and investigating the transition from the perturbative regime to the confinement regime.

Non-perturbative results
To compute the energies in a non-perturbative fashion we used the lattice formulation. We will first explain in very simple terms the essentials of the calculation and collect all the technical details in the appendix.
The SU(N) model is formulated on a finite lattice of spatial size L × L and twisted boundary conditions with flux k. The model depends on a single coupling b = 1/(aλ L ). The lattice energies E L are extracted from the exponential decay in time of correlation functions of Polyakov lines with the appropriate winding and projected to the corresponding minimum momentum value. One has still many possible operators to use, and that redundancy is used to maximize the coupling of the operator to the state whose mass one is trying to determine. The lattice quantities are all dimensionless because they are in units of the lattice spacing a. Dimensionless quantities are directly comparable. Hence, E L b = E and L/b = λl. Strictly speaking this identification has lattice artefact errors. These disappear when taking the lattice spacing to zero, implying b and L going to infinity with the JHEP07(2018)169  Table 1. Values of N , L and k used in the non-perturbative lattice simulations. ratio fixed. Comparing the results at different values of b and L one can quantify these errors. In our case, they are not expected to affect the main conclusions of the paper. We leave the detailed explanation of the procedure and other technical aspects to the appendix and in the following we will present the results.
Our study covers a wide range of values of N , k and x. In particular, we analyse all coprime values of k for 4 different gauge groups N = 5, 7, 11, 17 (cf. table 1) over a wide range of values of x -from the large-volume region to the small-volume perturbative one. Given that on the lattice one has x = N L/(4πb), we proceeded as follows. The results are divided in three groups of approximately constant value of N L. The coupling was then slightly adjusted to have the same set of x values within each group. In the second group L and b are doubled and in the third one multiplied by 4. This corresponds to dividing the lattice spacing by 2 and 4, as a way to measure the systematic lattice artefact errors. In addition, we also analysed N = 13 with k = −k = 5 for a subset of values of x, and N = 34 and N = 89 (with k = −k = 34 and k = −k = 13 respectively) for only a few values of x. These last additions correspond to values of N and k belonging to the Fibonacci sequence, which play a special role according to ref. [33] as will be explained later.
Our results extend those obtained in ref. [1] which concentrated mostly on the minimum momentum state | n| = 1. The most important properties of our result will now be enumerated.

Consistency with expectations for small and large volumes
At small x our energies agree with the values obtained by the perturbative calculation within errors. At large x they are consistent with the confinement behaviour. Indeed our calculation allows to obtain unprecedented information about the k-string spectrum: the dependence of the string tension on electric flux. This is encoded in the function φ(θ n/(2π)) appearing in eq. (3.4). One can extract this function from a fit to the data to be described later. The result is given in figure 1, where we display the values of the function σ φ( e/N )/λ 2 obtained from the fits compared to the function which is one of the characteristic dependences in simple models. We also obtain a value of the string tension √ σ /λ = 0.213(1), which is consistent with the value previously obtained in ref. [1]. This can also be compared with the analytic prediction of Nair [36] 1/ √ 8π and the value 0.19636(12) obtained from a recent lattice calculation [37]. We do not give much significance to this 7% difference, since the large x region is precisely where the lattice spacing corrections are expected to be larger. A detailed comparison would demand a continuum extrapolation and a thorough analysis of the systematic errors. This was not one of the main goals of this work, which covers such a wide range of sizes and values of N .

Transition from small to large volumes
Our data allow us to follow the dependence of the torelon energies at all values of x. The energies, that decrease as 1/x for small sizes, reach a minimum and then start to rise and eventually go linearly with x as predicted by confinement. Except in some cases to be mentioned later, this transition is smooth. This is illustrated in figure 2 and figure 3 for different values of N and k.
The remarkable fact is that the data are well described by the continuous lines which come from the following function    which is obtained by simply adding the leading formulas for large and small x. The function φ 0 , given in eq. (3.5), and the string tension value √ σ /λ = 0.213 were explained earlier. The Lüscher term function χ(z) which was expected to a have a slight electric flux dependence is fixed to a constant χ 0 = 0.6. If one tries to do better one can fit a different constant for each value of the argument z. The result is displayed in fig 4. In the range z ∈ (0.15, 0.45) it is consistent with a constant value χ(z) = 0.59(3) with a χ 2 per degree of freedom of 1.4. Hence, for our descriptive purposes it is much simpler, and hence better, to fix the value to 0.6.
It must be mentioned that eq. (3.6) is not supposed to be exact. The intermediate region might have a somewhat more complex structure. Indeed, in ref. [1] we added a term of the form With these two extra parameters we were able to fit the data points with good χ 2 . This holds also for the new data and has been used to obtain the results for the k-string spectrum presented in figure 1. Nonetheless, since our goal is more to understand the general dependence on the arguments than to obtain a precise determination of the parameters, we prefer the simple description provided by eq. (3.6). In terms of only two parameters σ /λ 2 and χ 0 we are able to approximate hundreds of measured values to within a few percent. This is clear in figures 2 and 3 but also holds for all the other values not displayed.

Continuity inθ
Our main hypothesis formulated in ref. [1] was that the energies depend on l, N and k only through the two combinations x andθ. However, since the latter quantity only takes discrete values for each finite N , the hypothesis necessarily implies continuity iñ θ. In principle, the validity of this assumption follows as a consequence of the capacity of our previous parameterization to fit the data. Nonetheless, thanks to the fact that we constructed our data at specific values of x we can also make a direct check of this continuity. In figure 5 we plot E When testing higher values of the momenta at small torus sizes a complication arises due to the existence of degenerate states at lowest order. The simplest situation occurs for p = (2p 0 , 0). The self-energy corrections are different for each state. A similar ambiguity appears at large volumes. We can have either a single string carrying the full electric flux or two strings with electric fluxes that sum up to the total one. The numerical method will just select the minimal energy of the two. Which one is smaller might depend on the value ofθ. Indeed, that seems to be happening according to our data. This is illustrated in figures 6a, 6b, where we show the dependence of E (n,0) onθ, for states with momentum | p| = 2p 0 and 3p 0 respectively, at several values of x. For clarity, the results at different values of x are displaced vertically by 0.4 starting from the lowest x-value. The continuous lines correspond to the predictions of our simple parameterization eq. (3.6) for E (n,0) and nE (1,0) . The data points indicate the presence of level crossings, but it is remarkable that our simple formula is able to predict where these crossings will appear. As we will see this simple exercise is a prelude of the full analysis done in the following subsection.

(Non-)existence of tachyonic instabilities
As mentioned in the introduction the negative value of the self-energy correction might lead to the condensation of some of the electric fluxes. This happens when the energies cross zero. In our previous publications [1,2]

JHEP07(2018)169
and cases in which it doesn't. Typically, the risk is higher for larger values of N and smaller values of the flux |k| or |k|. This fits nicely with the general proposal made in ref. [17] to avoid symmetry breaking in the 4 dimensional model. However, in our case it turns out that the approach of the energies to zero can be directly deduced from the approximate parameterization of eq. (3.6). The formula then provides a method to indicate whether there will be condensation of any electric flux torelon or not. This works perfectly for all of our data set, which covers a wide range of values of N , k and λ, but also allows a theoretical analysis of the presence or absence of tachyonic instabilities for arbitrary values of N and k. This was done in ref. [33]. Essentially, in order to enforce the absence of tachyonic instabilities one should guarantee that Thanks to our formula this translates into where the symbol e ⊥ N indicate that the two integers are coprime. Is it possible for any N to choose a flux k that guarantees this condition is met? In ref. [33] it was shown that this question translates into a still unproven conjecture in Number Theory, called the Zaremba conjecture [38]. Fortunately, it has been proven recently [39] that eq. (3.9) holds for almost all values of N . The analysis also shows that the largest values of Z min (N, k), and hence of the minimum energy, occur for N = F n , a member of the Fibonacci sequence, and k = |k| = F n−2 , for any value of n. Motivated by this result we included the simulation at N = 13 and k = −k = 5 which satisfies the criteria of maximal Z min (N, k). The corresponding energies are plotted in figure 3b. Next to it we have plotted the energies for N = 5 and k = 2, another Fibonacci optimal case. Notice the resemblance, and the ability of our simple parameterization to describe both data. Hence, we decided at the very late stage of this project to make some test runs at even larger values of N . In particular we studied N = 34 for L = 2 and k = |k| = 13 and N = 89, L = 1 (only a 1 point spatial lattice) and k = |k| = 34, which belong to the Fibonacci sequence. We only studied few values of x and with limited statistics. Furthermore, the code is different and we include less operators in the game. Nevertheless, the results came out to be rather good. In particular, we looked at the value x = 3.183 which was studied for all the other values of N . The time correlators of Polyakov lines show a clear exponential fall-off at moderate separations, which allowed us to obtain masses. In table 2 we show the values obtained compared with the expectations from our simple formula eq. 3.6. For completeness we also add the results for N = 13. Notice, the good agreement specially for the lighter states. The bigger value forn = 1 can be due to contamination with excited states due to the limitation in number of operators.
We will study some consequences of these findings in the next subsection.

Derived consequences
In this subsection we will extract several conclusions from the results presented in the previous subsection.

Conditions on the flux for the absence of phase transitions
As we saw in the previous section it seems possible (except perhaps for a set of exceptional N values) to choose the flux k such that the torelon energies are always bounded from below. This guarantees the possibility of taking the large-N limit avoiding any transitions as we move from small to large volumes of space. The condition is eq. (3.9) which is stronger than the previously defined |k| N , |k| N 0.1. Thus, our proposal made in ref. [17] might have to be modified accordingly. Nonetheless, the first prime number N in which the previous condition does not suffice is N = 61 (with k = 21) and the first non-prime N = 44 and k = 21. Notice that the pattern suggests to remain far from a simple fraction i.e. 1 2 , 1 3 , etc. Although our new condition has been obtained for the 2+1 dimensional system it might also affect the 4 dimensional case, but replacing N by √ N . Indeed, it can be shown that the value of Z min controls the size of the contribution of non-planar diagrams in perturbation theory [40]. However, if that also affects the existence of symmetry-breaking phase transitions in the non-perturbative region of the 3+1 system would only show up for N = 44 2 = 1936 or higher, which is hard to check. To avoid these problems, choosing k such that Z min (N, k) is large enough is recommended.

Implications for non-commutative field theory
As mentioned in the introduction, the system that we are studying is a particular case of a gauge theory on the non-commutative torus of size proportional to x and noncommutativity parameter proportional to x 2θ . The quantityθ/(2π) appears as a dimensionless non-commutativity parameter which takes rational values. Since the present theory is perfectly well defined both perturbatively and non-perturbatively, this would provide a window into the non-perturbative dynamics of non-commutative field theories. But is it possible to use this theory to investigate non-rational values of the dimensionless noncommutativity parameter? This is one of the questions raised in an interesting paper by Alvarez-Gaumé and Barbón [28] (see also [41]) where they also analysed in detail the possible limits that can be taken and their interpretation. They proposed to define the theory for irrational values ofθ by taking sequences of rational numbers converging to it: Using the new information that we have gathered in this work and in ref. [33], we posed ourselves the question if it is possible to take the limit avoiding the presence of tachyonic instability for all intermediate values where the so-called partial quotients a i are positive integers. Let us now call A max (N,k) = max i a i . In ref. [33] it is proven that Our condition on the absence of tachyonic instabilities (eq. (3.9)) then imposes A max < 10. The possible limiting irrationalsθ should have an infinite continued fraction built from an alphabet of only 9 digits. In this case it can be seen (see ref. [39]) that the set is uncountable, has zero measure and a Haussdorff dimension between 0 an 1. The conclusion is then quite striking. Only for a zero-measure set of values ofθ can one define a limiting theory without tachyonic instabilities. This is so since as N increases new light states appear in the spectrum that might cause a phase transition. 4 The glueball spectrum (zero electric flux)

General considerations
Yang-Mills theory is expected to have a mass gap in the large-volume limit and a whole spectrum of states, which are generically called glueballs. The masses of these states will be finite in λ units. In that limit the torelon states become infinitely massive and decouple from the theory. The value of the glueball masses are expected to remain finite when taking the large-N limit, defining the spectrum of the large N infinite volume theory.
If we now consider space to be a very large torus with twisted boundary conditions, these glueball states will still be there and with masses that are almost insensitive to the torus size l and the magnetic flux k. However, there will also be other states in the zero electric flux sector which are made of torelon pairs with masses proportional to the size of the torus l.
As we reduce the torus size the hierarchy between torelon pairs and glueballs gets reduced and there could even be mixing between them. Indeed, in the limit of very small volume the torelons are just single gluons and the lowest energy states in the zero electric flux sector are these torelon pairs. Making the volume slightly larger these masses get corrections which, within that domain, will continue to depend jointly on lN (actually on x) and also onθ. Our calculation at next to leading order of perturbation theory gives modifications to the masses of these two-gluon states coming from the self-energy of the individual gluons, plus an interaction term that splits the masses of the states into those that transform differently under the discrete rotation group. The state with the lowest energy is the rotationally invariant one and the interaction energy equals [1] written in terms of the twist angleθ. How does the transition from small to large volumes take place? Is there any type of volume independence present in this spectrum? In principle the answer to this last question seems to be no. As we saw in the previous section as the volume becomes larger JHEP07(2018)169 the lowest-lying torelon spectrum is far from being volume independent. When N grows new light states appear and these have energies that go with l, not lN . The corresponding zero-electric flux pairs will also have energies that grow with l. Thus, the statement that the spectrum of the theory at large N does not depend on the volume does not hold. At x fixed with low l and large N , some torelon pair states will remain in the low-lying spectrum as opposed to the infinite volume case. These considerations serve to prepare for the presentation of our non-perturbative results interpolating between both regimes.

Non-perturbative results
In this section we will present our results about the spectrum of the theory in the sector of vanishing electric flux. The methodology is similar to the one used for torelon states. Masses are extracted from the exponential decay in time of correlators between operators carrying no electric flux and projected to zero momentum. The spectrum is analysed using the GEVP method (see appendix) applied to the basis of operators described in section A. Let us focus first on the dependence for the mass gap. For x smaller than 1 or 2 it follows the characteristic 1/x behaviour predicted in perturbation theory. Indeed, the mass is quite close to twice the minimum torelon mass (| e| = |k|), indicating that this state is actually a torelon pair state. To see this, we display in the figure the lines obtained by multiplying by 2 the simple parameterization eq. (3.6) of the corresponding torelon. The difference between the double torelon mass and the mass gap is a measure of the torelon-antitorelon interaction energy. At higher values of x the mass gap has a minimum in the interval [1.5, 2] and then starts to rise. If we first focus on the N = 5 case (figure 7) we see that beyond x ∼ 4 the mass seems to level up and tends to a constant. The corresponding state in this larger x region would correspond to the true lowest mass glueball which is there in the infinite volume limit. Indeed, the value of the mass is rather compatible with the value m/λ = 0.873 (8) given in ref. [42] obtained for N = 5 on a 32 3 lattice at x = 7.76. This value is indicated by the horizontal red band in the figure. Hence, the N = 5 result is as expected and indicates a change of nature of the lowest mass state somewhere in the interval x ∈ [3.5, 4.5]. This change occurs when the mass of the torelon pair becomes higher than the mass of an initially more massive state which ultimately becomes the glueball. This is clear when looking at the energy of the next excited state, which beyond x ∼ 4 seems to extend nicely the behaviour of the lower x mass gap with the characteristic linear growth of a torelon-antitorelon pair. In the figure we also show the curve corresponding to twice the energy of the torelon with electric flux e = (|k|, |k|), which, consistently with expectations, describes the behaviour of the excited state at low values of x.
If we now focus on the N = 17 data (figure 8) we see that the behaviour of the mass gap and the next excited state for small x follow the same pattern as for N = 5, making clear its torelon-antitorelon nature. However, the behaviour of the mass gap for x > 4 is quite different to the N = 5 case, since the mass is actually decreasing in that region. The interpretation becomes obvious once we compare the data points with the blue line which is the predicted behaviour for a pair of torelons with opposite electric fluxes | e| = 1. This was the expected trouble for larger values of N arising due to the appearance of new light torelon pair states.
One may wonder if this result implies a failure of x-scaling in the glueball spectrum. We will argue this is not the case. For that purpose we have to look at the next excited state for the SU(17) glueball sector. We see on figure 8 that beyond x ∼ 3 its mass approaches the value corresponding to the infinite volume glueball, with very little x-dependence. These results indicate that, although the SU(17) glueball is not the lowest excited state in the zero-electric-flux sector, it is present in the spectrum for x 4 with a mass compatible to that of the SU(5) glueball -a blow-up of the comparison in figure 8 for the large x regime is presented in figure 9. For completeness, we also show in figure 8 that the predicted next excited torelon pair, depicted by the black curve, would have a higher mass than the bona fide glueball. Summarizing, the onset of the would-be large-volume glueball takes place in the same range of x-values for both N = 17 and N = 5. Note that x = 4 corresponds to λl = 0.8 for SU(5) and a much smaller physical volume λl = 0.24 for SU (17). We conclude hence that the quantity setting the onset scale for the appearance of the large-volume glueball is x and not the physical volume of the box.
Notice that this opens the door for the possibility of extracting a glueball spectra in the large-N limit at fixed values of x. This is in line with our x-scaling hypothesis translated to the glueball spectrum, which is a stronger statement than volume independence. Taking N to infinity at fixed x implies that the size of the torus goes to zero. This is one of the limits conceived in ref. [28] and called singular by their authors.
In section 4.2.2 we will discuss how to optimize the selection of the large-volume glueball by looking to the overlap of the zero-flux states onto two-torelon and Wilson loop operators. We will present evidence that the onset of the large-volume glueball has taken place for x 4 in all the cases we have analysed.
To close this section, let us mention that the Fibonacci sets are also optimal from the point of view of determining the true glueball mass. This is illustrated in figures 10 and 11 where, as before, we look at the x-dependence of the mass gap compared with Then by definition For large n the value quickly approaches 3− energy, and the minimum is always quite close to the value of the infinite volume glueball. It is tempting to think that this behaviour is not by chance.

Filtering out the glueball from the different torelon pairs
Our previous subsubsection shows, for x > 4 at least, the presence of a state in the spectrum with approximately the same mass as the infinite volume glueball. This happens for most values of N and k in our data set. The problem is that this state might not be the minimum energy state (above the vacuum). The question that we are trying to answer in this section is whether this state has some characteristics that allows to single it out from the remaining states in the low-lying spectra. A priori one would expect that the torelon pair states would couple more to the operators involving a product of the corresponding Polyakov lines. On the contrary, this glueball state would show preference for the Wilson loop operators. With this idea in mind we have analysed the spectra in the vanishing electric flux sector. Our data only allow to determine the mass of the two lowest mass states. For each state of our N = 5,7,11,17 data we computed the projection onto a certain restricted set of operators (the technical aspects of this procedure are explained in the appendix) and selected those states for which the projection onto Wilson loop operators is larger than 0.7. The resulting masses are displayed in figure 12  the relatively small spread centred around a mean value of 0.878. This value matches perfectly with the SU(5) lattice glueball mass from ref. [42] presented in previous plots, which appears as the yellow band in the figure. Our procedure does not work for our N = 13 data points, presumably because of the presence of almost degenerate torelon pair states. The N = 13 data points in the figure are just the value of the mass gap. The mass degeneracy is quite obvious. In units of the string tension, we obtain m/ √ σ = 4.10(4), compatible within errors with the large N glueball mass as estimated in ref. [37]. This provides a confirmation that beyond x ∼ 3.5 the glueball state is present in the spectrum.

Derived consequences
As we did in the previous section with the torelon spectra, we will explain in this subsection the consequences that one can extract from the results presented in the previous one.

x-scaling and volume independence in the e = 0 sector
Our results show a spectra in agreement with expectations. In the range of x values studied there are states with masses approximately given by the sum of two torelon ones, as given by our simple formula. They correspond to torelon pairs. It makes sense that for large N the mass of a torelon pair comes close to the sum of the masses. Interaction energies can be obtained from our data by taking the difference. Our data shows also a state coupled mostly to Wilson loops and with a mass which is close to the one of the infinite volume glueball. This state can be identified with the glueball. It appears in the low-lying spectrum at around a similar value of x ∼ 3 independently of the other parameters, despite the large difference in values of l and N involved. This implies that some sort of x-scaling is actually in play. However, the data does not allow to establish any pattern in the mass differences as a function of N orθ (see figure 12). It might be due to the fact that the differences are small and can be overtaken by all sorts of systematic and statistical uncertainties.
One can criticize the identification of the states as torelon pairs and glueballs. Since they have similar quantum numbers these states might mix and the mixing can affect the value of the masses too. Nonetheless, the mixing is expected to be bigger for states that are nearly degenerate. Most importantly the mixing is suppressed in the large-N limit and this makes the possibility of having a well defined glueball spectrum separated from a torelon pair spectrum viable even when the masses are not very different. Nevertheless, the criticism certainly applies at finite N where the states will be mixed with the torelon pair states which show a strong volume dependence. Still this might induce small corrections for N sufficiently large. What could happen in the large-N limit will be discussed in the next paragraph. Let us first clarify the possible options that occur when taking the large-N limit. The more standard way would be to take the limit at fixed value of λl. This implies x would go to infinity. According to volume independence the glueball spectrum should be independent JHEP07(2018)169 of l and of k. Although at finite l the torelon pair spectrum is still relatively light it is decoupled from the standard glueball spectrum. Our results show that the limit can be taken avoiding the presence of tachyonic instabilities paying the price of choosing the flux k appropriately for each N . This procedure is in line with the way in which one obtains results in 4 dimensions using the twisted Eguchi-Kawai model. Even with wrong choices of the flux the limit can still be taken for values of λl which are sufficiently large to lie beyond the region of instability. This would match nicely with the proposal of Narayanan and Neuberger [18] using zero-flux (periodic boundary conditions). However, using the right fluxes this minimum length can be decreased at will.
A much stronger type of large-N limit occurs when it is taken at fixed value of x, since in that case the size is driven to zero as N grows. This is the singular limit discussed in ref. [28] and also the one involved in x-scaling. Our results show that talking about a glueball spectrum in that limit might make sense at least for x beyond a certain threshold value. Furthermore, the mass gap has a very mild x dependence in that region. Whether this resulting spectra would depend onθ is somewhat more questionable. We saw that the torelon spectra only allows to take the singular limit for a zero-measure set ofθ values, without falling into a phase transition region at intermediate values of x. Still it could be possible that this does not affect the glueball states in question, but this is hard to defend this possibility without any type of solid argument.

Beating factorization for the glueball spectrum
Our paper has also implications in one of the main challenges of non-perturbative large-N gauge theories. This has to do with the difficulty of extracting the glueball spectrum due to factorization. If one computes the correlation of two Wilson loops at different times, the leading term will be the uncorrelated term, while the correlation carrying the signal of the time dependence is suppressed as 1/N 2 . Thus, it could be completely covered by the fluctuations of the uncorrelated part. One way to beat this problem is to consider large Wilson loops. In the confinement regime the expectation value of these large Wilson loop will become rather small, while this size does not in principle affect the coupling of the loop to the glueball. If we make a very naive estimate based on the area law, one might conclude that it would be enough to take the area to grow as log N to make the correlated and uncorrelated pieces of the same order. Even for reduced models in which the loop sizes are forced to be much smaller than N or √ N (which acts as an effective box size) this is acceptable. Our experience in the present case shows that this is indeed possible.

Conclusions
In this paper we have analysed the dependence of the spectrum of 2+1 dimensional Yang-Mills field theory on the volume of space. Ultimately, the goal is to have a complete understanding of the dynamics of this system, which is certainly a prerequisite before achieving the same goal for the more complex 3+1 dimensional theory. It would be nice to be able to test and substantiate analytic approaches as that followed in ref. [36]. Here, we use the spatial volume, and the boundary conditions on it as a probe, which allows us to better

JHEP07(2018)169
understand the different regimes present and the transition among them. Furthermore, we study a wide variety of values of N since we expect a simplification of that dynamics for large N . Our analysis has many implications for several sideline problems like that of volume independence and non-commutative field theories. The main implications of our results have been laid down at the end of the previous sections dealing with non-vanishing and vanishing values of the electric flux. We briefly summarize the main points below.
We have been able to synthesize the evolution of the torelon (non-zero electric flux) energies from the perturbative regime to the confinement regime into a simple formula. This has allowed us to predict the conditions for the avoidance of tachyonic instabilities in the system. This can be maintained at all stages when taking the large N limit, but only for a measure zero set of values of the dimensionless non-commutativity parameterθ. In the zero-electric flux sector we have been able to disentangle in the spectrum torelon pairs from genuine glueballs. The emergence of the latter with a relatively constant mass value takes place at a given value of the effective size lN and hence when volumes are still small enough for the torelon pairs not to become very massive. Thus, despite the complex volume dependence of torelon spectra a largely volume independent glueball state arises.
The data contained in this work represents a considerable effort given the large number of simulations involved. Our results cover a huge range of values of N and also of sizes implying not only an important computational method but also the necessity of dealing with different technical problems in different ranges of parameters. Nonetheless, this work is certainly improvable along various directions. In our opinion, the largest room for improvement is at the larger values of x. This would probably demand a more complete list of operators which could also allow a better determination of the spectra, a larger number of excited states and an exploration of different quantum numbers.

Acknowledgments
We have profitted with many conversations with colleagues during the development of this work. We signal very specially Carmelo Pérez Martín, Fernando Chamizo

A The lattice simulation
In this appendix we explain the methodology used in our lattice simulation.

A.1 Lattice model
We analyse a lattice SU(N ) model in 2+1 dimensions where the spatial components are defined on a torus with twisted boundary conditions and the temporal extent is taken periodic but always large enough to be able to neglect the effects of finite temperature, cf. section A.2.3. The action of the model is: The index n goes over the sites of the L × L × T lattice, with the physical size of the 2-torus given by l = La, with a the lattice spacing. The twist tensor z µν is equal to 1, except at one corner plaquette of each spatial plane where it is: The dimensionless coupling b is defined as b ≡ 1/(aλ L ), with λ L the 't Hooft coupling on the lattice, differing from the continuum λ by lattice artefacts. We will express most lattice computed quantities in units of λ L . Other choices are possible, for instance ref. [37] uses a mean-field improved coupling [43]. For the coarsest lattices that we are using, the difference amounts to about 10% and should go to zero in the continuum limit.

A.2 Determination of the spectrum
The spectrum of the non-zero electric flux states has been extracted directly from the exponential decay of the two-point correlation function of Polyakov loop operators with appropriate winding number and fixed non-zero minimal momentum allowed by the twisted boundary conditions. The Polyakov loops can be constructed in terms of single winding operators represented by the product of APE-smeared [44] link variables and analogously along the y direction. For the results presented in this paper we have taken a fixed number of APE smearing steps: s = 21. A generic Polyakov loop of winding w is given by: P w = Tr(P w 1 x P w 2 y ). 1 In most cases we have only considered Polyakov loops winding around one torus cycle but in a few cases we have also analysed correlation functions of loops of the form Tr(P e x P e y ). In order to improve the overlap onto the lowest mass state for the zero electric flux sector, we use the Generalized Eigenvalue Problem (GEVP) [45,46], that has become a standard tool to compute the spectrum in lattice gauge theories. We will briefly discuss below the particular implementation of the GEVP we have used to determine the glueball masses.

JHEP07(2018)169
The basis of observables used to extract the spectrum in the zero electric flux sector consists of rectangular Wilson loops W (n, m) and moduli of multiwinding spatial Polyakov loops |TrP n | 2 . 2 The latter capture the "torelon-antitorelon" states wrapping around the finite torus, while the former are most useful in the large-x regime where they couple to the large-volume glueballs.
We use three fixed levels of APE smearing [44] with 7, 14, 21 iterations (as well as the unsmeared operators, which are however too contaminated by excited states to be used in practice). Small Wilson loops can be dominated by UV effects, therefore to maximize the overlap to the physical states, we include in the basis the rectangular Wilson loops of large extent, ranging up to approximately 20-40 lattice units, depending on the value of the coupling. This can be seen as an alternative to the blocking procedure [47,48].
Note however, that this set-up is somewhat limited -it is best suited for the extraction of the lowest mass glueball, while the lack of decorated loops is expected to give at best a qualitative description of the excited glueball states.
Given the observables O i (t), the matrix of two-point connected correlation functions is given by: A typical size of the correlation matrix we used for the glueball mass extraction is between 15 and 30 operators.
To obtain a good overlap with the lowest energy states we perform the GEVP on the correlation-function matrix [45,46]: where the times are fixed to t 0 = a, t 1 = t 0 + a and α labels the energy states. The obtained eigenvectors v (α) (t 0 , t 1 ) are then used to find the improved basis of correlation functions [45]: for each value of t and then the diagonal valuesC (αα) (t) are used to find plateaux ranges and subsequently fit to the selected ranges. 3 Note that this is different from the approach of ref. [49] which scales t 0 proportionally to t (and therefore has better theoretical convergence properties). However, in our case increasing t 0 results in a very rapid growth of noise as a function of t and obtaining reliable GEVP plateaux with this method would require a huge increase of statistics.

A.2.1 Finding plateau ranges for the mass extraction
The masses of non-zero electric flux states are extracted in the usual way by looking for plateaux in the effective masses obtained from Polyakov-loop correlators. The plateau range is fixed by first determining the value of t where the effective mass is minimum among points with relative statistical error smaller than 4%. The fitting range around JHEP07(2018)169 that time is then adjusted to keep a χ 2 per degree of freedom smaller than 1. In most cases we take characteristic plateau ranges of 5 to 8 points with resulting χ 2 per degree of freedom below 0.5. The error is determined from the dispersion in fit parameters obtained by repeating this fitting procedure in jackknife bins.
The correlators used for the mass determination for vanishing electric flux are more noisy, also the energies are in some cases very closely spaced (level crossings), therefore a more systematic way of choosing the plateau ranges was introduced. We use singleexponential correlated fits, and the time ranges are selected using a variation of "fit histogram" method from ref. [50]. For each t min , t max we calculate the quantity measuring the "fit quality": where Γ is the cumulative distribution function of the appropriate χ 2 -distribution. Next we choose the fit interval which maximizes Q. We also analyse whether the mass gap fit result is stable among the fits with highest fit quality. In some cases, statistics had to be increased to fulfill that condition. On the other hand, for the next excited state this condition cannot always be fulfilled without a significant increase of the statistics, therefore we calculate the standard deviation of ten fit results with the highest Q and quote that as an estimate of the systematic error due to the choice of the fit range.

A.2.2 Selection of glueball states
As discussed in section 4.2.2 we have to face the problem that the glueball state is not always the lowest state over the vacuum in the zero electric flux sector due to the presence of light torelon-pair states. The identification of the glueball has required a procedure to determine the operator content of each state. We expect that torelon-pairs couple mostly to Polyakov loop operators while the glueball state would project more onto Wilson loop operators. In what follows we will briefly discuss how to implement this selection procedure.
Suppose we have a collection of normalized states |α , which have been identified as eigenstates of the Hamiltonian, and a set of operators O i that produce these states when acting on the vacuum. Now we can define a measure of the overlap of a given operator with a given eigenstate of the Hamiltonian as follows: The number in question is positive and smaller than 1. When this number becomes large we can say that we have a large overlap between the state (for example the lowest energy state) and the operator O i . The idea can be also applied in the case in which there are several operators having similar properties, like for example using different smearing levels. The philosophy is to gather all the operators of the same type into a given set. For example, suppose we have a set G as follows

JHEP07(2018)169
We can associate a vector space V G to the set G. This is given by the p dimensional subspace of the total Hilbert space generated by the vectors O i |0 . Now we can define a measure of the overlap of these set of operators with a given eigenstate of the Hamiltonian as follows: where the maximized quantity is the absolute value of the Hilbert space scalar product of the vector u with the state |α . Notice for a single operator it coincides with the previous formula.
The two lowest lying eigenstates of the Hamiltonian and their coupling to operators can be obtained from the results of the GEVP method described in section A.2. This approach seems to work remarkably well when working with a reduced family of operators including only one Wilson loop operator and one operator constructed from Polyakov loops for each torelon pair.

A.2.3 Suppression of finite T contributions in glueball correlation functions
One problem necessary to overcome when dealing with glueball correlation functions is the effect of the wrap-around states, coming from the finite time extent of the lattice. In perturbation theory one can show that the leading effect is a single wrap-around gluon, giving rise to a constant term in the correlators (even in the connected ones) proportional to exp(−m t T ), where m t is the mass of the lightest torelon state. While methods similar to the ones used in ref. [51], like midpoint subtraction of the correlators can be used to mitigate this problem, we find it most efficient to always use lattices with large enough T , that the contribution from those states is negligible.