Ising string beyond Nambu-Goto

A major result of the Effective String Theory (EST) description of confinement is the so called"low energy universality,"which states that the first few terms of the large distance expansion of any EST are universal and coincide with those of the Nambu-Goto action. Going beyond this approximation is one of the most interesting open problems in the EST. In the higher order terms beyond Nambu-Goto several important pieces of physical information are encoded, which could improve our understanding of the physical mechanisms behind confinement and of the physical degrees of freedoms which originate the EST. In this paper we evaluate numerically the first two of these corrections in the case of the three dimensional gauge Ising model. The first of them turns out to be negative: $\gamma_3=-0.00048(4)$, similar (but not equal) to the one recently measured in the $SU(2)$ Yang Mills theory in three dimensions and compatible with the bootstrap bound $\gamma_3 \geq -\frac{1}{768}$.


Introduction
One of the most promising approaches to understand and model the non-perturbative behaviour of confining Yang-Mills theories is the "Effective String Theory" (EST) description in which the confining flux tube joining together a quark-antiquark pair is modeled as a thin vibrating string [1,2,3,4,5] .
Recently, there has been a lot of progress in this context.In particular it has been realized that the EST enjoys the so called "low energy universality" [6,7,8,9,10,11,12,13] which states that, due to the peculiar features of the string action and to the symmetry constraints imposed by the Poincaré invariance in the target space, the first few terms of its large distance expansion are fixed and are thus universal.This implies that the EST is much more predictive than typical effective models and in fact its predictions have been be successfully compared in the past few years with lot of results on the interquark potential from Monte Carlo simulations in Lattice Gauge Theories (LGTs) (for recent reviews see for instance [13,14,15]).
At the same time it was recently realized that the simplest, Lorentz invariant, EST which is the well known Nambu-Goto model [1,2] is an exactly integrable, irrelevant, perturbation of the bidimensional free Gaussian model [10], driven by the T T operator of the D − 2 free bosons1 [16] which represent the transverse degree of freedom of the string.This observation stimulated lot of interesting results even beyond the original application to Yang-Mills theories [17,18,19,20,21,22,23,24]. In particular they are at the basis of a S-matrix bootstrap approach which can be used to constrain the EST action even beyond the Nambu-Goto approximation [25,26].Indeed, it is by now clear that the Nambu-Goto action should be considered only as a first order approximation of the actual EST describing the non-perturbative behavior of the Yang-Mills theories.Going beyond this approximation is one of the most interesting open problems in this context.In the higher order terms beyond Nambu-Goto several important pieces of physical information are encoded.Their study could be of great importance to understand the physical mechanisms behind confinement or the physical degrees of freedoms which originate the EST.
In particular, it is only by looking at these higher order terms that one may hope to find signatures, in the confining string, of the gauge group of the underlying LGT.For this reason there is an increasing interest in exploring these corrections in different LGTs [27].In this respect the three dimensional gauge Ising model that we shall study in this paper is a perfect laboratory to address this issue because, thanks to the duality transformation with the the 3d Ising model, one can use innovative, powerful, algorithms to estimate these corrections.Moreover, its gauge symmetry is very different from standard SU (N ) gauge groups and allows to test, for instance, which is the effect of the discrete vs continuous gauge symmetries on the confining string.A final important reason is that the scaling function of the string tension, which will play a central role in the following, is known (thanks again to duality) with very high precision.
Due to the low energy universality these corrections appear at a very high order in the large distance expansion and their evaluation is a delicate task.One must reach very precise estimates of the ground state energy of the string for a wide range of distances and, possibly, for a few different values of β to test the correct scaling behaviour.For these reasons we decided to evaluate them with two different approaches, using different algorithms (hierarchical Metropolis in one case and Swendsen-Wang in the other) choosing different observables (the ground state energy of the string in one case and its width in the other) and looking for an overall agreement of the final results between the two methods.
We were able to detect the first two corrections beyond Nambu-Goto which are described, as we shall see below, by the two parameters γ 3 and γ 5 [25].The value of γ 3 = −0.00048(4) agrees with the S-matrix bound found in [25] and is similar, but slightly more negative than the value γ 3 | SU (2) = −0.00037(6) found in [27] for the SU (2) LGT in (2+1) dimensions.The second correction is very small γ 5 = 3.0(4) × 10 −7 but not compatible with zero and its inclusion in the fit turned out to be mandatory to reach reasonable values of the reduced χ 2 .These values represent the first step toward a systematic study of γ 3 in LGTs This paper is organized as follows.In the second section we shall define the problem and set the notations, then in the third section we shall recall a few basic results of the Effective String Theory of confinement.In the fourth section we shall present the 3d gauge Ising model and discuss its properties.Then, in sections five and six, we shall discuss the two approaches that we used to evaluate γ 3 and γ 5 and present the main steps of the data analysis.Finally the last section is devoted to a summary of the results and a few concluding remarks.

General setting and notations
In the following we shall be mainly interested in finite temperature LGTs, which can be realized by imposing periodic boundary conditions in the time direction for the bosonic field (and antiperiodic for fermionic ones).In a finite temperature setting the compactified "time" direction does not have any longer the meaning of time (recall that we are describing a system at equilibrium in the canonical ensemble) but its size N t is instead a measure of the inverse temperature of the system.Thus a lattice of size N s a in the spatial directions and N t a in the timelike direction represents the regularized version of a system of finite volume V = (N s a) d at a finite temperature T = 1/(N t a).In the following we shall set the lattice spacing to a = 1, and systematically neglect it.
In a finite temperature setting one can define a new class of topologically non-trivial observables which are gauge invariant thanks to the periodic boundary conditions in the time direction: the Polyakov loops.If we define the link dynamical variables of the gauge model as U µ (⃗ x, z) (where µ denotes the direction of the link and (⃗ x, z) its coordinates in the the lattice), we may define the Polyakov loop P (⃗ x) as follows: In a pure LGT the Polyakov loop acquires a non-zero expectation value in the deconfined phase and is thus an order parameter of the finite temperature deconfinement transition.
The value β c (N t ) of this deconfinement transition in a lattice of size N t = 1/T can be used to define a new physical observable T c which is obtained by inverting β c (N t ).We obtain in this way, for each value of β, the lattice size in the time direction (which we shall call in the following N t,c (β)) at which the model undergoes the deconfinement transition and from this the critical temperature T c (β) ≡ 1/N t,c (β) as a function of β.We shall use this quantity to set the scale of our simulations.
In the following we shall be mainly interested in correlators of Polyakov loops: where R is the distance between the two Polyakov loops and the subscript N t in the expectation value reminds that the correlator was evaluated on a lattice at temperature T = 1/N t .G(R, N t ) can be used to define a finite temperature version of the interquark potential: In the confining phase we expect, for large values of R, a linearly rising behaviour for V (R, N t ), which implies the following behaviour for the correlator: From eq.( 4) we may extract the ground state energy E 0 (T ) of the confining flux tube joining the quark-antiquark pair.E 0 (T ) depends on the temperature of the model and vanishes at the deconfinement transition.The finite temperature behaviour of E 0 (T ) will play a major role in the rest of the paper.It is interesting to notice that the observable (2) has the topology of a cylinder whose circumference is fixed by the (inverse) temperature N t and the length by the interquark distance R.
3 Effective String Theory

The Nambu-Goto action
The behaviour of the flux tube in a confining LGT is well described by the Effective String Theory which models the flux tube as a thin vibrating string and allows to evaluate the contribution to the Polyakov loop correlator of the quantum fluctuations of this flux tube.The simplest possible EST fulfilling the constraints imposed by the Lorentz invariance in the target space is the Nambu-Goto action [1,2] defined as follows: where g ≡ det g αβ and is the metric induced on the reference world-sheet surface Σ by the mapping X µ (ξ) of the world sheet in the target space, and ξ ≡ (ξ 0 , ξ 1 ) denote the worldsheet coordinates.This term has a simple geometric interpretation: it measures the area of the surface spanned by the string in the target space and has only one free parameter2 : the string tension σ 0 .
In order to perform calculations with the Nambu-Goto action one has first to fix its reparametrization invariance.The standard choice is the so called "physical gauge".In this gauge the two worldsheet coordinates are identified with the longitudinal degrees of freedom of the string: ξ 0 = X 0 , ξ 1 = X 1 , so that the string action can be expressed as a function only of the (D − 2) degrees of freedom corresponding to the transverse displacements, X i , with i = 2, . . ., (D − 1) which are assumed to be single-valued functions of the worldsheet coordinates.We shall comment below on the problems of this gauge fixing choice.
With this gauge choice the determinant of the metric can be written as and the Nambu-Goto action can then be written as a low-energy expansion in the number of derivatives of the transverse degrees of freedom of the string which, by a suitable redefinition of the fields, can be rephrased as a large distance expansion.
The first term of this expansion is exactly the gaussian action, i.e. a two dimensional Conformal Field Theory (CFT) of D − 2 free bosons which represent the transverse degrees of freedom of the string, and the remaining terms combine themselves so as to give an exactly integrable, irrelevant perturbation of this CFT [10], driven by the T T operator of the D − 2 free bosons [16].
Thanks to this exact integrability, the partition function of the model can be written explicitely [7,28].For the Polyakov loop correlator in which we are interested here3 , the expression in D space-time dimensions is, using the notations of [7,28]: where K 1 2 (D−3) is the modified Bessel function of order D−3 2 , R denotes, as above, the distance between the two Polyakov loops, N t the size of the lattice in the compactified direction, w n is the multiplicity of the closed string states which propagate from one Polyakov loop to the other, and E n their energies: At large distances the correlator is dominated by the lowest state whose multiplicity is w 0 = 1 and the Polyakov loop correlator, neglecting irrelevant constants, reduces to which, in the D = 3 case in which we are interested simplifies to with Thanks to the exponential term in the asymptotic expansion of the Bessel function, we find at large distance, as expected, a linearly rising behaviour of the interquark potential controlled by the ground state energy E 0 .On top of this we have a set of subleading corrections, encoded in the asymptotic expansion of K 0 , which represent a specific, unique, signature of the Nambu-Goto action and must be taken into account when fitting the results of Monte Carlo simulations.
An important side consequence of this result is that we can extract from the tachyonic singularity of E 0 an estimate for the critical temperature T c,N G measured in units of the square root of the string tension √ σ 0 [31,32] which is, for generic values of D, and corresponds to the value of the ratio √ σ 0 for which the ground state energy E 0 vanishes.Using this result we can rewrite the ground state energy as where we use the notation | N G to stress the fact that this is only the Nambu-Goto estimate for the ground state energy of the string, which we may expect to be modified by other terms in the EST action.
The estimate quoted above for the critical temperature turns out to be in remarkable (but not exact!)agreement with the results obtained from Monte Carlo simulations, both for non-abelian LGTs and for the three dimensional gauge Ising model.However, the remaining small deviations of T c,N G from the Monte Carlo results, together with the fact that the Nambu-Goto EST predicts (as can be seen looking at eq.( 17)) a deconfinement transition of the second order, with a mean field value for the critical index (which is in disagreement with the Monte Carlo results for all known LGTs) suggest that the Nambu-Goto action cannot be the end of the story and that some correcting terms beyond Nambu-Goto should be present in the EST.

Beyond Nambu-Goto
It is by now clear that in the actual EST of the confining string the Nambu-Goto action is only the first term of an infinite series of contributions.Indeed, there are several reasons why the Nambu-Goto action is unsatisfactory and must be completed with some higher order correction.Besides the above mentioned disagreement at the deconfinement transition, a major problem of the Nambu Goto action is that it is, so to speak, too universal.It predicts the same behaviour for all LGTs, with no dependence on the gauge group.Moreover, it is well known that the physical gauge fixing discussed above is anomalous and it is widely expected that this anomaly could be cured by higher order terms in the EST action.
The requirement of Poincarè invariance in the target space strongly constrains the terms which can be included in the EST beyond Nambu Goto [6,7,8,9,10,11,12,13].In D = 3 the first few allowed terms can be written as follows: where the γ i are new coupling constants, R denotes the Ricci scalar constructed from the induced metric, and K is the extrinsic curvature, defined as K = ∆(g)X, with the Laplacian in the space with metric g αβ .In principle the new coupling constants γ i , should be fixed, as we do for σ 0 , by comparing with experiments (or more likely, with computer simulations).However this process is simplified by the observation that K 2 vanishes on shell and that the term proportional to R is a topological invariant in two dimensions.Since in the long-string limit in which we are interested one does not expect topology-changing fluctuations, both these terms can be neglected and the first non-trivial contribution appears only at higher orders [13].This result is known as "low energy universality" [10] and strongly constrains the form of the EST.It implies that the first correction beyond the Nambu-Goto action can only appear at the order 1/R 7 (or 1/N 7 t in the finite temperature setting in which we are interested in this paper).This explains why the Nambu-Goto model has been so successfull over these last forty years to describe the infrared behaviour of confining gauge theories despite its simplicity and why the deconfinement temperature predicted by Nambu-Goto is so close to the one obtained in Monte Carlo simulations.At the same time this explains why identifying these corrections is so difficult and only in the last few years it was possible to unambiguously detect them in Monte Carlo simulations [27,33,19,34,35] The first non vanishing term γ 3 K 4 is only the first of an infinite sequence of terms, obtained combining higher order curvature invariants.It turns out that the best way to organize these terms is to study the 2 → 2 scattering amplitude of the string excitations [10].It can be shown that in D = 3 this S-matrix S = e 2iδ can be written in a particularly simple and elegant form: where the the first term s/4 leads to the energy spectrum of the Nambu-Goto action, the fact that the term proportional to s 2 is missing is the way in which the low energy universality is realized in this S-matrix approach and the next non trivial term is exactly the S-matrix description of the K 4 correction mentioned above.By using the analiticity properties of this S-matrix and requiring the UV completion of the underlying theory it is possible to obtain a set of important results on the ground state energy of the string [10,25].
• both the 1/N 7 t and the 1/N 9 t terms are controlled only by γ 3 and the next independent parameter γ 5 only appears at the order 1/N 11 t • It is possible to set bounds on these parameters.In particular, defining γn =γ n + (−1) (n+1)/2 1 n2 3n−1 one finds [25] γ3 ≥ 0 γ5 ≥ 4γ From the S-matrix, by using the so called Thermodynamic Bethe Ansatz one can obtain [25] the following expression for the non universal corrections up to the order 1/N 11 This is the expression which we shall compare with the results of our simulations.

The boundary corrections problem and how to deal with it.
It is clear from the previous section that measuring the γ i coeffcients on the lattice is a highly non-trivial task.In particular it is essentially impossible in the standard "zero temperature" scenario, in which the contribution of the effective string to the interquark potential manifests itself as an expansion in powers of 1/R and the corrections in which we are interested, which appear in this expansion at the order 1/R 7 , are shadowed by the boundary corrections which are proportional to 1/R4 [36,11,37,38,39].
There are in principle two ways to avoid this problem.
• The first is to study observables without boundaries.This can be done for abelian gauge models using duality and looking at the finite size effects of the interface free energy (choosing interfaces with periodic boundary conditions) [40,41] .In the case of the Ising model this approach was recently discussed in [42] where it was shown that corrections beyond Nambu-Goto certainly exist in the gauge Ising model and are rather large.However it turned out to be difficult to quantify these corrections, most probably due to the interactions between nearby interfaces 4 .
• The second is to study the model in the finite temperature regime (just below the deconfinement transition) in the limit of very large separation of the two Polyakov loops (R >> N t ).
It can be shown that in this regime the boundary corrections become subleading and can be neglected [43,15].Moreover, this is exactly the limit discussed in the previous section, in which the results obtained with the S-matrix approach and TBA hold.Thus, by choosing this geometry in our simulations we shall be able to make immediate contact with eq. ( 23), with no interference from the boundary terms and extract from the data reliable estimates for the γ i coefficents.This was the approach recently used to study these corrections in the SU (2) gauge model in three dimensions in [27].
Once the geometry of the observables in which we are intersted is fixed, the remaining task is to obtain estimates for these observables precise enough to detect and quantify the tiny corrections in which we are interested.In this geometry this requires studying the system at large interquark distances and standard algorithms are affected by an exponentially decreasing signal to noise ratio in this limit.The main advantage of studying abelian models is that, thanks to duality, it is possible to avoid this limitation and to study (using for instance non-local cluster algorithms as in [44] or the so called "snake algorithm" [45]) Polyakov loops correlators at any interquark distance R with the same signal to noise ratio.This is the main reason behind the choice of the Ising model as a laboratory to study the EST.
There is indeed a long track record of applications of this kind of methods to the 3d gauge Ising model to study subtle features of EST.In particular, in [44,46,47,48] which may be considered as precursors of the present work, the 1/R 3 correction to the interquark potential was precisely measured for the first time and shown to be exactly the one predicted by the Nambu-Goto action (in agreemeent with the low energy universality).Later the same approach was adopted for the 3d U (1) model in [49] and allowed to unambiguously detect corrections to the Nambu-Goto actions as the continuum limit was approached.
We shall discuss in sect.5 below the result of a study performed with the same snake algorithm used in [46,47].This algorithm allowed us to obtain a first estimate of γ 3 which however turned out to be affected by a rather large statistical uncertainty.In order to improve this uncertainty and to test the robustness of the result we decided to evaluate the same quantity with a completely different method, which was proposed a few years ago in [50] and that we describe in detail in sect.6.
This second estimate agrees within the errors with the previous one, is more precise and allows to quantify with rather good precision even the next to leading correction γ 5 .
The agreement between the two estimates strongly supports the reliability of our analysis.
In the following sections we shall first describe the main fetaures of the 3d gauge Ising model and then discuss the two approaches that we used to evaluate the γ i coefficients.

The 3d Gauge Ising Model
The three dimensional Gauge Ising Model (also known as Z Z 2 gauge model) was proposed over fifty years ago by Wegner [51] as a tool for understanding the properties of lattice models with gauge symmetries.This model exhibits a local Z Z 2 symmetry, which is realized by the choice of σ l ∈ {1, −1} as the dynamical Z Z 2 link variables.The plaquette action, derived from the familiar Wilson action, is tailored specifically for the case of Z Z 2 link variables and can be defined as follows: The action S Z 2 is a sum over all the plaquettes of a cubic lattice, Despite its apparent simplicity the 3d Gauge Ising Model shares with more complex LGTs several important properties.It is charaterized by a confining string with a non-trivial spectrum of string excitations [52,53] and has a glueball spectrum very similar to the one of more complex 3d LGTs [54] .For these reasons it is a perfect laboratory to test, with high precision, non trivial properties of the confining strings in LGTs.
This model is known to have a bulk (i.e. at zero temperature) deconfinement transition at β c = 0.76141330(6) (this value is obtained via duality from the critical temperature of the 3d Ising model quoted in [55], see below).For values of the coupling β < β c the model is in the confining phase, while for β > β c it is deconfined.The transition at β = β c is of second order and, due to the duality relation (see below), it belongs to the same universality class of the standard magnetization transition of the 3d Ising model.This model also possesses an (infinite-order) "roughening transition" at β r = 0.47542(1) [56] (in the confined phase), which separates the strong coupling regime (for β < β r ) from the so-called "rough phase" (for β r < β < β c ).
In the following we shall be interested in the behaviour of the model in the confining phase, in the scaling region near the critical point.In particular we shall study the model at three different values of the coupling constant β (see tab.1), for which the (finite) deconfinement temperature is known with high precision [57] so as to be able to precisely set the scale for the distances between Polyakov loops and for the lattice size in the time direction (i.e. the inverse of the finite temperature of the model).These values are all located in the rough phase and close enough to β c so as to be within the scaling region.
The major reason of interest of this model is that it is related to the ordinary three dimensional Ising spin model by an exact duality mapping à la Kramers and Wannier (see [58] for a general review on duality transformations): where the Ising spin model is defined by the usual action: where, as usual, the s i ∈ {1, −1} are Z Z 2 spin variables, the i and j indices denote the sites of the dual lattice, the notation ⟨i,j⟩ means that spin variables interact with their nearest-neighbours only and {s i } denotes the sum over spin configurations.
The critical temperature of the model is known with remarkable precision β c = 0.22165462(2) [55], moreover thanks to the recent advances in the bootstrap approach [59,60] also the critical indices of the two relevant operators, the magnetization M and the energy ϵ, are known with high precision: ∆ M = 0.5181489 (10) and ∆ ϵ = 1.412625 (10) respectively [59,60].From ∆ ϵ we can extract the critical index ν = 1  3−∆ϵ = 0.6299708.... Thanks to duality this is the same critical index which drives the critical behaviour of the string tension in the 3d gauge Ising model.More precisely σ(β) ∼ σ c (β c − β) 2ν .We shall further discuss the scaling behaviour of the string tension in the appendix.
The main reason of interest of this mapping is that a similar construction can be performed also in presence of external source terms for the gauge model (for instance, a pair of Polyakov loops).This can be easily realized by introducing sets of topological defects in the spin system.As a result, it is possible to show [46] that, for instance, the Polyakov loop correlator in which we are interested is mapped into the partition function of the spin system with anti-ferromagnetic coupling on a well defined set of links: with: where the value of the J ⟨i,j⟩ coupling is +1 everywhere, except on a set of bonds, which pierce a surface (in the direct lattice) having the source worldlines as its boundary: for such a set of bonds, J ⟨i,j⟩ = −1.
Similar mappings can be constructed essentially for any observable of interest in the gauge model, from Wilson loops to glueball correlators 5 .
Both the numerical algorithms that we shall use in the following exploit this duality of the model, simulating the Ising spin system, and measuring (the first algorithm) ratios of partition functions associated with different stacks of defects (which we shall use to express the expectation values of Polyakov loops pairs in the original gauge model) or expectation values of spins in presence of the defect surface (the second algorithm).
A particularly useful advantage of numerical simulations in the dual setting is the fact that this method overcomes the problem of exponential signal-to-noise ratio decay, which is usually found when studying the interquark potential V (r) at larger and larger distances.
Another important feature is that, thanks to duality, the valus of the Polyakov loop correlators that we obtain in this way are not affected by corrections due to the periodic boundary conditions in the spacelike directions.This greatly simplifies the study of these correlators (no additional terms must be included in the fits to keep into account these corrections) and allow to study the system for (relatively) small lattice size in the spacelike directions.
We perfomed simulations for the three values of β quoted in tab.1: β = 0.751800, 0.756427, 0.758266 which were chosen because for these values the deconfinement temperature is known with very high precision and coincides with 1/T c = L c = 8, 12, 16 respectively [57].1: Some information on the three values of β, listed in the first column, which we chose for the simulations.In the second column we report the corresponding values of β for the (dual) 3d Ising model.In the last three colums columns we report respectively the (inverse of) the deconfinement temperature, the string tension (taken from ref. [41]) and the value of α (see below for the definition of α(β) and its evaluation from the scaling function of the model).

Analysis with the snake algorithm
The first method that we used to estimate γ 3 is the Ising implementation of the snake algorithm [45] discussed in detail in [46].The main feature of the algorithm is the hierarchical organization of the updates.In our particular implementation we chose five sublattice levels with size {6, 13, 17, 21, 24} lattice spacings respectively.
We studied the first two values of β reported in tab.1 corresponding to a deconfinement temperature of N t,c = 8 and N t,c = 12 respectively.Details on the simulations are reported in tab.2.For each β we studied seven values of N t just above the critical value N t,c (i.e.just below the deconfinement temperature).Then for each value of N t we simulated 8 different values of R as reported in the table.The value of the lattice size in the space direction was chosen to be ten times the value of N t,c to avoid finite size effects (which in any case are strongly suppressed thanks to the duality transformation).The values of R were chosen so as to make higher order energy levels in eq.( 9) negligible within the errors.Thus we could fit the R dependence of our Polyakov loop correlators using eq.( 13).From the snake algorithm we directly obtain the ratio of two nearby correlators F (R, N t ) = G(R + 1, N t )/G(R, N t ) which we thus fitted with the following expression: with E 0 (N t ) as the only free parameter of the fit.For all values of N t we found good values of the reduced χ 2 .We report in tab.3, as an example of the results obtained with the snake algorithm, the values obtained for the largest Polyakov loops correlators that we studied i.e. those for β = 0.756427 and N t = 24.As anticipated there is no increase in the signal to noise ratio as R increases and we could estimate the ratio of the two Polyakov loop correlators with less than 1% error for areas as large as 24 × 84 lattice spacings.9,10,11,12,14,16,18 8, 12, 16, 20, 24, 32, 40, 48 80 0.756427 12 13,14,15,16,18,20,24 18, 24, 30, 36, 48, 60, 72, 84 120 Table 2: Some information on the simulations.
We report in tab.4 and 5 the results of these fits.These are the values that we compared with the expectation of eq.( 23).We performed two types of fit.In the first we kept as free parameters only σ 0 and γ 3 .Accordingly we truncated the square root of the Nambu-Goto action to the same order O(N 9 t ) at which the γ 3 terms appears.Results of these fits are reported in tab.6 In the second we included also γ 5 and, accordingly, truncated the square root to the order O(N 11 t ).Results of these fits are reported in tab.7.A few observations on this result: • In the first type of fits the values of σ 0 that we obtain are definitely larger (more than three standard deviations) than the expected ones.Accordingly, if we try to fit the data keeping σ 0 to the expected value we also found very large values of χ 2 .Moreover, the fact that the  30) for β = 0.756427.In the last column the reduced χ 2 of the fits.reduced χ 2 is larger than 1 suggests that the inclusion of the next order term in the fits could lead to non negligible corrections and in fact the values of γ 5 in the fits truncated at O(N 11 t ) are different from zero within the errors.We also tried to fit the data truncating the expansion at the order O(N 7 t ) and keeping only the first correction proportional to γ 3 , but we found values of σ 0 even further away from the expected values.Table 6: Results of the fits of E 0 (N t ) according to eq.( 23) truncated at the order O(N 9 t ) .Table 7: Results of the fits of E 0 (N t ) according to eq.( 23) truncated at the order O(N 11 t ) .
• Including the O(N 11 t ) term in the fit the values of σ 0 that we obtain move toward the expected values and for both values of β are within two standard deviations from the values reported in tab.16 .
• Accordingly, in all these fits, if we force σ to the known values we always find unacceptably high values of the χ 2 which were associated to large deviations in the best fit values of γ 3 and γ 5 .
• The values of γ 3 and γ 5 that we find for the two values of β are compatible with each other within the errors.
Looking at the results we see that the estimates of γ 3 and γ 5 are affected by rather large errors and are strongly influenced by the value of σ.This is due to the fact that the fit is dominated by the Nambu-Goto part of the fitting function and in particular by the σN t term and by the Lüscher correction.It seems difficult to improve the overall precision of the result in the framework discussed in this section since the simulations (due to the peculiar structure of the snake algorithm and the need to increase the size of the lattice in the inverse temperature direction) become more and more expensive as β moves toward the continuum limit.For this reasons we decided to approach the task following a different strategy which could partially overcome these problems.
6 Using the mean flux density in presence of the Polyakov loops to estimate the EST ground state energy.
To avoid the above problems we tried a completely different approach.Following [50] instead of looking at the interquark potential, we studied the changes induced in the flux configuration by the presence of the Polyakov loops.We shall show below that as a consequence of this choice the explicit dependence on σ 0 and on the Lüscher term vanish.This makes this observable an unique tool to explore higher order corrections.Another reason of interest of this approach is that it is deeply related to another important issue of the effective string description of LGTs, i.e. the study of the flux tube thickness.It can be shown that in the high temperature regime in which we are presently interested the width of the flux tube increases linearly with the interquark distance and not logarithmicaly as one would naively expect [62,63,64].This linear increase is related to the linear increase in the flux energy that we observe here.In both cases the slope is temperature dependent and contains information on the higher order effective string corrections in which we are interested.
The lattice operator which measures the flux through a plaquette p in presence of two Polyakov loops P , P ′ for a generic LGT is: where U p is the trace of the ordered product of the link variables along the plaquette and in our case coincides with the product σ 2 introduced above.This is the quantity which is typically used to study the profile of the flux tube.In our analysis we are actually interested in a much simpler observable, the mean flux density, i.e. the sum of ϕ(p; P, P ′ ) over all the positions and orientations of the plaquettes, normalized to the number of plaquettes of the lattice.Due to translational invariance this quantity will depend only on the distance R between the two Polyakov loops and on the inverse temperature N t .Let us define where N p = 3N 2 s L denotes the number of plaquettes of the lattice.It is easy to see from the definition of G(R, N t ): that the mean flux density ⟨Φ(R, N t )⟩ can be written as: Since β appears in the observable only through the scale σ 0 (β) the above equation can be rewritten as This choice has two important consequences: • In the term proportional to R the string tension, which was the major source of systematic uncertainty in the previous calculation, is substituted by its derivative with respect to β, which can be evaluated with high confidence thanks to the very precise knowledge we have of the scaling behaviour of σ 0 (β) (see below).
• The first string correction (the "Lüscher term") which is universal and does not depend on σ 0 disappears.This makes the above observable a perfect tool to explore higher order corrections of the effective string.
From eq.( 13) we have where K ′ 0 denotes the derivative of the K 0 Bessel function, and α is defined as Using the identity K ′ 0 (z) = −K 1 (z) the logarithmic derivative K ′ 0 (z)/K 0 (z) can be expanded in powers of 1/z as follows: which gives: where the three functions A, B and C are defined as follows A crucial role in the analysis is played by α(β), a precise estimate of this quantity allows to strongly constrain the results of the fits.α can be extracted from the scaling function of the model and in its determination we leverage the very precise knowledge we have of this scaling function, thanks to the bootstrap results for the critical indices of the 3d Ising model.A detailed derivation can be found in [50].We report for completeness the main steps of the derivation in the appendix A. The values we used in the fit are listed in tab. 1.Once the value of α is fixed, we can use the values we obtain for A(N t ) to estimate the corrections in which we are interested 7 .By setting x = π 3σ 0 N 2 t we see that we can express the Nambu-Goto expectation for A (see eq.( 14)) as As anticipated the expansion starts at the order x −2 i.e.N −4 t , this makes this observable particularly suited to evaluate higher order corrections.
If we introduce the corrections beyond Nambu-Goto (see eq.( 23)) we find where the expansion of A(N t ) N G is truncated at the same order at which the additional terms proportional to γ 3 and γ 5 appear.This is the function which we shall use to fit the results of our simulations.

Numerical simulations
To estimate the function A(N t ) we used again duality and mapped the Polyakov loops correlator into the partition function of a 3d Ising spin model in which we changed the sign of the coupling of all the links dual to the surface bordered by the two Polyakov loops.This is the same method which was used in [65,63] to estimate the width of the flux tube.
We then estimated ⟨Φ(R, N t )⟩ by simply evaluating the mean value of the plaquette in presence of these frustrated links.We chose periodic boundary conditions in the original gauge Ising model.These b.c. are mapped by duality into a mixture of periodic and antiperiodic b.c. in the dual spin model.However we always chose N s large enough to eliminate any contribution from the antiperiodic sectors (which are depressed by terms proportional to e −σ 0 NsNt or e −σ 0 N 2 s ).Since, as discussed above, we are interested in the following only in the term proportional to R in ⟨Φ(R, N t )⟩, we may neglect the disconnected component ⟨U p ⟩ in the evaluation of ⟨Φ(R, N t )⟩.
Details on the simulations can be found in tab.8 We perfomed simulations for all the three values of β reported in tab.1.For each value of β and N t we simulated several values of the distance R between the two Polyakov loops.
For each simulation we used 10 5 iterations to thermalize the lattice and then performed 10 6 measures using a Swendsen-Wang algorithm.The values of R were chosen large enough so as to make the last term in eq.( 39) negligible, thus allowing to perform a simple linear fit to extract the values of A(N t ).We report in tab.9 an example of the data we obtained from the simulations (to allow a comparison, we chose the same values of β and N t reported in tab.3) and in tab.s 10,11 and 12 the values of A(N t ) extracted from these linear fits.
We then fitted these values with eq.( 44) keeping as only free parameters γ 3 and γ 5 .Results are reported in tab.13.Thanks to the high precision in the determination of α, the systematic uncertainty on γ 3 and γ 5 due to the uncertainty on α is negligible and we quote in tab.13 10: Results of the linear fits of the first two terms of eq.( 39) for β = 0.751800.
Looking at these results we see that there is a remarkable agreement between the values of γ 3 and γ 5 obtained with this method and those obtained in the previous section with the snake algorithm.We also see, as anticipated, that with this method there is a significative decrease of the uncertainty on the determination of γ 3 .We also see that the values obtained (with both methods) for β = 0.751800 do not agree within the errors with those obtained with the other two values of β.This suggests that, within the precision of our analysis, β = 0.751800 is still slightly outside the scaling window, while the data for β = 0.756427 and β = 0.758266 agree between them thus showing a good scaling behaviour.
We quote as our final result, obtained combining together the values obtained at β = 0.756427 and β = 0.758266 with the present approach.These values are compatible within the errors with the value obtained with the snake algorithm at β = 0.756427.

Concluding remarks
We studied, using two different methods and different algorithms the correction beyond Nambu-Goto for the confining potential in the three dimensional gauge Ising model.We found a good agreement between the two approaches.The two largest values of β that we studied show a good scaling behaviour and lead to values for γ 3 and γ 5 compatible within the errors.Our final estimate for these parameters is γ 3 = −0.00048(4)and γ 5 = 3.0(4) × 10 −7 .
The value that we obtained for γ 3 agrees with the bound of eq.( 22) , while γ 5 is slightly below the bound of eq.( 21) which, inserting the value of γ 3 and keeping into account the uncertainty in the determination of γ 3 , becomes γ 5 > 1.6 × 10 −6 .This difference is most probably due to the truncation in the perturbative expansion.Keeping into account higher order terms, and including also γ 7 might fill this gap.
It is interesting to compare our result with the existing estimates for γ 3 for other LGTs.In [27]  Table 13: Results of the fits of A(N t ) according to eq.( 44) truncated at the order O(x 6 ) i.e.O(N 12 t ) . ) and the red one is the prediction of eq.( 44) with the best fit values of γ 3 and γ 5 for β = 0.758266.
it was shown that also for the SU (2) LGT in three dimension γ 3 is negative.The value quoted in [27] is γ 3 | SU (2) = −0.00037(6)which is similar, but not compatible within the errors, with the one we obtained here for the Ising model.On the contrary for SU (6) a positive value γ 3 ≈ 3 × 10 −4 was found [33,19,34].These values represent the first steps toward a classification of EST models for LGTs.Indeed, in the past years, one of the main problems of the EST description of Yang-Mills theories was its universality, i.e. the fact that it predicted essentially the same behavior (with only a mild dependence on the number of spacetime dimensions), for models as different as the three-dimensional Z 2 gauge model and the four-dimensional SU(3) Yang-Mills theory.This feature is now understood as a universality that manifests itself only at low-energy (or, equivalently, a side-effect of the high accuracy of the Nambu-Goto approximation of EST), while the details related to the gauge group (and, possibly, to the confining mechanism into play) are instead encoded in the γ i corrections, which are not expected to be universal.In particular, from the values quoted above, it seems that γ 3 for ordinary LGTs takes very small values, which seem to increase with the complexity and size of the gauge group.This should be contrasted with the case of the 3d U (1) model where sizeable deviations from the Nambu-Goto prediction were observed in several quantities [49,66,67] which most likely point to a much larger, positive, value of γ 3 .Finally, let us add a final comment on the numerical side of our analysis.As we have seen, duality plays a crucial role in our analysis and for this reason the approach discussed in this paper is particularly suited for abelian gauge theories, however, apart from the numerical convenience, there is no obstruction to extend the flux based method discussed in sect.6,given enough computational power, also to non-abelian models.Moreover we have seen that with this approach the error in the determination of A(N t ) is dominated by the statistical uncertainty while the systematic error due to σ and α is fully negligible.This means that there would be in principle no obstruction to improve the estimates of γ 3 and γ 5 with a larger statistics.

A Appendix: Evaluation of α(β)
A first approximation for α can be obtained assuming the known scaling behaviour for the string tension in the 3d gauge Ising model which leads to a simple and elegant expression for α α(β) = 2νσ 3N 2 s (β c − β) However this is not enough for our purposes.In order to estimate higher order effective string corrections we need to evaluate the flux density with an uncertainty smaller than 1% and thus it is mandatory to include in the expression the next to leading terms of the scaling function.Following [50] and using the results of [41] we can approximate the scaling function as σ(β) = σ c t 2ν × (1 + at θ + bt) , where t = β − βc is the dual of the reduced temperature, θ = 0.5241 (33) is the next to leading scaling exponent and the constants take the following values: σ c = 10.083(8), a = −0.479(26) and b = −2.12(19).
Inserting this correction in the definition of α and approximating for simplicity 2θ ∼ 1 we finally obtain This is the expression that we used to obtain the values listed in tab.8.

Figure 1 :
Figure 1: Values of A(N t ) for the three different values of β.The horizontal line corresponds to the first order approximation in which the EST is truncated to the Lüscher term only, the violet curve is the Nambu-Goto prediction (truncated at the order x 6) and the red one is the prediction of eq.(44) with the best fit values of γ 3 and γ 5 for β = 0.758266.

Table 3 :
Results of the snake algorithm for β = 0.756427 and N t = 24.

Table 4 :
Results of the fit with eq.(30) for β = 0.751800.In the last column the reduced χ 2 of the fits.

Table 5 :
Results of the fit with eq.(

Table 8 :
only the statistical errors of the fits.Some information on the simulations.

Table 9 :
Results of the algorithm for β = 0.756427 and N t = 24.