Simulating the All-Order Strong Coupling Expansion V: Ising Gauge Theory

We exactly rewrite the Z(2) lattice gauge theory with standard plaquette action as a random surface model equivalent to the untruncated set of its strong coupling graphs. By extending the worm approach applied to spin models we simulate such surfaces including Polyakov line defects that randomly walk over the lattice. Our Monte Carlo algorithms for the graph ensemble are reasonably efficient but not free of critical slowing down. Polyakov line correlators can be measured in this approach with small relative errors that are independent of the separation. As a first application our results are confronted with effective string theory predictions. In addition, the excess free energy due to twisted boundary conditions becomes an easily accessible observable. Our numerical experiments are in three dimensions, but the method is expected to work in any dimension.


Introduction
This paper continues our series where we try to build on and generalize an idea by Prokof'ev and Svistunov [1] on alternate ways of formulating and simulating simple statistical systems and Euclidean lattice field theories. In our context, their idea has two aspects. One is that, instead of summing over the field configurations in the original form, one exactly rewrites the models as a finite or infinite sum over associated strong coupling or high (classical) temperature graphs. The second aspect for spin models is, that enlarging the class of such graphs from those of the partition function to a larger set associated with the two point correlation function makes simulations ('diagrammatic Monte Carlo') both simpler and much more efficient with respect to critical slowing down. In [2] it was shown that we may in addition exploit some freedom to design the enlarged ensemble to our advantage and thus achieve an excellent signal to noise ratio for interesting observables. It is this aspect that we will extend in this paper to Polyakov line correlators in Abelian gauge theories while critical slowing down is unfortunately not eliminated at the same time.
Some remarks are in order here. The standard usage of strong coupling expansions is to evaluate a truncated series for some suitable observables and then take the thermodynamic limit within this approximation. These expansions in powers of some β usually have finite radii of convergence related to singularities or phase transitions. In systems with a large but finite number of compact degrees of freedom however, these expansions converge for many quantities like for instance the partition function itself. Therefore, for such systems -and we never simulate anything else -the untruncated set of strong coupling graphs furnishes an exact reformulation of the theory in which we can attempt to apply stochastic summation methods. For a given action or Hamiltonian the graphs derived from it reproduce finite lattice observables exactly and not just universal features. For spin models we thus encounter representations in terms of loops drawn on the lattice and the obvious generalization to gauge models studied here consists of surfaces. The admissible graphs are restricted by constraining rules, for instance on the number of lines that may touch at a site, and we sum over such constrained objects. Solving the constraints in terms of independent new variables would complete our reformulation to a duality transformation, which is however not done in [1]. What the reformulations have in common is that rather different observables may be advantageously computed in one or the other. In models similar to Ising models duality [3] relates the strong coupling expansions in one model to the weak coupling expansion in another. It is trivial, but amusing to remark, that a diagrammatic Monte Carlo of the weak coupling (low temperature) expansion of an Ising model is nothing but a standard simulation, as spin-flip related pairs of configurations correspond to one naive weak coupling graph. In the systematic weak coupling approximation one just truncates these contributions guided by the size of their Boltzmann weights.
Our previous work in this series extended the method from the Ising model [2] to the O(N ), CP(N − 1) models and to fermions, see [4] for an overview and more references. For fermions an unsolved sign problem hampers simulations in more than two dimensions while in the other cases the technique is expected to be efficient in arbitrary dimension. Studies of O(N ) as a special case of certain loop graph models on particular lattices are also found in [5], [6].
There have been other recent attempts to apply the 'worm' method to Abelian gauge models [7], [8], [9], [10], [11]. The important new feature in the present study is in our mind however, that we succeed in generalizing the low noise estimators for fundamental correlations from spin to gauge systems. Also the rôle of twisted boundary conditions on the torus will be discussed which leads to interesting topological universal finite size observables.
As we get access to precise results for Polyakov line correlators, a comparison with low energy effective string models suggests itself as a first interesting application of our simulation technique. Such studies form an active research area and we here offer a demonstration of the usefulness of our method in the context. A very profound such study is beyond the scope of this publication. The field has a long history with a boost, as far as we see, after the papers by Lüscher and Weisz [12], [13] which interestingly were also triggered by algorithmic improvement. The low energy description of large Wilson and Polyakov loops in confining gauge theories was pushed to the two loop level in the spirit of a general effective field theory approach reminiscent of chiral perturbation theory vis a vis QCD. It was found that to this order and in three dimensions and taking into account all symmetries, no free low energy coupling constants enter. This investigation was later even extended to the three loop level in [14] so that we have remarkable absolute predictions to check with our low noise long distance results for the Polyakov loop correlator.
For the Z(2) Ising gauge theory in three dimensions that we adopt here as a test case, there actually exists an enormous body of precise results in the literature, see [15], [16] for recent papers with further references. Here the duality of the gauge theory with the spin model is exploited and the simulations are conducted there. The method is described in detail in [17]. A Wilson loop in the gauge theory translates into a ratio of partition functions in the spin model differing by flipped bonds. Such a ratio is factorized into many ratios differing in single bonds, which are evaluated in separate simulations. While this is an intrinsically three dimensional method, the same is not true for our more direct approach, although we here test it in D = 3 as well.
The paper is organized as follows. In Section 2. we introduce the model, its boundary conditions, rewrite it as a surface ensemble and collect some formulae of the effective string description. This is followed by the development of our Monte Carlo method for the surface ensemble in Section 3. and a report of our numerical experiments in Section 4. We end on conclusions in Section 5. and an Appendix listing numerical data.

Z(2) gauge theory
In this section we setup our model on a hypercubic lattice in D Euclidean dimensions. We include the definition of fluctuating twisted boundary conditions. Then we introduce the one-to-one reformulation as a surface model, which generalizes the loop (re)formulation of spin models.

Twisted boundary conditions
We consider a gauge field σ(x, µ) ≡ σ µ (x) = ±1 ∈ Z(2) defined on the links of a D-dimensional hypercubic periodic lattice of extent L µ in directions µ = 0, 1, . . . , D − 1. The standard Wilson action is defined on plaquettes (x, µ < ν) by The plaquette dependent background field P (x, µ, ν) ∈ Z(2) will be useful later and may be first imagined to be unity until further notice. Generalized periodic boundary conditions 1 demand that the gauge field is periodic up to gauge transformations 2 [18] σ(x + L νν , µ) = σ(x, µ)τ ν (x) τ ν (x +μ) (2) specified by fixed 'external' transition functions τ ν (x). Because the shifts form an Abelian group it is necessary that shifts x → x + L ν → x + L ν + L λ come with the same gauge transformation as the double shift in the opposite order, which however still leaves the possibility The case γ νλ = −1 is possible because the action of τ ν on gauge fields is independent of the global sign of the transition function. Then we have twisted boundary 1 A fruitful point of view here is to consider all fields on the infinite lattice with a finite subset of independent variables as all others are 'locked' by periodicity. 2 Gauge invariant densities will thus be periodic in the usual sense. This generalizes antiperiodic boundary conditions in the Ising model which leave Z(2) symmetric composites periodic.
conditions in the νλ plane while we call planes with γ νλ = +1 untwisted or just periodic below.
If we gauge transform in the ordinary sense, this implies the following change of the transition functions they transform like parallel transporters on a set of 'super-lattices' with spacings L µ . The γ µν represent the gauge invariant content of the transition functions. For a given set of twists γ µν we now define reference transition functions where we round upwards to an integer ('ceil') in the exponent. Then the product ν (x) has trivial twist and can be gauged to unity. Hence we may assume that the transition functions have the form τ (γ) and then the gauge field is periodic except (possibly) for By a further change of variables, we may arrive at fully periodic σ(x, µ) again if we absorb signs into P γ (x, µ, ν) = (γ µν ) δ xµ,0 δ xν ,0 with periodic δ symbols (period L µ , L ν respectively). Thus for each twisted plane we have a D − 2 dimensional set ('stack') of negative plaquettes that represent a background flux or disorder ('vortices') concentrated on a line, sheet, ... for D = 3, 4, . . .. This step is analogous to absorbing antiperiodic boundary conditions into a background gauge field, see [19] for example.

Gauge theory as a surface model
Next we introduce the partition function with current insertions In this formula we sum over strictly periodic σ(x, µ) = ±1 independently at all 0 x µ < L µ and µ, and j(x, µ) ∈ {0, 1} is a periodic external field. The number of independent sites, links, plaquettes and 3-cubes (for later use) is Performing a local gauge change of variables (4) with periodic ρ (x) one derives with the divergence being the sum over the 2D links surrounding x. As ρ can be arbitrary this shows that only divergence free currents (in the Z(2) sense), for which holds at all sites, yield nonzeroZ γ [j]. In addition we may flip the gauge field σ µ on any D − 1 dimensional hyperplane x µ = z orthogonal to the µ-direction ('layers').
A similar argument as before yields the requirement that the layer sums must be even for all layers, which are hence pierced by an even number of current quanta. Equivalently we may say that the link field j (x, µ) must have vanishing Z(2) winding number with respect to all torus directions. Using e βσ = cosh(β) n=0,1 t n σ n , t = tanh(β) (15) for each plaquette, we introduce a field n(x, µ, ν) ∈ {0, 1}. Then we may average over the original σ(x, µ) which leave behind constraints only. We arrive at The sign is given by where for each n (x, µ, ν) configuration we have introduced the 'wrapping' numbers They are topological quantities 3 that count (modulo two) how many times the surfaces made of plaquettes with n (x, µ, ν) = 1 wind around planes and represent two dimensional generalizations of the winding numbers of Ising loops. To define the divergence of the plaquette field, we extend 4 The constraint (mod 2) combines the 2(D − 1) plaquettes surrounding each link We note that consistently this constraint can only be satisfied by j that obey (13) and (14). The proof looks slightly funny, for example An easy to interpret case is a configuration j ≡ 0. Then the plaquettes n(x, µ, ν) = 1 form a surface which may branch and consist of disconnected components. The zero divergence condition demands that they are closed and have no boundaries, each link is surrounded by an even number of surface elements. Boundaries arise if j is nonzero and, due to (13) and (14) they may be written as a superposition of contractable closed current loops. It will be helpful for the reader to work out the closely analogous but geometrically simpler loop formulation of the Ising spin model in this language. In analogy to the steps taken for spin models we now form the current ensemble where the non-negative weight R −1 [j] will be specified later and Z without subscript stands for trivial twist Z = Z γ≡1 . Note that we here include all wrapping numbers in the sum over n with no extra signs, which corresponds to trivial twist in all planes. Expectation values in this ensemble are given by In addition we define 'vacuum' expectation values on the subset of defect-free configurations. Such quantities do not depend on the choice of R. It is now obvious that ratios of partition functions of different twist, which lead to interesting observables, are given by topological observables We shall later see that there are Monte Carlo algorithms that are ergodic only in the sector of trivial wrapping in some or even all planes. It is clear now that the corresponding ensembles can be considered as arising by dynamically averaging over twisted and untwisted boundary conditions for these planes (fluctuating boundary conditions).
By differentiating Z δ[∂ * µ n µν ] with respect to β we obtain the relation between the average plaquette E of the original theory and the total surface area in vacuum configurations Below we shall find that close to the critical point in D = 3 we have the rather high values n 0 ≈ 1/3 and E ≈ 0.95.

Polyakov loop ensemble
The source j (x, µ) allows to place a large variety of defect configurations like arbitrary Wilson loops. We here define however a highly restricted framework involving only two Polyakov lines in the 0-direction that are located at u = (u 1 , . . . , u D−1 ) and an analogous v. The corresponding conserved current is For coinciding u = v it vanishes and there is no defect. We consider the ensemble Here the weight R[j ( u, v) ] has been specialized to 0 The two point function (in the original gauge theory) of the Polyakov loop operator is now given by a ratio of expectation values of counters if we normalize ρ( 0) = 1.

Contact with transfer matrices, effective string theory
The Polyakov line correlator is the ratio of partition functions with and without two static charges and thus given by where e −Vn( x) labels the corresponding eigenvalues of the transfer matrix in the 0-direction with charges separated by x and integer weights w n account for degeneracies [13].
Alternatively we may consider a transfer matrix in a spatial direction, k = 1 for example. Then the Polyakov line operator excites a flux state whose energy strongly depends on its length L 0 (now a 'transverse' direction) and we may consider the correlation where we have projected to zero momentum for the k > 1 directions and v n is a nontrivial matrix element in this case.
In effective string theories one can compute the energiesẼ n as an asymptotic expansion in 1/L 0 . On the basis of the Nambu-Goto action for example one predicts for D = 3 the relation [20] [13] for the ground state energy which is expected to be relevant for z, s → ∞, and σ is the (zero temperature) string tension. Formula (32) implies asymptotic expansions and, expanding s in z −1 , It is interesting to note that in general the Nambu Goto picture is not expected to hold to all orders in the asymptotic expansion, but the terms exhibited above have been uniquely derived for general effective low energy actions just restricted by symmetries [13], [14].
In addition it is shown in [13] that, up to rotational symmetry breaking cutoff effects, the original correlation is given in terms of theẼ n by the expansion with r = | x | and the modified Bessel function K ν . Note that here an infinite volume is assumed.
The estimator for C (y) in our extended ensemble (27) follows from (29) and is given by the simple observable that simply follows from the statistics of the defect-line separations. Note that this result is independent of the choice of ρ which we shall hence optimize with regard to the numerical simulation of (27). In practice one will symmetrize over spatial directions if they are symmetric with respect to extent and boundary conditions.

Simulation algorithms
We now introduce several algorithms to simulate the surface representation of Z(2) lattice gauge theory.

Defect conserving update
A very simple update step consists of the following sequence CF ('cube flip') of operations 5 • take a 3-cube c ≡ x, µ < ν < λ, • propose to flip all 6 plaquettes associated with this cube which is the set i.e. n (p) → 1 − n (p) for all p ∈ P, • accept this proposal with the Metropolis probability min (1, q) with Such steps may be iterated either with randomly chosen c or as a systematic sweep covering each c once in some order. This can be shown to lead to an ergodic algorithm at fixed j (vacuum graphs with j ≡ 0 for example) and within a fixed wrapping number sector. We expect (and confirm) however dynamical exponents close to two for such local updates that preserve the vacuum graph constraints. This is similar to just flipping links around plaquettes in the strong coupling loop representation of Ising spin models. A natural generalization of the worm algorithm [1] would be to allow to 'open up' vacuum graphs and allow excursions to the enlarged state space of (allowed) j = 0 (possibly traversing 'useful' configurations regarding observables of interest) and eventually returning to a vacuum graph. We have a long record of not really successful experiments of this type, see [9] for experiments with U(1) gauge theory. In Z(2) we have simulated correct ergodic algorithms where we included a chemical potential per defect link j (x, µ) = 1 to control their number. We could find values that have led to ensembles of randomly shaped (typically irregular) Wilson loops where also vacuum configurations (re-)appeared at a reasonable rate. None of these attempts has so-far led however either to fast dynamics or to easily interpretable observables. Therefore we come up with the hybrid below which combines cube flips with Polyakov line pair defects only where we have to tolerate however some critical slowing down.

Polyakov line moving update
We now describe worm-type updates referring to the ensemble (27). If (say) u is moved to a neighboring site u in the spatial direction i the corresponding line defect moves and the whole 'ladder' of L 0 plaquettes in the 0i plane 'above' this spatial link is flipped if the move is accepted. Such a move clearly preserves the constraint. To specify the procedure in detail it is helpful to define the auxiliary spatial link field Now we execute the following sequence PS ('Polyakov shift') of (standard worm) steps: • If u = v holds, we randomly re-locate both together to a new site on the D − 1 dimensional lattice.
• Pick one of the 2 (D − 1) spatial neighbors u of u, such that u = u ±î.
• Accept the proposal with the probability (pre-tabulated, of course) If this happens, all plaquettes in the ladder are flipped and u is changed to u . In the case of rejection the previous configuration is kept.

Complete update sequence
In this subsection we specialize to D = 3 for simplicity. One complete iteration of algorithm A1 consists of N x /L 0 repetitions of the following steps • Apply CF to 4L 0 cubes around the temporal line at u. We use a helical order and first consider the 4 cubes around the link (u, 0) with u 0 = 0, then at u 0 = 1 and so on.
• This is followed by n ps applications of PS. In all our simulations given below we took n ps = 16.
The cost of one iteration A1 is of the order volume like a conventional sweep. We here try to mimic the successful worm algorithm in so far as we attach our update moves to the defect, which now is a (straight) line instead of a point. A difference in the case at hand is that defect moves alone -within our restricted class as discussed before -are by far not sufficient for ergodicity. We have mimicked this conglomerate also in the Ising model by proposing flips of links around only those plaquettes of which a defect forms a corner. This works but offers no advantage in this case. In the algorithm just described the defects at u and v are not treated on the same footing -we only shift u -although, due to the relocation step in PS they can both reach all positions. This is in contrast to (27) where they play a symmetric role. Analogous options appear already for the point defects in spin models, where we found the asymmetric algorithm easier to program and equally efficient as a manifestly symmetric variant. We make no such attempt here.
As described before only temporal defect line pairs are included. Their migrations over the torus allow to change wrapping numbers in the 0k planes, but not in the purely spatial planes (only 12 in D = 3). Therefore A1 simulates an ensemble where one dynamically sums over γ 12 = ±1 with the other twists trivial.
A second version A2 of the algorithm results, if we allow a pair of Polyakov lines in any of the 3 directions, but always only one pair at a time. We denote by j (µ, u, v) the current with lines in the µ-direction with u, v locating these lines. In (27) an additional summation over µ is included, and D µ is the D − 1 dimensional sublattice of sites x with x µ = 0. In addition the formulae in subsection 3.2 have to be modified in the obvious way. The practical difference between A2 and A1 is simply that in PS, when u = v is encountered, a new value µ ∈ {0, 1, 2} is proposed together with a random position of u = v in D µ . The proposal is accepted with the probability min (1, L µ /L µ ). This results in an ensemble as in (21) where all wrapping numbers contribute. If all extensions L µ are the same the acceptance step is trivial and the resulting histograms of occurring u − v values can be combined. If this is not the case, they are collected separately and yield simultaneous results about Polyakov loop correlations in several geometries.

Rejection free Polyakov shifts
A potential problem with algorithm A1 (and also A2) can be small acceptance rates for the shifts in PS, if L 0 is large. After all we propose nonlocal albeit only one-dimensional changes to the configuration. Simulating L 3 lattices at the critical point, we indeed find acceptances falling with L which reach values of about 1.6% only for L = 64.
In [5] a so-called rejection-free version for general worm-type algorithms was proposed. If we consider a PS step in A1, then the total probability that one of the four conceivable moves takes place is given by Note that 1 − A > 0 is the probability of keeping u unchanged. It is easy to show [5] that if we replace PS by a step PSnr • if u = v holds, the same steps as in PS are taken, • for u = v one of the four moves discussed before is chosen and executed with the normalized probabilities then PSnr has the Boltzmann weight in as a fixed point. To construct a complete simulation algorithm A1nr we need to replace beside PS → PSnr also the CF step by one, where the Boltzmann ratio corresponding to (43) is formed for the acceptance step. Finally also a version A2nr is programmed, where A ( u; k) is suitably generalized to depend on the direction of the defect lines. Since in any case the relative weights of vacuum graphs are the same in (27) and (43) all expectation values (23) are as before.
We see that the term 'rejection free' refers to the Polyakov shifts in non-vacuum configurations only. Acceptance rates in CF are however at efficient levels in all our simulations. If now the remaining acceptance rate of PSnr out of vacuum gets too small, we may exploit our freedom of choosing ρ to improve this. For example the choice with c < 1 reduces the relative weight ρ −1 of vacuum configurations and thus enhances the probability to leave them. A few brief experiments allowed us to arrange for reasonable acceptance rates for these moves, too. We find this simple choice for ρ to be a good one at the critical point, where we will implement A2nr, while in the confined phase in connection with A1 other choices will be more favorable.

Simulations close to the critical point
We first report on a number of simulations at which corresponds to the dual of the estimate given for the spin model in 6 [22]. On each lattice of sizes L µ ≡ L = 6, 8, 12, 16, 24, 32, 48 we have executed 8 × 10 6 iterations of A2nr with ρ from (44) and c = (8/L) 2 after a few tests. All observed acceptance rates of CF are around 50% and the remaining acceptance of PS out of vacuum configurations drops from 34% at L = 6 to 23% at L = 48 with these choices.
One of the simplest observables is the average fraction of plaquettes participating in the random surface of vacuum configurations which is analogous (and equivalent, see (25)) to the average plaquette. We measure values like for example Θ = 0.33751 (2) at L = 6 and Θ = 0.33493 (4) at L = 48. In Fig. 1 we visualize a typical graph occurring in the L = 6 simulation, larger systems look even more cluttered. Histograms counting the frequencies of separations u − v are accumulated after each individual PS step and turn out to be rather flat, as expected in the critical theory, where the Polyakov correlation decays only slowly. During each iteration of A2nr we separately accumulate contributions to observables from the subset of vacuum configurations encountered during this iteration, like for example for Θ. In the end we perform an autocorrelation analysis of this time series to estimate errors as described in [24]. The such defined integrated autocorrelation time of Θ is shown in Fig. 2.
It is in units of A2nr iterations which cost proportional to L 3 . We see that while the absolute values are not too large on our lattices the rate of growth exhibits critical slowing down τ int ∝ L z with an effective dynamical exponent for our range of lattice sizes close to two. Note that we cannot make statements on the truly asymptotic dynamical exponent on the basis of these data.
Of greater interest than Θ are the topological observables (18). Taking into account the symmetries between all planes we here measure the set of observables R n = δ n,w 01 +w 02 +w 11 0 , n = 0, 1, 2, 3.
(47) Figure 1: A typical configuration on a 6 3 lattice at t = t c . The Polyakov line defects are shown as fat red lines.
Due to (24) each R n can be simply related to ratios of partition functions with twisted and untwisted boundary conditions and is hence expected to possess finite continuum or scaling limits. In particular, at the exact critical point there are definite finite values R * n for each R n in the limit L → ∞ at t = t c . In fact, this property may be used to determine t c . Note that the values R n are not independent, of course, but have to sum to unity. Our raw results on these quantities are compiled in Table 3 in the Appendix. We here consider the expected finite size scaling behavior close to the critical point in the form where ω is the exponent of the leading scaling violations. In our simulations at t = t c we find the results shown in Fig. 3. In the variable L −ω a linear approach to R * n is thus expected. Recent determinations of the exponent ω can be found in

Simulations in the confined phase
In the confined phase β < β c we simulate with algorithm A1, such that the distribution of u − v has a direct relation to the temporal Polyakov loop correlator (29). We do not generate surfaces wrapping around the 12 plane corresponding to dynamically summing over twisted and untwisted boundary conditions for this plane. The Polyakov correlator is however periodic in space in this situation. For our test we took the coupling β = 0.73107 with L 1 = L 2 = L = 64 and a series of inverse temperatures L 0 = 6, 8, . . . , 26, 28. Such values (but smaller L 0 ) have been adopted in [17] and a string tension in lattice units ofσ = 0.0440 (3) has been estimated based on earlier literature. In our simulations in the confined phase we want to use a ρ ( x) in (27) that essentially captures the decay of the two The massM is estimated fromσ by using the first three terms on the right hand side of (34). The sum over integer k, which makes ρ periodic, is actually truncated after a few terms as the remainder would be insignificant compared to roundoff due to the exponential decay of the Bessel function. We emphasize that any imperfection in ρ is related to the efficiency of the simulations but does not entail any systematic error. lines in our run for L 0 = 24. With the exception of very short torus distances we find a rather flat behavior (note the fine scale of this logarithmic plot) which leads to very uniform bin heights, which, up to correlation effects, would naively imply constant relative errors. Via (29) this allows to measure the correlator which for the same run is plotted in Fig. 5. To determine errors here, we had to store the time series for these separations rather than just the histogram. We see that the exponential decay can be followed over the whole lattice with relative errors staying small. The plot jointly exhibits 'on axis' and 'diagonal' separations which form a very smooth curve and thus exhibit small violations of rotation invariance. The spatial correlation length is of order one and therefore the periodicity is just barely visible close to r = L/2 and r = L/ √ 2. The closed string mass gapẼ 0 is estimated from (36) and our rather precise results are compiled in Table 1. A typical case of the determination of the mass gap is shown in Fig. 6. Effective masses m (y + 1/2) are determined by solving C (y) /C (y + 1) = cosh (m (y − L/2)) / cosh (m (y + 1 − L/2)). The horizontal lines show the error band from a fit deeply in the plateau region which leads to the entries forẼ 0 in the table. Observed autocorrelation times τ int for these derived quantities [24] range from 0.508(4) (L 0 = 6) to 0.598(6) ( L 0 = 28). The remarkable feature here is, that the errors in the effective masses do not grow with separation. This is due to our judicious adaptation of the simulated ensemble by choosing ρ. The situation is analogous to the one in the Ising spin model [2] where the reasons for this success are discussed in more detail.
What is the implication of our fluctuating boundary conditions in the 12 plane for the correlation (36)? The falloff in the direction x 1 is associated with the transfer matrix in this direction which acts on wave-'functionals' ψ [σ 0 (x 0 , x 2 ) , σ 2 (x 0 , x 2 )]. The Polyakov operator (28) acts on them by multiplication. The integration over links σ 1 in the Euclidean theory implies the inclusion of a projector on gauge invariant states in the thermal trace given by the path 'integral'. This refers to gauge transformations that are periodic in both directions. One may however in addition consider transformations [18] that are topologically nontrivial, here antiperiodic in x 2 , and physical states may be even or odd under them. By summing over twisted and untwisted boundary conditions in the 12 plane we include a projector also with respect to this quantum number. Under the assumption that the ground state is even, such a projection is no disadvantage. To analyze our data we perform a fit of the form To have visible errors at all we immediately plot the difference between the left and right hand side of such a fit against L −2 0 in Fig. 7. The fit here was derived from the subset L 0 12 of the data. The point with L 0 = 8 is much higher while L 0 = 6 is far off the panel. At these time extents the low temperature expansion has clearly and rather abruptly broken down. In fact, for our β the phase transition is close to L 0 ≈ 4 according to [17].
We note that (50) is exactly equivalent to the Nambu Goto form (32) if we identify Nambu Goto : In Table 2 we list a number of fits (50) with free c 0 , c 1 together with the resulting ratio r for whose error the correlation between c 0 and c 1 is taken into account.  We see a small but significant deviation of r from one by only about 1%. The most likely explanation for this are cut off effects in our opinion, given that the string tension σ in lattice units, which may be taken as indicative for their size, is about 0.044. It remains to be discussed to which terms in the expansion (33) our fits are actually sensitive. The errors of ourẼ 0 are 7%, 24% , 68% of the last term in (33) for L 0 = 10, 12, 14. The analogous numbers for the next term are 50%, 236%, 897%. This next term reads 5π 4 / (10368σ 3 L 7 0 ) and this is only taken as a model of conceivable next order terms. We thus have to conclude that in spite of our rather high precision we are just still sensitive to the last (known) universal term of order L −5 0 and cannot confirm or disprove the coefficient of L −7 0 here. In future simulations it will be desirable to simulate smaller lattice spacings to check the size of cutoff effects and possibly to enlarge the statistics to the point of seriously probing another term in (34).

Conclusions
We have attempted to generalize the worm algorithm for the Ising model to gauge theories with Ising link variables. The labeling of all-order (in tanh β) strong coupling graphs as configurations of constrained plaquette fields was straight forward. The essential method of [1] to achieve efficient updates of the graphs was to allow for a pair of point defects. The possible defects in the gauge case form a very large set of generalized loop networks. We could not identify among them a suit-  Table 1 and the second fit of the form (50) listed in Table 2.
able subset that strongly reduces critical slowing down while staying 'close to' the vacuum or simple Wilson loop defects. The defects in the spin model allowed at the same time for a very precise estimation of the fundamental correlation at large distance [2]. We successfully have generalized this aspect to the gauge model and could compute the Polyakov line correlator with similar precision. At the same time, the algorithmic efficiency is rather high on large lattices in three dimensions in spite of critical slowing down.
Other Abelian gauge theories have been considered in the surface representation, see for example [7], [9], [10]. The techniques for correlators in this paper can clearly be generalized to these cases, i.e. to the groups Z(N ) or U(1). Theories with SU(N ) variables, both gauge and the principal chiral spin models 8 are open problems in this context to which we hope to come back in the future.

Appendix
In this appendix we list the raw data on which Fig. 3 is based. Each line represents 8×10 6 iterations of algorithm A2nr as described in subsection 3.4 except for L = 48 where the number has been doubled. The error of these quantities seems to grow roughly proportionally to L 1.2 with a fixed number of iterations whose costs scale like L 3 .  Table 3: Results for the observables (47) from simulations on L 3 lattices with periodic boundary conditions in all directions at the estimated critical β given in (45).