SU(2) lattice gauge theory in 2+1 dimensions: critical couplings from twisted boundary conditions and universality

We present a precision determination of the critical coupling beta_c for the deconfinement transition in pure SU(2) gauge theory in 2+1 dimensions. This is possible from universality, by intersecting the center vortex free energy as a function of the lattice coupling beta with the exactly known value of the interface free energy in the 2D Ising model at criticality. Results for lattices with different numbers of sites N_t along the Euclidean time direction are used to determine how beta varies with temperature for a given N_t around the deconfinement transition.


Introduction
It is widely accepted today that the deconfinement transition in pure SU(N) gauge theories at finite temperature is driven by the dynamics of center vortices [1]. In the vortex picture of confinement, Wilson loops acquire a disordering phase factor from every vortex that they link with. The area law for timelike Wilson loops in the pure gauge theory comes from the percolation of spacelike vortex sheets in the confined phase. Their free energies have been measured over the deconfinement phase transition at finite temperature in the 3 + 1 dimensional pure SU (2) gauge theory from ratios of partition functions with 't Hooft's twisted boundary conditions in temporal planes, forcing odd numbers of Z 2 center vortices through those planes, over the periodic ensemble with even numbers [2,3]. A Kramers-Wannier duality is then observed by comparing the behaviour of these center vortices with that of 't Hooft's electric fluxes which yield the free energies of static charges in a well-defined (UV-regular) way [4], with boundary conditions to mimic the presence of 'mirror' (anti)charges in neighbouring volumes. This duality follows that between the Wilson loops of the 3-dimensional Z 2gauge theory and the 3D-Ising spins, reflecting the universality of the center symmetry breaking transition. Here we study the vortex free energies of pure SU(2) in 2 + 1 dimensions over the deconfinement transition because the relevant interface free energies of its universal partner, the self-dual Ising model in 2 dimensions, are all known analytically. Moreover, the vortex free energies in 2 + 1 dimensions are much cheaper to simulate and discretisation effects vanish more rapidly than in 3 + 1 dimensions. Together these reasons allow for numerical studies of much higher precision.
In Sections 2 and 3 we briefly outline the basic concepts and our numerical procedure. Our results are presented in Sec. 4. These include the precise determination of the critical coupling for various lattices with up to N t = 9 sites in the Euclidean time direction, an analysis of the finite volume corrections and a brief comparison with the corrections to scaling in the 2D Ising model. We then determine how the critical coupling depends on N t , including 1/N t corrections, and use this result to derive how the lattice coupling β varies with temperature for a given fixed N t around the deconfinement transition. This is needed, for example, for a detailed finite-size scaling analysis of the vortex free energies in 2+1 dimensional SU(2) [5].

Concepts and Methods
For pure SU(2) gauge theory, 't Hooft's twisted boundary conditions fix the total number of Z 2 vortices modulo 2 through each plane of a finite box [6]. Twist in a plane corresponds to an ensemble with an odd number of Z 2 vortices through that plane.
In a L 2 × 1/T Euclidean box, we can distinguish between two types of twist. Magnetic twist is defined in the purely spatial plane and forces vortices that run along the temporal direction. They may spread independently of the temperature T and play no role in the deconfinement transition: their free energy is expected to vanish in the thermodynamic limit at all temperatures which has been demonstrated explicitly in 3 + 1 dimensions [7]. On the other hand, vortices from twist in the two temporal planes are squeezed more and more as T is increased. They may no longer spread arbitrarily and this is what drives the phase transition.
In this paper we're interested in configurations with a twist in one of the temporal planes. If we denote the partition function of this ensemble by Z tw (L, T ), its free energy per T is defined via the ratio where Z 0 is the partition function of the periodic ensemble. As we approach the thermodynamic limit, Z tw /Z 0 converges to a non-trivial fixed point at the critical temperature. It can therefore be used as a phenomenological coupling in the same manner as the Binder cumulant. Typically, one uses the pairwise intersections of Binder cumulant curves from different lattice sizes to estimate the critical coupling or temperature [8,9]. Using a ratio of partition functions in this way was proposed by Hasenbusch for the 3D Ising model [10].
Here we can go a step further. Since SU(2) in 2+1 dimensions is in the same universality class as the 2D Ising model, the the value of Z tw /Z 0 at the fixed point is exactly known. A single temporal twist in SU(2) corresponds to a single anti-periodic direction on a N × N Ising lattice. This forces an odd number of spin interfaces perpendicular to the anti periodic direction. The corresponding ratio of partition functions Z ap /Z pp in the Ising model gives the free energy of the system in the same way as equation (1). In both cases, configurations of minimum action/energy dominate the partition function. Thus, in the thermodynamic limit, the free energy of a single spacelike vortex in SU(2) is identified with the free energy of a single interface in the 2D square Ising model at the respective critical points. This universal value is given by [11] lim N→∞ Z ap (T c )/Z pp (T c ) = 1/(1 + 2 3/4 ). ( For SU(2) on the lattice, 1/T = N t a where N t is the number of sites in the time direction and the lattice spacing a ≡ a(β) depends on the coupling. So the critical temperature T c corresponds to a critical lattice coupling β c,∞ for each N t . The subscript ∞ reminds us that we only have a strict phase transition and hence critical coupling in the limit of infinite spatial volume. Still, the intersection of Z tw /Z 0 with the universal value (2) gives a reliable estimate of β c,∞ (N t ) provided that the spatial length L = N s a is large enough. This will be more precise in general than the estimates obtained via pairwise intersections. 1 We assume a finite size scaling (FSS) behaviour of the form where ω is a correction to scaling exponent that should be approximately independent of N t , and ν = 1 is exactly known from the 2D Ising model. Furthermore, β c ≡ β c,∞ (N t ), and we define a 'pseudo-critical coupling' in a finite volume, β c (N t , N s ), by the requirement that the corrections to the universal value in  (3) vanish. These estimates then converge to the infinite volume critical coupling β c,∞ (N t ) as Here, the coeffiecient d = c/b is an N t dependent fit parameter. For notational simplicity we will hereafter drop the subscript ∞ on the critical coupling. Unless N s is explicitly mentioned, β c refers to the infinite volume limit. We have used a variety of spatial lattice sizes and the above scaling formula (3) to first find β c (N t , N s ) for N t = 4, 5, 6, 7, 8 and 9 with high precision, in order to then determine the corresponding β c (N t ) from fitting the volume dependence of β c (N t , N s ) by (4).

Numerical recipe
Twist was implemented in the usual way [13]. We kept periodic boundary conditions but flipped the coupling β → −β of a stack of plaquettes perpendicular to the plane of the twist. In other words, we introduced a Z 2 Dirac string that is closed by lattice periodicity. The result is a transformation of the usual Wilson action.
In practice, the overlap of Z tw and Z 0 is poor. To overcome this we interpolated in the number of flipped plaquettes using the snake algorithm of Ref. [14]. This was combined with the other variance reduction tricks therein.
For each combination of N t and N s we performed simulations for ∼ 10 values of β around the intersection of Z tw /Z 0 with the exact Ising value (2). Random errors were estimated via the boostrap method. Since SU(2) lattice gauge theory is less computationally expensive in 2+1 dimensions than in 3+1 dimensions, we were able to perform a very large number of measurements. 1−30 million configurations were used for each β, depending on the lattice size.
It turns out that the free energy F tw (β) = − ln Z tw /Z 0 has less curvature than Z tw /Z 0 near the critical point, so it's a better candidate for linear approximation. Therefore we performed least squares fits, with parameters f 1 and β c , of the form  to obtain estimates β c = β c (N t , N s ) for the critical lattice coupling. In each case the reduced χ 2 was ∼ 1, which justifies the linear ansatz (5). See Fig. 2 for some representative fits. Note, however, that a small reduced χ 2 does not exclude the existence of significant systematic errors in β c (N t , N s ). On finite lattices, F tw (β) has positive curvature near the critical coupling, so β c tends to be underestimated. To control this, we carefully chose the size of our fitting windows. For each N t , we performed precise measurements of F tw in a quadratic fitting window for one or more of our smallest lattices. By writing F tw as a function of the finite size scaling variable where t = T/T c − 1 is the reduced temperature and ξ = ξ 0 ± |t| −ν are the correlation lengths for T < > T c with ν = 1 for the 2D Ising model, we were able to translate this data to larger N s . See Fig. 3 for an example of FSS data collapse in a large window for several lattice sizes. 2 The relationship between temperature and coupling is described in Section 4.4. It requires a paramerisation of β c vs N t = 1/T c a c , which we roughly obtained using literature values of the critical coupling and our own preliminary results. The translated data was used to estimate the slope and curvature of F tw (β) near the critical point for large lattices without performing additional simulations. We then adjusted the linear fitting windows such that the systematic errors in β c (N t , N s ) should be less than one quarter of the quoted random error for each extrapolated β c (N t ).

Determining β c (N t )
We have obtained estimates of β c (N t ) for every N t between of β c (N t , N s ) vs N s for N t = 4, 5 and 6. The data was fitted according to the equation (4) with all three parameters free. The plots for N t = 7, 8 and 9 look very similar. Note the rapid convergence to the infinite volume values, which Hasenbusch also observed for the pairwise intersection of Z ap /Z pp curves in the 3D Ising model [10]. He found much more rapid convergence than for the intersections of Binder cumulants. We summarize our results in Table 1. For each N t we include the estimates β c (N t , N s ) from our two largest lattices as well as the fitted values of β c (N t ). It's clear from the reduced χ 2 s that the data is very well described by the FSS ansatz (4).
For N t = 4, 5 and 6 we were able to surpass the precision of current literature values of the critical couplings by two orders of magnitude. These lattices also gave us good precision for the correction to scaling exponent ω. Our results for N t = 7, 8, 9 are somewhat less precise, especially the fitted values for ω. This is because the lattices used had smaller aspect ratios and we collected fewer statistics. Still, ω is consistent between each of the fits. In all, we obtain a weighted average for ω of 1.61 (9), which is also consistent with the value of 1.64 obtained by Engels et al. [15] in early study of Polyakov loop averages.
For reference, we have included extrapolations of the crit- Figure 5: Critical beta estimates vs N s relative to their infinite volume limits for N t = 4, 5 and 6, and the corresponding fits to (4).  ical coupling with ω fixed at 1.61. Due to correlations, the quoted errors for β c (N t )| ω=1.61 should be taken with a grain of salt. Nevertheless, fixing ω may lead to more accurate values of β c (N t ) for our large N t results.
We furthermore obtain the fit parameter d(N t ) in (4) for each of the six N t values. This in turn allows us to fit its N t dependence. Using a two parameter form d fit (N t ) = d γ N γ fit t with a single effective exponent γ fit , a fairly good description is obtained with γ fit = 3.97(7) (and d γ = 0.06 (1)). This fit works best for the smaller N t but it deteriorates somewhat towards N t = 9. The N t dependence of d(N t ) might well be determined by several competing terms with nearby exponents and is thus difficult to extract reliably from the data. Alternatively, using 3 parameter fits of the form we obtain γ fit = 3.8(2) for δ = 0, γ fit = 3.7(4) for δ = ω, or γ fit = 3.5(6) for δ = ω + 1, for example, where we have used ω = 1.61 corresponding to our global fit for the correction to scaling exponent ω. Because ω + 2 = 3.61, this suggests that one might also try a form with ω fixed at 1.61. This form is particularly interesting because it will allow us below to parametrise the finite-size corrections entirely in terms of the aspect ratio in the form N t /N s . Unfortunately, the coefficients d 0 and d 1 of the two subleading powers in ω are much too correlated to determine all the 3 parameters in this fit reliably from the available data. But we can get good fits with either d 1 = 0 or d 0 = 0, see Fig. 6. Keeping the resulting two parameters fixed afterwards, a further fit to the remaining one (d 1 or d 0 ) then reproduces a zero result with very high accuracy (i.e., values smaller than 10 −4 in either case). It's interesting to check how much was gained by using the exact universal value (2) instead of pairwise intersections.  Figure 6: The parameter d(N t ) from the fits of β c (N t , N s ) to Eq. (4) together with 2 parameter fits of its N t dependence via the form in Eq. (8) with d 1 = 0 (long-dashed) and d 0 = 0 (short-dashed), respectively.
To this end we performed the pairwise intersection method for N t = 4. See Table 2 for the results. The left side of the table shows the intersection coupling for lattices with N s and N ′ s = N s /2. On the right hand side are comparative results from intersections with the universal Ising reference line.
Since the simulations were not catered for pairwise intersections, it may be unfair to directly compare the extrapolations in Table 2. Still, due to the rapid convergence of β c (N t , N s ), our data was quite well centered around the intersection points except for the smallest lattices. And knowing what value of Z tw /Z 0 to concentrate the efforts around was a big advantage of our method. Anticipating the location of pairwise intersection points is much more troublesome. In all, it's clear that the ex-N ′ s − N s intersection coupling N s β c (4, N s ) from (5)   ploitation of the universal number F tw (β c ) = ln(1 + 2 3/4 ) gave a significant boost to the precision of our results.

2D Ising model
For the sake of comparison we repeated our procedure for the square 2D Ising model. Taking the exact solutions for Z ap and Z pp on N × N lattices from Ref. [19], we used Mathematica [20] to find θ = k B T/J at the intersection of Z ap /Z pp with the universal value (2). Here J is the coupling of nearest neighbour spins. In Figure 7 we plot the results for N ∈ [100, 640]. The error bars are representative only (the errors are limited only by the working precision and should be smaller than 10 −16 ).
Again we fit the data with a FSS ansatz of the form In this case, we hold θ c,∞ fixed at the known value of 2/ ln(1 + √ 2). We obtain for the correction to scaling exponent ω = 2 + δ with δ → 0 + as we increase the lower bound N min used in the fit, i.e., ω tends towards 2 from above (δ starts at around 2 · 10 −4 for the full range of N shown in Fig. 7, and it falls below 10 −5 at N min around 400). Since the exponent of the leading irrelevant operator that breaks rotational invariance is predicted to be exactly 2 [21], our result is consistent with the conjecture of Ref. [22] that the only irrelevant operators that appear in the 2D nearest neighbour Ising model are those due to the lattice breaking of rotational symmetry. This may be tested on a triangular 2D lattice where the leading irrelevant operator to break rotational invariance leads to ω = 4 instead, while the leading rotationally invariant operator would give an isotropic correction to scaling with ω = 2 in either case [22].
On the other hand, our correction exponent for SU(2) is clearly at odds with ω = 2. In this case, it's possible that there exists an irrelevant operator that is not present in the 2D Ising model or the corresponding conformal theory. It may be more likely, however, that our exponent is really an effective exponent. When there are several nearby competing exponents it is extremely difficult to extract the smallest one from simulations. 3 3 It was noted in [15] that the observed correction to scaling exponent of 2+1-

β c vs N t
In 2+1 dimensions the coupling g 2 3 has the dimension of mass and sets the scale for the theory. The bare lattice coupling is then given by [17] Substituting this expansion into β gives Note now that T = 1/ (N t a). So at criticality we have The first correction to the bare lattice coupling gives the yintercept of β c (N t ) vs N t while the second is a correction to linearity for small N t . In Fig. 8 we plot our results for β c (N t ), using the fitted values with ω free. We include also the estimate β c (3) = 4.978(35), which is the average of the literature values 4.943 (13) [16] and 5.013 (15) [18]. For the uncertainty we use their standard error, which is simply half their difference. The actual error used here has very little influence on the results. The data is fitted with both a straight line and with a function of the form (12). To the naked eye, it would appear that the linear approximation is good all the way down to N t = 3. However, our data is precise enough to pick up the deviation from linearity. In fact, the reduced χ 2 of the linear fit is ∼ 70. With, the 1/N t correction the dimensional SU(2), ω = 1.64 in their case or ω = 1.61 (9) in ours, agreed well with some predictions for the universality class of the 2D Ising model which included 1.6 [23]. At the time it was discussed whether such non-integral correction exponents could arise in other ferromagentic models of this class, and whether the corresponding correction amplitudes happened to vanish identically in the exactly solvable pure Ising model, see [24]. This seems to be ruled out by a conformal field theory analysis: there is no irrelevant operator with ω < 2 in any unitary model of the 2D Ising class, see [25,22]. reduced χ 2 drops to 1.5. Thus, the data is very well described by Eq. (12). We obtain the parameterisation Consequently, from Eq. (12) for N c = 2, we obtain This is compatible with the estimate T c /g 2 3 ≃ 0.385 with an error of ±0.010 as quoted in Ref. [18]. Note that it is the 1/N t corrections included in (12) that lead to a somewhat lower value here as compared to [18] which is however within their estimate of the systematic uncertainties. In order to translate our estimate into units of the zero-temperature string tension σ we use the mean value with standard error of the four N = 2 values for √ σ/g 2 3 listed in [26] as pairs of upper and lower bounds including some systematic uncertainty, which gives √ σ/g 2 3 = 0.3347 (5). With uncorrelated error propagation, this together with our estimate (15) then corresponds to which agrees very well with the corresponding result of [16], T c / √ σ = 1.1224(90). From Eq. (13), the other constants in Eq. (12) are analogously determined as c 1 = −0.176(5) , c 2 = 0.0675 (5) .
Eq. (13) can be used to obtain accurate estimates of the critical coupling at large values of N t . We can furthermore include the finite-size corrections in our 'pseudo-critical coupling' (4) for the intersection, at finite N s , of the vortex free energy F tw in (5) with the universal value ln(1 + 2 3/4 ). If we assume a form as given in Eq. (8), these corrections can be expressed entirely in terms of the aspect ratio A ≡ N t /N s . In particular, the leading corrections for small A can then conveniently be combined with the large N t expansion (12) as follows, where we again used N c = 2 and ρ = ω + 1/ν. Since ν = 1 in the 2D Ising model, our global fit to the correction to scaling exponent ω amounts to ρ = ω + 1 = 2.61(9). If we assume d 1 = 0, our fits yield and a reduced χ 2 of 2.6. The corresponding fits are shown in Fig. 6. Apart from this systematic uncertainty in the otherwise elegant parametrisation (8) of the finite-size corrections, we have therefore determined all parameters in (17).