Revisiting the flux tube spectrum of 3d SU(2) lattice gauge theory

We perform a high precision measurement of the spectrum of the QCD flux tube in three-dimensional $\SU(2)$ gauge theory at multiple lattice spacings. We compare the results at large $q\bar{q}$ separations $R$ to the spectrum predicted by the effective string theory, including the leading order boundary term with a non-universal coefficient. We find qualitative agreement with the predictions from the leading order Nambu-Goto string theory down to small values of $R$, while, at the same time, observing the predicted splitting of the second excited state due to the boundary term. On fine lattices and at large $R$ we observe slight deviations from the EST predictions for the first excited state.


Introduction
While the microscopic origin of confinement in Quantum Chromodynamics (QCD) remains elusive due to the lack of analytic methods to solve QCD at small energies, the formation of a region of strong chromomagnetic flux, a flux tube, between quark and antiquark provides a heuristic explanation for quark confinement. Strong evidence for flux tube formation has been found in numerous simulations of lattice QCD, both in the quenched approximation, e.g. [1], and in simulations with dynamical fermions, e.g. [2] (for reviews and more detailed lists of references see [3,4]).
For large quark-antiquark distances R the flux tube is expected to resemble a thin energy string, so that its dynamics will be governed by an effective (bosonic) string theory (EST). Here and in the following we neglect the effects of dynamical quarks, which enable the formation of a light quark-antiquark pair from the vacuum and a breaking of the string. Following two decades of tremendous progress, a number of features of the EST, including the spectrum up to O(R −5 ) [5,6,7], 1 are by now rather well understood. See Refs. [8,4] for recent reviews. 2 Most of the parameters in the EST are constrained by Lorentz symmetry and take universal values. The first non-universal parameter in the EST, denoted asb 2 (in its dimensionless version), appears at order O(R −4 ) and is the coefficient of a boundary term.
The spectrum of the flux tube excitations can be computed in numerical simulations of pure gauge theories, where dynamical fermions are absent. Typically good agreement between the EST predictions and the lattice results has been observed down to qq separations where the EST is expected to break down (for a compilation of results see [4]). Since the non-universal coefficients appear in subleading terms, their extraction requires very accurate results for the energy levels up to large values of R. So far sufficient accuracy has only been achieved for the coefficientb 2 in three dimensions (3d) [15,16,17,18,19], where the simulations are at least an order of magnitude less expensive than in the 4d case. First hints for a non-vanishingb 2 in 4d SU(3) gauge theory, however, have recently been obtained both from the groundstate energy -the static potential -and the flux tube profile at non-zero temperature [20,21]. Despite a parametric suppression, the most precise computation ofb 2 comes from the static potential, which can be obtained with much higher accuracy than the excitation spectrum. The consistency between the boundary coefficient extracted from the static potential [18] and the excited states has so far only been checked for a comparably coarse lattice spacing [15,19]. In this article, we will extend and improve on our initial studies of the excitation spectrum in 3d SU (2) gauge theory, Refs [22,23,24,15]. In particular, we aim at investigating the continuum approach of the excited states at large R, which has not been possible with sufficient accuracy previously, and compare our results to the EST predictions using the EST parameters extracted in Ref. [18].
Currently, the main challenge from the theory side concerns the inclusion of corrections due to the vorticity of the flux tube, possibly showing up as massive modes on the worldsheet (e.g. [25,26,27,28]; see also the more detailed discussion and list of references in Ref. [18]).Candidate states with contributions from massive modes have been seen in 4d SU(N ) gauge theories [29,30,31,32]. For closed strings the results for such an anomalous state are in good agreement with the contribution of a massive pseudoscalar particle on the worldsheet, known as the worldsheet axion [33,34] (see also [35]). Similar anomalous states do not appear in 3d, which might be related to the absence of the topological coupling term and, consequently, less sensitivity of the energy levels to massive modes. In 3d, the presence of massive modes leads to an additional term which mixes with the boundary correction term at O(R −4 ) and, consequently, impacts the numerical result forb 2 . Since the 3d massive mode contributions to the excited state energies are unknown, we cannot perform a direct comparison of the cases with and without massive modes as done for the static potential in [18] and leave this for future studies.

Excited states of the QCD flux tube
In this section we will focus on the extraction of the excited states of the flux tube in 3d SU(2) gauge theory for multiple lattice spacings. The comparison with the EST is left for the next section.

Computation of Wilson loops for excited states
While for the extraction of the static potential the ideal observables are Polyakov loop correlation functions (e.g. [36,37,18]), spatio-temporal Wilson loops of extents R and T , including operators coupling to particular quantum number channels on the spatial lines, are suitable observables for the extraction of excited states. The main difficulty is the need for large loops for which the signal-to-noise ratio decreases exponentially with the area of the loops. At a given value of the qq distance R, the reliable high-precision extraction of the energy state demands the sufficient suppression and control of the contributions of excited states to the Wilson loop expectation value (cf. eq. (2)). In the first comprehensive study, large correlation matrices including highly optimized operators in combination with smearing and anisotropic lattices have been used to extract the excitation spectrum in different gauge theories and three and four dimensions [30,38,31]. After careful diagonalisation, excited state contributions are suppressed in the individual eigenvalue correlators, 3 so that the convergence to the T → ∞ limit is enhanced for the low-lying states. An alternative strategy has been pursued in [22,23,24], using an improved version of the Lüscher-Weisz multilevel algorithm [40] for error reduction to reliably extract large loops for a small number of spatial operators and temporal extents. The residual excited state contaminations have been removed using suitable fit functions.
Both strategies give reliable results for intermediate values of R, but for large values of R the contamination due to excited states becomes more and more severe, since the energy gaps to the excited states decrease. Consequently, there is doubt about the sufficient suppression of the excited contributions in this regime. A combination of the two methods, using the improved multilevel algorithm in combination with a large set of operators, correlation matrices and fits including excited state contributions, has been used in Ref. [15], leading to accurate results for the excited states in 3d SU(2) gauge theory at large values of R, albeit at a single and comparably large lattice spacing. In this study we use this strategy, further improving on the analysis by employing improved fits and more temporal extents in the analysis, to extract the spectrum at smaller lattice spacings.
In 3d the string energy levels can be classified by the quantum numbers of charge conjugation and parity (C, P ). The individual combinations are denoted as channels. Using suitable sets of spatial operators S j i , the correlation matrices with respect to these quantum numbers can be transformed to a block-diagonal form using the linear combinations To extract excited states in a given (C, P ) channel, we use the 8 different operator sets introduced in Ref. [15]. Setting up a generalized eigenvalue problem for these correlation matrices for the reliable extraction of the eigenvalues in the limit T → ∞ is problematic, since the correlation matrices might be illconditioned already for the smallest temporal extent available. As in Ref. [15], we thus diagonalize the correlation matrices for each value of T separately utilizing the QR reduction method. The resulting eigenvalues are denoted as λ CP n (R, T ) with n = 0, . . . , 7. To check that we are identifying the right eigenvalues with increasing T , we look at the associated eigenvectors. Those belonging to the same states for different T should be almost parallel, whereas those belonging to different states perpendicular. In the following we will mostly focus on the groundstates in the individual channels, for which the identification of the eigenvalues is unambiguous. The only excited state which we consider is the first excited state in the (+, +)-channel, for which the identification becomes more difficult at smaller lattice spacings. In contrast to the eigenvalues obtained from a generalized eigenvalue problem, the resulting eigenvalues might include contaminations from other states in the channel, albeit with overlaps which are strongly suppressed for T → ∞.

Removal of excited state contaminations
The eigenvalues of the correlation matrices generically obey the spectral representation [41] Here E CP n are the energies in the (C, P )-channel, are the energy differences, β CP n the overlap of the eigenstate with energy state n and α CP k,n (R) the ratio of overlaps of the eigenstate with energy states k and n. At finite T , it depends on the size of the energy gaps and the reduction of overlaps which state dominantes the sum on the right hand side of eq. (2).
At large T and with sufficient suppression of contributions from other states, the energies can be extracted using the asymptotic T → ∞ formula In practice, however, contaminations from other states are not negligible due to the high accuracy for the eigenvalues, so that eq. (4) cannot be used to extract the energies reliably. To remove contaminations from other states we perform a simultaneous fit to the results for the eigenvalues for all available T a and T b to the leading order formula (see also [22,24,15]) Here T a < T b , α CP n (R) ≡ α CP n+1,n (R) and δ CP n (R) is the energy gap to the closest state (or the one with the dominant contribution to the sum on the right-hand-side of eq. (2)) in the channel. In most of the cases these fits lead to accurate results with acceptable control over the systematic effects. In some cases, however, these fits become unstable with respect to statistical fluctuations of the eigenvalues. This is particularly true for large values of R and higher excited states. To further constrain the fits, we include the data for the individual eigenvalues in the global fit, leading to the additional relations where γ CP n (R) = ln β CP n (R) is an additional fit parameter. Despite this additional parameter, the joint fits using eqs. (5) and (6) are generically more stable. A particular example for such a fit is shown in Fig. 5 in appendix B. To control the systematics in these fits we implement further checks which are also described in appendix B.
The energy differences ∆E CP ;CP nm (R) from eq. (3) can be extracted independently from the total energies, so that they might serve as independent crosschecks. Using Eq. (2), one can derive the analogue to eq. (5) for the energy differences where α = α CP ;CP nm (R) is a suitable combination of the overlaps and δ = δ CP ;CP nm (R) corresponds to the energy gap to the closest state in the CP channel for eigenvalue m (which is equivalent to the gap in the CP channel for eigenvalue n to leading order in 1/R). In analogy to the fits for the total energies we can improve and stabilize the fits by including the analogue of eq. (6) for the energy differences in the global fit, . In the analysis we considered fits including all loops, fits for which we left out the data from loops with the smallest temporal extents, which was beneficial when for the smallest temporal extents the contaminations from higher excited states were still sizeable, and, in very few cases, fits where the data from the largest temporal extent has been excluded. The latter was beneficial when the results from the largest temporal extents showed large fluctuations, but has only been considered when none of the other fits obeyed the quality criteria of appendix B. In the following we will only use results for which the removal of the excited states has been successful, i.e. those energies and differences for which at least one of the aforementioned fits passed the checks of appendix B. We estimated the systematic uncertainty associated with the removal of the excited states from the difference of the best fit result -being the one with data from the maximal number of temporal extents included for which the fit passed the checks -with the result from a fit where the data from the loop with the smallest temporal extent of the best fit has been excluded, whenever the latter fit was considered to be reliable.

Results for the flux tube spectrum
We perform simulations of 3d SU (2) gauge theory employing the Wilson plaquette action with the common mixture of heatbath [42] and three overrelaxation [43] updates. The simulation parameters together with the parameters of the multilevel algorithm are collected in Tab. 1. We set the scale using the Sommer parameter r 0 [44], which has been determined for the present parameters with high accuracy in Ref. [18]. The energies include a lattice spacing dependent additive normalization. We get rid of this normalization by subtracting the normalization constant V 0 obtained in Ref. [18] from fitting the potential to the EST prediction. Both r 0 and V 0 are listed together with other EST parameters in Tab. 2. Note, that the ensembles at β = 5.0 have already been analysed in Ref. [15]. Here we reanalyse the data using the improved fits described above.
We show the results for the energies in Fig. 1. Ancillary files including the jackknife bins of the energies and energy differences are provided along with this paper. Unfortunately, the removal of contaminations from other states Table 1 Parameters of the simulations. t sub s,t are the temporal extents of the sublattices in the LW algorithm and N sub s,t the number of sublattice updates (see appendix A). N skip indicates every which point has been used in a timeslice for the evaluation of Wilson loops, for the purpose of memory reduction in the simulations. The parameter set at β = 5.0 is the one already discussed in Ref. [15]. 192 × 96 2 6 36000 1200 Table 2 Results for the Sommer parameter r 0 , the string tension σ, the normalization constant V 0 and the boundary coefficientb 2 for different β values and in the continuum limit from Ref. [18]. Forb 2 the first uncertainty is purely statistical, the following systematic uncertainties are associated with the unknown higher order correction terms in the EST, the choice of fitrange in the extraction ofb 2 and, for the continuum value, the continuum extrapolation.

Comparison to the EST spectrum
The energy levels and differences discussed in the previous section can now be compared to the predictions of the EST.

EST predictions for the spectrum
We first discuss the EST predictions for the energy levels of the flux tube relevant for this study. For an in-depth discussion of the EST properties we refer the reader to the reviews [8,4] and the extended discussion of section 2 in [18]. Using the action up to 6 derivatives order and the constraints for the coefficients [45,46,16,47] from Lorentz invariance, the spectrum of a flux tube (open string) of length R to O(R −5 ) is given by [5,6] Table 3 String states in the lowest four energy levels, their representation in terms of creation operators αm in the Fock space and the values for the coefficients B l n and C l n . Note, that in three dimensions C l n always vanishes.
Following the arguments from [32,15,48,34], the first term on the right-handside is the full light-cone spectrum [49] b 2 is the dimensionless leading-order boundary coefficient and B l n and C l n are dimensionless coefficients tabulated in table 3 for the lowest few string states. The B and C coefficients depend on the representation of the state with respect to rotations around the string axis and lift the degeneracies of the light-cone spectrum. The lowest order correction term to eq. (9) is expected to appear with an exponent ξ = 6 if the next correction originates from another boundary term.
The effective string theory is expected to break down for √ σR 1, where the energy of the degrees of freedom reaches the QCD scale. The EST does not account for several QCD processes, which are allowed generically in the microscopic theory. Among them are glueball emission and (virtual) exchange, as well as inner excitations of the flux tube. The latter are expected to appear as massive excitations on the worldsheet and are not included in the standard form of the EST. For the static potential, rigidity or basic massive mode contributions can be included in the EST analysis. They contaminate the extraction of the boundary coefficientb 2 and also add an additional term to the potential (see sections 3 and 5 in Ref. [18]). For excited states, the explicit form of such rigidity or massive mode corrections has not been computed so far, so that we cannot test the presence of such corrections in this analysis.

Comparison of the data to the EST predictions
To enable the visibility of small differences at large R, we from now on plot rescaled energies and differences following (see also [24]) for which the expansion of the LC spectrum, eq. (10), to O(R −1 ) yields E rsc n = n and ∆E rsc nm = n − m. In Figs. 2 and 3 we show the rescaled energies in comparison to the predictions of the EST including continuum parameters. Note, that differences in the LC curves for the different lattice spacings would not be visible in the figures. For the full EST prediction, containing also the boundary term, we also include the curve with the parameters of β = 5.0. For the other β-values the full EST curves lie between the β = 5.0 and continuum curves, but are very close to the latter. Generically, both for n = 1 and 2 the data qualitatively follows the LC curves down to very small values of R, where the EST is no longer expected to describe the flux tube dynamics. Quantitatively, small deviations are visible, which, however, appear to remain more or less constant with R. Note, that this is to some extend an artifact of the rescaling, eq. (11), which includes subtractions and a multiplication with R. Consequently, a decreasing deviation, like the one in Fig. 1, might appear constant in Figs. 2 to 4.
For n = 1, both the β = 5.0 and 7.5 data become consistent with the full EST prediction around R ≈ 2.5r 0 (please note the difference in normalization and the reanalysis for β = 5.0 as compared to the results presented in Ref. [19]), even though the β = 7.5 data shows a slight tendency to overshoot the curve for R 3r 0 . For β = 10.0 this trend continues and the data lies above the continuum curve for R 2.5r 0 . This might hint to deviations from the curve in the continuum limit, but it could also be an artifact of insufficient removal of the excited states contributions which become more severe in this region. In case of the former, it could be a sign for a higher order corrections, or massive mode contributions, needed to describe the energies accurately at While conclusions concerning the agreement with the full EST from the first excited state in the (+, +) channel are difficult, the difference between this state and the groundstate in the (−, −) channel might be more sensitive to the boundary term. Generically ths difference might be useful to investigate the form of correction terms in the EST as long as they lift the Nambu-Goto degeneracies, since all universal terms belonging to the n = 2 states cancel. Unfortunately, this difference is also very difficult to compute with control over the systematic effects and so far we have only been able to obtain reliable results at β = 5.0. Those results, obtained from the individual analysis of the energy differences, are shown in Fig. 4. In this normalization the difference appears almost constant with R while in fact it decreases with R −1 . It is unclear whether it eventually approaches the EST prediction for larger R due to the large uncertainties in this region. We note, that the R the LC prediction could well be a sign for a massive mode being responsible for this difference to be non-vanishing. At the same time it could be a combined contribution of higher order corrections which mimics such a R −1 correction.

Conclusions
We have extracted the spectrum of the open flux tube in 3d SU(2) pure gauge theory up to the second excited state for multiple lattice spacings. The combination of the multilevel algorithm and a variational method allowed for precise results up to qq distances of about 4r 0 . Excited state contaminations have been removed using a sophisticated fitting procedure with several checks for systematic effects. The results qualitatively follow the energy levels of the Nambu-Goto string theory in the light cone quantization, eq. (10), down to small values of R where the EST is not expected to provide a valid description of the flux tube excitations. Quantitatively, however, we observe deviations which we compare to the predictions of the full EST, including a boundary term on top of the Nambu-Goto action with coefficients computed in Ref. [18]. We observe that lattice artifacts are small in general, confirming the findings from Ref. [24]. However, some lattice artifacts might be visible for the first excited state at large values of R.
While the results tend to agree with the EST predictions, in particular the results for the second exited state show the expected splitting predicted by the EST, we observe some deviations for the first excited state at smaller lattice spacings. This could be a sign for higher order or massive mode corrections becoming important in the continuum limit, a generic disagreement with the predictions at large R, or uncontrolled systematic effects. To verify either of these scenarios, further and more accurate results at large R are needed. For n = 2 the data apparently tends to approach the full EST predictions asymptotically for all lattice spacings, even though the approach is slower on the finer lattices. A particularly interesting quantity with respect to correction terms to Nambu-Goto energy levels in the EST is the difference between the first excited state in the (+, +) and the groundstate in the (−, −) channel, which is vanishing for the LC spectrum. Results with the current precisionwe were only able to extract results on our coarsest lattice -decrease with R −1 . It is, however, difficult to judge whether the results will first converge to the boundary correction at large R or quantitatively disagree with this term.
Despite the drastic increase in precision, results for larger values of R and better precision for the higher excited states are needed to fully confirm or falsify the agreement between spectrum and EST predictions for the excited states. So far we could not observe any unambiguous discrepancy between data and EST, despite the fact that some deviations become apparent on the finer lattices. However, these could still be remnants of systematic effects, which, generically, become harder to control for larger values of R. Of particular importance for future studies is the inclusion of possible corrections due to massive modes in the EST predictions, which in 3d so far have not been computed within the EST framework.
Acknowledgements This paper is dedicated to the memory of Pushan Majumdar, a dear colleague and friend. He introduced me to the world of programming and high-perfomance computing in science, QCD in particular, and was the co-supervisor of my diploma thesis on QCD strings. In the years after I visited him several times in India and whenever we met we had lively discussions on all possible topics connected with physics. QCD strings were one of his favorite topics and the extension and improvement of our initial study [24], reported in this article, was always something which was on his mind.
The simulations have been done in parts on the Athene cluster at the University of Regensburg and the FUCHS cluster at the Center for Scientific Computing, University of Frankfurt. I am indebted to the institutes for offering these facilities. During this work I have received support from DFG via SFB/TRR 55 and the Emmy Noether Programme EN 1064/2-1.

A Error reduction for Wilson loops
The extraction of the flux tube spectrum requires the accurate computation of large Wilson loops including non-trivial spatial gluonic operators. For error reduction we use the variant of the Lüscher-Weisz multilevel algorithm [40] discussed in [23,24] (see also [50]). The multilevel algorithm exploits the locality properties of the Wilson plaquette action to perform intermediate averages of parts of the operators located in the interior of sublattices, separated by time-slices with fixed spatial links. While in the original application of the algorithm to Wilson loops the spatial operators have been put on the boundaries of the sublattices [40,22], in the improved algorithm the spatial operators are located in the middle of a sublattice.
The Wilson loop expectation value can be split into spatial and temporal sublattice operators, L α i (x 0 ) and T(t). The spatial sublattice operator consists of the spatial operator S α i ( x, x 0 ) from eq. (1), located in the middle of the spatial sublattice of extent t sub s , and two-link operators T (x 0 ) of spatial extent R in direction with unit vector i, connecting the spatial operator with the upper boundary of the sublattice. With L † we denote the operator which includes the spatial operator [S α i ( x, x 0 )] † , but connects with the lower boundary of the sublattice through two-link operators (in an abuse of the notation †). The product of two-link operators is defined by where we have used the sum convention for indices appearing twice on one side of the equation. The spatial sublattice operator for the spatial sublattice with lower boundary at time coordinate x 0 is then given by The temporal sublattice operator for a temporal sublattice of extent t sub t with lower boundary at time coordinate x 0 consists of the multiplication of two-link operators from the lower to the upper boundary, Using these two types of sublattice operators and denoting sublattice averages with {·}, we can decompose a Wilson loop of extents T and R as Here the temporal extents of the sublattices have to fulfill T = t sub s + k · t sub t with k ∈ N and we have made use of the fact that the temporal sublattice operators obey the two-link operator multiplication law, eq. (13).
The algorithm contains several parameters that can be tuned to achieve optimal error reduction. For the sublattices including the temporal sublattice operators the parameters are the sublattice extent t sub t and the number of updates N sub t . Both can be tuned following the lines of Ref. [51]. The optimal number of sublattice updates increases with the size of the loops and with decreasing lattice spacing we found it beneficial to increase t sub t . For the sublattices including the spatial sublattice operators we similarly have to tune t sub s and N sub s . As for the other sublattices we found it beneficial to increase the former with decreasing lattice spacing. For the excited states, large values of N sub s are beneficial, as described in Ref. [24]. However, when considering 8 different operator sets including longer contours the computational cost for the computation of the operators the associate sublattice averaging typically contributes more than 90% of the overall computational cost, so that N sub s cannot become overly large. To make efficient use of the full temporal extent of the lattice, we vary the temporal lattice extent for different loops. Since all temporal extents are comparably large, we do not expect to see relicts of this in the data. Since the algorithm is inherently memory consuming, we compute the Wilson loops only on a fraction of the points on a given timeslice for the larger lattices at smaller lattice spacings. When going to the finer lattices we expect this to not affect uncertainties significantly, since neighbouring points become more and more correlated.

B Control of the excited state fits
The extraction of the spectrum heavily relies on the fits used to remove the contaminations due to excited states, discussed in Sec. 2.2. A particular example for such a fit (employing the quality criteria discussed below) is shown in Fig. 5, where we plot the input effective energies E eff defined by the left-hand-side of eq. (5), together with the curves obtained from the fit to the central values. Good control of the systematics of these fits is essential to obtain β = 7.5; (+, −); R = 15a reliable results. χ 2 /dof typically shows acceptable values even if the fit misses some of the points at large temporal extents, which have large uncertainties but are the most important ones concerning the extrapolation. In addition to χ 2 /dof we thus install additional quality criteria and constraints. We first constrain the fitparameters so that their values will be in the physically relevant regime. For the fits to extract the energies, eqs. (5) and (6), the relevant parameters are the energies E, the overlap ratio α, the gap to the first excited state δ and the logarithm of the overlap γ (here all indices are suppressed). Both E and δ should be positive and the latter of the order of the energy differences between the lowest energy levels. If δ was an order of magnitude larger we discarded the fit. The ratio of overlaps α is expected to be a number of O(0 − 100) and not much larger. For the extraction of the energy differences ∆E, eqs. (7) and (8), the relevant parameters are α, δ and γ. For α and δ similar criteria as for α and δ apply, whereas the energy difference ∆E can not be expected to be positive. In most of the cases they should be, but we also consider differences, in particular the difference ∆E ++,−− 10 , which are expected to be negative. In addition to these constraints, we also apply the additional quality criteria introduced in Ref. [24]. In particular, we compare the excited state contribution from the fit parameters and for eqs. (5) and (6), respectively (similar for the energy differences with α → α, δ → δ and γ → γ), to the actual difference of the asymptotic energy with the effective energy, and for eqs. (5) and (6), respectively (for the energy differences E → ∆E and the second terms on the r.h.s. are replaced by the terms of the l.h.s. of eqs. (7) and (8)). For each fit, we plot the results for ∆ versus ∆ together with the expectation ∆ = ∆. While we allow deviations for the two to three largest values of ∆ and ∆, we only keep those fits for which the other values agree with the ∆ = ∆ line within uncertainties and do not show a systematic trend away from this line. Example plots for acceptable and non-acceptable fits have been shown in Ref. [24].