Spectrum of very excited $\Sigma_g^+$ flux tubes in SU(3) gauge theory

Spectra with full towers of levels are expected due to the quantization of the string vibrations, however different theoretical models exist for the excitation spectra. First principle computations are important to test the different models and to search for novel phenomena, but so far only a few excited states of QCD flux tubes have been studied with pure gauge SU(3) lattice QCD in 3+1 dimensions. We thus aim to study a spectrum of flux tubes with static quark and antiquark sources up to a significant number of excitations. We specialize on the spectrum of the most symmetric case, namely $\Sigma_g^+$, where up to two levels are already published in the literature. To achieve the highest possible excitation level, we construct a large set of operators with the correct symmetry, solve the generalized eigenvalue problem and compare the results of different lattice QCD gauge actions with different lattice spacings and anisotropies.

Spectra with full towers of levels are expected due to the quantization of the string vibrations, however different theoretical models exist for the excitation spectra. First principle computations are important to test the different models and to search for novel phenomena, but so far only a few excited states of QCD flux tubes have been studied with pure gauge SU(3) lattice QCD in 3+1 dimensions. We thus aim to study a spectrum of flux tubes with static quark and antiquark sources up to a significant number of excitations. We specialize on the spectrum of the most symmetric case, namely Σ + g , where up to two levels are already published in the literature. To achieve the highest possible excitation level, we construct a large set of operators with the correct symmetry, solve the generalized eigenvalue problem and compare the results of different lattice QCD gauge actions with different lattice spacings and anisotropies.

I. INTRODUCTION
Understanding the confinement of colour remains a main theoretical problem of modern physics. Its solution could also open the door to other unsolved theoretical problems. An important evidence of confinement, where we may search for relevant details to understand it, is in the QCD flux tubes [1] computed in lattice QCD. Experimentally, flux tubes are suggested by the Regge trajectories [2,3] in the hadron spectrum. String models can account for Regge trajectories, and they also imply the linear confining quark-antiquark potential similar to the one used in quark models [4,5].
Presently, an approximate understanding of flux tubes is quite developed, for example with string models. The dominant behaviour of the string-like flux tubes has a single scale: the string tension σ. The main analytical string model utilised in the literature to explain the behaviour of the QCD flux tubes is the Nambu-Goto bosonic [6,7] string model [8]. It assumes infinitely thin strings, with transverse quantum fluctuations only. The quantum fluctuations predict not only a finite profile width of the groundstate flux tube, increasing with distance [9], but also an infinite tower of quantum excitations [10,11]. However, the string-like behaviour obscures the details of confinement and of other possible hadronic phenomena.
Clearly, at short quark-antiquark distances, the flux tube deviates from the string model. The Nambu-Goto model in four space-time dimensions has an imaginary tachyon [11] at short distances for the groundstate, whereas the QCD flux tube has a real Coulomb potential [11]. At really short distances, lattice QCD has recently shown the potential becomes dominated by perturbative QCD [12].
Another instance where the groundstate flux tube devi- * bicudo@tecnico.ulisboa.pt † nuno.cardoso@tecnico.ulisboa.pt ‡ alireza.sharifian@tecnico.ulisboa.pt ates from the string model is in the flux tube profile. Recently, our lattice QCD collaboration PtQCD [13] studied the zero temperature groundstate flux tube of pure gauge QCD, and found evidence for a penetration length λ [14], as a second scale other than the string tension σ, contributing to the colour fields density profile of the flux tube. This intrinsic width of the flux tube may be the reason why the QCD flux tube is stable in four dimensions, unlike the Nambu-Goto string which is rotational invariant only in 26 dimensions.
There is also an ongoing puzzle in the excited spectrum of mesons, as reported [15] in measurements by the Cristal Barrel detectors [16]: the Regge slope for radial excitations is similar to the one for angular excitations, and this cannot be explained with a quark model [2]. A large degeneracy, larger than the chiral restoration symmetry [17], has been analysed [18][19][20][21]. Possibly there is a new principal quantum number [3]. Notice a principal quantum number is already present in the Nambu-Goto spectrum.
Besides, an infinite tower [22] of these excitations can also be obtained with constituent gluon models, in the denominated hybrid three body quark-gluon-antiquark systems [23,24]. This may possibly be explained by [3] a constituent gluon with an effective mass [25,26], and an effective quark gluon-gluon potential [27,28].
Moreover, in there is evidence for another particle, the string worldsheet axion [29], in the lattice gauge theory spectra of closed strings.
Thus precise first principle theoretical studies are necessary to go beyond the string models and clarify the details of QCD flux tubes, or hybrid. In lattice QCD it is straightforward to study groundstate open flux tubes, using the technique of Wilson loops [1], the simplest gauge invariant correlations in lattice QCD.
There are two classes of flux tubes, the closed flux tubes and the open ones. In lattice QCD, they are respectively studied with torelons and Wilson loops. There has been a serious effort of extracting higher excitations in the closed flux-tube channel [30,31]. The spectrum of the closed flux-tube can be partially approximated by arXiv:2105.12159v2 [hep-lat] 13 Sep 2021 the Nambu-Goto model.
In this paper we study the excited open flux tubes in lattice QCD, in particular we opt to address the puzzling sector of radial excitations. Some excited states have indeed been observed by lattice QCD computations [32][33][34][35][36][37][38], so far compatible with the string dominance of the QCD flux tube. Moreover the structure of excited flux tubes started to be studied as well [39][40][41], where some evidence for a constituent gluon is present in some states. However only the first excitations have so far been studied.
The spectrum of radial excitations of the groundstate is the Σ + g spectrum. We now review the symmetry group of our flux tubes. With two static sources, it is equivalent to the one of the molecular orbitals of homonuclear diatomic molecules. It is the point group denominated D ∞h . We thus utilise the standard quantum number notation of molecular physics, already adopted in the previous studies of QCD flux tube excitations [32-38, 40, 41]. D ∞h has three symmetry sub-groups, and they determine three quantum numbers.
The two-dimensional rotation about the charge axis corresponds to the quantum angular number, projected in the unit vector of the charge axis Λ = |J g ·ê z |. The capital Greek letters Σ, Π, ∆, Φ, Γ . . . indicate as usually states with Λ = 0, 1, 2, 3, 4 . . . , respectively. The notation is reminiscent of the s, p, d · · · waves in atomic physics. In the case of two-dimensional rotations there are only two projections J g ·ê z = ±Λ.
The permutation of the quark and the antiquark static charges is equivalent to a combined operations of charge conjugation and spatial inversion about the origin. Its eigenvalue is denoted by η CP . States with η CP = 1(−1) are denoted by the subscripts g (u), short notation for gerade ( ungerade).
Moreover, there is a third quantum number, different from the phase corresponding to a two-dimensional pwave. Due to the planar, and not three-dimensional, angular momentum there is an additional label for the s-wave Σ states only. Σ states which are even (odd) under the reflection about a plane containing the molecular axis are denoted by a superscript + (−).
With these quantum numbers, the energy levels of the flux tubes are labeled as The states we opt to study are most symmetric system ones, corresponding to the Σ + g spectrum. This is the spectrum where up to two levels have so far been reported in the lattice QCD simulations literature [36,37]. For the other spectra only the groundstate level has been reported in the literature [32][33][34][35][36][37][38]. In principle, the Σ + g is the most amenable for a first study of very excited states.
In Section II we thus detail the different state of the art lattice QCD techniques adequate to study the excited spectrum of the gluonic flux tubes [32][33][34][35][36][37] produced by a static quark-antiquark pair. Working with pure gauge SU(3) fields discretised in a lattice, we utilise Wilson loops with a large basis of gluonic spacelike Wilson lines to include different excitations of the quark-antiquark flux tube. Moreover, we apply different smearing techniques, an improved action, and an anisotropic lattice to improve the signal over noise ratio. We combine our operators, to block diagonalize our basis the angular momentum and parity quantum numbers of the D ∞h point group and project on Σ + g states. We numerically solve the correlation matrix via the Generalized Eigenvalue Problem and compute the corresponding effective masses. We also discuss our computational efficiency. The number of gluonic operators combined with the space points where we compute the flux tube densities turn out to be very large, and we resort to computations in GPUs and to CUDA codes.
In Section III we show, for the excited levels where the signal is clear, the results of our computations for the spectrum. We evaluate how the different techniques improve the results. Finally in Section IV, we analyse our results for the first excitations of the flux tube and search for signals of novel phenomena beyond the Nambu-Goto string model; in Section V we conclude our work.

II. LATTICE QCD FRAMEWORK TO COMPUTE THE Σ + g FLUX TUBES
We discuss our strategy to compute as many radially excited states as we can for the Σ + g flux tubes. We consider a large basis of operators, use the standard action and an improved ones, apply different smearing techniques in the configurations, and use different lattices with space-time anisotropy. A. Our operator basis for the Σ + g excited states Our gauge invariant operators are a collection of Wilson loops, a closed path of Wilson lines. Since we have static charges, our temporal Wilson lines are quite simple, they are straight lines.
We want our spatial operators to have the same symmetry as the flux tubes we are studying. Moreover, since we are using a correlation matrix, we want any linear combination of our operators to also have the same symmetry. Otherwise we could be producing unwanted states. The symmetry is provided by the spacial Wilson lines, who close the Wilson loop and make it gauge invariant.
In the study of the flux tubes, we utilise a basis of spacial Wilson line operators, with the necessary and sufficient operators to produce the Σ + g spectrum and avoid producing the other spectra. As usual we choose our frame such that the charge axis is the z axis, and the origin is set at the midpoint between the quark and the antiquark, with distance R. The x and y axis are in the two perpendicular directions.
We first must choose a basis of operators, consisting of a linear combination of a set of spatial Wilson curves connecting the two static charges, with the same symme-tries as those of Σ + g . Each curve we consider first departs from the antiquark charge in the xy plane perpendicular to the charge. In Fig. 2 we show different such curves. Then the path of the operator is continued with a straight Wilson line in the z direction, and finally it is completed with an opposite curve in the x, y plane to join the quark charge. The operators must be invariant for rotations of angles multiple of π/2 around the charge axis, thus we must have a sum of n op different paths in order to achieve this invariance.
In Fig. 2 we show the denomination of the operator O(l 1 , l 2 ), the distance l between the straight section of the Wilson line and the charge axis, and the number of terms n op in the operator to make it symmetric. For instance the first operator O(0, 0) is the simple strait Wilson line directly connecting the two charges, it has l = 0 and just n op = 1 term. The next operators are deformations of the strait Wilson line. Then the operator O(1, 0) is a staple with l = 1, but we need to add N op = 4 of these staples in different directions to have and invariant operator for rotations of angles multiple of π/2 around the charge axis. The subsequent operator O(1, 1) has its straight z direction line at distance l = √ 2 from the charge axis, but for it to be invariant for the two-dimensional rotation around the charge axis and invariant for the parity inversion about the median point we need to have a sum of all n op = 8 of these possible Wilson curves to get a fully symmetric operator. Continuing along this way we can construct an infinite tower of symmetric operators.
Obviously, we must truncate this sum up to some distance l max of the order of half the spatial length of the lattice since our lattice has periodic boundary conditions. Ideally, we should have l max larger than the finite width of the flux tubes.
However, in contradistinction with the continuum, in a cubic lattice there is a mixing between particular different angular momenta Λ operators. Let us consider for instance the operators such as the O(l 1 , l 2 ), which are a linear combination of sub-operators pointing in four directions separated by right or flat angles of 0, π/2, π, 3π/2. An operator with planar angular momentum Λ must have [42], for each of its spatial links as described for instance in Fig. 2, the phase where ϕ if the azimuthal angle corresponding to the respective spatial link. In the four directions of our operator, any angular momentum multiple of four, has the same phase in all four directions. Thus we must be aware of this problem and do our best to mitigate it. This is further complicated by the phenomenology of flux tubes: in the Nambu-Goto model there is as well a degeneracy between the excited states of these operators.
In particular, we must be careful to prevent the combinations of the symmetric operators O(l 1 , l 2 ) to generate other symmetries. These operators are discrete, they include in general Wilson lines in four to eight directions and all these directions have the same phase.
Having the same phase in all directions is necessary for a fully symmetric Σ + g operator. However, when combining different operators with different directions, a nonsymmetric operator may be generated. For instance O(1, 0) which azimuthal directions in the xy plane are ϕ = 0, π/2, π, 3π/2 and O(1, 1) which azimuthal directions in the xy plane are ϕ = π/4, 3π/4, 5π/4, 7π/4 may be combined with an opposite phase, for instance O(1, 0) − O(1, 1) in the case they are both properly normalized, which would correspond to a Γ g state.
Indeed if these O(l 1 , l 2 ) operators are all included in a correlation matrix, its diagonalization will pick up not only Σ + g states but also states with large angular momentum Λ = 4 about the charge axis. We explicitly verified that, when using a wider basis of operators O(l 1 , l 2 ), we would get more energy levels, some of them nearly degenerate, than expected. Besides, the signal of the excited potentials would be less clear.
We thus restrict our basis of operators to avoid as much as possible states other than Σ + g . We only consider operators which set of Wilson curves have the same directions. It is convenient to use the operators on the directions x and y of the lattice axes. For completeness, we show in Fig. 3 the four sub-operators, spatial Wilson line paths from the antiquark to the quark, used to construct the gauge field operators at Euclidean time t 0 . The inverse Wilson lines are used for the operators at time t. Our basis to construct the correlation matrix consists of the operators of type O(l, 0) only, but considering several different l. It is interesting that using a smaller operator basis and less dense in the space of the flux tube lead to clearer results for the spectrum of excited states as we found out in our computation.
B. Solving the GEVP for the correlation matrix to compute the excited spectra We utilise the correlation matrix C kl (t) to compute the energy levels of the excited states, as done previously in the literature [43][44][45][46]. Now the sub-indices k and l stand for the spacial operators in the operator basis defined in Section II A, denoted O k . The spacial operators, defined in Figs. 1, 2, 3, are connected by temporal Wilson lines L, Notice each matrix element corresponds to an evolution operator in Euclidean space, where all energy levels E i contribute, with coefficients depending on how close the operator is to the actual physical states, with the Euclidean damping factor exp(−E i t).
The first step to compute the energy levels, is to diagonalise the Generalized Eigenvalue Problem for the correlation matrix for each time extent t of the Wilson loop, and get a set of time dependent eigenvalues λ i (t). With the time dependence, we study the effective mass plot and search for clear plateaux consistent with a constant energy E i in intervals t ∈ [t iini , t ifin ] between the initial and final time of the plateau. We have the option of choosing the initial time and we choose t 0 = 0 because it produces the clearest results. Moreover, the results with t 0 = 0 are compatible with the ones obtained with other small values of t 0 . The different energies levels E i , should correspond to the groundstate and excited states of the flux tube. If our operator basis is good enough, then E 0 is extremely close to the groundstate energy, E 1 is very close to the the first excited state, etc.
Moreover, with the diagonalization we also obtain the eigenvector operators [46] corresponding to the groundstate, first excitation, etc. We get a linear combination of our initial operators,  Table I: Our ensembles, for the isotropic Wilson action and the improved anisotropic S II action. ξ is the bare anisotropy in the Lagrangian and ξ R is the renormalized anisotropy. The renormalized anisotropy and the lattice spacings are computed with the prescription of Section II D.
Notice this result must be interpreted with a grain of salt. The eigenvector operators O i do not exactly correspond to the respective state as in quantum mechanics, but they get the clearest possible signal to noise ratio, among our operator basis. The eigenvector operators O i and the respective correlation matrix can be used in the same time interval t ∈ [t iini , t ifin ] ideal for the effective mass plateaux of the energy spectrum.
The number of gluonic operators turns out to be large, requiring a large computer power. We thus write all our codes in CUDA and run them in computer servers with NVIDIA GPUs. Due to the GPU limited memory this requires an intensive use of atomic memory operations. For instance, to compute a 13 × 13 correlation matrix, per configuration, our CUDA code takes approximately 380s for a lattice volume of 24 3 ×96 using a GeForce RTX 2080 Ti with 7.5cc.

C. Gluon actions and configuration ensembles
We compute our results using six different ensembles, defined in Table I.
We use the anisotropic Wilson action [1] computed with plaquettes, where W c = c 1 3 Re Tr(1 − P c ) where s, s runs over spatial links in different positive directions, P s,s denotes the spatial plaquette, P s,t the spatial-temporal plaquettes and ξ is the (unrenormalized) anisotropy.
Moreover, to improve our signal we also resort to the improved anisotropic action S II developed in Ref. [47], with u s = 1 3 Re TrP ss 1/4 , u t = 1. W ss,s and W ss,t , instead of plaquettes, include 2 × 1 rectangles.
So far, the results with more excited states shown in the literature, two states [36,37] and up to three in an unpublished work [48], have been obtained with this action.
We generate five different ensembles of configurations, with the parameters defined of Table I.
The anisotropy is used in order to have a smaller temporal lattice spacing a t . This enables a better estimation on the extraction of excited states as well as a more precise result since we have more time slices for the same time intervals.
Notice an anisotropic action enables us to use larger distances with the same number of spatial lattice points. However, in order to have a good spatial resolution, it is then convenient to use an improved action, and S II uses plaquettes extended up to 1 × 2 rectangular shapes. S II is specially designed to eliminate spurious high-energy states from the gluon spectrum.
Moreover, in order to improve the signal over noise ratio, for the Wilson ensemble, we use the multihit technique in the temporal Wilson lines and the APE smearing in the spatial Wilson lines [14]. The multihit technique, [49,50], replaces each temporal link by its thermal average, Here it is not possible to utilise the extended Multihit technique as defined in Ref. [14], because our operators in the spatial Wilson line have a broader structure. In particular, we use MultiHit with 100 iterations in time followed by APE smearing [51] with α = 0.4 and 20 iterations for the Wilson action without anisotropy. For the S II ensemble and for the Wilson ensembles with anisotropy, we use MultiHit with 100 iterations in time followed by Stout smearing [52] in space with α = 0.15 and 20 iterations. It turns out it is more economical, using GPUs, to perform all our computations on the fly, rather than saving configurations. In general, we first generate a configuration, then apply smearing, then compute the correlation matrix with our full operator basis. What we save to disk is the correlation matrix. Since we generate the computations on the fly, we also list the ensemble Table I the sets of operators used. Notice we may as well turn off smearing, to check the importance of smearing, and we study the results of ensembles W 1 , W 2 , W 4 and S 4 with and without smearing.
In what concerns the efficiency of our codes, for instance using a lattice volume of 24 3 × 96 and the Wilson action, it takes 5s to generate a new configuration. To decorrelate the configurations, we run 50 iterations between the used configurations. Each iteration is composed by 4 heatbath steps followed by 7 overrelaxation steps. It takes 0.8s for 100 iterations of Multihit and 0.014s for 20 steps of Stout smearing in space.

D. Computing the lattice spacing and the renormalized anisotropy
In the case of isotropic actions, there is only one independent scale, arising from dimensional transmutation, since the action is conformal invariant. The physical scale is the string tension σ present in the linear term of the quark-antiquark potential. From its value, we determine the lattice spacing a. In the case of anisotropic actions we have two different lattice spacings, the spatial a s and the temporal a t .
To set the scale of the lattice spacing a of the isotropic Wilson action (with ξ = 1) in physical units, corresponding to the ensemble W 1 , we use the equations fitted in SU(3) pure gauge lattice QCD by Ref. [53], where g is the coupling constant of the Wilson action, and [53] the parameters are b 0 = 0.01364, b 1 = 0.2731, b 2 = −0.01545 and b 3 = 0.01975, valid in the region 5.6 ≤ β ≤ 6.5,â (g) = f (g 2 ) f (g 2 (β = 6.0)) , As for anisotropic actions, the renormalized anisotropy ξ R = a s /a t can be determined as a function of the bare anisotropy ξ. The groundstate potential is computed with the Wilson loop, two different directions for the time direction, either the usual direction 4 now anisotropic, or using one of the isotropic distances, say 3. Comparing the short distance potential with both time directions, which is well determined without any smearing and is fitted with a simple function, the ratio a s /a t is determined.
Once the physical anisotropy ratio ξ R is determined, we determine the physical scale from the string tension. We fit the data for the groundstate potential at large interquark distance R with the linear plus constant ansatz c 0 + c 1 R. We did not find any improvement in the results fitting the data with the ansatz c 0 + c 1 r + c 2 /r. We utilise smearing to be able to fit the large distance potential. The linear term is obtained from fitting the difference for different distances R of the ground-state potential, with a plateau. The linear term depends on the spatial lattice spacing a s , the temporal lattice spacing a t and the string tension σ as, therefore Note that, using the above formulas, all the results presented in this work are in units of √ σ for the potentials and in units of 1/ √ σ for the distances. This is the standard two step technique to determine the scales of our lattice. Besides, there is another option to fit both lattice spacings a s and a t in one step, when we fit the whole excited flux tube spectrum. This will be discussed in Section IV when we will analyse the spectrum.  Here we are interested in checking the effect of using a large basis of operators. If this improves the results, it can be used together with other techniques. Notice in general in the literature smearing is mandatory to improve the results, already for the groundstate quarkantiquark fluxtube. We thus compare the effect of using a large basis of operators to using smearing.
In Fig. 4 we compare the groundstate potential for different ensembles, using just one operator (the standard Wilson loop with O(0, 0)), or 13 sets of operators of the type O(l, 0) up to l = 12. It turns out the use of 13 operators improves the groundstate signal. But to have the best signal for a distance large enough to be physically interesting, we need smearing. According to the plots, smearing maintains the results unchanged for R > 1a s , in particular the string tension error bar is reduced with smearing but it remains consistent with the string tension with no smearing.
In Fig. 5 we use the 13 operators to search for at least one excited state, and we compare the results without and with smearing. We also include the Nambu-Goto spectrum in dashed lines to guide the eye. Indeed the operators are necessary to get the excited sates, but in most cases smearing is necessary as well and produce similar results to Nambu-Goto. Only in ensemble W 4 , the results for the smaller distances R < 4a s can be obtained without any smearing for the first two excited states. Having an excited state with no smearing, similar to the smearing ones and to Nambu-Goto confirms that smearing is not distorting our spectrum. This W 4 result also shows that an anisotropic action may be more effective to study excited states.

B. Degeneracies with O(l1, l2) operators
We now study the results obtained with a more complete set of operators, using the ensembles O 1 and O 2 of Table I. In these simulations we include all 11 possible operators O(l 1 , l 2 ) with l < 4.5, corresponding to the operators with l 1 l 2 = 0, 0; 1, 0; 2, 0; 3, 0; 4, 0; 1, 1; 2, 1; 3, 1; 2, 2; 3, 2; 3, 3. As discussed in Subsection II A this may produce states with different symmetries from Σ + g . For instance a phase difference between operators O(l, 0) and operators O(l 1 , l 2 ) may be selected by the GEVP, and this could produce operators with angular momentum four, of symmetry Γ + g . This is verified in Fig. 6 where we find more states than in Fig. 7, with just the sets of operators O(l, 0) which are less prone to generate unwanted symmetries other than Σ + g . Indeed we find approximately degeneracies among the N = 2, N = 3 states, the N = 4, N = 5, N = 6 and the N = 7, N = 8 states. This makes sense in the Nambu-Goto perspective where there is a principal quantum n = 2 * n r + l, similar to the principal quantum number for an harmonic oscillator vibrating in a twodimensional x − y space.
From this perspective we expect to have the degeneracies, n = 0, 1 state : n r = 0, n = 2, 1 state : n r = 1, n = 4, 3 states: n r = 2 or l = 4, n = 6, 3 states: n r = 3 or n r = 1 and l = 4, n = 8, 5 states: n r = 4 or n r = 2 and l = 4 or l = 8, . . . Indeed we find in Fig. 6 nearly degenerate states N = 2 and N = 3, compatible with n r = 2, l = 0 and n r = 0, l = 4 states. They only differ for small distances, where possibly there are Coulomb contributions beyond the Nambu-Goto model. We also find approximate degeneracies for the higher states.
This confirms as expected that to get clearer and unambiguous Σ + g signals it is preferable to have a smaller class of operators, using only the O(l, 0) sets of operators. Using 13 sets of on-axis operators O(l, 0) and the Wilson action, with no anisotropy, we find already several energy levels, clearly ordered. However, the level N = 3 is a bit out of the remaining patterns, with larger error bars a bit unexpectedly closer to the N = 2 level. This is clear in the top of Fig. 7, ensemble W 1 .
The energy levels also have a similar dependence on R as in the Nambu-Goto model, constant at short distances and linear at large distances.
However, the higher levels seem to be shifted vertically when compared with the Nambu-Goto levels de- picted with dotted lines.

D. Results with the Wilson action and anisotropic lattice
To get potentials at large distances, we first use the anisotropic Wilson action. The results are shown in Fig.  7 centre, both for ξ = 2 and ξ = 4. Then we are able to go not only up to large distances but also to clearly see one more level, going up to N = 6.
The pattern of the separation of the levels is also more striking even than in the isotropic case.
However, the values of the higher potentials get distorted at shorter distances, perhaps because the anisotropy somehow enables the degeneracy already discussed in Sections II A and III B to set in.
Nevertheless, the higher levels remain undistorted for R > 4a s .

E. Results with the improved anisotropic action
To solve the problem at smaller distances for the anisotropic lattices, we resort to the improved action. Indeed the short distance is improved although it remains distorted, as seen in Fig. 7 bottom, ensemble S 4 .
Nevertheless we are able to get two more levels, N = 7 and 8. This is the highest level we are able to get.
We also compare the results of all actions in Fig. 8, where we study one potential level in each panel. The excited levels tend to be higher than the Nambu-Goto model.
To conclude on the different spectra, we are able to get several levels, with on-axis operators, with smearing, and with anisotropic lattices. But we should discard the first four smaller distances for the anisotropic action because the levels above N = 2 are distorted and tend to get degenerate.

IV. ANALYSIS OF OUR RESULTS
The simplest description of the QCD excited flux tubes -with charges in the triplet representation of SU(3) -is given by the bosonic string model, based on the Nambu-Goto action [6,7], which spectrum for an open string with ends fixed with Dirichlet boundary conditions is the Arvis potential [11], In Eq. (24), D is the dimension of the space time D = 4, and n = 2 n r +l is the principal quantum number. In our case of Σ + g , l = 0 and the only quantum number we have is n = 2 n r where n r is the order of the radial excitation. Because of this simple analytical form, we opt to fit the excited spectrum with the Arvis potential. Any deviation may indicate other phenomena, say a constituent gluon.
Note that, for the groundstate, the Arvis potential is tachyonic at small distances since the argument of the square root is negative, moreover rotational invariance is only achieved for D = 26. Nevertheless the first two terms in the 1/R expansion, including the Coulomb term, are more general that the Arvis potential, since they fit the D = 3 and D = 4 lattice data quite well beyond the tachyonic distance. The Coulomb term is independent of the string tension σ and for the physical D = 3+1 has the value − π 12 1 R . This is the Lüscher term [57]. The energy spectrum of a static quark-antiquark and of its flux tube is certainly well defined (not tachionic) and this was the first evidence of flux tube vibrations found in lattice field theory.
In what concerns the groundstate, is well known, as detailed in Section I that the Nambu-Goto fails, because we have no tachyon at short distances. However, the leading terms in a large distance expansion are accurate not only at long distances, but also at medium range short distances, where the tachyon is replaced by the Lüscher [57] coulombic potential, which confirms the factor n − D−2

24
in the Arvis potential. At very short distances the correct potential matches perturbative QCD [12].
Nevertheless, for the excited states the intrinsic width of the flux tube [14] should become negligible compared with the quantum vibrations of the string. The Nambu-Goto model should then be adequate to analyse the potentials we compute with lattice QCD. Notice the Arvis potential produces a tower of levels, as we observe in our lattice QCD data. We now start analysing our data, shown in Figs. 7 and 8, where they are compared explicitly with the Nambu-Goto model.
If we exclude the smaller four distances showing some energy degeneracy, evident in Fig. 7 in the compression of the higher levels with N > 2, we see no evidence for the Γ g states, discussed in Subsection II A, degenerate with our Σ + g states. Notice the degeneracy in these smaller four distances is possibly due to higher harmonics of the string vibrations. The harmonics with an even number of nodes have the same quantum numbers Σ + g as the fundamental harmonic. Since our operators only have two nodes, at the open ends of the string, they should not have a large overlap with the higher harmonics, except at shorter distances where the square shape of our operators may be too crude to only overlap with the fundamental harmonic. Since we are only interested in studying the radial excitations, From now on, we will exclude the four shortest distances from the fits of our data.
Moreover, we see no evidence for string breaking [58][59][60], noticing with the anisotropic actions we are able to go to distances much larger than the string breaking distance. Thus we are confident to be analysing the pure hybrid Σ + g states. The first two excited states seem to be in general compatible with the Nambu-Goto model within error bars. However, the next states depart from it, apparently the energy levels are in general higher than in Nambu-Goto.
Thus, we analyse quantitatively the difference to the Nambu-Goto potential, fitting the data with an extra parameter. We multiply the constant term in the square root by an independent prefactor. The simplest way we find to do this is to change σ in to σ 2 in the denominator or the constant term depending on n r . The modified Nambu-Goto model is, In particular we fit the results of Fig.  7 for each excited state separately with E 2N (R)/(Rσ) to 1 + 2π σ2R 2 2N − 1 12 and extract σ 2 . The resulting fits of all the energy levels for the different ensembles are shown in Fig. 9, and detailed in Table II.  Besides, we also use another approach with a gobal fit of all levels. Instead of fitting each excited state separately, we perform a non-linear multifit for all the excited states in each ensemble in Fig. 10. For each emsemble, we display the resulted σ 2 and χ 2 /dof. We find that for ensembles W 1 , W 4 and S 4 , σ 2 is up to 10 % smaller than σ (meaning that the energy levels are higher than in the Nambu-Goto model), but for W 2 and σ 2 are similar.
Furthermore, since the results in Fig. 7 are in units of the string tension obtained from the groundstate fit, we also perform an even more global fit, fitting all parameters to the excited spectrum (and not fitting the groundstate). In Fig. 11 and in Table III we depict the non-linear multifit for all the excited states in lattice In this case what we have to compare with the previous σ 2 /σ is now σ 2 /ξ R /σ 1 . Again this goes down, up to 10% in ensembles W 4 and S 4 , while in ensembles W 1 and W 2 this is of order of 1. This departure of up to 10 % from the Nambu-Goto model, may possibly be interpreted as the string being non-homogeneous, or be interpreted with the existence of a constituent gluon in the more excited states.

V. CONCLUSIONS AND OUTLOOK
We computed the potentials for several new excitations of the pure SU (3) flux tubes produced by two static 3 and3 sources, specializing in the radial excitations of the representation Σ + g . We succeeded in obtaining the spectrum of several new excited flux tubes, using a large basis of operators , employing the computational techniques with GPUs of Ref. [40] (improved to be able to study field densities), and utilizing different actions with smearing and anisotropy. Previously in the literature and websites, states up to N = 2 have been shown, and we go up to N = 8.
In general the excited states of the the Σ + g flux tubes are comparable to the Nambu-Goto string model with transverse modes, only depending on the string tension σ and the radial quantum number n r .
Nevertheless we find evidence that the Σ + g flux tubes cannot be exactly modelled by the Nambu-Goto string model. A second parameter, we parametrize as a second σ 2 up to 10 % smaller than σ, essentially corresponding to a larger energy splitting between levels than in the Nambu-Goto model leads to better fits. Interesting for the theoretical studies of the QCD flux tubes, this may indicate the formation of a non-homogenous string with a lump or a constituent gluon, or of an extra symmetry higher in the spectrum as referred in Section I.
To clarify in more detail this small tension with the Nambu-Goto model, we would need to have more computational power, to be able to use larger lattices and more operators. Another direction of research is on investigating higher angular excitations of the flux tube (here we investigated radial excitations), but a new approach would be necessary to overcome the limitations of a cubic lattice. Searching in different quantum numbers for explicit evidences of a constituent gluon [41] or of an axion [29] is also motivating. We leave this for future works.   Fig. 11 extracted from a nonlinear multifit of all the energy levels, in the interval [5a s − 12a s ]. The groundstate was excluded from the multifit. σ is the linear term of the groundstate fitted alone.