Closed flux tubes and their string description in D=3+1 SU(N) gauge theories

We calculate the energy spectrum of a confining flux tube that is closed around a spatial torus, as a function of its length l. We do so for various SU(N) gauge theories in 3+1 dimensions, and for various values of spin, parity and longitudinal momentum. We are able to present usefully accurate results for about 20 of the lightest such states, for a range of l that begins close to the (finite volume) deconfining phase transition, and extends up to l.sqrt(K)~6 (where K is the string tension). We find that most of these low-lying states are well described by the spectrum of the Nambu-Goto free string theory in flat space-time. Remarkably, this is so not only at the larger values of l, where the gap between the ground state energy and the low-lying excitations becomes small compared to the mass gap, but also down to much shorter lengths where these excitation energies become large compared to sqrt(K), the flux-tube no longer `looks' anything like a thin string, and an expansion of the effective string action in powers of 1/l no longer converges. All this is for flux in the fundamental representation. We also calculate the k=2 (anti)symmetric ground states and these show larger corrections at small l. So far all this closely resembles our earlier findings in 2+1 dimensions. However, and in contrast to the situation in D=2+1, we also find that there are some states, with J,P=0,- quantum numbers, that show large deviations from the Nambu-Goto spectrum. We investigate the possibility that (some of) these states may encode the massive modes associated with the internal structure of the flux tube, and we discuss how the precocious free string behaviour of most states constrains the effective string action, on which much interesting theoretical progress has recently been made.


Introduction
The idea that the strong interactions may be described by a string theory is even older than Quantum Chromodynamics, e.g. [1], and the idea that gauge theories (at least in the planar limit) may have such a description is only a little younger [2]. The more recent and radical version of this idea is gauge-string duality [3]. To learn something about this string theory it is natural to start by focusing upon any degrees of freedom that are manifestly string-like and, in linearly confining theories such as SU(N) gauge theories in D = 2 + 1 and D = 3 + 1, these are long confining flux tubes. One can ask what effective string theory describes their dynamics. Recently there has been substantial analytic progress towards answering this (old) question which, roughly speaking, tells us that the dynamics governing very long flux tubes is, to a certain approximation in powers of 1/l, that of a Nambu-Goto free bosonic string theory [4,5,6,7,8]. (We shall review this in more detail below.) At the same time our numerical calculations of the the low-lying excitations of closed flux tubes in D = 2 + 1 [9] have shown that all the calculated energies are remarkably close to those of the free bosonic string theory even when the flux tube length l is not much greater than its width so that an expansion in powers of 1/σl 2 (where σ is the string tension) no longer converges. Moreover the fact that we have not observed additional massive string modes associated with the non-trivial structure of the flux-tube, suggests that these interact very weakly with the usual massless stringy modes. These complementary numerical and analytic results strongly suggest that the effective string action in D = 2 + 1 SU(N) gauge theories is some small perturbation about the Nambu-Goto action (in flat space-time), so that the latter, rather than the traditional classical string configuration, should be the appropriate starting point for any analytic investigation.
In this paper we extend our calculations to 3 + 1 dimensions. As in 2 + 1 dimensions, we close our flux tubes around a spatial torus of length l so as to stabilise them at a chosen length. Such flux tubes have non-trivial rotational properties in 3 spatial dimensions and hence a richer spectrum than in 2 spatial dimensions. Thus, despite the fact that our calculations turn out to be less accurate than in D = 2 + 1, we are able to obtain usefully accurate results on a substantial number of low-lying eigenstates. What we shall show is that most of these states are well-described by the Nambu-Goto model down to very short flux tube lengths, just as in D = 2 + 1. However we shall also find that a few states, with specific quantum numbers, are very far from the Nambu-Goto prediction, leaving it unclear whether they approach the free bosonic string value as l increases or not. We discuss in some detail the possibility that this is a signal of the massive modes that had eluded us in our earlier D = 2 + 1 calculations.
We shall begin, in Section 2, by summarising the theoretical background and briefly describing where the analytic calculations of effective string theories now stand. In Section 3 we turn to the technicalities of our numerical calculation. We describe the lattice Monte Carlo set-up, how to calculate the eigenspectrum of closed flux tubes, the specific lattice operators that we use, and the quantum numbers of the states. We finish Section 3 with some brief remarks designed to place our calculations in the context of earlier lattice work. We then turn to our results in Section 4. We begin, in Section 4.1, with an analysis of how the (absolute) ground state energy of such a closed flux tube varies with its length l. We compare our values of the energy to the free string prediction and fit the observed (small) deviation to some theoretically motivated correction terms. In Section 4.2 we consider non-zero momentum along the flux tube and calculate the energies of the lightest states in each such sector. In Section 4.3 we calculate a number of excited state energies. It is here that we shall identify and discuss a number of states with energies very far from the free string prediction. In Section 4.4 we briefly deviate from the main theme of this paper and consider the properties of the k = 2 (anti)symmetric ground states. Finally, we discuss the implication of these results in Section 4.5 Section 5 contains our conclusions.
An important aim of our work is to provide 'data' that can be used to inform theoretical analyses of the effective string action. While we do attempt to provide comparisons with those analyses available to us at the time of writing, this is an area where there is rapid current progress, and other comparisons may soon be of interest [10]. So to maximise the usefulness of our work, we have deliberatively provided in an Appendix a comprehensive tabulation of our numerically determined eigenenergies, in a form that is intended to be readily usable.
We remark that a brief summary of some of this work has appeared previously in [11].

Effective string theory
We are in the confining phase of an SU(N) gauge theory on a 4-torus, with partition function Z. We label the coordinates by (x,y,z,t). Since space-time is Euclidean we may in fact choose any co-ordinate as 'time' and write Z as a sum over eigenstates of the Hamiltonian defined on the orthogonal 3-dimensional 'space'. The same is true for partition functions including sources, which we define below, and the resulting 'dual' descriptions play an important role both in constraining the effective string theory and in the set-up of our numerical calculations.

General considerations
Consider a static fundamental source at x = (0, y, z), with a conjugate source at x = (r, y, z). Let τ be the the Euclidean time extent of our 4-torus (with all other tori large). The partition function for the system with these two sources is Z ss (r, τ ) = dA l † (x = r, y, z)l(x = 0, y, z) exp{−S[A]} (1) where l( x) is a Polyakov loop. (The traced path-ordered exponential of the gauge potential along a path that encircles the t-torus at x = (x, y, z), which is the phase factor arising from the minimal coupling, j µ A µ = j 0 A 0 , of our static source.) If r ≫ 1/ √ σ there will be a flux tube between the sources which, as it evolves in time, sweeps out a sheet bounded by the periodic sources. This sheet clearly has the topology of a cylinder. The partition function can be written as a sum over energy eigenstates 1 Z Z ss (r, τ ) = n e −En(r)τ (2) where E n (r) is an energy of two sources separated by r. The states are (excited) flux tubes that begin and end on a source and which evolves around the t-torus. There is another way to interpret Z ss (r, τ ). Take x to label the 'time', so l is now a Wilson line that winds around what is now a spatial torus of length τ . What Z ss (r, τ ) represents, in this point of view, is a correlation function whose intermediate states consist of flux tubes that wind around this same 'spatial' torus of length τ and propagate the distance r between the two Wilson lines. This partition function can therefore be written as a sum over these energy eigenstates 1 Z Z ss (r, τ ) = n,p c n (p, τ )e −Ẽn(p,τ )r whereẼ n is an energy of an (excited) flux tube that winds around a spatial torus of length τ . E n and E n are different functions because the flux tubes have different boundary conditions (although later on we shall use E n for both). Since the operators l are localised at y, z, the winding flux tubes are delocalised in transverse momentum, and we have made explicit in eqn(3) the sum over that momentum. The c n are the wave-function factors for the overlap of a state |n, p on the state obtained by applying the operator l to the vacuum state |vac : c n = | vac|l † |n, p | 2 . Suppose that we have an effective string theory, governed by an effective action S ef f , which reproduces the long distance physics of flux tubes. Then the string partition function, taken over the r × τ cylinder considered above, should reproduce the flux tube contribution to the field-theoretic partition function in eqn ( 1): where we integrate over all surfaces S spanning the cylinder. From eqn (4) we see that Z cyl (r, τ ) can be written as a sum over open flux tube states as in eqn (2). This is nothing but a Laplace transform in τ of Z cyl (r, τ ). So if we have a candidate string action, S ef f [S], we can perform this Laplace transform and extract a prediction for the the open flux tube energies E n (r). We can also express Z cyl (r, τ ) as a sum over closed flux tube states using eqn (3). This is a Laplace transform in r of Z cyl (r, τ ), and so a candidate S ef f [S] will imply a prediction for the the closed flux tube energies,Ẽ n (p, τ ). However in this case the energies are not all independent: the relationship between energy eigenstates that differ only by their transverse momentum p is determined by Lorentz invariance. (In fact one can show that doing the integral over p turns eqn(3) into a sum over Bessel functions [4,12].) There is no reason that, for some arbitrary choice of S ef f [S], the Laplace transform of Z cyl (r, τ ) in r should reproduce an energy spectrum that satisfies this relationship. As first realised in [4], this provides a powerful constraint on the permitted form of the effective string action. It was then observed in [7] that the constraints arising from the above open-closed duality associated with a cylinder can be extended to other geometries that are less obvious from a field-theoretic perspective. In particular one can extend the analysis to an r × τ torus [7]. A world sheet spanning such a torus corresponds to a closed flux tube of length r propagating over a time τ or one of length τ propagating over a time r. Thus we have a closed-closed string duality: dSe −S ef f [S] = n,p e −Ẽn(p,τ )r = n,p e −Ẽn(p,r)τ (5) which turns out to provide useful new constraint on S ef f [S] [7]. It may perhaps be that we have not exhausted all the possibilities and that other boundary conditions and set-ups may provide further useful constraints. This general discussion needs some qualifications. Our effective string theory is in 2 + 1 or 3 + 1 dimensions and far from the critical dimension where we can consistently define a string theory. However the resulting anomalies, which show up in different ways depending on how one 'gauge-fixes' the diffeomorphism invariance in one's calculation, typically die off at long distances [13] (although the implications of this are not unambiguous), and when one considers smooth fluctuations around a long string [14]. Thus it can make sense, at least technically, to consider a string path integral over single large surfaces, in an effective string theory approach outside the critical dimension [14]. Since this represents the world sheet swept out by a single long fluctuating string, this means that an effective string theory is limited to describing the dynamics of a single long fluctuating flux tube. This is clearly an important physical limitation. In reality, a sufficiently excited flux tube can decay into a flux tube of lower energy and a glueball, and such states inevitably appear in the field-theoretic sum over states in eqn (2) and eqn (3). In the string picture a glueball is a contractible closed loop of string whose length is O(1/ √ σ) (for light glueballs). There is no guarantee that an effective string theory can consistently describe such extra small surfaces. One can partially circumvent this by only considering low-lying string states which are too light to decay: where m G is the energy of the lightest glueball. (Or of the lightest state composed of a flux loop and its conjugate, if that is lighter.) However even such states will be affected by mixing through virtual glueball emission, which corresponds to small handles on our large surfaceagain something that would be problematic for the string theory. There is of course a limit in which mixings and decays do vanish, and that is the N → ∞ limit. In the SU(∞) theory one can consistently discuss a partition function over states containing a single flux-tube, as in eqns (2,3), and ask if it is given by an integration over single surfaces with some effective string action. It is then also plausible that as we move continuously away from that limit, to finite N, the corrections will be small, so that it makes sense to analyse all SU(N) gauge theories within such an effective string framework [7]. The fact that there is very little dependence on N in the low-lying flux tube spectrum, both in D = 2 + 1 [9] and in D = 3 + 1 (this paper), lends credibility to this assumption.
Consider a flux tube that winds around a spatial torus of length l. The excited states of this flux tube are obtained from the ground state, E 0 (l), by exciting some of the modes living on the flux tube. If the excited mode is massive we expect the energy to be If the excited mode is massless, we would expect the extra energy to equal the (longitudinal) momentum which, for bosons, is quantised to be p = 2πk/l; k = ±1, ... by periodicity. So we then expect So once our flux tube is long enough, with π/l ≪ √ σ, its lightest excitations are given solely by the excitations of the massless modes. One should therefore be able to construct an effective string theory solely out of its massless modes if what we want is a description of the low-lying spectrum of long flux tubes.
In general we expect modes to be massless for symmetry reasons. In the case of a flux tube there are D − 2 obvious massless modes. These are the Goldstone modes that arise from the fact that once we have specified the location of our flux tube, we have broken spontaneously the translation invariance in the D −2 directions transverse to the flux tube. So one usually starts with an effective bosonic string theory involving just these Goldstone fields. By comparing with the numerical spectrum at large-l one can see if the presence of other, less obvious, massless modes is indicated.
In practice most numerical work has involved the ground state energy, E 0 (l). The massless modes, when quantised, contribute O(1/l) zero-point energies to E 0 (l). Summing over frequencies, gives the 'universal Lüscher correction' [15] where c is proportional to the central charge of the string theory and will differ from the bosonic value shown if other massless modes are present. For example, one expects it to be zero in the case of N = 1 SUSY 1 . (While this is for closed flux tubes, there is a similar expression for open flux tubes.) There have been increasingly accurate lattice determinations of c over the last 25 years that leave little doubt that in pure gauge theories the only massless modes on the flux tube are indeed the transverse Goldstone modes.
To proceed with explicit computations, one needs to fix convenient coordinates to describe the surface in the path integral. This 'gauge-fixing' of the fundamental diffeomorphism invariance typically makes the constraints that follow from this string symmetry less obvious. Here we follow the 'static gauge' approach of [15,4,7] and do not discuss the general effective string approach of [14], which has been used in [5,6,8] to obtain comparable results (in 'conformal gauge').
Suppose we are integrating over the surfaces of the cylinder discussed above. There is a minimal surface which can be parameterised by x ∈ [0, r] and t ∈ [0, τ ]. Other surfaces are specified by a transverse displacement vector h(x, t) that has two components in the (y, z) directions (but only one in D = 2 + 1). This way of parameterising a surface is called 'static gauge'. We can now write the effective string action in terms of this field h and the integral over surfaces becomes an integral over h(x, t) at each value of (x, t) ∈ [0, r] × [0, τ ]. Since the field h is an integration variable in (0, ∞), we can take it to be dimensionless. Moreover, since the action cannot depend on the position of the flux tube (translation invariance), it cannot depend on h but can only depend on ∂ α h where α = x, t. That is to say, schematically, (10) and we can perform a derivative expansion of S ef f in powers of derivatives of h: (very) schematically where the derivatives are with respect to x and t and indices are appropriately contracted. The coefficients c n have dimensions [length] (2n−2) to keep the terms dimensionless. So we can expect that for the long wavelength fluctuations of a long string, such a higher order term will make a contribution to the energy that is down by O(1/(σl 2 ) n ) compared to the leading linear piece and so the importance of these terms is naturally ordered by the number of derivatives. All this is entirely analogous to the familiar way chiral Lagrangians depend on their Goldstone fields. And just as the applicability of chiral Lagrangians is typically bounded by the energy scale of the lowest resonance, this derivative expansion is designed to capture the physics on energy scales smaller than the O( √ σ) dynamical mass scale. For a surface with a boundary, like a cylinder, there can also be operators localised on the boundary, which contribute odd powers of 1/r, and these we have ignored. Since we are interested in the spectrum of winding strings, that can be obtained from the torus geometry which has no boundaries, it is safe for us to do so.
Note that our chosen 'static-gauge' parameterisation does not work for general surfaces. To describe a string with an 'overhang' or any kind of 'back-tracking', the field h(x, t) would have to multivalued, which is something the standard treatments do not allow. That is to say, we arbitrarily exclude such rough surfaces from the path integral. For a flux tube, its finite width provides a physical lower distance cutoff on such fluctuations: any overhang that is within a distance 1/ √ σ will in effect be a fluctuation in the intrinsic width of the flux tube i.e. a massive mode excitation. Any backtracking/overhang that is larger will increase the length by ∆l > 1/ √ σ and hence the energy by ∆E ∼ σ∆l > √ σ. In both cases the associated excitation energies will be much greater than the O(1/l) gap to the stringy modes, once l is large enough. Thus this should not be a significant issue for the long wavelength massless oscillations we have discussed above. But it needs to be addressed in any analytic treatment that aims to be more ambitious than that.

Gaussian action
The first non-trivial term in our effective string action in eqn(11) is the Gaussian piece: S gauss It has the fewest derivatives and so will provide the leading large l correction to the linear piece of the string energy. Since the action is Gaussian one can evaluate Z cyl exactly, giving Z cyl (r, τ ) = e −σrτ |η(q)| −(D−2) : q = e −πτ /r (13) in terms of the Dedekind eta function (See e.g. [4] whose notation we borrow.) If we expand the product in eqn (14) we have a sum of powers of q, which, using q = e −πτ /r , becomes a sum of exponentials in τ . This is precisely of the form given in eqn (2). If we match this sum to eqn(2), we obtain the prediction : Gaussian (15) for the energy levels (as well as a prediction for their degeneracies). This is the exact result for the energy levels of a string with ends fixed to static sources, a distance r apart, for the case of a Gaussian S ef f . We note that the excitation energies display an O(1/r) gap as expected from eqn (8).
To obtain a corresponding prediction for closed strings, we use the Dedekind eta function's modular invariance under q = exp{−πτ /r} →q = exp{−4πr/τ }. This allows us to rewrite the expression for Z cyl in eqn(13) as a sum of exponentials in r rather than in τ . However this sum turns out not to be precisely of the form shown in eqn (3), where the eigenstates fall into subsets that are related by Lorentz invariance (so that, for example, their energies satisfy the usual energy-momentum dispersion relation). Thus a Gaussian S ef f does not encode the open-closed string duality exactly and cannot be considered as a candidate for an exact description of strings on a cylinder. However if, instead, we view the Gaussian S ef f as an approximation for long closed strings, we can expand the momentum dependence in eqn(3) in powers of p/στ , and then match to this approximate form. We thus obtain a prediction for the closed string energy levels, and also for the overlaps c n [4].
The O(1/r) correction to the leading linear term in E 0 (r) in eqn (15) is the famous Lüscher correction [15] for a flux tube with ends fixed on static sources. Physically it arises from the regularised sum of the zero-point energies of all the quantised oscillators on the string. It depends only on the long wavelength massless modes and so is universal: any bosonic string theory in which the only massless modes are the transverse oscillations will have precisely this leading correction. The same applies to the O(1/τ ) correction to the leading linear term iñ E 0 (τ ) in eqn (16).
Although the above results for E n (r) are obtained in the Gaussian approximation to S ef f [h], this approximation becomes exact as r → ∞, and these predictions for the leading O(1/r) correction are exact and universal. And similarly forẼ n (τ ) as τ → ∞.

Nambu-Goto action
A string theory that is simple enough to have a calculable energy spectrum [16] is the freestring Nambu-Goto model in flat space time: where we integrate over all surfaces, with the action proportional to the invariant area. A free string theory is not necessarily unrealistic as an effective string theory: after all, we recall that at N = ∞ confining flux tubes do not interact. Moreover we have found that the flux tube spectrum in D = 2 + 1 is remarkably close to the Nambu-Goto spectrum [9]. For these reasons we shall use this model as our main point of comparison. Consider a string winding once around the x-torus. Working in static gauge, perform a Fourier decomposition of h(x, t). Upon quantisation the coefficients become creation operators for massless 'phonons' on the string with momenta ±2πk/l along the string and energy 2πk/l. (The k = 0 mode is not included since it corresponds to a shift to a different transverse position of the whole string i.e. to another vacuum of the spontaneously broken symmetry.) There are 2 transverse directions, so h is a 2-component vector and the phonons carry angular momentum of ±1. We call positive momenta left-moving (L) and the negative ones rightmoving (R). Let n ± L(R) (k) be the number of left(right) moving phonons of momentum 2πk/l and angular momentum ±1. Define the total energy and momentum of the left(right) moving phonons as 2πN L(R) /l, then: If p = 2πq/l is the total longitudinal momentum of the string then, since the phonons provide that momentum, we must have The angular momentum around the string is One can now write down the expression for the energy levels of the Nambu-Goto string: where the degeneracies corresponding to particular values of N L and N R will depend on the number of different states that can be formed from combinations of values of n ± L and n ± R that give the same values of N L + N R and q in eqn (18). In discussing the states, we shall often refer to the left and right moving phonon creation operators of (absolute) momentum 2πk/l and angular momentum ±1 as a ± k and a ± −k respectively, and the unexcited string ground state as |0 .
Let us specialise, for simplicity, to q = 0, i.e. N L = N R = n, and make some general comments.
• The energy E n (l) can be expanded for large l in inverse powers of 1/σl 2 : We note that the first correction to the linear piece is exactly as in eqn (16) -as it should be if it is 'universal'.
• The ground state energy becomes tachyonic at small l: One can regard this as the Hagedorn/deconfining transition in the Nambu-Goto model. We note that the length at which this transition occurs, l H √ σ = 2π/3 ≃ 1.45, is not accessible in pure gauge theories at large N, because of a first-order 'deconfining' transition that occurs at l c √ σ = √ σ/T c ≃ 1.65 [17].
• The expansion of the square root expression for the energy E n in eqn (22) only converges for σl 2 > 8πn (ignoring the negligible D − 2 term). So the more excited is the state, the larger is the value of l to which we must go to be able to use such an expansion. The derivative expansion of the action is not uniform in frequency: it is formal.
• One can show (see Appendix C of [4]) that the Nambu-Goto model satisfies open-closed duality exactly. This is in contrast to the Gaussian string action. We assume (although we have been unable to locate the proof in the literature) that it also satisfies the closed-closed duality associated with a world sheet on a 2-torus. Thus, if these dualities are our only constraints, the Nambu-Goto model is a viable candidate for providing a consistent effective string action. An important consequence of this is that where imposing these constraints allows us to completely fix the expansion coefficients of E n (l) up to some order in 1/l, these coefficients will have to be precisely the same as those obtained by expanding the Nambu-Goto expression in eqn (22) to that order. (Or the corresponding expression for strings with fixed ends.) Finally we note that a comparison between our lattice results for closed flux-tubes and the Nambu-Goto spectrum of eqn (21) can be seen as analogous to comparing glueball spectra to Regge trajectories [18]. Specifically, if we focus on excited states for which the second term in eqn (21) is dominant, and if for simplicity we take q = 0, or N R = N L ≡ N, then we can write This linear relation between the energy squared and the quantum number of the state is very reminiscent of the M 2 ∼ J + n Regge relation that one might expect for glueballs (the 'Pomeron' trajectory) and indeed one immediately sees from eqns (18)(19)(20) that the maximum value of J increases linearly with the value of N.

Recent theoretical progress
The seminal work in analysing flux tubes in a string description in static gauge [15] (reviewed above) and the later more general Polchinski-Strominger approach using conformal gauge [14] (not reviewed here) led to an understanding of the universality of the leading O(1/l) Lüscher correction to the linear growth of the flux tube energy. Until recently there was, however, very little further analytic progress along these lines. The situation changed in 2004 when major progress took place within both approaches. 1) In [4] it was shown that the open-closed duality (discussed above) could be used to provide useful constraints on the higher order terms in the expansion of the effective string action. In particular it was shown that in D = 2 + 1 the next, O(1/l 3 ), term is also universal and that the coefficient is precisely what you get by expanding the Nambu-Goto square-root expression to that order (as must be the case). In D = 3 + 1 the coefficient is not fixed but there is a relationship predicted between the coefficients of the two terms in the effective action that contribute at that order. 2) Simultaneously, the next O(1/l 3 ) term was calculated within the Polchinski-Strominger framework in [5] (and also, later, in [6]). The same conclusion was reached as in [4] for D = 2 + 1, but a stronger conclusion was obtained in D = 3 + 1, where the O(1/l 3 ) term in the action was shown to be universal (and equal to the value in the Nambu-Goto expansion).
In 2009 further substantial progress was achieved [7] using the static gauge approach and including the new constraints that arise from the less obvious closed-closed (torus) duality. In D = 2 + 1 the terms up to and including O(1/l 5 ) have been shown to be universal (and equal to the Nambu-Goto values) [7]. In D = 3 + 1 the O(1/l 3 ) contribution has been shown to be universal [7], in agreement with [5], and the O(1/l 5 ) contribution to the string ground state (and partition function) has also been shown to be universal. Moreover, during the writing of this paper, some papers have appeared [8] extending the Drummond-Polchinski-Strominger approach [14,5] and claiming that the terms up to O(1/l 5 ) are universal in all the eigenstates.
It is now clear that the effective string theory of long flux tubes is that of the Nambu-Goto free string theory to quite high order in the derivative expansion.

Lattice methodology
In Section 3.1 we outline the lattice framework of our calculations. In Section 3.2 we discuss in some detail how we extract the energies of (excited) flux tubes. For the reader who wishes to skip this sub-section, the important message for interpreting our subsequent results is that the larger the energy of the flux tube, whether because it is longer or because it is more highly excited, the larger are the errors. (Particularly the systematic errors that can be hard to estimate.) In Section 3.3 we describe our labelling of the quantum numbers of the flux tubes and tabulate the corresponding Nambu-Goto description -all of which is important for understanding our later results. We present the details of our lattice operators in Section 3.4, which the general reader may again wish to omit. Finally in Section 3.5 we provide a very brief and incomplete sketch of relevant earlier work, to provide some context for our calculations.

Setup
Our numerical calculations are entirely conventional. In order to remove the troublesome phase factor in the Minkowski path integral, we perform our calculations in the Euclidean theory, and to make the problem finite we work on a hypercubic lattice on a finite hypertorus. Schematically, where we express the discrete space-time points as x = na or x µ = n µ a, where {n µ } is a D-tuple of integers, a is the lattice spacing, and φ L is a dimensionless lattice field variable, with action S L , chosen so that where dim(φ) is the length dimension of the field φ. Because the theory has a finite mass gap, the leading finite size corrections are exponentially small, and because it is asymptotically free one can show that the leading finite lattice spacing corrections for the action we shall use are O(a 2 ). The gauge degrees of freedom are SU(N) group elements, U l , that are assigned to the forward going links, l, of the lattice. (With U † l for backward-going links.) For a link that joins the site x to the site x+aμ, the link matrix U µ (x) will transform as U µ (x) −→ V (x)U µ (x)V † (x+ aμ) under a gauge transformation V (x). To construct a gauge-invariant action we note that the trace of the product of link-matrices around any closed path c, Tr l∈c U l , is gauge invariant. So the simplest choice for the lattice gauge action is to use the path that is an elementary square on the lattice, called a plaquette: where U p is the path-ordered product of link matrices around the plaquette p. The p ensures that the action is translation and rotation invariant. Taking the real part of the trace ensures that it has C = P = +. So our lattice path integral is where β is a constant whose value will determine the lattice spacing a. Since we see that β = 2N/g 2 a→0 −→ ∞. In D = 2 + 1 where g 2 has dimensions of mass, the dimensionless bare lattice coupling is ag 2 and β = 2N/ag 2 . Since a smooth large N limit requires keeping g 2 N fixed, this means that we must vary β ∝ N 2 in order to keep a approximately fixed as we vary N.
To calculate the expectation value of some functional Φ[U] of the gauge fields, we generate a set of n g gauge fields {U I }; I = 1, ..., n g distributed with the Boltzmann-like action factor included in the probability measure, i.e. dP ∝ l dU l exp{−βS[U]}. We use a Cabbibo-Marinari algorithm applied to the N(N − 1)/2 SU(2) subgroups of the SU(N) matrix, with a mix of standard heat-bath and over-relaxation steps. We can now calculate the expectation value of some functional Φ[U] of the gauge fields as follows: where the last term is the statistical error.

Calculating energies of closed flux tubes
We begin by recalling the standard decomposition of a Euclidean correlator of some operator φ(t) in terms of the energy eigenstates of the Hamiltonian H: where the energies are ordered, E i+1 ≥ E i , with E 0 that of the ground state. The only states that contribute are those that have overlaps c j = vac|φ † |j = 0, so we need to match the quantum numbers of the operator φ to those of the state we are interested in, which here is a flux tube winding around the x-torus, of length l = l x = La. The simplest such operator is the elementary Polyakov loop: where we take a product of the link matrices in the x-direction once around the x-torus, and the sum over translations projects onto zero transverse momentum (p y , p z ) = (0, 0). The operator is invariant under translations in x and therefore has zero longitudinal momentum p x = 0. It is invariant under rotations about its axis and so has J = 0. It is clearly invariant under parity (accompanied where necessary by charge conjugation C to reverse the flux). Clearly we will need to use deformed versions of the Polyakov loop to access non-trivial quantum numbers for the flux tube (as discussed in detail later on). An operator like the Polyakov loop which winds once around a torus has zero overlap onto states such as glueballs that are localised and which are described by operators constructed around closed loops that can be continuously contracted to a point. This orthogonality is enforced by a center symmetry that corresponds to SU(N) gauge transformations that are periodic -on the torus -up to an element z ∈ Z N . These leave any contractible loop, l c , invariant but not ones with unit winding, l p → zl p . Thus l † c l p = z l † c l p = 0 but only as long as this center symmetry is not spontaneously broken. As we decrease l there is a critical length l c = 1/T c , where T c is the usual deconfining temperature, below which this center symmetry is spontaneously broken and we are in a 'deconfining' phase in which we no longer have confining flux tubes around the x-torus. (Although large Wilson loops orthogonal to x continue to display a confining area law.) In this paper we always restrict l to be > l c .
So, having chosen a suitable operator φ, we calculate its correlator using eqn(30) with Φ = φ † (t)φ(0). This will produce an estimate of φ † (t)φ(0) with a finite statistical error. To extract a value of aE 0 using eqn (31), we need to have significant evidence for the exponential behaviour ∝ e −aE 0 nt , over some range of t = an t and this range needs to begin at small enough n t that the decreasing exponential is still clearly visible above the statistical errors. (In a pure gauge theory one can show that the error tends to a non-zero value at large t, so that the error/signal ratio increases exponentially with t.) This is obviously harder to achieve for larger E 0 , so the corresponding systematic error will be larger for heavier states. However, even for the lightest state we need the normalised overlap to be close to unity, |c 0 | 2 ≃ 1, to obtain an accurate energy estimate: our operator needs to be a good wavefunctional for the state whose energy we are interested in, so that its correlator is dominated by this state even at small n t . For the ground state of the flux tube this can be achieved by using the Polyakov loop operators in eqn(32), but with 'smeared' [19] and 'blocked' [20] link matrices replacing the elementary ones. These produce winding operators that are smooth on various transverse size scales including ones that are comparable to that of the physical flux tube, and, not surprisingly, this works well for the ground state which presumably has a smooth wavefunctional. Starting with this basis of operators, {φ i ; i = 1, ..., n}, we consider the vector space V φ generated by them, calculating all the correlation functions φ † i (t)φ j (0) , and determine the linear combination that is the 'best' operator ψ 0 by a variational criterion where t 0 is some convenient value of t. Then ψ 0 is our best variational estimate for the true eigenfunctional of the ground state. We now use the correlator ψ 0 † (t)ψ 0 (0) in eqn(31) to obtain our best estimate of the ground state energy. This generalises in an obvious way to calculating excited state energies. One constructs from V φ the vector space orthogonal to ψ 0 , repeats the above within this reduced vector space, and obtains ψ 1 which is our best variational estimate for the true eigenfunctional of the first excited state. And so on.
For the ground state the main systematic error is to overestimate aE 0 by extracting it at a value of t where there is still a significant contribution from excited states. This becomes an increasing problem as we increase l and hence aE 0 (l), and so decrease the range of t where the statistical errors are small compared to the 'signal'. At the same time the gap to the excited states decreases ∝ 1/l, increasing their relative contribution to the correlator. This is the main obstacle to accurate calculations of flux tube energies for very large l and can obviously undermine an attempt to calculate corrections to the large-l behaviour of E 0 (l). We therefore restrict our calculations of the ground state energy to values of l where we estimate this systematic error to be smaller than our quoted statistical errors. For excited states this systematic error is less easily avoided because aE i (l) is often quite large for all values of l. One therefore needs to be aware of this in any comparison between our results and the predictions of an effective string model in this paper. In addition there is now an additional error that arises from the fact that our best variational wavefunctional, ψ i , may contain some admixture of lighter eigenstates. At sufficiently large t these will dominate and give an energy that is lower than E i . Note, however, that if the state is the lightest one for some given quantum numbers, then this will not occur and we will often refer to these states as 'ground states', at least when there is no risk of confusion with the absolute ground state. Such ground states will be of particular interest to us because their calculation is typically much more reliable.
Both the systematic errors described above can be controlled if ψ i is a very good wavefunctional for the corresponding excited state. In that case there will be a useful range t ∈ [t 1 , t 2 ] where the desired state dominates because the small admixture of higher excited states has died away by a small value of t 1 and any small admixture of lower states only becomes significant beyond t 2 . That we are able to obtain usefully large overlaps onto quite a large number of states with different quantum numbers, is the main technical feature of the present work. It requires the construction of a large basis of winding operators, a strategy that proved successful in our earlier calculations of the closed flux tube spectrum in D = 2 + 1 [9] and which we describe in detail in Section 3.4. The detailed construction is partly motivated by the range of quantum numbers we wish to encode, and so we turn to this now.

Quantum numbers
Consider a confining flux tube, with the flux in the fundamental representation and let it wind around the x-torus. Its quantum numbers can be chosen as follows.
• The length l = aL of the x-torus around which the flux tube winds.
• The number w of times (modulo N) that the flux tube winds around this torus. In this paper we shall only consider w = 1, except for the discussion of k-strings in Section 4.4.
• The momentum along the flux tube which is quantised by periodicity to be p = 2πq/l where q is an integer. Since the string ground state is (presumably) invariant under longitudinal translations, p = 0 will require some non-trivial excitation along the string with an additional excitation energy. By contrast, transverse momentum p ⊥ should simply lead to the usual dispersion relation and so is uninteresting to us, and we only consider p ⊥ = 0. Thus p will always refer to the longitudinal momentum although we shall sometimes label it by p , where this increases clarity. Since the energy does not depend on the sign of q, we will restrict ourselves to q ≥ 0.
• The flux tube can rotate about its symmetry axis and can have a corresponding angular momentum of J = 0, ±1, ... along that axis. • There is a 'transverse parity', P t , in the plane that is transverse to the symmetry axis and which is just like parity in D = 2 + 1: P t : (y, z) → (y, −z). It is plausible that the absolute ground state, with energy E 0 (l), is invariant under reflections, and so will have P t = +. Thus the lightest P t = − state requires a non-trivial excitation of the flux tube and the difference between the P t = ± ground state energies will measure the energy of that excitation. This parity does not commute with the angular momentum projection since under P t : J → −J. So when we choose to use this parity to label our states, we have to use |J| rather than J as the spin label (although we shall continue to refer to it as J where there is no ambiguity). States with |J| = 0 come in degenerate P t = ± pairs. • There is also a reflection symmetry across the string midpoint, which defines a corresponding 'longitudinal parity' that we call P l ; it reverses the momenta of the individual phonons. (The reversal of the direction of flux is compensated for by the simultaneous application of charge conjugation, C, which is to be understood from now on.) Since P l reverses the momentum it is only a useful quantum number if p = 0. So when p = 0 we label states by P t , |J| and |p|.
• Under charge conjugation, C, the direction of the flux is reversed. Thus a flux tube has zero overlap onto its charge-conjugated homologue (the center symmetry argument again) and so linear combinations of definite C are trivially degenerate. (Except for SU(2), where C is trivial, and for k-strings when k = N/2.) Thus C is, by itself, not interesting in general.
Our comparison to the effective string theory will be in static gauge (as described earlier) and so it is useful to discuss these symmetries and quantum numbers in terms of the transverse deformation of the string, h(x), and the associated phonon creation operators, a ± k , that arise from quantising the coefficients in the Fourier decomposition of h(x), and which carry momentum 2πk/l and spin (around the flux-tube axis) of ±1. The longitudinal and angular momenta of the flux tube are simply the sum of all the phonon momenta and spins. Under P t we have a + k → a − k 2 while under P l we have a ± k → a ± −k . For later reference we give in Table 1 the explicit phonon content of the lowest-lying states in the free string Nambu-Goto model.
The above discussion has been for a flux tube in continuous space. The cubic lattice has implications for angular momentum and the two-dimensional parity, P t . The relevant irreducible representations are those of the two-dimensional lattice cubic symmetry [21] which we denote by C 4 . It is a subgroup of O(2) that corresponds to rotations by integer multiples of π/2 around the tube axis. This makes values of angular momenta that differ by an integer multiple of four indistinguishable on the lattice, and factorizes the Hilbert space into four sectors: J mod 4 = 0, J mod 4 = ±1, J mod 4 = 2. In the continuum states of nonzero J are P t parity degenerate, but on the lattice this is true only for the odd J sector. In the even J = 0 sector the P t = ± partners will in general experience O(a 2 ) splittings. All this means we can denote our states by the 5 irreducible representations A 1,2 , E, B 1,2 of C 4 whose J and P t assignments are: one-dimensional except for E which is two-dimensional.

Operators
To calculate excited states efficiently we need a large basis of operators that explicitly includes ones that are less symmetric than the simple Polyakov loops. In our D = 2 + 1 calculations [9] we contructed a basis of ∼ 200 operators which provided us with good overlaps onto the ∼ 10 flux tube states that were light enough to be accessible within our statistics. In D = 3 + 1 we have two rather than one transverse direction so, very naively, the corresponding number of operators would be ∼ 200 2 = 40000. This is much too large. Instead we limit ourselves to ∼ 1000 operators which are partly chosen by looking at which operators turned out to be most important in the D = 2 + 1 calculations. By projecting these operators to the different quantum numbers of interest, we have matrices of correlation functions that are no larger than O(100), which is manageable. The operators we construct have shapes that lead to certain values of J, P t , P l , and p. This is achieved by choosing a linear combination of generalised Polyakov loops whose paths consist of various transverse deformations of simple Polyakov loops, at various smearing and blocking levels. All the paths used for the construction of the operators are presented in Table 2 and all together they lead to a basis of around 700 operators. (The blocking levels for different parts of such an operator need not be the same.) To construct, for example, an operator with a certain value of the spin J we proceed as follows. Begin with the operator φ α (n µ ) that has a deformation that we label as being at an angle α in the plane transverse to x (assuming the flux tube winds around the x-torus) and that we label as being located at the site x µ = an µ . We can construct an operator φ J (n µ ) that belongs to a specific representation of C 4 by using the formula: This construction assumes that we are in the rest frame or in a frame where the momentum and spin are aligned. In our case, where the spin is aligned along the flux tube, this means fixing the transverse momentum to zero, which we can achieve by summing over all transverse spatial translations: φ J (n t , p , p ⊥ = 0) ≡ ny,nz φ J (n t , n x , n y , n z ). It is straight-forward to show that φ J=0 belongs to either A 1 or A 2 (depending on its value of P t ), that φ J=±1 belongs to E, and that φ J=±2 belongs to B 1 or B 2 . The projection onto certain values of P t and P l is demonstrated pictorially in eqns(35,36,37) for an operator of J = 0, 2, 1 respectively.
If k = +1 then the operator φ E projects onto {E, P l = +} and if k = −1 it projects onto {E, P l = −}. It is important to note that the above construction does not distinguish between J and J + 4n where n is an integer. Thus when we label a state as J = 0, what we mean is J = 0, ±4, ±8, ..., and similarly for our J = ±1, ±2 labels. Usually one can safely assume that the lightest state will have the lowest value of J, but excited states may well have higher spins. From now on we shall use the simple J = 0, ±1, 2 label, but the reader should always be aware for what this is a shorthand.  So if we wish to project onto a flux tube that not only has spin J but in addition has non-zero momentum p = 2πq/l along its axis, while still having zero transverse momentum p ⊥ = 0, we can perform the sum Since the spin component along the axis of a boost is invariant under boosts, such an operator produces states of the desired spin even though it has not been defined in the rest frame. Note that since the longitudinal parity flips p → −p , it is not a useful quantum number when p = 0, and in that case we set i = k = 0 in sums like those in eqn(35) and we simply label our states by J or by |J|, P t . Since J = 2 and J = −2 differ by ∆J = 4, the operators for them are identical. However we can construct lattice operators that have transverse parity P t = ± as described above. Hence for even J we shall choose to label the states by |J| and P t . For odd J, on the other hand, ∆J = 4n and we can use eqn (34) to construct orthogonal ±|J| operators. Since these must be degenerate, we choose to calculate only one of these two states; but it should be understood that every such odd-J state is automatically accompanied by a degenerate −J state that we have chosen not to show explicitly.
In terms of our choice of operators, the present work should be regarded as exploratory. Not surprisingly we will find that the overlaps are not as good as in our earlier D = 2 + 1 study. However for a number of states we do have overlaps |c| 2 ≥ 0.9, at which point we can regard the difficult-to-estimate systematic errors as being under control. Of course, for heavier states our identification of an 'asymptotic' exponential behaviour becomes less reliable, and therefore so does our estimate of the overlap. On the other hand, in D = 3 + 1 one has more channels with different quantum numbers than in D = 2+1, and each of these channels has its own ground state(s) which corresponds to some excitation of the absolute ground state. Since these ground states are afflicted by fewer systematic errors, this allows us to obtain relatively reliable results for a number of string excitations even with our poorer overlaps.

Earlier lattice calculations in D = 3 + 1
Numerical exploration of the stringy nature of open and closed flux tubes, at the level of testing the Lüscher correction to the ground states, dates back to the mid-1980's. However the pioneering calculations for excited string states date to the early-90's, e.g. [22]. The interest here was both theoretical and phenomenological: the excited string energy can be used as a potential in a Schrödinger equation to get predictions for the masses of hybrid mesons where some of the quantum numbers are carried by excited glue. More or less simultaneously, there were calculations of Wilson loop expectation values investigating the match between string theory predictions and what one obtains in various gauge and spin models by, in particular, the Torino group, e.g. [23]. Expectations of Wilson loops are transforms of eigenspectraalthough care has to be taken with the self-energies associated with the boundaries -and provide an alternative way to test string models. Indeed it is in this body of work that one first sees a prolonged and serious focus on matching the Nambu-Goto model to numerical results.
Much of the work on flux tubes has focused on open flux tubes i.e. the static potential and its excitations. In this case there is a smooth transition between short-distance perturbative physics (the Coulomb potential) and the long-distance confining physics (the flux tube). While this transition is potentially of great interest (especially at N → ∞) it can introduce extra ambiguities in the extraction of the flux tube spectrum. The spectrum of closed flux tubes, stabilised by closing them around a spatial torus, is a particularly clean way to investigate the properties of flux tubes, and that is why we have focused on that approach in both our D = 2 + 1 [9] and D = 3 + 1 calculations.
As far as the eigenspectrum of closed flux tubes in D = 3 + 1 gauge theories is concerned, the most ambitious calculation we are aware of is the one in [24] which was part of a larger calculation largely focused on the spectrum of open flux tubes [25]. This calculation was in SU (3) and was performed at one value of the lattice spacing and for three values of the flux tube length, l. The timelike and spacelike lattice spacings were chosen to be different, with the former quite small, while the spatial lattice spacing was chosen quite large, a √ σ ∼ 0.5.
The energies of 12 excited states, of various quantum numbers (some being ground states within their channels) were calculated. The three flux tube lengths chosen were all quite long, l √ σ ∈ [4,8], the intention being to compare with the leading O(1/l) correction to the linear piece, σl, at large l. The results show an initial tendency to approach the theoretical Nambu-Goto expectation but then, typically, to overshoot at the largest value of l. (See Fig.2 of [24].) At the larger two values of l a rough agreement is observed with the theoretically expected level ordering and degeneracies, although there is an additional fine structure that is still prominent even at the largest value of l. Shortly after this work a calculation for SU (4) and SU(6) appeared, [26], that focused on the ground state of the closed flux tube, but also pointed to the qualitative agreement between the first excited state and the full prediction of the Nambu-Goto model, even at modest values of l.
Our calculations attempt to improve upon these earlier calculations in several ways. First, our point of view is that the simplest picture should emerge at N = ∞ and we therefore perform calculations not only for SU (3), but also for SU (5) and SU (6). We find that the 1/N 2 corrections are negligible for N ≥ 3, at our level of statistical accuracy, which means we can treat our extensive SU(3) calculations as being valid for larger N. Secondly, what we are interested in is continuum rather than lattice physics. Indeed, the (spatial) lattice spacing used in [24] was large, and so carried the risk of significant finite-a corrections. To control these in our calculations we mainly work at a much smaller value of a, and then explicitly check for O(a 2 ) corrections by performing a calculation at an even smaller value of a. We find that any dependence on a is negligible and so we can treat our results as being for the continuum limit. Thirdly, motivated by our D = 2 + 1 calculations [9], which show that the approximate agreement with the Nambu-Goto spectrum begins at remarkably small values of l, we focus within our calculations mainly on small to medium values of l, i.e. l √ σ ∈ [1.8, 4.8], where we can be reasonably confident that we have the systematic errors under control. As we remarked in Section 3.2, as one makes l larger, the energies of interest become larger, and one runs a rapidly increasing risk of a systematic error that typically leads to an overestimate of the energy. It appears plausible to us that this is the source of the large-l overshoot observed in [24]. Since we find that we can obtain interesting results on the spectrum without venturing to such large l, we choose not to do so in this study, except for one calculation at l √ σ ∼ 6.3 which we use both to check the qualitative behaviour of some states and to ensure that we have some appreciable overlap with the range of l where analytic expansions in 1/l should be applicable (for the lightest excited states). In addition, and perhaps most importantly from the technical point of view, we incorporate a a very large basis of operators in our calculations. Although the benefits are not as striking as in our earlier D = 2 + 1 calculations, it is primarily this technical improvement that allows us to maintain usefully good overlaps onto a large number of excited states at the small lattice spacings where we can confirm that we are close to the continuum limit.

Results
Our calculations of the spectrum of closed flux tubes that wind around a spatial torus of length l have been performed for the SU(N) groups and inverse bare couplings, β = 2N/g 2 , listed in Table 3. For orientation we also list some basic physical properties of the lattice gauge theories at those values of β. First there is the string tension a 2 σ which we extract from our calculation of the ground state string energy E 0 (l). (We fit E 0 (l) to the Nambu-Goto expression plus a O(1/l 7 ) correction. The value of a 2 σ extracted is not sensitive to any reasonable changes to the fitting function.) This tells us what a is in physical units. To express a in more intuitive 'fermi' units one can set σ to its real world value, √ σ ∼ 440MeV, although in all our field theories such units are, of course, fictitious. In these units we see that most of our calculations are for a ∼ 0.09fm while the SU(3) calculation at β = 6.338 corresponds to a ∼ 0.06fm. These values of a are small enough that the lattice corrections are, in general, small. We also show the critical value of the flux tube length, l = l c , below which the flux tube no longer exists because we are in a finite volume (anisotropic) deconfined phase that is a precise analogue of the finite temperature deconfined phase. (These values have been obtained by extrapolating/interpolating in β, and for SU (5) in N, the values in [17].) Our calculations will go down to values of l that are very close to l c . Here the flux tube should be about as wide as it is long and, naively, no longer looks anything like a thin string. Finally we list the value of the mass gap am G , which here is the mass of the lightest scalar glueball. (These values are obtained from the values in [29] and, for N = 5, in [30].) N β l/a ∈ a √ σ l c /a am G 3 6.338 [16,24] 0.12878(69) 11.99(9) 0.448(11) 3 6.0625 [9,20] Table 3: Parameters of our calculations with some corresponding properties of the gauge theories: the string tension, σ, the deconfining length, l c , and the mass gap, m G .
We have performed high statistics calculations with a small number of operators in SU(6) at β = 25.55 and in SU(3) at β = 6.0625. This allows us to obtain very accurate values for the ground state energy, E 0 (l), and to analyse its l-dependence in some detail. However to calculate the energies of excited flux tubes we need much better overlaps on the excited states than these calculations provide, and this we attempt to achieve with our much larger basis of operators, as discussed above, and with correspondingly lower statistics (because of the expense of such calculations). We perform the latter calculations in SU (3) at both values of β, so as to allow an explicit control of lattice spacing corrections. These corrections turn out to be small, so we only need to perform a similar SU(5) calculation, which allows us to control finite-N corrections, at the coarser value of a, as listed in Table 3. We shall see that the finite-N corrections are also small. This means we can focus our calculations of the energy spectrum on one N and on one lattice spacing, i.e. SU(3) at β = 6.0625, without loss of generality. Here we cover a much larger range of l and have much higher statistics, leading to a more reliable extraction of flux tube energies. We list separately the value of a √ σ extracted from this large operator calculation in Table 3, since this is what we shall use in our later analysis of the excited states. We see that this value is considerably more accurate than the other calculations with a large number of operators. We begin with a detailed analysis of our high statistics calculations of the ground state energy, E 0 (l). We then move on to the lightest states with different non-zero values of longitudinal momentum, p = 2πq/l, for q = 1, 2. These, as we have seen, must possess non-trivial excitations, and we compare the excitation energy to the predictions of the free-string Nambu-Goto model. Next we consider those states which correspond to the first, and sometimes second, excited energy levels of the p = 0 and p = 2π/l sectors, and whose energies we are able to calculate with some reliability. In one particular case we also look at excited states in the p = 4π/l sector. We follow this with a brief aside on the ground state of the k = 2 flux tube in SU(6). Finally we compare our results to the theoretical predictions described in Section 2.

Ground state
In this section we analyse the l-dependence of the ground state energy, E 0 (l), of a flux tube that winds once around a spatial torus. Our aim is to see what we can say about the universal corrections to the leading linear piece, σl. Previous lattice calculations, both for open and for closed flux tubes, e.g. [25,27,28], have provided good evidence that the coefficient of the O(1/l) Lüscher correction has the value appropriate to a simple bosonic string theory where the only massless modes are those of the transverse fluctuations, and, moreover, that this continues to be the case for larger N [26]. Motivated by the recent theoretical developments summarised above [4,5,7], we shall see if we can obtain any information on the coefficients of the higher order terms that we now believe to be universal. In addition, motivated by what we have observed in D = 2 + 1 [9], we shall test how well the Nambu-Goto free string theory describes the ground state energy.
The calculations in this subsection differ from our main calculations in that we use a small basis of operators which we know, from previous experience, will provide a very good overlap onto the ground state. Such calculations are much faster and this enables us to achieve a greater statistical accuracy for E 0 (l). We will present two calculations, one in SU(3) and one in SU (6). The SU(3) calculation is at β = 6.0625. Here one finds a √ σ ≃ 0.2, which translates to a ≃ 0.09fm if one wants to use such semi-fictitious units. Thus this is a small value of a, where one expects lattice spacing corrections to be small. We perform calculations for l/a ∈ [9,20]. The smallest length is only just above the first-order finite-volume deconfining transition which, at β = 6.0625, occurs at l c /a ≃ 8 [17]. The SU(6) calculation is at β = 25.55 which, as we see in Table 3, corresponds to roughly the same value of a in units of √ σ.
The value of l c is slightly larger than in SU(3), and the range of flux tube lengths studied is l/a ∈ [10,18]. We shall begin by using the SU(3) calculations to quantify some of the systematic errors, so as to ensure that they are small in our later calculations of the flux tube spectrum. The first part of this section is devoted to this issue, which can be skipped by the reader who is primarily interested in our results.

systematic errors
The systematic errors we focus upon are those associated with finite volume effects and with the extraction of an energy from a given correlation function.
The finite volume corrections are of two kinds: those that arise from the finite spatial extent and those that arise from the finite temporal extent. An example of the former is the emission by a flux loop of a virtual glueball that propagates around one of the orthogonal spatial tori, of length l s , before being reabsorbed. This will alter the flux tube energy by ∆E ∝g 2 e −m G ls where m G is the mass of the glueball. If l s is small enough then the lightest state will be a combination of a flux loop and a conjugate loop and in that case the leading large-l s contribution will come from its propagation around an orthogonal torus, giving ∆E ∼ g 2 e −2E 0 (l)ls . So we expect that as we reduce l we will have to have a larger spatial volume to minimise these corrections. Since we expect thatg 2 ∝ 1/N 2 all these finite volume corrections should disappear as N → ∞ (as long as the volume is large enough that we remain in the usual confining phase) but will be present at finite N.
The finite temporal extent l t means that the partition function Z in the denominator of eqn(30) receives contributions not just from the vacuum but from other states propagating around the time torus. (It is a 'finite temperature' partition function.) These same states will contribute to the path integral that appears in the numerator of eqn (30). At N = ∞ colour singlet states do not interact, and so this extra contribution will simply cancel between numerator and denominator. (In this sense there is no T dependence at N = ∞ in the confining phase.) However at finite N a state propagating around the time torus will interact with the flux tube that propagates between the Wilson line operators in the numerator and this will imply an incomplete cancellation and a shift in the flux tube energy. Typically such corrections will be ∆E ∼ O(e −mlt ) where m is the lightest state. For small l the lightest state is the winding flux tube, so that m = E 0 (l), and we must make l t larger as l decreases, so as to maintain E 0 (l)l t ≫ 1.
To monitor these finite volume corrections one can perform calculations for various transverse and temporal tori at each l. Such a detailed study can be found for SU (6) in [26]. Since it is not too expensive to make just one lattice torus very large, we choose to do so for the temporal torus and then we vary the spatial tori to monitor the corresponding finite spatial volume corrections. The values of l and lattice sizes of our SU(3) calculation are listed in Table 4. We calculate the ground and first excited state energies as described in Section 3.2. We show the ranges of t = n t a which we fit with the single exponential in order to obtain an energy. For the ground state we typically have a normalised overlap |c 2 0 | ∼ 0.98 − 0.99 allowing an accurate and reliable extraction of its energy. For the first excited state we have a much worse overlap, |c 2 1 | ∼ 0.8, and so we have to fit at larger t and the statistical error is about ten times larger. Higher excited states are even worse: this is why we need to use a much larger basis of operators, as we shall do in later sections.
In Table 4 we can see that the finite volume corrections on E 0 (l) are not large: no more than ∼ 3 − 4% even for the smallest volume. However since we are interested in identifying leading and subleading corrections to the linear dependence of E 0 (l), these are important. Since the leading large volume correction decreases exponentially with the spatial size, we can assume that the corrections to our largest lattices are much smaller than any variation we observe between our second largest and largest spatial volumes. We therefore use the values of E 0 (l) calculated on the largest volumes listed in Table 4, in the expectation that any finite volume effects are smaller than the quoted statistical errors. We also see from Table 4 that the finite volume corrections to the energy of the first excited state, E 1 (l), are at most comparable to the statistical errors. We therefore assume that even with our later more accurate calculations, we are safe in using the largest, or even second largest, volumes in Table 4 for our eigenspectrum calculations. SU(3) ; β = 6.0625 l/a L y × L z × L t sweeps×10 6 n t ∈ aE 0 (l) n t ∈ aE 1 (l) 9    We now consider two of the main systematic errors in extracting a ground state energy from a correlation function C(t). 1) We typically perform the fit using a single exponential to a range t ∈ [t 1 , t 2 ] as described in Section 3.2. However this cannot be entirely correct: the fact that there is a visible contribution from excited states for t < t 1 means that it also exists for t ≥ t 1 , even if small. This is an error which leads to our estimates of E 0 (l) being too high, albeit by an amount which hopefully is no greater than the statistical error. However since this error is systematic, its contribution after a combined analysis of several calculations, may well be significantly larger than the final statistical error. 2) Secondly, when performing the exponential fit we typically treat the errors on C(t) at different values of t as being independent, which is also incorrect. The calculations at different values of t are based on the same lattice field configurations and so may be quite highly correlated (particularly if the energy in lattice units is small). The value of χ 2 calculated with independent errors will be smaller than the true χ 2 and we thus run the danger of extracting E 0 (l) from a range [t 1 , t 2 ] where the fit appears acceptable but is in reality poor. (The error estimates on E 0 (l) are obtained by a jack-knife method and should be less affected.) To account for this, our heuristic procedure is typically to accept a fit only if the χ 2 per degree of freedom is ≪ 1. In practice we check that the χ 2 is not significantly decreased if we translate our fitting interval to larger t.
We will now explicitly estimate the effect of both the above errors on the SU(3) ground state calculations carried out on the largest spatial volumes listed in Table 4.
In the first column of Table 5 we display the ground state energies that we obtain with our default method, which uses a single exponential fit to the correlator of the best (variational) operator, and which assumes uncorrelated errors ('nocortt'), with a heuristically reduced estimate of what constitutes an acceptable χ 2 . We also show the range in t used for each fit.
In the second column of Table 5 we show the result of incorporating the error correlations in the statistical analysis ('cortt'). This is done in a standard way: one calculates the correlation beween the fluctuations of the correlator at t and the fluctuations at t ′ for all t, t ′ that are of interest. This gives a correlation matrix whose inverse then appears in the definition of χ 2 . What one is effectively doing is determining the linear combinations t c t φ † (t)φ(0) whose fluctuations are independent, and then using these to determine the χ 2 . In practice this procedure can produce anomalies if all the elements of the correlation matrix, whose inverse we take, are not determined very accurately. This is one reason we do not apply it systematically to our later calculations but rather test what difference it makes in this high statistics calculation. Since this is the correct definition of χ 2 , the criterion for an acceptable fit is now simply the normal one: unity per degree of freedom. Comparing the values of E 0 (l) in the first and second columns of Table 5, we observe that for larger values of l the results are, within statistical errors, the same. It is only for the smallest l that there is any sign of what might be a small difference.
We now try to estimate the maximum possible effect of an excited state contribution, by assuming that it is the first excited state that is providing this. (Since it is plausible that in choosing the lightest excited state we maximise the shift in E 0 .) We take the energy of the excited state, E 1 , from Table 4. We then perform a 2-exponential fit, with E 0 and the two overlaps, |c 2 0 | and |c 2 1 |, as free parameters. There is a partially subjective choice to make about the start of the fitting range t ∈ [t 1 , t 2 ]. (The value of t 2 is usually chosen for convenience and plays a minor role.) In practice we use the lowest value of t 1 for which we can get an acceptable fit for some values of the parameters, to determine the possible range of these parameters. We show the results of this in the third and fourth columns of Table 5, using uncorrelated and correlated errors respectively. Comparing these to the values in the first two columns we see that the shift in E 0 is only noticeable at the largest values of l. And, as we remarked, this estimate is intended to provide something like an upper bound on the effect of excited state contributions.

ground state analysis
In Tables 5 and 6 we present our high statistics results for the ground state energy in SU (3) and SU (6) respectively. SU (6) ; β = 25.55 l/a L y × L z × L t sweeps×10 6 n t ∈ aE k=1 0 (l) 10   What do these results for E 0 (l) tell us about the coefficient of the universal 1/l correction? In particular, do we find evidence for which would correspond to the simplest bosonic string universality class where the only massless modes on the flux tube are the transverse oscillations? There are many ways to address this question, and we proceed as follows. We have calculated E 0 (l) for the lengths l = l 1 , l 2 , ... where l i+1 ≥ l i . We define effective Lüscher coefficients and string tensions by and determine the parameters c ef f and σ ef f for each pair of values E 0 (l i ) and E 0 (l i+1 ). We expect that at small l the value of c ef f may be very different from its large-l limit, because of additional corrections that are higher powers in 1/l. So the relevant question is whether c ef f → 1 as l → ∞. Motivated by what we have observed in D = 2 + 1 [9] we perform a similar exercise for the Nambu-Goto expression, writing We first consider SU(3). We take the values of E 0 (l) from the second column of Table 5. (The other sets of values give very similar results.) We obtain the values of c ef f (l i , l i+1 ) shown in Fig. 1. Note that while the vertical bar on each point indicates the statistical error, the horizontal bar indicates the range, l i √ σ to l i+1 √ σ, from which the value of c ef f has been extracted. Given that additional massless modes typically contribute a multiple of ±0.5 to the Lüscher term in eqn (40), what we see in Fig. 1 provides quite convincing evidence for the validity of eqn(39). Equally interesting is the fact that when we use eqn(41) we obtain c ef f close to unity for all values of l. Thus any corrections to the Nambu-Goto expression for E 0 (l) must be small, even at our smallest values of l where the flux 'tube' is presumably a fat blob, not much longer than it is wide. This recalls the situation in D = 2 + 1, with the difference that here the deviations from c ef f = 1, although small, are large enough to be visible. The corresponding results for SU (6) are shown in Fig. 2. This is very similar to what we have just seen in Fig. 1, so we may assume that for all N ≥ 3 the central charge corresponds to a bosonic string theory where the only massless modes on the flux tube are the transverse oscillations. All this corroborates the claims from earlier calculations that the effective string theory of confining fux tubes in D = 3 + 1 SU(N) gauge theories is in the universality class of the simple bosonic string theory.
The theoretical analysis of [7] goes much further than this and tells us that when we expand the ground state energy E 0 (l) in powers of 1/l not only is the O(1/l) Lüscher correction universal, but so also are the O(1/l 3 ) and O(1/l 5 ) terms, and that their coefficients are precisely what one gets by expanding the Nambu-Goto free string expression in eqn (21) to that order. To test this prediction we fit to our results for E 0 , where the c N G n are the appropriate Nambu-Goto coefficients, and on theoretical grounds we expect γ ≥ 7. For a fixed γ the fitted parameters are σ and the coefficientc of the leading unknown correction. Beginning with SU(3) we take (as above) the values of E 0 (l) listed in the second column of Table 5. We fit these values with eqn(42) for various powers of γ and show in Fig 3 the χ 2 per degree of freedom of the best fit as a function of γ. We repeat the exercise for SU(6) using the values in Table 6. We see that the SU(3) calculation favours a value γ < 7. In fact the acceptable values are γ = 5 and 3. Since our two lowest values of l are very close to l c (see Table 3) and since the SU(3) transition at l = l c is weakly first order, it might be that these fits are influenced by this transition. In that case one would expect a more reliable result for SU (6) where the transition is robustly first order. As we see in

Ground states with p = 0
We now consider the lightest states with non-zero momentum along the flux tube axis, p = 2πq/l = 0. To have p = 0 a flux tube cannot be invariant under longitudinal translations: it must have transverse deformations, corresponding to the excitation of some modes along the tube, and it is the energies of these that interest us here. So to project onto these we need to use our extended basis of operators. Of course, the energy not only increases from the inclusion of this excitation energy but also from the p 2 term. So to avoid the systematic errors associated with large energies, we only present results for the lightest energy levels with q = 1 and, with caveats, for those with q = 2. Of course we also simultaneously calculate the q = 0 ground state energy, in the same basis of operators, to obtain a value for the string tension a 2 σ, as described above. This value is then inserted into eqns (19,20,21) to provide parameter free predictions for the energies of states with q = 0 in the Nambu-Goto model. The parameters of our calculations are listed in Table 8 with some of the physical properties listed in Table 3.
When presenting our results showing how the flux tube energies E n vary with the flux tube length l, it is clearly preferable to do so using physical units rather than the lattice units, (aE n and l/a) in which these quantities are actually calculated. We choose to use the calculated value of a √ σ as our unit and typically plot E n / √ σ versus l √ σ. Since, for this purpose, any variation of the extracted value of a √ σ with the choice of how one fits the q = 0 ground state is completely negligible, the reader should not be concerned that there might be some hidden circularity in our plots. Using physical units makes it easier to assess the significance of the plots, and since we find that the O(a 2 ) lattice spacing corrections to our energies are not visible, there is nothing misleading about this use of physical units. Before continuing, a technical caveat. Our q = 2 energies are amongst the largest of any energies calculated in this paper. So the corresponding correlation functions decrease rapidly and it is frequently the case that one 'identifies' an effective energy plateau starting at some t = t 0 where the error to signal ratio for t > t 0 is too large to be useful. This means that the plateau identification is of low significance and the possibility of the energy being a significant overestimate is serious. In this case it makes sense to be guided, in identifying the energy plateaux, by those calculations that have the greatest statistical accuracy. As we see from Table 8, that is the SU(3) calculation at β = 6.0625. So in our SU(5) calculation, which is at roughly the same value of a √ σ, we use the same t-ranges as in these SU(3) fits in order to extract the heavy q = 2 energies. For orientation we begin by recalling the spectrum in the Nambu-Goto model. For q = 1 the lightest state has one phonon with minimal non-zero momentum: p = 2πk/l with k = 1. In the continuum limit this will produce two degenerate states, with J = ±1, or alternatively two |J| = 1 states with P t = ±. The lightest q = 2 states can either be formed from a single phonon with k = 2 or from two phonons, each with k = 1. The former gives two states with |J| = 1 and P t = ±, while the latter gives states with J = 0 and J = ±2 (the latter we organise into states with |J| = 2 and P t = ±). There is only one J = 0 state with P t = +, because the phonon operators commute so that a + 1 a − 1 |0 > and a − 1 a + 1 |0 > are identical states. On a cubic lattice the two |J| = 1 states belong to a two-dimensional representation and are exactly degenerate even at finite a. So in this section we choose to show the energy of only one of these states (whether for q = 1 or for q = 2). With this convention, we expect to find one J = 1 state occupying the lowest energy level of the q = 1 sector, and four states (J = 0, P t = +; |J| = 2, P t = +; |J| = 2, P t = −; J = 1) in the lowest energy level of the q = 2 sector. These low-lying states of the Nambu-Goto model are displayed in Table 1.
We begin with our high statistics calculation of the lightest q = 1 and q = 2 energy levels, in SU(3) at β = 6.0625. In Fig. 4 we plot the energies against the flux tube length, all in units of σ, and we also include the q = 0 ground state energies from which we extract the value of a 2 σ. We observe that the q = 1 ground state energy agrees extremely well with the parameter-free Nambu-Goto prediction, as do the nearly degenerate q = 2 states, albeit within significantly larger statistical errors. The number and quantum numbers of these q = 2 states are also exactly as predicted by Nambu-Goto. It is striking that the quantitative agreement persists down to the shortest flux tubes where l √ σ ∼ 2 and where the flux tube width, presumably ∼ 1/ √ σ, is comparable to its length.
In Fig. 5 we display the analogous results for SU (5). Again the q = 1 ground state energy agrees very well with the Nambu-Goto prediction. However for the q = 2 states the errors are now visibly larger, and while it is clear that they are roughly consistent with Nambu-Goto, it is hard to make a stronger statement than that.
As we remarked above, a substantial part of the extra energy of the q = 0 states may arise from their non-zero momentum, i.e. from the (2πq/l) 2 term in E 2 , and we are not interested in that. Rather we want to isolate the flux tube excitation energy and see how that compares with, for example, the Nambu-Goto model prediction. To do so we follow [31] and define where the Nambu-Goto value follows from eqn(21), E(q; l) is our result for the ground state energy with momentum p = 2πq/l, and E 0 (l) is the energy of the (absolute) ground state with p = 0. We choose to subtract the calculated value of E 0 (l) rather than the Nambu-Goto prediction for it, since the former already has some small corrections which we would regard as belonging to the Casimir energy rather than to the excitation energy that we are trying to isolate here. Note also that we assume any lattice corrections to the continuum energymomentum dispersion relation to be negligible at this β: an uncontroversial assumption, but one which it would be useful to check explicitly. In Fig. 6 we plot some of the q = 1 and q = 2 excitation energies, obtained after processing the energies plotted in Fig. 4 and Fig. 5 through eqn(43). We omit the SU(5) q = 2 values, because the errors are too large for them to add much here. By removing the momentum contribution, which will presumably arise in any string model, this analysis gives us a more precise appreciation of the significance of the apparent agreement in Figs 4 and 5 between our results and the predictions of the free string Nambu-Goto model. For q = 1 the errors remain small and the agreement down to the smallest values of l √ σ is as striking as before. For our lightest, nearly-degenerate, q = 2 states, which we find to have precisely the quantum numbers predicted by the Nambu-Goto model, the relative errors are now larger. Nonetheless, while one clearly sees significant deviations from the free string prediction at the smaller values of l, one also clearly sees that these states converge rapidly to the Nambu-Goto value as l increases, and are in agreement with it within errors for l √ σ ∼ 4. So, just as in D = 2 + 1 [9], the level of agreement with the free string model, for these modest flux tube lengths, is surprising. It is clear from Fig. 6 that the lightest q = 1 excitation energy shows no dependence on N within the small statistical errors. What can we say about the q = 2 excitation energies? Since in this case the statistical errors are too large for a direct comparison to be informative, what we can do instead is to compare the effective energies extracted from t = a to t = 2a, where our statistical errors are small enough for a significant comparison. Of course these energy estimates will contain a significant admixture of higher excited states, since our SU(3) effective energy plateaux typically set in for t ≥ 2a, so we are not looking for agreement with the Nambu-Goto prediction. Rather we simply want to compare the SU(3) and SU (5) values as a measure of the N-dependence of the q = 2 spectrum. The result of processing these effective energies through eqn(43) shows that beyond the very shortest flux tube length, the SU(3) and SU(5) excitation energies show no significant N-dependence within reasonably small errors.
The above SU(3) and SU(5) results have been obtained at a single common value of the lattice spacing, a √ σ ∼ 0.2. Do these results survive the continuum limit? We address this question with an SU(3) calculation at β = 6.338 where the lattice spacing is significantly smaller, a √ σ ∼ 0.13. The lightest q = 0, 1, 2 energies are plotted in Fig. 7 and if we compare to Fig. 4 we see that there appears to be no significant a-dependence within the errors. Reassured by all this we choose not to attempt a comparable small-a SU(5) calculation. We can assume from these checks that our calculations, although at fixed lattice spacings and fixed values of N, do in fact describe the continuum limit of the gauge theory, at arbitrary N ≥ 3, within the errors of our calculations.

Excited states
Having examined the ground state energy levels for q = 0, 1, 2, we now turn to the first, and sometimes second, excited energy levels. We shall mostly restrict ourselves to flux tubes with q = 0 and q = 1, since only here are the errors small enough for a precise analysis.
For each value of p = 2πq/l, we shall begin by recalling the excitation spectrum in the Nambu-Goto model. The string tension σ is calculated from the ground state energy so there are no free parameters. We shall find that not only are the quantum numbers of the lightest states precisely as predicted by this model, but that the calculated energies of most of these states are close to the Nambu-Goto predictions down to remarkably small values of the flux tube length l. However we shall also find that, for each value of q, the lightest excited state has an energy that is much lower than predicted, and that this energy increases only weakly as l grows. In each case the anomalous state has J = 0, P t = − quantum numbers and, partly because it is light, it is accurately determined in our calculations.

q = 0 excited states
In the Nambu-Goto model, the states in the first excited q = 0 energy level have two phonons of opposite momentum, with this momentum taking its minimal absolute value |p| = 2π/l. Each phonon carries unit spin, and these can be either parallel or anti-parallel, leading to two states with J = 0 and two with J = ±2. (See Table 1.) These can be re-organised into four states with quantum numbers (|J|, P t , P l ) = (0, +, +), (0, −, −), (2, +, +), (2, −, +). In the Nambu-Goto model all these states are degenerate.
In our numerical calculations we are able to obtain, with reasonable accuracy, the first four excited states above the absolute ground state, and we find that these have precisely the quantum numbers predicted by Nambu-Goto. We then process the calculated energies through eqn(43) so as to obtain the excitation energies, ∆E 2 (q = 0, l). We display the results in Fig 8, together with the Nambu-Goto prediction. We observe that 3 of the 4 states show the by now familiar rapid convergence to the Nambu-Goto prediction as l √ σ increases. However one state, the one with 0 −− quantum numbers, displays a very different behaviour: it is much lower than the predicted value and approaches that value only very slowly as l √ σ increases.
The fact that the energy is so much lower is robust against all our systematic errors: that is to say, this anomalous behaviour is a reliable result.
Since it is only in the N → ∞ limit that we should expect a simple stringy description to be possible, it is natural to ask if this anomalous behaviour is a finite N effect. To investigate this we repeat the calculation in SU(5) at a value of a that is similar to that of SU(3) at β = 6.0625. Since the leading large N corrections are expected to be O(1/N 2 ), if the anomaly is a finite-N correction it should be smaller in SU(5) by a factor of ∼ 3. We plot the excitation energies obtained in SU(5) in Fig 9. We see that the anomalous behaviour of the lightest 0 − state becomes more rather than less pronounced in SU (5). It is clear that this anomaly is a prominent feature of the gauge theory at N = ∞.
Does this anomalous behaviour survive the continuum limit? To address this question we repeat the calculation at β = 6.338 where a 2 σ is smaller by a factor of ∼ 2.3. The latter is the relevant measure for lattice spacing corrections, since these are expected to be O(a 2 ) for our plaquette action. We plot the corresponding excitation energies of the 0 −− state in Fig 9. We see that it continues to be anomalous, much as it was at the coarser value of a. This confirms that what we are seeing is continuum physics.

q = 1 and q = 2 excited states
For q = 1 the first excited energy level in the Nambu-Goto model is composed of states that have either two or three phonons. (See Table 1.) In the states with two phonons, one has momentum p = 2πk/l with k = 2, and the other has k = −1. The two unit spins add to give 2 states with J = 0 and two with J = ±2 and these can, as usual, be re-organised into states with P t = ± for each of |J| = 0, 2. In the three phonon case, two of the phonons have k = 1 and one has k = −1. There are now three unit spins to be added, giving three sets of states: two pairs with J = ±1 and another pair with J = ±3. On a cubic lattice the states that become J = 1 and J = 3 in the continuum limit reside in the same E representation and cannot be distinguished on straightforward symmetry grounds. Moreover the P t = ± partners are in the same multiplet, and are exactly degenerate, and so it suffices to focus on the P t = + states. Thus we are looking for 4 states with even J and 3 with odd J. We shall refer to the latter as J = 1, although of course the state may in fact have any odd value of J -in particular J = 3.
In our numerical calculations we are able to obtain, with reasonable accuracy, a number of q = 1 excited states above the ground state. In the even-J sector we find that the lightest four states have the quantum numbers J Pt = 0 ± , 2 ± , which is precisely what one expects for the 2 phonon component of the first excited energy level. Moreover, in the odd-J sector, we find three states that are close in energy, which is, again, precisely what one expects for the 3 phonon component of the first excited energy level.
We begin by processing the even-J energies through eqn(43) to obtain the excitation energies, ∆E 2 (q = 1, l). We first do this for our most accurate and extensive calculation, SU(3) at β = 6.0625, using the energies listed in Table 11. The results, for lengths l/a ≤ 24, are displayed in Fig 10, where they are compared to the Nambu-Goto prediction. For comparison we also show in this plot the q = 1 ground state, which has J = 1 and was shown earlier in Fig. 6. We see, just as we saw for the lowest 4 excited states with q = 0, that 3 of the 4 states show a rapid convergence to the Nambu-Goto prediction as l increases, but one state, again with 0 − quantum numbers, displays a very different behaviour: it is much lower than the predicted value and only appears to approach that value very slowly as l increases. As before, the fact that the energy is lower means that this anomalous result is robust against essentially all our systematic errors.
As in the case of q = 0, we repeat the calculation at the smaller value of a. The results for the 0 − state are plotted in Fig 11, and we see that there is no significant dependence on a. So once again, we are able to claim that this anomalous behaviour is clearly part of the physics of the continuum theory. Similarly a calculation in SU(5), whose results are also displayed in Fig 11, shows that any large N corrections to this anomalous behaviour are small. Once again we can claim that this represents the continuum physics of all theories with N ≥ 3.
We turn now to the J = odd states with q = 1. In Fig. 12 we plot the energies of the ground state and the first four excited states. We see that above the ground state there are three nearly degenerate excited states, consistent with the degeneracy predicted by Nambu-Goto for the first excited energy level. Moreover their energies are close to the Nambu-Goto prediction for that energy level. We may therefore assume that what we are seeing is in fact two J = 1 states and one J = 3 state. The fourth state is well separated from the other three and is very close in energy to what Nambu-Goto predicts for the second excited energy level. There should of course be several J = odd states in this energy level, but we do not attempt to identify these since their energies are quite large, so that the errors will be large as well. However we can say that there do not appear to be any anomalously light odd-J states belonging to this energy level. So we can conclude that in the q = 1 odd-J sector we see no sign of any anomaly up to the second excited energy level.
We end this section by extending our analysis to the q = 2 sector. In Fig. 13 we show the energy of the lightest 0 − state in our two SU(3) calculations. For comparison we show the Nambu-Goto q = 2 ground state energy, to which the lightest 0 + , 2 ± states rapidly converge, and the energy of the first excited Nambu-Goto level, to which the lightest 0 − should converge as l → ∞. This plot clearly confirms that the lightest 0 − state is anomalous for q = 2, just as it is for q = 0 and q = 1.

The anomalous 0 − states
We have seen that nearly all the lightest states of a closed flux tube converge very rapidly to the free-string Nambu-Goto prediction as l increases. The striking and sole exceptions have been the lightest J Pt = 0 − states, and this is so in every sector of longitudinal momentum that we have investigated. These 0 − states show a very large deviation from the free-string prediction, with what is at best a very slow approach to the latter. This deviation persists, perhaps becoming even slightly larger, at larger N. This makes these states very interesting: they are our first candidates, whether in 2+1 or 3+1 dimensions, for states that might be reflecting the non-stringy massive modes that one expects to be present for a confining flux tube.
One possibility is that this 0 − state is a mixture of the massive and massless modes, with this mixing becoming large at small l. The coupling will inevitably involve derivatives of the massless Goldstone field and since the momenta of the phonons in the first excited Nambu-Goto state is |p| = 2π/l → 0 as l → ∞, the mixing will decrease as we increase l. Moreover as l increases the gap between the massless modes decreases ∝ 1/l and so the number of these modes that are close enough in energy to the massive mode to mix with it grows ∝ l. So the mixing with any individual mode must eventually decrease as l grows. In such a scenario we expect this 0 − state to approach the lightest 0 − Nambu-Goto energy level as l increases.
A different possibility is that the lightest 0 − state is simply a massive excitation of the ground state string. Because of the underlying flux tube, the energy of this state should grow with l, and because its energy maintains a finite gap with respect to the ground state, it will initially appear to approach the first excited Nambu-Goto level, but at some higher l it will cross that level, and at even larger l it will cross higher excited Nambu-Goto energy levels.
In such a scenario there should be an additional 0 − state that is the first massless excitation and which approaches the lightest 0 − Nambu-Goto energy level as l grows. Although such an excited 0 − state will lie at a higher energy at small and moderate values of l, and so will be harder to identify accurately, it is clear that its identification is important if we are to make any progress with understanding the nature of our anomalous states.
We therefore begin by trying to identify the next excited 0 − states, to see if they include plausible candidates for stringy states. In Fig. 14 we plot the excitation energies, obtained using eqn(43), of the lightest and first excited p = 0 0 − states, taken from all our three calculations. We observe that in the range of smaller l, l √ σ ≤ 3.2, where we have results from all three calculations, the excited 0 − state appears to be approaching the second excited Nambu-Goto energy level as we increase l, much as the lowest 0 − state appears to be approaching the first excited Nambu-Goto level, and with a similar anomalously large deviation. So this would suggest that there is no higher excited 0 − state to approach the first Nambu-Goto 0 − level, and therefore that the lightest 0 − will have to do so: and this would point to the first scenario discussed above. However the larger l calculations in Fig. 14, which we have performed only in SU(3) at β = 6.0625, radically alter this picture: as l grows the excited 0 − state turns over and begins to approach the first excited energy level. We include here, for the first time, some results from our l = 32 (l √ σ ∼ 6.3) calculation. As warned, the errors are large, but the result is qualitatively striking: the lightest 0 − appears to have crossed the first excited Nambu-Goto level and is approximately degenerate with the first excited 0 − , which is consistent with approaching this Nambu-Goto level. This now suggests that the physics is in fact closer to the second scenario described above: the state that is the lightest 0 − at smaller l is (mainly) a massive excitation that crosses the stringy energy levels, while the first excited 0 − is (mainly) the stringy state associated with the first excited Nambu-Goto energy level. We now turn to the q = 1 0 − states where we are able to obtain some results for the second as well as for the first excited states. We plot the corresponding excitation energies in Fig. 15. For l √ σ ≤ 3.2, what we see is similar to what we saw for q = 0 in this range of l: the first excited 0 − state also appears to be anomalous, approaching only very slowly the corresponding free-string prediction (the upper horizontal line). However this behaviour continues at our larger values of l, in contrast to what we have just seen for the q = 0 states in Fig. 14: there appears to be no clear signal that the 0 − excited state turns over, so as to approach the Nambu-Goto 0 − ground state (although the errors are large enough to allow that at larger l), or that the 0 − ground state crosses this energy level. This difference between the q = 0 and q = 1 results is highlighted if one considers the difference of the excitation energies of the lightest and first excited energy levels (labelled E 0 , Looking at the q = 1 results in Fig. 16 we see that the difference remains consistent with the Nambu-Goto expectation for all values of l. The large anomalous contribution in the lightest and first excited q = 1 0 − states appears to be a common additive contribution that simply cancels in the difference leaving just the excitation energy of the massless stringy modes. But with the important caveat that the errors are quite large at larger l. By contrast, the q = 0 results in Fig. 17 show that after an initial behaviour similar to that seen for q = 1, the gap in the excitation energies unambiguously decreases towards zero. Since the q = 0 results for the turnover appear unambiguous, and the q = 1 results at larger l have much larger errors, the simplest hypothesis consistent with these errors is that the q = 1 excited 0 − has in fact begun to turn over at the largest values of l, and so is paralleling the behaviour of the q = 0 state, albeit at somewhat larger values of l. This latter behaviour is to be expected because of the phonon contents of the q = 1 and q = 0 Nambu-Goto states, (a + 2 a − −1 − a − 2 a + −1 )|0 and (a + 1 a − −1 − a − 1 a + −1 )|0 respectively. As remarked earlier, one would expect higher momentum phonons to interact more strongly with the massive modes, because the interaction involves derivatives of the massless fields. Since the phonon momentum is 2πk/l, one expects the k = 2 phonon for a length l to have the same interaction as a momentum k = 1 phonon at length l/2. This heuristic argument would roughly suggest that the interaction energy of the q = 1 state of length l will have corrections that are a mean of those of the q = 0 states of length l and l/2, so that the turn-over (indicating a decrease in this interaction energy) will be delayed for the q = 1 state as compared to the q = 0 state. All this would seem to be consistent with what we see for the first excited 0 − states, although without more accurate calculations one cannot say anything stronger than that.
To decide whether the q = 0 and q = 1 0 − ground states are in fact excitations of massive modes, a quantitative analysis is needed. Since we are not in a position to provide that at present, the following simple heuristic analysis may be useful in deciding whether such a scenario is plausible 3 . Let us make the simple assumption that the lightest state with a massive mode excitation of the flux tube has an energy where m is the mass scale of this mode, p is the momentum it carries, and E 0 (l) is the energy of the flux tube ground state. One then finds We plot this function in Figs 14 and 15 with a value for the mass scale m √ σ = 1.85.
chosen to provide a qualitative fit. We note that this is an entirely reasonable value for the mass scale, being roughly half the mass gap, m ≃ M G /2, where G is the lightest scalar glueball. However, while this simple ansatz readily fits both the q = 0 and q = 1 0 − ground states and the differences between them, it does not have anything specific to say about the behaviour of the higher 0 − states. If on the other hand we wish to include both the ground and first excited 0 − states in a single analysis, a simple way to do so is to assume a simple 2-body mixing between a Nambu-Goto stringy state and a new, possibly massive state, with a mixing parameter that eventually decreases as l increases (for the reasons given earlier). Such a mixing leads to an equal splitting in ∆E 2 around the average of the un-mixed energies and since, as we see from Fig. 14, this average is close to the Nambu-Goto energy (taking seriously the large errors on our largest value of l) and since the unmixed stringy state already has this energy, this means that our putative massive mode is also close in energy to that of the Nambu-Goto state and so looks like a massless mode -an unexpected conclusion. Moreover the same analysis does not then describe the q = 1 energies in Fig. 16. So such a simple mixing picture does not appear to be tenable.
Turning now to the second excited state in the q = 1 0 − sector, we see from Fig. 15 that it does not seem to be anomalous, but is close to the appropriate Nambu-Goto prediction. Recall that while the lightest q = 1 0 − Nambu-Goto state is constructed from two phonons, the first excited 0 − energy level involves states with both 2 and 4 phonons, with both sectors producing 0 − excited states. Thus there is room in the Nambu-Goto first excited energy level for this 0 − state, irrespective of whether the anomalous first excited state turns out to belong to the lightest or to the first excited string energy level.

Closed loops of k=2 flux
In this paper our main focus is on confining flux tubes that carry flux in the fundamental representation. However, one can also study flux tubes carrying flux in higher representations. In particular consider fluxes that emanate from sources that transform as the product of k fundamentals and any number of adjoints, where 0 ≤ k ≤ N/2. The corresponding flux tubes are called k-strings. At finite N such a flux tube will be unstable if it can decay via gluon screening to a flux tube with lower energy per unit length. However the latter will also be in the same k-sector since gluons are adjoint particles. So the k-string with the smallest string tension, σ k , will be absolutely stable. It turns out [32,28] that σ k < kσ, so the lightest k-string is not merely k independent fundamental flux tubes, but may be thought of as a bound state of these. Such a flux tube clearly has additional internal structure to that possessed by a fundamental flux tube and it would be interesting to see how this affects the energy spectrum of closed k-strings. This would provide us with some insight into how such extra massive modes are encoded in the effective string theory that describes the k-strings. Precisely such a study has been carried out in D = 2 + 1 gauge theories [31]. In this paper we do not attempt such a study but initially limit ourselves to the k = 2 ground state, as calculated in our high statistics SU(6) calculation (where our operator basis is too small to attempt a full spectrum calculation). Let ψ 0 be the best operator that our variational calculation produces for the absolute ground state, using our full basis of operators. We list in Table 7 the energy of this ground state as a function of its length l. Now, the operator basis we use contains equal numbers of totally symmetric, k = 2S, and totally anti-symmetric, k = 2A, operators. (We refer to [31] and [28] for a more detailed discussion of the appropriate operators.) What we find is that the absolute ground state, ψ 0 , is almost entirely k = 2A, corroborating the observation in [29]. To be specific, if we denote by Φ 2S the space of k = 2S operators, φ 2S , spanned by our basis, then the overlap of ψ 0 on Φ 2S turns out to be as follows: : l = 16 = 0.00003 (1) : This is extremely small and tells us that to a very good approximation the k = 2 absolute ground state is in the totally antisymmetric representation. (And indeed, if we restrict our operator basis to k = 2A, the ground state ψ 0 that we obtain differs insignificantly from the one obtained using the full basis.) This strongly suggests that, just as we showed in D = 2 + 1 [31], the low-lying eigenspectrum consists of states that are, to a good approximation, either wholly k = 2S or k = 2A. While we cannot attempt here an analysis as complete as the one we performed in D = 2 + 1, we are able to calculate the ground state in the k = 2S sector, whose energies we list in Table 7. This is again a state that changes negligibly if we use the full k = 2 basis rather than the restriction to k = 2S.
In Fig. 18 we show how the energies of the k = 2A and k = 2S ground states vary with l. On the plot we also show the best fits using the free string Nambu-Goto expression, where the parameters being fitted are the respective string tensions, σ k=2A and σ k=2S . On such a plot, where both E and l are rescaled by the appropriate string tensions, the two Nambu-Goto curves will of course coincide. What is interesting is that the two sets of energies also appear to coincide: the (slightly) unstable k = 2S flux tube behaves just like the stable k = 2A flux tube. And both show the approximately linear rise with l associated with the formation of a confining flux tube.
It is clear from Fig. 18 that at smaller l the corrections to the Nambu-Goto fit become quite substantial. To resolve this more clearly we perform fits using eqns(40) and (41) just as we did in the case of the fundamental flux tube. The results are shown in Fig. 19 for the k = 2A flux tube. (The errors on our k = 2S energies are too large to allow a useful analysis of this kind.) Comparing to the comparable plots for fundamental flux tubes, in Figs 1 and 2, we see that the corrections to Nambu-Goto are very much larger here. This is just what we observed in our D = 2 + 1 [31] calculations. (And has previously been observed for the effective Lüscher correction in the D = 3 + 1 SU(6) gauge theory in [26].) Nonetheless this plot does provide some additional significant evidence that the k = 2 flux tube is also in the bosonic string universality class.

Effective string theory: discussion
Our above results show that most of the low-lying excited states of the fundamental flux tube are remarkably close to the predictions of the free-string Nambu-Goto model, even down to values of l where the flux tube is not much longer than it is wide (∼ 1/ √ σ). Now, we recall that the effective string analysis in [7] has shown that if we expand the excited string state energies E n (l) in inverse powers of l then they are universal to O(1/l 5 ) in D = 2 + 1 and to O(1/l 3 ) in D = 3 + 1, and that they must therefore agree with Nambu-Goto to that order. (If they are in the bosonic string universality class as is strongly implied by the numerical results for the ground state.) It is therefore natural to ask to what extent our numerical results are merely reflecting this theoretical result.
To address this question we take, as an example, the q = 0 excited states in SU(3) at β = 6.0625, whose excitation energies were plotted in Fig. 8, but omitting for the moment the 'anomalous' J Pt = 0 − state. We plot their energies as a function of length l in Fig. 20. We also plot the theoretical prediction to O(1/l 3 ), which includes all the terms that are currently known to be universal. The value of σ comes from fitting the ground state energy so this prediction is in fact parameter-free. (The variation of σ with the precise form of the ground state fit is completely negligible in the present context.) For comparison we also plot the prediction to just O(1/l), i.e. just the 'Lüscher correction', and also the full Nambu-Goto prediction. We see from Fig. 20 that the agreement of our numerical results with Nambu-Goto is, in fact, far too good to be a mere corollary of the known theoretical universality results. Does this perhaps indicate that the next one or two terms, in powers of 1/σl 2 , are also universal, and would this then suffice to explain our numerical results? The answer to the latter question is no. As indicated by the way that the O(1/l) and O(1/l 3 ) curves in Fig. 20 oscillate around the Nambu-Goto curve, the expansion of the Nambu-Goto expression in powers of 1/l is divergent in the range of l where most of our numerical results are located. Indeed, as we see from eqn (21), the expansion of the energy in powers of 1/σl 2 , for these N R = N L = 1 states, ceases to be convergent for So the analytical derivative expansion for this state only holds for l √ σ 4.8. Although we have one or two values in this range, and therefore make contact with the region to which the theoretical analysis applies, the bulk of our calculations lie at smaller values of l. Thus in reality they complement the theoretical analysis, while maintaining enough overlap to ensure some continuity between the numerical and analytical results.
Our results for the range of smaller l where the derivative expansion no longer converges, are telling us something simple about that expansion: it should, to all orders, remain close to that of Nambu-Goto. In effect our numerical results are telling us that the appropriate starting point for an effective string picture of closed flux tubes is neither the classical background string nor the Gaussian action, but rather the full free string Nambu-Goto action. That is to say, if we write the energy of the the n ′ th eigenstate as where E N G n (l) is the Nambu-Goto energy, then Fig. 20 tells us that for the 0 + and 2 ± states that are shown there this correction term, δE n (l), is 'small' over the whole range l √ σ ≥ 2.5: Since the theoretical analysis in [7] tells us that the expansion of δE n (l) in powers of 1/l cannot start with a power smaller than 1/l 5 , we know that at large l it must have the form: The different possibilities can be usefully illustrated by the following two scenarios.
• The fact that the corrections are small down to small l would follow most simply if the expansion coefficients in eqn(53) were to be small so that the first term dominated even at small l: Indeed if we fit the observed deviations of the states in Fig. 20 from Nambu-Goto, we find that they can be roughly described by simple O(1/l 5 ) corrections. This is demonstrated in Fig. 21, where the dotted curves incorporate such a correction to the Nambu-Goto prediction.
The magnitude of the coefficients in the two fits are c n ≃ −16, +9 respectively. These may appear quite large, but they are in fact ∼ 50 times smaller than the coefficient of the O(1/l 5 ) term in the expansion of E N G (l), and so should be regarded as 'very small' in natural units that include appropriate powers of 4π. Nonetheless, while this qualitative fit demonstrates the possibility of such a scenario, our calculations are not accurate enough to encourage a more precise analysis.
• Alternatively, it may well be that as we decrease l from large values, the series expansion for δE n (l) diverges at some value of l well above l c , in the same way as the Nambu-Goto expansion does. For lower l it should be resummable into some finite expression, since we know the flux tube energy does not diverge for any finite l > l c . This means that at small l √ σ the apparent power of l may well be very different from the asymptotic ∝ 1/l 5 behaviour.
This can be illustrated by the following simple example of a resummed series: (We assume that d n /σl 2 ≫ 1 over most of our range of l, as this is dictated by the radius of convergence we have assumed.) At small l this can decrease faster or slower than 1/l 5 depending on the sign of γ n (which in turn depends on how the coefficients in eqn(53) oscillate in sign). Let us assume, for example, that this series expansion diverges around the same length as Nambu-Goto, i.e. at l √ σ ∼ 5 in the context of the states in Fig. 20. This implies that d n ∼ 25 in eqn(55). Given that the deviation from Nambu-Goto in Fig. 21 is δE/ √ σ ∼ 0.5 at l √ σ ∼ 2, we see from eqn(55) that, very roughly, we need c n ∼ 16/6 γn . If γ n > 0 this leads to a value for c n that seems unnaturally small when compared to the corresponding Nambu-Goto coefficient. To allow a larger value of c n one needs to have γ n < 0. In that case, as we see from eqn(55), the effective power of 1/l that one sees at small l will be smaller than 5. One sees from this example the linkage between the coefficient c n of the leading asymptotic correction, the value of l √ σ at which the expansion of δE n diverges, and the effective power of the correction at smaller l √ σ. We illustrate this in Fig. 22 with a fit to the excitation energies of these states, using a simple c/(l √ σ) 4 correction term in ∆E 2 . It is straightforward to see that this corresponds to a correction term in eqn(51) of δE = (E 2 N G + ∆) 1/2 − E N G where ∆ = 4πσc/(l √ σ) 4 , which is a resummed expression that diverges at the same place as Nambu-Goto and whose leading piece is of the form in eqn(55) with γ n = −1/2. Clearly this illustrates that such a fit is possible.
So far this is very similar to what we found in 2+1 dimensions [9]. However, the important difference is that here we have also identified some additional states, the ground and first excited 0 − states with longitudinal momenta q = 0 and q = 1 (as well as the ground state with q = 2), which are far from the Nambu-Goto predictions over much of our range of l.
At moderate values of l √ σ the deviation from Nambu-Goto decreases much more slowly with increasing l than O(1/l 5 ): the dotted line in Fig. 22 represents a very slowly varying ∝ 1/l 1.7 correction to ∆E 2 . Such a behaviour could, in principle, be accommodated by the small-l behaviour of a correction of the form shown in eqn(55). In addition, we have also seen that at larger l the p = 0 0 − state does appear to break away from this slow variation, approaching the Nambu-Goto energy level more rapidly. This can also be naturally accommodated by a correction as in eqn (55), since at larger l the term d n /σl 2 becomes small enough that the power behaviour of the correction returns to its asymptotic value (∝ 1/l 5 in this case). Such an interpretation is however less plausible if we take seriously the apparent crossing of the Nambu-Goto level at our largest value of l. And it is not consistent with the apparent behaviour of the p = 0 excited 0 − state which at larger l appears to turn over and approach the lightest 0 − Nambu-Goto energy level. As discussed in some detail in Section 4.3.3, a much more natural explanation of these states is that they reflect massive mode excitations, and these cannot be easily described in the framework of a simple effective string action.

Conclusions
In this paper we have calculated the energies of ∼ 20 of the lightest excitations of a confining flux tube that is closed around a spatial torus. We have varied the size l of the torus from the very short, close to the critical length l √ σ ∼ 1.6 where the theory undergoes a finitevolume deconfining transition, up to moderately large lengths, l √ σ ∼ 6.3. Although our most extensive calculations have been in the SU(3) gauge theory at a fixed value of the lattice spacing a, we have also performed calculations in SU(5) and SU(6), at the same a, so as to have control over the large-N limit, and in SU(3) at a smaller value of a, so as to have control over the continuum limit. We find no significant a or N dependence of our energy spectra, within our statistical errors, so that our results can be read as being for the continuum gauge theory, and for all N ≥ 3. Our goal has been to learn more about the effective string action that describes confining flux tubes and which can hopefully teach us something about the string theory that may describe the physics of a confining gauge theory, at least at large N. In the last few years there has been startling analytic progress in this direction [4,5,6,7], and one has learned that if one expands the effective string action in powers of 1/l, then the first few terms are universal. This also implies that, to this approximation, the effective action is the same as Nambu-Goto in flat space-time, i.e. a free string theory.
Our calculations are largely complementary to this analytic work. The latter focuses on long flux tubes, where 1/l √ σ is small, and on excited states whose gap to the ground state is small: ∆E(l) ≪ m G where m G is the mass gap of the theory. By contrast our calculations concentrate on small to medium values of l √ σ. This is largely dictated by the fact that most of the systematic errors in a numerical estimate of an energy increase as that energy increases (as discussed in Section 3.2). To control these errors in our calculation, which makes important use of a variational approach, we use a very large basis of lattice operators, as described in Section 3.4. This is the main technical innovation of our calculations. This strategy works reasonably well, but is significantly less successful than in our similar calculations in 2+1 dimensions [9]. For that reason the present calculations should be regarded as exploratory.
A better choice of operators could increase the range and accuracy of such calculations quite dramatically. This is an important direction for future work. One of our main systematic errors does not affect the lightest state in a sector with particular quantum numbers, and most of our calculated states fall into this category. These include ground states with momenta along the flux tube p = 0, 2π/l and 4π/l, and with spin around the flux tube axis of |J| = 0, 1, 2 (a continuum and so not entirely reliable labelling of the content of the relevant representations of the cubic group), and with various parities, as described in Section 3.3. Almost all these quantum numbers entail non-trivial excitations of the flux tube and so are dynamically interesting. (In contrast to transverse momentum which would teach us nothing new and which we therefore ignore.) Our high statistics calculation of the absolute ground state, described in Section 4.1.2, demonstrates quite accurately that the O(1/l) universal Lüscher correction is in the simple bosonic string universality class, and so confirms earlier work. (In a brief aside in Section 4.4, we also confirmed that k = 2 flux tubes appear to belong to the same universality class.) In the case of SU(6), but not SU(3), we obtained some weak numerical evidence that the leading non-universal O(1/l γ ) term is consistent with γ ∈ [3,7]. This is to be compared to the recent analytic result that γ ≥ 7 [7]. It will be hard to do much better than this because one is analysing very small 'Casimir energy' corrections to the ground state energy, ideally on long flux tubes that have a large energy. By contrast, the excited states receive a much larger contribution from the excitations of the string and thus provide a much better arena within which to test ideas about the effective string action.
The most striking feature of our spectrum of excited states is how closely most of these states track the predictions of the Nambu-Goto free string theory, and that they do so all the way down to lengths l √ σ ∼ 2 −3 where the flux tube is presumably not much longer than it is wide. This is evident from the various plots in Section 4.2 and Section 4.3, where we compare our calculated energies to the parameter-free predictions of the free string theory. (The only parameter is σ and that is determined by fitting the absolute ground state.) It is of course natural to ask whether this might not simply be a consequence of the fact that we now know [7] on general grounds that an expansion of the energy E(l) in powers of 1/l must coincide with the corresponding expansion of the Nambu-Goto energy E N G (l) up to at least O(1/l 3 ).
The answer to that question, as illustrated in Fig. 20, is 'no'. Indeed, the agreement with Nambu-Goto persists down to values of l √ σ where the expansion of E N G (l) in powers of 1/l is no longer convergent and one has to use its full, resummed expression. Note that our range of flux tube lengths, while focused on smaller values, does in fact extend into the region where the Nambu-Goto expansion converges (for these states), so that we overlap with the range of l to which the analytic calculations apply. So, just as in our earlier 2+1 dimensional calculations [9], the implication of this result is very clear: the sensible starting point for analysing the effective string action is the Nambu-Goto free string theory and not the conventional Gaussian string action (in static gauge). If we do the latter then we will need to calculate corrections to all orders in 1/l in order to access the shorter values of l where we find that the spectrum still maintains a simple stringy character. If, on the other hand, we do the former, then it appears that something like a modest O(1/l 5 ) correction might be sufficient, as illustrated in Fig. 21, or some resummed version, with a modest leading coefficient, as illustrated in Fig. 22. This is, however, not the whole story. We have also found a few states for which the corrections to Nambu-Goto are very large. These corrections reduce the energy, making the state relatively light -incidentally ensuring that our calculations of these energies are accurate and reliable. These states are of particular interest because they potentially encode the excitation of the massive modes that we expect to be associated with the non-trivial intrinsic structure of the flux tubes. (Something we were unable to detect in our more accurate analysis of fundamental flux tubes in D = 2 + 1.) All of these 'anomalous' states, or at least all those that we have identified so far, turn out to have J Pt = 0 − quantum numbers. (Here P t is the two dimensional parity in the plane transverse to the flux tube axis.) The lightest 0 − flux tube with the lowest non-zero longitudinal momentum, p = 2π/l, remains far from the Nambu-Goto prediction over our whole range of lengths l √ σ ∈ [1.8, 6.3] and, as we see in Fig. 15, it appears to approach the stringy prediction only very slowly. By contrast, the lightest 0 − flux tube in the p = 0 sector shows a very similar behaviour till about l √ σ ∼ 4, but then turns upwards and appears to have crossed the Nambu-Goto energy level in our largest-l calculation at l √ σ ≃ 6.3. At the same time, as we see in Fig. 14, the first excited 0 − , after initially appearing to increase towards the first excited Nambu-Goto level, turns downwards, appearing to approach the Nambu-Goto 0 − ground state as l increases. A natural interpretation of this is that the lightest state is (largely) a massive excitation, while what we took to be the first excited state becomes the 0 − stringy ground state at large l. To go beyond such a qualitative observation requires a good understanding of how to simultaneously describe the massive and massless modes of a flux tube, and this we do not have at present. However, as we saw in Section 4.3.3, a simple ansatz for the massive modes naturally describes the behaviour of the lightest 0 − states for both p = 0 and p = 2π/l with a mass scale that is about one half of the lightest glueball mass. While such a simple picture does not say anything specific about the observed large corrections to the next, heavier 0 − states, it is natural, if these states are (approximately) the lightest massless mode excitations, that they should show the somewhat different behaviour for p = 0 and p = 2π/l that is observed, as a results of the derivative couplings such Goldstone fields must possess. In summary, our result is that the lightest closed flux tube states appear to fall into two distinct classes. The first class includes most of the states and here the spectrum is remarkably well described by the free bosonic string theory. This complements and extends to shorter lengths recent analytic results for very long flux tubes. The second class includes only a few states, all of which (so far) possess 0 − quantum numbers, and here the energies are remarkably far from the predictions of the simple Nambu-Goto model. The lightest of these states, for both zero and non-zero longitudinal momenta, display a behaviour that is consistent with being primarily a massive excitation that does not couple too strongly with the massless modes. It is clear that a significant, but realistic, improvement to our calculations would have the potential of providing convincing answers to some very interesting questions.

Note added:
As the revised version of this paper was being prepared, two further papers have appeared containing interesting new analytic results on the effective string action [33]. A summary of these and other results may be found in a recent talk [34]. From the point of view of our work in this paper, a very interesting result is that the coefficients of all those operators in the full effective action that appear in the expansion of the Nambu-Goto action (the 'scaling 0' operators) are in fact universal and equal to their Nambu-Goto values [34]. This result provides a conceptually clean separation of the full effective action into a piece that is Nambu-Goto and a qualitatively different piece that is the correction to it. Another interesting result is that the leading correction to Nambu-Goto is in fact universal [34]; this arises not from the constraint of open-closed duality described earlier in this paper, but from the constraint that world-sheet conformal symmetry should be maintained at the quantum level -a constraint that is not obvious in static-gauge, but can be seen in the conformal gauge approach of Polchinski-Strominger, and which the Nambu-Goto model fails to satisfy. N β L x = l/a L y × L z × L t sweeps×10 6 num ops 3 6  aE gs (q, l) q |J| P t l/a N = 3 N = 5 l/a N = 3f 0 0 + 10 0.2448 (27)  aE(q = 0, l) |J| P t P l l/a N = 3 N = 5 l/a N = 3f 0 + + 10 0.905(12) 0.945 (13) , l), of the four lightest excited flux tube states with zero longitudinal momentum p = 2πq/l = 0. Quantum numbers as shown. Calculations are labelled as in Table 9.
Calculations are labelled as in Table 9.
aE(q = 1, l) , J =odd l/a N = 3 N = 5 l/a N = 3f 10 1.417 (24) Table 12: Energies of the four lightest excited flux tube states with longitudinal momentum p = 2πq/l = 2π/l and with J =odd. See Table 9 for the lightest odd-J state. Calculations are labelled as in Table 9.
aE(q, l) , J Pt = 0 − q l/a N = 3 N = 5 l/a N = 3f 0 10 1.144 (12) Table 13: Energies of the lightest excited 0 − flux tube states with longitudinal momentum, p = 2πq/l, apart from those states in Table 10 and Table 11 Calculations are labelled as in Table 9.