Spectrum of the open QCD flux tube and its effective string description I: 3d static potential in SU(N=2,3)

We perform a high precision measurement of the static $q\bar{q}$ potential in three-dimensional SU($N$) gauge theory with $N=2,3$ and compare the results to the potential obtained from the effective string theory. In particular, we show that the exponent of the leading order correction in $1/R$ is 4, as predicted, and obtain accurate results for the continuum limits of the string tension and the non-universal boundary coefficient $\bar{b}_2$, including an extensive analysis of all types of systematic uncertainties. We find that the magnitude of $\bar{b}_2$ decreases with increasing $N$, leading to the possibility of a vanishing $\bar{b}_2$ in the large $N$ limit. In the standard form of the effective string theory possible massive modes and the presence of a rigidity term are usually not considered, even though they might give a contribution to the energy levels. To investigate the effect of these terms, we perform a second analysis, including these contributions. We find that the associated expression for the potential also provides a good description of the data. The resulting continuum values for $\bar{b}_2$ are about a factor of 2 smaller than in the standard analysis, due to contaminations from an additional $1/R^4$ term. However, $\bar{b}_2$ shows a similar decrease in magnitude with increasing $N$. In the course of this extended analysis we also obtain continuum results for the masses appearing in the additional terms and we find that they are around twice as large as the square root of the string tension in the continuum and compatible between SU(2) and SU(3) gauge theory. In the follow up papers we will extend our investigations to the large $N$ limit and excited states of the open flux tube.


Introduction
Flux tube formation between quark q and antiquarkq provides a possible mechanism for quark confinement and has been verified in simulations of lattice QCD (e.g. [1]). While the microscopic origin of the formation of the flux tube is still debated, it is common consensus that the long distance properties, in particular the spectrum, are well described by an effective string theory (EST). The EST is a two dimensional effective field theory for the Goldstone bosons associated with the breaking of translational symmetry, the quantised transversal oscillation modes. While the basic properties of the theory are known for a long time [2][3][4][5][6] a number of features have only been elucidated in the past decade [7][8][9][10][11][12][13][14][15][16][17][18]. In particular, it has been clarified [14] that the leading order spectrum agrees with the light cone quantisation [19] of the Nambu-Goto string theory (LC spectrum, eq. (2.7)) and the corrections to the LC spectrum have been computed up to O(R −5 ) [9,11,16], where R is the distance between quark and antiquark. For more details and a recent review see [20]. The predictions of the EST can be compared to lattice results for the excitation spectrum of the flux tube and good agreement has typically been found, proceeding down to qq separations where the EST is not expected to be reliable (for a compilation of results see [20]). In particular, the EST predicts a boundary correction of O(R −4 ) with a free and, presumably, non-universal coefficientb 2 . This coefficient has first been extracted from the excited states in 3d SU (2) [21] and later from the groundstate in 3d Z 2 [12] and SU(N = 2, 3) [22] gauge theory for a number of different lattice spacings. Indeed, strong evidence for a non-universal behaviour has been found. The good agreement down to small values of R is surprising, given that the flux tube profile only agrees with the EST starting from around 1.0 fm, at least [1,[23][24][25]. A possible explanation could be that the flux tube consists of a solid, vortex-like core whose fluctuations are governed by the EST. This would explain the good agreement of the profile with the associated exponential decay, e.g. [25][26][27], rather then with the Gaussian profile predicted by the EST [28,29], while leaving the spectrum untouched to leading order. 1 In fact, it has been found [24] that the profile can be analysed using a convolution of the EST and vortex profiles, which would indicate that the flux tube shares features of both.
One can expect that a particular class of corrections to the EST might show up as massive modes on the worldsheet. Indeed, candidates for states receiving contributions from massive modes have been seen in 4d SU(N ) gauge theories [36][37][38][39] and the results for closed strings [39] are in good agreement with the energy levels obtained by including a massive pseudoscalar particle on the worldsheet known as the worldsheet axion [17,40]. Recently, the behaviour of the mass of the worldsheet axion for N → ∞ has been investigated [41] and a finite result has been found in this limit. It is interesting to note that possible massive modes are only seen in 4d, which might be due to the fact that the topological coupling term associated with the worldsheet axion only exists for d > 3 (see section 2.3). Consequentially, any massive modes in 3d can only couple via terms contributing to higher orders in the derivative expansion or via vertices including more than one massive field. In this case results at intermediate distances are potentially less affected and the massive modes appear as quasi-free modes on the worldsheet.
Apart from the contribution of massive modes, there could also be contributions from rigidity. The associated term satisfies the symmetry constraints of the theory and naturally turns up in cases where the EST can be derived from a vortex solution of the underlying field theory [30][31][32][33][34][35]. While the contributions due to this term start at higher orders in the derivative expansion of the EST (and thus should be negligible up to O(R −7 ), at least), it has been found, using zeta function regularisation, that its non-perturbative (in the sense of the 1/R expansion) contribution cannot be expanded systematically in R −1 . It thus leads to contamination at all orders which are, however, exponentially suppressed for large values of R [18,[42][43][44]. In fact, the leading order contribution is formally equivalent to the contribution of a free massive particle, as discussed in section 2.3.
In this series of papers we will investigate in detail the agreement of the spectrum of the open flux tube in SU(N ) gauge theories with the predictions from the EST. Particular emphasis lies on the reliable extraction of the boundary coefficientb 2 from the static potential and its continuum extrapolation. In this context it is also important to control the higher order corrections. For the case that the other corrections are regular corrections that are part of the derivative expansion, the next term would be of O(R −γ ) with γ ≥ 6. Such corrections are comparably easy to disentangle from the boundary term. If, on the other hand, there are contributions from the rigidity term, or equivalently from massive modes, those can contaminate the extraction ofb 2 [18,43,44]. We will see that we cannot exclude the latter possibility, but we can show that in a certain range of distances the exponent of the correction term with respect to the LC spectrum is indeed 4 and it is not well described by the additional terms alone. Since we cannot exclude the contamination ofb 2 , we carry out two independent analyses, excluding and including the rigidity/massive mode contributions. Eventually we are interested in the large N limit of the non-universal contributions. In particular, there is hope that some of the non-universal parameters can be used to constrain the possible holographic backgrounds and eventually help to find the holographic dual of large N Yang-Mills theory.
In this first paper of the series we will introduce the methods that we use to analyse the potential in SU(N ) gauge theory. However, we focus only on results for N = 2 and 3 in three dimensions and leave the extraction of the results for N > 3 and the extrapolation N → ∞ for the next paper in the series. In follow up papers we plan to consider excited states, extending the studies from [21,45], and the four dimensional theory. The paper is organised as follows: In the next section we briefly discuss the relevant predictions and properties of the EST and its limitations. The consequent section 3 is devoted to the basic analysis of the lattice data for the static potential. In particular, we perform a reliable extraction of the Sommer parameter and the string tension and show the leading correction to the LC spectrum is indeed of O(R −4 ). In sections 4 and 5 we extractb 2 excluding/including massive mode contributions, respectively. Section 6 provides the summary and discussion of the main results of the two analyses and we conclude in section 7. The details of the lattice simulations are presented in appendix A and we discuss the control of systematic uncertainties in appendix B.

Effective string theory and its limitations
The EST is the low energy effective field theory (EFT) describing a single, stable, noninteracting flux tube. Here we will focus on the case of a flux tube stretched between two sources in the fundamental representation (quark and antiquark), but one can also study the hypothetical case of a closed flux tube wrapping around a compactified dimension. The presence of the flux tube breaks translational invariance in the transverse directions, leading to d − 2 Goldstone bosons (GBs), the quantised transverse oscillation modes. For brevity we will consider the string of length R to reside in the (x 0 , x 1 )-plane (x 0 being the temporal direction) with endpoints located at x 1 = 0 and R. Then one can parametrise the fluctuation field by X µ (x α ) = (x α , X i (x α )), where α = 0, 1, i = 2, . . . , d and the X i are the GBs of the EST. The action up to 6 derivatives order has been derived in [8,10,12] together with the constraints for the coefficients (the result for the coefficient c 4 has been clarified in [14]). Here we will write down the action in the static gauge (instead of the diffeomorphism invariant form discussed in [15,20]) and in Minkowski space for simplicity. Since we are considering an open string the action consists of two parts, where S c is the bulk action, where the integral runs over the worldsheet M of the flux tube, and S b the boundary action where the integral runs over the woldsheet boundary ∂M. The coefficients in eqs. (2.2) and (2.3) are constrained by the residual symmetries of the effective field theory, which for the EST is Lorentz symmetry. The results for the coefficients are [8,10,12,14] while b 2 remains unconstrained and thus represents the first subleading low energy constant not fixed in terms of σ. In the next section we will discuss the spectrum which follows from the effective action. Note, that b 2 is a dimensionfull quantity, so that, for the purpose of lattice simulations and the associated continuum limit, it makes sense to introduce the dimensionless couplingb which we will use from now on.
On top of the terms in eq. (2.1), there is one more class of terms allowed within the EST, constructed from powers of the extrinsic curvature [15,18]. The leading order term is known as the rigidity term and has first been proposed by Polyakov [46]. It also turns up in those cases where vortex solutions could be constructed from the underlying microscopic theory [30][31][32][33][34][35]. In this sense, the presence of rigidity terms can be seen as evidence of vortex contributions to the EST energies. Rigidity contributions to energy levels start at higher orders in the perturbative expansion, but it has been found that it leads to a nonperturbative contributions to the energy levels in zeta-function regularisation [18,47,48]. However, as pointed out in [13], this regularisation scheme does not preserve the non-linear Lorentz symmetry of the EST, leading to possible counterterms that need to be taken into account. It is thus mandatory to crosscheck this result using a regularisation scheme which preserves Lorentz invariance. In this paper we will use the result as presented in [18], which we summarise and discuss in more detail in section 2.3.
The EST is expected to break down at the point where the energy of the fluctuation modes reaches the QCD scale Λ QCD ≈ √ σ. The energy of the modes is of the order 1/R, meaning that the EST is expected to break down for Owing to the derivative expansion in eqs. (2.2) and (2.3) this makes sense, since each derivative is of the order of one unit of momentum of the degrees of freedom, ∂ ∼ p ∼ 1/R, so that the EST corresponds to an expansion in ( √ σR) −1 which is expected to break down when eq. (2.6) is fulfilled. On top of this, there are also several processes that are allowed in the microscopic theory, but are not accounted for within the EST. Among them is the emission of glueballs. If those are on-shell, the emitting state has to be an excited state with E ≥ E + m G , where E is the state it decays to and m G is the mass of the lightest glueball with appropriate quantum numbers. Consequently, such a process can only appear for excited states with E n > E 0 + m G , where E 0 is the groundstate. On top of this on-shell emission, there can also be virtual glueball exchange. This process will always be present (at finite N ) and leads to corrections of unknown size. These meson-glueball interactions and mixings are suppressed with increasing N and vanish in the N = ∞ limit (e.g. [49]). Furthermore, within the EST the string is not allowed to develop knots or to intersect with itself, which would correspond to handles on the worldsheet. In addition, the flux tube is likely to have an intrinsic width, allowing for inner excitations, contributing in the form of massive excitations on the worldsheet which have to be added to the EST. We will discuss these additional modes and their possible contributions in section 2.3. In general, the EST only describes a single non-interacting and stable string. We would like to emphasise that this is not the case realised for the flux tube in QCD, which can break due to the creation of a light qq pair from the vacuum. Furthermore, the quarks are taken to be static, meaning that one takes the limit of infinitely heavy quarks. The effect of finite quark masses can be included using non-relativistic effective field theories of QCD [50][51][52][53], for instance (see [54] for a recent evaluation of the 1/m 2 corrections to the potential using the EST predictions). In this paper, however, we will focus on pure gauge theory, where effects owing to finite quark masses are absent and the EST holds up to the limitations discussed above.
We close this section by noting that one possibility for an extension of the standard EST to include the effect of possible internal excitation or external influences on the flux tube is to view them as 'defects' on the worldsheet [55,56]. Here a defect is a place on the string where the derivatives of the coordinates are discontinuous. The associated effect has been worked out for the potential and found to be in reasonable agreement with observed anomalous states in four dimensions [56].

Spectrum of the open string and boundary corrections
We will now discuss the spectrum that follows from the action (2.1). The constraints in eq. (2.4) for the coefficients set c 2 to c 6 to the values they obtain in the NG theory. Consequently, one can expect that the leading order spectrum is equivalent to the NG one, up to the corrections due to the boundary term proportional tob 2 . This expectation is corroborated by the formulation of the EST in diffeomorphism invariant form [15], where the full NG action appears as the leading order term in the action. The main question thus concerns the spectrum following from the NG action in d dimensions. For d = 26 and d = 3 the exact spectrum is known from the light cone (LC) quantisation and takes the form [19] Note, that we only consider open strings, so that n denotes the number of phonon excitations and as such labels the excitation level. The first few excited states in terms of phonon creation operators α i −m (for phonon momentum m) are listed in table 1. However, the LC quantisation breaks Lorentz invariance explicitly so that, away from d = 26 and d = 3, counterterms are necessary for a consistent quantisation. The first counterterm is proportional to the c 4 term in eq. (2.2) with a coefficient [13,17] Thus, the first correction in the bulk action to the LC spectrum appears at order R −5 , but, since the NG action only contains one free parameter, the correction is universal. One might thus wonder whether the square root formula in eq. (2.7) should be taken as the starting point of the expansion, or whether it is rather its expansion in orders of R −1 . In fact, this has been one of the biggest puzzles associated with the EST in the past decades. Lattice data (see [20] for a compilation) show excellent agreement with the full square root formula even down to values of R where its expansion breaks down. It thus has been conjectured that it should be the full formula which provides a reasonable starting point (e.g. [21,39]). The discussion above agrees with this conjecture in the sense that the LC spectrum can schematically be written as [13] Table 1. String states in the lowest three energy levels and their representation in terms of creation operators α m in the Fock space together with their representation with respect to SO(d − 2). We also list the associated values for the numbers B l n and C l n appearing in eq. (2.10). 'sym. tracel.' for the state |2, 3 stands for 'symmetric traceless'.
Consequently, it will always be the full square root formula which appears when we solve the above equation for E NG n . Furthermore, an analysis using the machinery of the Thermodynamic Bethe Ansatz (TBA) shows, that the leading order S-matrix is integrable, leading to energies given by the full square root formula [13,17]. Note, that the boundary term can also be included in this TBA analysis [16].
Following the above discussion and the corrections to the LC spectrum computed in [9,11] the spectrum up to O(R −5 ) is thus given by where we have inserted the dimensionless couplingb 2 introduced in the previous section. B l n and C l n are dimensionless coefficients tabulated in table 1. They depend on the representation of the state with respect to rotations around the string axis, transformations of X i with elements of SO(d − 2), and thus lift the degeneracies of the LC spectrum. Note once more, that C l n vanishes identically for d = 3 since the state |2, 3 does not exist in this case. This can also be seen from the associated term in the EST action, the c 4 term in eq. (2.2), which is trivial for d = 3. The next correction term to the LC energies can be expected to appear with an exponent ξ = 6 or ξ = 7, depending on whether the next correction originates from another boundary term (which will generically be the case if the associated coefficient does not vanish identically due to symmetry) or a bulk term.

Beyond the standard EST: massive modes and rigidity term
As mentioned already in the introduction and in section 2.1, there is also another class of terms allowed within the EST containing the extrinsic curvature. The leading order term of this type is known as the rigidity term, whose presence has first been found in [42]. In terms of the derivative expansion the contributions of the rigidity term start at 8th derivative order, but the presence of this term in the action can give a non-perturbative contribution to the potential [18,57,58]. Including the first two terms from the bulk action, eq. (2.2), together with the leading order contribution from the rigidity term and evaluating the resulting Gaussian integral leads to a Euclidean potential of the form (for the details see [18]) where T is a finite temporal extent of the spacetime, ∆ is the 2d Laplace operator, and we have introduced the mass m = σ/2α .
The first determinant is the one resulting from a free massless boson field, leading to the Coulombic term within the EST, first computed by Lüscher [59], while the second is reminiscent of the one of a free boson field with mass parameter m. At this point it is important to stress that the contribution of the rigidity term resembles that of a massive excitation on the world sheet in Euclidean spacetime. However, the above result is only the leading order result, originating from the Gaussian part of the rigidity action. Higher order contributions will potentially spoil this equivalence. Nonetheless, at this order the two types of contributions cannot be disentangled.
Using zeta-function regularisation, the determinant of the massive boson in eq. (2.11) leads to a leading order term in the potential of the form [18,57,58] which appears on top of the leading order Coulomb term originating from the massless determinant. Note, that the result has been obtained in zeta-function regularisation which breaks Lorentz invariance explicitly. Owing to the discussion in section 2.1 this means that one has to care for the inclusion of potential higher order counterterms in the action to obtain the correct result. The higher order terms in the action can be included perturbatively [18,43], leading to the additional term [43,47] V rig which contaminates the R −4 boundary term from eq. (2.10). Naively we expect a similar term in the presence of a free boson on the worldsheet.
To assess the effect of the term V rig 0 it is instructive to consider the two different limits with respect to the EST. First, let us consider the large R limit, so that mR 1. In this case the dominating contribution comes from the term in the sum with k = 1, which gives a leading order contribution of the form [18] Consequently, V rig 0 will be exponentially suppressed with respect to all other terms in the EST, so that it can only give a relevant contribution at intermediate or small values of R for a given mass m. Next, let us consider the region of small mR < π. In this case one obtains [18] (keeping only the terms relevant for R → 0) Thus, for small distances the rigidity term leads to another Coulombic term, doubling the standard Lüscher term within the EST. What is particularly important is the fact that in some intermediate regime (which, for not too large values of m, will still be within the validity region of mR < π expansion) the Coulombic term will dominate over the R −4 correction (while the logarithm is already negative). Consequently, we expect a negative Coulombic correction to the LC energy levels in some regime before seeing an R −4 increase for R → 0 if theb 2 term is absent.
As a remark, we would like to stress that the above computation only applies to the potential. The modification of the excited energy levels due to the rigidity term are still unknown. For the purpose of this paper we will only need the potential, but a distinction between the boundary corrections and the rigidity term is notoriously difficult in this case. It would be very interesting to get a result for the associated modification of excited states. The hope is, that for the excited states the predictions including the massive modes is incompatible with the splittings between the energy levels, so that one can distinguish between the cases with or without massive mode contributions. In particular, it has been found in [21] that the splitting of the first few excitations is well described by the boundary term, which could be different when massive modes are included.
In the presence of additional massive bosonic degrees of freedom the presence of terms proportional to the extrinsic curvature also allows for a topological coupling term between the boson field and the GBs, proportional to the mode number [17,40] 2 . Consequently, since the boson couples to the worldsheet analogue of the θ-term, the boson can naturally be referred to as the worldsheet axion [40]. The action for the worldsheet axion can then be written as [17] where ∇ is the covariant derivative with respect to the induced metric h αβ = ∂ α X∂ β X.
To leading order, i.e. up to coupling terms of the form φ∂ α ∂ β φh αβ and higher orders, ∇ α can be replaced by ∂ α . The associated leading order spectrum has been computed with the TBA method for closed strings in [17,40] and has been found to be consistent with states showing an anomalously slow approach to the LC energies in the spectrum of the closed flux tube for d = 4 [39]. Note that the topological coupling term in eq. (2.17) is only present for d > 3, which could be a reason why no anomalous states have been observed in the 3d flux tube spectrum. In 3d massive modes appear as free bosons up to the coupling terms mentioned above. It would be interesting to include those coupling terms in a computation of the energy levels to check their influence on the spectrum. In the following we will always denote the contributions discussed in this section as the contributions from "massive modes" for brevity. We would like to emphasise, however, that strictly speaking we can only compare our results to those obtained from the rigidity term, or from the leading order contribution, neglecting any direct couplings, of a massive boson on the worldsheet. Whether direct coupling terms are indeed negligible is an open question, which we cannot answer in the course of this study.

EST and gauge/gravity duality
We close the discussion of the EST with the remark that the EST action can potentially be computed (e.g. [8] and references therein) from the original 10d fundamental string theory appearing in a generalisation of the AdS/CFT correspondence for large N gauge theories [60][61][62]. For recent computations of properties of the flux tube in this framework see [63][64][65][66]. Consequently, constraining the possible terms and their couplings in the EST action could allow to constrain the fundamental string theory relevant for Yang-Mills theory.
In particular, in such cases where the background in the fundamental string theory is weakly curved, the additional bosonic and fermionic degrees of freedom (e.g. the additional coordinate fields in the directions transverse to the gauge theory plane and the supersymmetry partners of the bosonic fields) can be integrated out perturbatively. This has been done for closed strings and a special class of backgrounds in [8] and the resulting terms have been found to agree with the terms appearing in the EST bulk action, eq. (2.2). Furthermore, in [10] it has been shown that the same class of backgrounds considered for an open string ending on two infinitely stretched D-branes leads to the boundary term proportional tob 2 in eq. (2.3). In terms of the masses of the additional bosons its contribution is given by [10] where ξ labels the transverse directions to the gauge theory plane, m b ξ is the mass of the boson field associated with the direction ξ, b f 2 is the contribution from the fermionic degrees of freedom and BC(ξ) depends on the boundary conditions for the fields in this direction. In the case of Dirichlet boundary conditions BC = 0, while for Neumann boundary conditions one has BC = 1. The ellipses in eq. (2.18) stand for terms originating from possible other fields present on the gravity side of the duality. We would like to emphasize that, even though cancellations can appear, there is no reason why the terms in eq. (2.18) should add up to zero, except for certain fine-tuned situations (see also the discussion in [10]), so that, generically, one can expectb 2 to be non-vanishing when the EST originates from a duality with this types of string background. It would be interesting to see whether the duality can also account for the rigidity term in the EST action and make a statement about its coefficient.
The aforementioned computations have neglected the appearance of possible additional massless scalar modes on the worldsheet. 3 If present, they would contribute to the Lüscher (Coulomb) term and thus inherently change the large R behaviour of the confining string. Since the Lüscher term is well reproduced by lattice data, this situation is basically ruled out. This issue can be resolved if the action for these additional degrees of freedom nonperturbatively develops a mass gap, so that these additional fields contribute as massive degrees of freedom on the worldsheet.  Listed are the range of qq separations, the temporal extent of the LW sublattices t s , both in units of the lattice spacing, the number of sublattice updates n t and the number of measurements. We also list the temporal lattice extent in units of the Sommer scale r 0 .

Results for the potential
To extract the static potential we have used the spatial Polyakov loop correlation function. The details of the simulations and the extraction of the potential from the correlation function are discussed in appendix A, where we also discuss the Lüscher-Weisz (LW) multilevel algorithm [67], which is mandatory to achieve the precision needed for the extraction of the subleading coefficients in the EST. For the suppression of excited state contaminations and finite size effects our simulations have been done on large lattices. The parameters of the simulations are tabulated in table 2. To demonstrate that, even for the high precision achieved, the aforementioned effects are indeed negligible, we report on extensive checks in appendix B.

Scale setting and string tension
For scale setting purposes we use the Sommer parameter r 0 [68] (see appendix A). For SU(3) and d = 4 the associated continuum value is r 0 = 0.5 fm which can be used to convert to "physical units". We extract r 0 /a from the force using four different methods: (a) a numerical polynomial interpolation; (b) a numerical rational interpolation; (c) a parameterisation of the form [68] for the values of R corresponding to the four nearest neighbours of r 0 (motivated by the LO EST); (d) the parameterization of (3.1) with f 2 = 0 for the two nearest neighbours of r 0 .
As final estimate for r 0 /a we will use method (d). The systematic uncertainty associated with the interpolation to obtain r 0 /a can be estimated from the maximal deviation of r 0 /a obtained from methods (a), (b), and (c) compared to method (d). The results for r 0 /a are tabulated in table 3.
Another option for scale setting is to use the string tension σ, associated with the R → ∞ asymptotics of the potential. We extract the string tension in two different ways: (i) we fit the data for the force to the form following eq. (A.4); (ii) we fit the potential to the leading order EST prediction, eq. (2.7) for n = 0, adding a normalisation constant V 0 .
Both ansätze are correct only up to a certain order in the 1/R expansion, so that, in the region where higher orders become important, we will get incorrect results for the string tension. To isolate the asymptotic linear behaviour of the potential we investigate the dependence on the minimal value of R included in the fit, R min . The strength of using these particular two ansätze is, that the corrections to the fit formula are different, so that we can determine the region where higher order terms are negligible by comparing the results. In the region where both sets of results, in dependence on R min , show the onset of a plateau and agree within errors, the estimate for σ will be reliable (with the present accuracy). The results obtained from the two methods, denoted as σ (i) and σ (ii) , respectively, are shown in figure 1 for SU(2) (left) and SU(3) (right). The plots indicate that in most of the cases the extraction of σ is reliable. The most critical case is β = 11.0 for SU(3), whose value of σ we exclude from the following analysis. Another critical case is β = 16.0 for SU (2), where the two plateaus for σ (i) and σ (ii) do not agree within errors. In that case we use the extent where the two results are closest (the discrepancy is of order 1σ). The resulting values for σ are listed in table 3 together with the other fit parameters and are indicated by the bands in the figures. The plots indicate that within the region where the two fits agree the particular choice for R min does not matter within the given uncertainties. In the following analysis we will use σ (ii) , whose value is shown as the red band in figure 1.

The static potential
In figure 3 we show the results for the static potential, rescaled via In this rescaled form the leading order potential to O(1/R) in the 1/R expansion is normalised to 0, so that small differences become visible. We have rescaled the energies using the string tension for each individual value of β. The solid lines in the plot correspond to the LC spectrum, eq. (2.7), and are evaluated using the continuum extrapolated string tension. For SU(2) gauge theory we have not plotted the potential for the two largest lattice spacings for which the results look similar.
Corrections to the LC spectrum start at around R/r 0 ≈ 2 (i.e. around 1 fm in 4d physical units) and are positive, except for one case, namely β = 11.0 for SU (3). This means that the dominant corrections are not expected to be due to the rigidity term from section 2.3. In that case one would expect to obtain a negative correction for intermediate values of R/r 0 . However, this does not rule out the presence of the rigidity term in general. It can still be a subleading correction for the lattice spacings considered. The only exception is the SU(3) β = 11.0 lattice (with a ≈ 0.15 fm -in physical units defined via r 0 ≡ 0.5 fm), where the negative correction could be due to a dominant rigidity term. In all of the cases we observe that the magnitude of the corrections becomes stronger when we approach the continuum limit. In particular, they are larger for SU(2) gauge theory than for the SU(3) case.

Isolating the leading order correction terms
On top of the leading order behaviour (the LC spectrum from eq. (2.7)), the EST predicts corrections starting at O(R −4 ). We would now like to test whether this prediction is true. The first step is to isolate the leading order correction to eq. (2.7). If eq. (2.7) is incorrect, corrections will appear with an exponent 0 < m < 4, if the EST predictions are correct andb 2 = 0, we expect m = 4, if corrections only appear at higher orders we will obtain m > 4. If eq. (2.7) is correct to all orders we should obtain m ≈ 0.
To investigate the power of the leading order correction we fit the data to the form where E LC 0 (R) is the potential obtained from the LC energy levels (2.7) and η and m, together with σ and V 0 , are fit parameters. The results for the exponent m versus the   (2), we see that the result for m is typically a bit smaller than 4 (around 3.6). This could indicate that we observe a mixing of two types of corrections where the one of O(R −4 ) is the dominant one. For R min /r 0 > 1 we cannot resolve the correction terms reliably, so that m is not determined sufficiently.

Analysis of boundary corrections without massive modes
In this section we will now turn towards the extraction of the boundary coefficientb 2 . In this first part of the analysis we will neglect contributions from massive modes (or the rigidity term) and extract the value ofb 2 in this setup. We will see howb 2 changes in the presence of these terms in the next section.

Extraction of the boundary coefficient
To extractb 2 we use the groundstate energy from eq. (2.10). To check for the impact of higher order correction terms our general fit formula is of the form where we have included appropriate powers of σ multiplying the higher order terms to keep the coefficients dimensionless. In practice, we perform five different types of fits: A we use the string tension and V 0 extracted from method (ii) in the determination of the string tension from section 3.1 and use eq. (4.1) withb 2 , γ From the analysis in the previous section, we expect the correction terms to be relevant starting with R/r 0 ≈ 1.2. To test the region in which the data is well described by the higher order terms we perform the fits for several values of the lower cut in R and check at which value for R min we get a good description, indicated by acceptable values for χ 2 /dof. For the final result, i.e. for the R min in the final fit, we pick the second smallest distance for which the fit gave an acceptable χ 2 /dof< 1.5. We use the results with minimal R value of R min ± 1a to estimate the systematic uncertainty of the extraction of the fit parameters for this particular fit. For the fits including terms of O(R −6 ) or O(R −7 ) we would expect that these work well down to even smaller values of R min than the fits including only the R −4 term, since these fits include higher order correction terms which should improve the agreement with the data.
The results of the fits are listed in table 5 for SU(2) gauge theory and in table 6 for the SU(3) case. Let us discuss the individual fits before moving on with the analysis. In general, the fits lead to small values of χ 2 /dof (since we picked the second possible value for R min ), so that it is difficult to make statements about the agreement of the EST predictions with the data based on these numbers. Accordingly, our main argument to judge the agreement will be based on the values of R min which can be used in the fits. At coarse lattice spacings all fits work equally well, leading to similar of R min . When going to smaller lattice spacings, however, the picture changes slightly. In these cases the worst fits are typically fit B (which is not surprising since higher order terms have been neglected) and fits A and E. For fit A this might be due to the fact that the fit does not allow σ and V 0 to change, which, apparently, is to restrictive, even though both do not change significantly for the other fits. For fit E R min needs to be larger even though higher order terms are included in the fit, indicating less agreement with the data. This implies that the scenario withb 2 = 0 is disfavoured, in agreement with the analysis from section 3.3. For the following analysis we will thus include the results from fits C and D together with the one from fit B, for which the larger value of R min has been expected. We note, however, that we cannot fully exclude the possibility that the coefficient of the R −4 term is 0 and that the result from section 3.3 is due to a fine-tuned combination of higher order terms.
We determine the associated result forb 2 on a single lattice spacing using the average over the results from fits B to D, weighted with the associated uncertainties. To determine the systematic error due to the particular choice for R min we repeat the same procedure with R min ± 1a and take the maximal deviation as an estimate. The results are listed in table 7. In the next section we discuss the associated continuum extrapolation. It is interesting to note that in SU(3) gauge theory the systematic errors are typically an order of magnitude smaller than in the SU(2) gauge theory. A possible reason could be, that the   Table 6. Results of the fits for the extraction ofb 2 for SU(3) gauge theory. agreement between the effective string theory predictions and the energy levels becomes better with increasing N . This appears natural in the light of possible corrections to the EST discussed in section 2.1, which are expected to be further suppressed with increasing N .
The results forb 2 are plotted versus a 2 in figure 5. The plot displays the rather smooth behaviour towards the continuum (a = 0). The exception is the data point at β = 11.0 for gauge group SU(3), which shows a rather strong upwards trend. We have excluded this data point from the following analysis, since it appears to lie outside of the scaling region.  Table 7. Final results forb 2 for the individual lattice spacings. The first error is the statistical uncertainty, the second the systematic one due to the unknown correction terms, estimated by calculating the maximal deviations of the results forb 2 from fits B to C, and the third is the systematic one associated with the choice for R min . Using this parameterisation we now perform three different fits:

Continuum limit of boundary corrections
(1) a fit including all terms in eq. (4.2) and all lattice spacings (except for β = 11.0 for SU(3)); (2) fit (1) but with bb 2 ,2 = 0; As discussed above, we have excluded the data set with β = 11.0. In all cases we have included systematic errors due to the functional form and the particular value for R min   (1) to (3) (see text). The first error is the statistical uncertainty, the second one the one associated with the unknown higher order correction terms in the potential and the third one the one due to the choice for R min .
used for the extraction ofb 2 by performing these fits for the results from fits B to D and the fits with a minimal R value of R min ± 1a individually.
To determine the final result, we have averaged the results weighted with the individual uncertainties of the extrapolations for fits B to C and estimated the systematic uncertainty due to the choice for R min by doing the same for R min ± 1a. The procedure is the same as the one used for the averaging at the individual lattice spacings. In this way we obtain a final result for the three different fits (1) to (3). The continuum results forb 2 for the different fits are given in table 8. As can be seen from figure 5, the data is consistent with a straight line, which is also indicated by good values for χ 2 /dof, below but close to 1, for fit (2), even though they still might show a slight curvature. As our final result we thus use the result from fit (2), including the spread of the results from fits (1) and (3) as a systematic uncertainty associated with the continuum extrapolation. The associated curves from fit (2) with the main value for R min are shown in figure 6.
The final continuum results forb 2 , given in eq. (6.3) in section 6, are already shown in figure 5. The figure indicates that the results forb 2 in SU(3) gauge theory are significantly larger than the results for SU (2), indicating the non-universality of the boundary corrections in the EST. In particular, this difference remains in the continuum limit. From the tendency between SU(2) and SU(3), one might think thatb 2 tends to zero for N → ∞. This will be investigated further in the next publication of the series.

Analysis of boundary corrections with massive modes
Up to now we have ignored possible contributions from massive modes. Our aim in this section is to see how the presence of these modes changes the result forb 2 from the previous section.

Testing the consistency with the potential
To include the massive modes in the analysis we include the terms from eqs. (2.13) and (2.14) in the fit function eq. (4.1), The main difficulty when fitting to eq. (5.1) is the presence of the infinite sum over the modified Bessel functions of the second kind. In practice, however, the sum is always completely dominated by the first few terms, since K 1 (nc) decays exponentially with increasing The last fit constitutes a check whether the presence of the R −4 term due to the massive mode (i.e. the term from eq. (2.14)) is already sufficient to describe the R −4 correction found in section 3.3.
The fit results are tabulated in tables 9 and 10. The first thing we note is that fit J needs a value of R min which is a factor between 1.5 and 2 larger than R min for the other fits. In particular, R min is typically larger than those values of R where we found the R −4 correction term to be dominant in section 3.3. This basically rules out the possibility that the boundary term is absent, in which case the full R −4 correction would be given Fit √ σr 0 aV 0b2 · 10 2 r 0 m γ (1)/(2) 0 · 10 3 χ 2 /dof R min /r 0 β = 5.0 F 1.2320 (2) Table 9. Results of the fits for the extraction ofb 2 and m for SU(2) gauge theory.
by the correction term from the massive modes. For the other fits the values of R min are comparable, but typically a bit smaller than those of the fits in the previous section. This is not surprising given the fact that each fit has an additional free parameter compared to fits B to D from the previous section. When looking at the result forb 2 and m we see that the fits F to H always agree within uncertainties, since the uncertainties of the fits G and H are large. The large uncertainties are most likely due to the fact that the available data does not allow to constrain the additional fit parameter in these fits compared to fit  Table 10. Results of the fits for the extraction ofb 2 and m for SU(3) gauge theory.
F. In the following analysis we will thus only use the results from fit F. The uncertainties of the parameters from the other fits do not allow for any conclusions about the values of the parameters and continuum extrapolations. In this case our results for the individual lattice spacings are given by the values for fit F in tables 9 and 10 and they only have two uncertainties, namely the statistical one and the systematic uncertainty associated with the particular choice for R min (estimated in the usual way). Unfortunately we cannot investigate the systematic uncertainty due to unknown higher order correction terms, but we believe that the presence of these terms is only a minor effect if massive modes are present.
The results forb 2 and r 0 m are plotted versus a 2 in figure 7. In comparison to the results forb 2 from the previous section (figure 5), the results for SU (2) show stronger fluctuations in the approach to the continuum while the results for SU(3) obtain larger uncertainties. In general, the results forb 2 are closer to zero when massive modes are included. The reason for this is obviously given by the presence of the additional R −4 term in eq. (5.1). This term contaminates the boundary correction term and thus reduces the magnitude ofb 2 . The results for r 0 m show a rather smooth approach to the continuum with the exception of the point at β = 25.0 for SU(3) gauge theory, which, however, also has a large uncertainty.

Continuum extrapolations
Ultimately we are interested in the comparison between the two sets of results in the continuum. To this end we perform a continuum extrapolation similarly to the one in section 4.2 by using fits of the type (1) to (3)  and we perform fits of the types (1) to (3) where bb 2 ,2 is replaced by b m,2 in the fit definitions. The difference to section 4.2 is, that we only perform continuum extrapolations for fit F, so that no averaging in the continuum limit is necessary. Once more the fits have also been done for the results from R min ± 1a to assess the uncertainty associated with the choice for R min .
The results for the continuum extrapolations with fits (1) to (3) are given in table 11. We usually obtain a rather good continuum extrapolation, even though forb 2 in SU(2) gauge theory the extrapolation leads to rather large values of χ 2 /dof around 5 to 10, due to the fluctuations ofb 2 for the small values of a. Since the continuum extrapolation is a bit more problematic in this case we will use the results from fit (3) for our final results. This fit appears to be a bit more stable than the continuum extrapolation linear in a 2 using all data points. Once more the uncertainty concerning the continuum extrapolation is estimated via the maximal difference of the result of fit (3)

Summary and discussion of the final results
The previous three sections contained a number of interesting results, which are, however, difficult to extract from the somewhat lengthy discussion. Before we move on to the conclusions let us thus summarise and discuss the main findings in some more detail.

Summary of results
In section 3 we have extracted the Sommer parameter and the string tension with high accuracy from our results for the static potential. The results are given in Here the string tension has been extracted by parameterising the potential with the LC spectrum, eq. (2.7), up to the point where the result for σ did agree with the NLO expansion of eq. (2.7). The first uncertainty is the statistical one and the second uncertainty originates from the continuum extrapolation. In three dimensions any comparison to physical units is meaningless. Nonetheless, identifying r 0 with 0.5 fm [68] (to define the unit 'fm' in three dimensions) and using the standard conversion factor c = 197.3(. . .) fm MeV we obtain √ σ cont = 487.6(2)(1) MeV for SU (2) 486.4(2)(1) MeV for SU(3) , (6.2) respectively. A comparison to the latest results for the continuum string tension in three dimensions [69,70] is only possible in terms of the Karabili-Kim-Nair prediction [71], i.e.
in terms of the continuum extrapolated coupling. We leave this type of comparison to the next publication. Our results at finite lattice spacing are fully in agreement, but more precise, than the values given in [72,73]. With a constant continuum extrapolation the SU(2) results from [73] give a continuum result around 1.234, which is a bit smaller but comparable to our result. After the extraction of the leading order (linear) behaviour we have compared the results for the potential to the leading order EST prediction, namely the LC spectrum, and extracted the exponent of the leading order correction in 1/R. The results, shown in figure 4 are consistent with an exponent of 4, as expected from the EST, where the leading order correction is the boundary term proportional tob 2 . We then extracted the value for 3) The first error is purely statistical, the second is the systematical error due to the unknown correction terms to the potential, the third is the one associated with the particular choice for the minimal value of R included in the fit for the extraction ofb 2 , R min , and the fourth systematic uncertainty is the one due to the continuum extrapolation (for details see section 4). These values indicate that the magnitude ofb 2 decreases with increasing N , meaning that it could potentially vanish in the limit N → ∞. The result also shows that b 2 , indeed, is non-universal, as expected since its value is not constrained in the EST. Concerning the extraction ofb 2 the main uncertainty comes from the fact that massive modes, or, equivalently, a possible rigidity term in the EST, lead to a contaminating additional R −4 correction. The presence of this term is expected to change the value ofb 2 . At present it is unclear whether such a contamination will ultimately be present or not, so that we have to take this possibility into account. On the basis of our simulations we also cannot exclude this possibility, since the fits to the potential reported in section 5.1 work equally well compared to the results from section 4.1. We can exclude, however, the possibility that the R −4 correction is fully due to the correction associated with the massive modes, since fit J (see section 5.1) leads to a much larger value for R min than the other fits. Repeating the whole analysis, we see that the accuracy is only sufficient for fit F from tables 9 and 10. Using these results we obtain the continuum results listed in table 11  Here the first error is purely statistical, the second is the one associated with the particular choice for R min and the third systematic uncertainty is the one due to the continuum extrapolation. In this analysis we have only been able to use one of the fits for a reliable extraction ofb 2 (fit F from tables 9 and 10), which is why one of the systematic uncertainties cannot be estimated. The comparison of eqs. (6.3) and (6.4) reveals that the inclusion of the contamination reduces the magnitude ofb 2 but does not alter the qualitative feature of a decrease in magnitude ofb 2 with increasing N . From the extraction ofb 2 when massive modes are included we also obtain an estimate for the mass of the massive modes. The continuum results are listed in table 11, and as our final continuum results we get r 0 m cont = 2.61 (6)(28)( 5) for SU (2) 2.71 (9)(29)(11) for SU(3) , (6.5) where the uncertainties are as in eq. (6.4). The interpretation of this "mass" is an open question. The first possibility is that it represents a massive mode on the string worldsheet, which could indicate that we observe the three-dimensional analogue to the worldsheet axion (cf. section 2.3). Recall, that the topological coupling term does not exist in 3d, so that the massive mode will appear as a quasi-free mode on the worldsheet up to coupling terms including two massive boson fields. On the other hand, the contribution could also be due to the rigidity term, in which case the results from eq. (6.5) can be translated into results for the coupling α via eq. (2.12). It is intriguing to note that the results for m, dividing by √ σ from eq. (6.1), are very similar to the results for the mass of the worldsheet axion from [40,41], at least for the case of SU(2) gauge theory. For SU(3) gauge theory our result is somewhat larger than the result from [40,41], however, in contrast to those results our result is extrapolated to the continuum. In fact, our results around β = 25.0 are already fully compatible with the results from [40,41]. This opens up the possibility that, indeed, we are looking at a massive mode on the three dimensional worldsheet which is similar in nature to the worldsheet axion. It is also important to ask whether the value for m extracted in our analysis does still comply with the framework of the EST. From eq. (2.6) we expect modes of masses down to a few times √ σ to be integrated out. The masses we have obtained here are around twice as large as √ σ. In our fits we could typically go down to R/r 0 ≈ 0.7. When we assume that this is a sign for the scale where the EST is bound to break down, this leads to a cut-off scale of 1.2 × √ σ. However, this estimate could well underestimate the true cut-off scale, so that a mass of 2 × √ σ can potentially be below that bound.

Comparison between continuum EST and data
We will now compare the continuum predictions for the potential from the EST to the results for the potential on our finest lattice spacings. A suitable quantity to visualise the subleading contributions to the potential is the curvature, associated with the second derivative (the uncertainties are smaller than the line), and the red band includes the boundary term from eq. (2.10) with the continuum boundary coefficientb 2 from eq. (6.3). The inclusion of the boundary term obviously enhances the agreement with the data significantly down to small values of R. We also note that the two curves are basically indistinguishable at the scale of the plot (they are not identical, however). For SU(2) gauge theory the data at finite lattice spacing already agrees rather well with the continuum EST, while for SU(3) gauge theory the data lies below the curve. We have not shown the curve including massive modes in the plot since it overlaps with the 'LC+O(R −4 )' curve and compares very similarly to the data.

Conclusions
In this article we have performed a high precision study of the static qq potential in threedimensional SU(N ) gauge theory with N = 2 and 3 and compared the results to the potential obtained from the effective string theory. In particular, we obtained accurate results with full control over the systematic effects for the continuum string tension, the non-universal boundary coefficientb 2 and, for the extended analysis, the "mass" of the possible massive mode within the EST. The results are summarised and discussed in detail in section 6.1. In particular, we could show that the leading order correction to the light cone spectrum is of O(R −4 ) (cf. section 3.3). If massive modes are present, the results for b 2 change (see eqs. (6.3) and (6.4)) due to another correction of O(R −4 ), contaminating the result forb 2 . However, the contribution from the massive modes only is not enough to describe the data, so that we can conclude thatb 2 = 0. 4 The data forb 2 shows an interesting trend towards zero with increasing N , leading to the possibility thatb 2 could vanish in the limit N → ∞. This result is also independent of whether or not we include massive modes in the analysis. In the context of generalisations of the AdS/CFT correspondence to large N gauge theories, this is an interesting result since it imposes constraints on the fundamental string theory for the dual version of Yang-Mills theories at large N . We will investigate this issue further in the next publication in this series of papers.
The main obstacle of our analysis forb 2 is the fact, that it is impossible to judge whether or not massive modes, or, equivalently, the rigidity term (or even both) are present. While some properties (like the decrease in magnitude for N → ∞) are independent of this issue, the exact value ofb 2 can only be extracted once the issue is resolved. To this end it would be desirable to have an expression for the excited state energies in the presence of a massive mode, or, alternatively, the rigidity term. This could potentially help to discriminate between the existence or non-existence of such terms and, in addition, it could help to discriminate between the two types of additional contributions.
It is intriguing to see that the results for the mass of the possible massive modes is in good agreement with the masses found in [40,41]. It is thus possible that we are seeing a massive mode on the worldsheet which is similar in nature to the worldsheet axion. Note, that in 3d the topological coupling term is absent, so that the analogue to the worldsheet axion appears as a quasi-free mode on the worldsheet coupling to the GBs via the covariant derivative and possible higher order coupling terms. It will be interesting to see whether the mass of the massive modes remains consistent with the one of the worldsheet axion for N → ∞.

A Simulation setup
This appendix contains the details of the numerical simulations. For the configuration updates we have used the, Cabbibo-Marinari heatbath algorithm [74] for SU(N ) gauge theories with all SU(2) subgroups, in combination with the improved SU(2) heatbath algorithm from [75]. For each update we perform three overrelaxation steps [76].
To obtain a clean result for the groundstate potential in terms of contaminations from excited states the spatial correlation function of two Polyakov loops is the most promising observable. Its spectral representation is given by for R < L/2, where T = a N t and L = a N s are the temporal and spatial lengths of the lattice, b i is the overlap between the operator and the associated energy eigenstate and E i (R) is the i'th energy level. We always take energies to be ordered in ascending order, i.e. E 0 < E 1 < E 2 < . . . . In the limit of T → ∞ the leading contribution to eq. (A.1) is given by the groundstate E 0 (R) = V (R) and excited states are suppressed with exp{−[E i (R) − E 0 (R)] T }. For large values of T excited states can thus be neglected to a good approximation and we check in appendix B whether this is true for the temporal extents used in the present study. If we can neglect contaminations from excited states we can extract V (R) via Note, that the potential is determined from the effective string theory up to a constant, which we denote as V 0 (it can also be related to the constant term µ in the boundary action (2.3); here, however, we are not really interested in this quantity). The string tension can either be extracted from a fit to the potential, or via the force Note, that the determination via the force does not demand the determination of V 0 . Here we will use both methods to extract the string tension and compare the results. The string tension can also be used to fix the lattice spacing in physical units and observables can then be expressed in units of a √ σ. However, the determination of a √ σ demands to control the asymptotic R → ∞ behaviour. A more suitable observable to set the reference scale is the Sommer parameter r 0 [68], which is defined implicitly by r 0 can be obtained directly by interpolation of the results for R 2 F (R).
The reliable extraction of the potential demands the control of the contaminations from excited states and thus large temporal extents. Since the signal-to-noise ratio of such Polyakov loop correlation functions decays exponentially with T , the reliable extraction of the expectation values demands the use of a suitable error reduction algorithm. Our choice is to use the multilevel algorithm introduced by Lüscher and Weisz [67]. In particular, we use one level of averaging for all lattices and the parameters such as the number of sublattice updates, n t , and the temporal extent of the sublattices, t s , are given in table 2.

B Control of systematic uncertainties
For the measurements of the potential there are several systematic effects that need to be controlled for a reliable extraction of the potential at a given lattice spacing. In particular, this concerns the contaminations from excited states and finite size effects. Checks concerning these two effects are reported in the following. There is yet another effect at very fine lattice spacings, owing to the "freezing" of the topological charge in certain sectors. We control this effect by using a completely new run (including full thermalisation starting from a random configuration) for at least each two measurements. This can be done since the thermalisation (we take a very conservatively long thermalisation with at least 1000 update steps) needs relatively little time compared to one measurement with the multilevel algorithm.

Excited state contributions
Following the spectral representation (A.1) for the Polyakov loop correlation function, the correction to eq. (A.2) owing to excited states is given by (see also [77]) The contamination from the second term on the right hand side can be checked by performing simulations for two different temporal extents T of the lattice. Here we have done simulations for SU(3) with β = 14.0 on 48 and 64 × 48 2 lattices and compare the results for the potential. The results for the difference of the two energies are shown in figure 9.
As can be seen from the plot, the two sets of energies agree perfectly within errors for the whole range of qq separations considered, meaning that excited states are negligible. Since for the different ensembles T changes only little in physical units, we expect this to hold for all lattices considered. Note, that the energy levels (and the overlaps) of the flux tube, in general, change only on the subleading level for different values of N . Consequently, we would expect this result to be valid for different values of N as well. Since E 0 (R) growths linearly with R, we expect the winding contribution to be strongly suppressed when L − R > R.
Assuming the leading order linear relation σR for the bare energy in lattice units and a relative error of 10 −3 for the Polyakov loop correlator (which is a lower bound for what we get for the larger R values), we see that the second term is negligible if where in the second step we have used that for our lattices T = L and √ σL 12.3. This condition is always fulfilled in our study.
As an additional check we computed the potential on a 48 × 64 2 lattice for SU(3) with β = 14.0 and compared the results for the potential with the 48 3 run used in the analysis. The results are shown in figure 10. The plot indicates that the effects due to the finite lattice extent are negligible. Note, however, that finite size effects can be present for these values of R if the lattice is chosen to be smaller. This can be seen from the comparison to the results for a 48 × 32 2 lattice, also shown in figure 10.
We also note that another finite size effect could be due to a glueball exchange between the two Polyakov loops. Generically we expect this effect to be suppressed by a factor of exp(−m G R), or exp(−m G (L − R)). Clearly, as long as m G > √ σ (which is always the case) this contribution appears as an excited state effect, which we have already ruled out in the check for excited state contaminations.