Induced QCD II: Numerical results

We numerically explore an alternative discretization of continuum $\text{SU}(N_c)$ Yang-Mills theory on a Euclidean spacetime lattice, originally introduced by Budzcies and Zirnbauer for gauge group $\text{U}(N_c)$. This discretization can be reformulated such that the self-interactions of the gauge field are induced by a path integral over $N_b$ auxiliary bosonic fields, which couple linearly to the gauge field. In the first paper of the series we have shown that the theory reproduces continuum $\text{SU}(N_c)$ Yang-Mills theory in $d=2$ dimensions if $N_b$ is larger than $N_c-\frac{3}{4}$ and conjectured, following the argument of Budzcies and Zirnbauer, that this remains true for $d>2$. In the present paper, we test this conjecture by performing lattice simulations of the simplest nontrivial case, i.e., gauge group $\text{SU}(2)$ in three dimensions. We show that observables computed in the induced theory, such as the static $q\bar q$ potential and the deconfinement transition temperature, agree with the same observables computed from the ordinary plaquette action up to lattice artifacts. We also find that the bound for $N_b$ can be relaxed to $N_c-\frac{5}{4}$ as conjectured in our earlier paper. Studies of how the new discretization can be used to change the order of integration in the path integral to arrive at dual formulations of QCD are left for future work.


Introduction
In the strong-coupling limit of lattice gauge theories, gauge fields do not interact directly with each other, leading to a factorization of the link integrals in the path integral. This allows both for analytical investigations and for the construction of new simulation algorithms (e.g., [1][2][3][4][5]). Away from the strong coupling limit the self-interactions of the gauge field need to be taken into account and, with the standard actions, gauge integrals no longer factorize, spoiling the applicability of these strong-coupling methods. A particular way to overcome this problem is to reformulate the gauge action in terms of auxiliary degrees of freedom so that the gauge fields only couple to these unphysical degrees of freedom rather than among themselves (e.g., [6][7][8][9]). In this approach the gauge action is "induced" in a well-defined limit only after the auxiliary degrees of freedom have been integrated out. Typically this involves taking the limit to an infinite number of fields, rendering the resulting theories impractical for numerical simulations.
In [10], Budczies and Zirnbauer (BZ) developed a method which induces the pure gauge dynamics already for a fixed and small number of auxiliary bosonic fields. The key idea is to give up on the exact reproduction of the Wilson gauge action at finite lattice spacing in favor of an alternative lattice discretization of Yang-Mills theory which allows for a formulation in terms of auxiliary bosons. The vital ingredients for this idea to work are (a) the existence of a continuum limit and (b) its equivalence with continuum Yang-Mills theory. For the "designer action" (or weight factor) introduced in [10], BZ could show these properties for gauge group U(N c ) as long as the number of auxiliary bosonic fields, N b is larger or equal to N c . In QCD we are interested in gauge group SU(N c ) and in [11] (the first paper of this series) we adapted the BZ approach to this case, avoiding a spurious sign problem in the bosonization of the original formulation by a slight reformulation of the original weight factor. In particular, we could show the existence of the continuum limit as long as N b is larger than or equal to N c − 1 (or N c − 5/4 if we allow N b to be non-integer) and the equivalence with SU(N c ) Yang-Mills theory in the continuum limit. As in the original BZ paper, the latter could be shown in d = 2 dimensions if N b ≥ N c (or N b ≥ N c − 3/4 if we take N b to be non-integer), but it is a conjecture for d > 2. A brief review of the main findings from [11] which are of importance in this article is included in section 2.
In the present article we will investigate the conjecture numerically for the simplest nontrivial case, namely SU(2) gauge theory in three dimensions. We simulate both the standard theory with the Wilson plaquette action, which involves a single parameter β, and the induced theory with a fixed number of bosonic fields N b , which involves a single parameter α. We set the scale for both theories by computing the Sommer scale r 0 [12] from the static quark-antiquark potential. Matching r 0 from both theories gives us a relation between α and β. Using this relation we can compare other observables, such as quantities connected to the static qq potential and the finite-temperature phase transition. We find that the results agree very well already away from the continuum limit and that the agreement improves as the continuum limit is approached. A preliminary analysis of data from SU (3) gauge theory in four dimensions shows that the modified BZ method also works as expected, supporting the universality argument given in [10]. Therefore the modified BZ method, combined with a suitable strong-coupling approach, can be used to reformulate lattice gauge theories in a number of different ways. It will be very interesting to explore such reformulations in the future to see whether they may have advantages over the traditional formulation (and perhaps even solve the sign problem afflicting lattice QCD at nonzero density). We remark that an exact rewriting of the pure gauge action in terms of auxiliary fields has recently been achieved with Hubbard-Stratonovich transformations [13]. This leads to a qualitatively similar reformulation of the theory, even though the auxiliary fields are rather different. These two types of reformulations can thus be seen as complementary approaches with different properties concerning possible reformulations in terms of dual variables. This paper is structured as follows. In section 2 we review the results from the first paper in the series [11]. The matching between the lattice couplings in Wilson's pure gauge theory and in the induced gauge theory is discussed in section 3. In sections 4 and 5 we compare the results for the static qq potential and the deconfinement phase transition, respectively, before we conclude in section 6. The details concerning the simulation algorithms and the extraction of the observables as well as the details of the simulations for comparison between perturbation theory and numerical results presented in [11] are collected in the appendixes. First reports of our study have been published in [14,15].
2 Induced pure gauge theory for gauge group SU(N c ) We start by reviewing the results from [11] which are relevant for this paper. The weight factor of the alternative discretization of Yang-Mills theory can be written in the form where we have adopted the notation used in [11]. Here, 0 ≤ α < 1 is the lattice coupling, the index p labels (unoriented) plaquettes, U p is the product of link variables around the plaquette p, and N b is the number of bosonic fields in the bosonized version of the weight factor. The weight factor (2.1) is a reformulation of the original weight factor introduced in [10] so that one obtains a real action after bosonization (cf. section 2.2 in [11]). Note that in this formulation of the theory N b is just a (positive) number and thus can take any non-integer value, while the bosonization is only possible for integer N b . The weight factor is designed in such a way that, for a given value of N b above the bound in eq. (2.2), the continuum limit is obtained when α → 1.
In [11] we have shown for gauge group SU(N c ) that the theory associated with the weight factor (2.1), which we denote as "induced pure gauge theory" (IPG), approaches a continuum limit for α → 1 as long as Furthermore, in two dimensions the theory in the continuum limit is equivalent to continuum Yang-Mills (YM) theory if For d > 2 the equivalence is a conjecture based on a universality argument made in [10]. As shown in [11], another way of approaching the continuum limit is to send N b → ∞ while keeping α constant. In particular, the lattice theory becomes equivalent to pure gauge theory with the Wilson plaquette action [16] S

PT bound
No continuum limit divergent Figure 1. The conjectured phase diagram of SU(N c ) induced pure gauge theory in the parameter space of the number of bosons N b and the coupling α. The darker shading in red indicates the approach to the continuum theory, indicated by the red lines. The dashed black line in the lower right corner indicates the bound α = 1/3 below which perturbation theory for N b → ∞ is valid. The blue arrows indicate the interesting region for simulations. at lattice coupling β when sending α → 0 and N b → ∞ while keeping β ∼ N b α fixed. We will refer to the lattice theory with action (2.4) as "Wilson pure gauge theory" (WPG). The phase diagram of IPG theory is shown in figure 1 in the parameter space of N b and α.
In standard discretizations of Yang-Mills theory one can investigate the nature of the continuum limit by using lattice perturbation theory, making use of the fact that the bare lattice coupling g goes to zero in the continuum limit. As noted in [11, sec. 4.1], such an expansion is not possible in IPG theory around α = 1 due to the absence of a Gaussian saddle point. An alternative is to expand the theory around the saddle point for fixed α < 1/3 and N b → ∞ and to analytically continue the results to the region where α → 1. The associated bound for perturbation theory is shown as the black dashed line in figure 1. Using this expansion one finds that a suitable definition of the coupling in IPG theory in the region α → 1 is given by where d 0 is a perturbative coefficient. The relation between the couplings g W and g I in WPG and IPG, respectively, is then given by with another perturbative coefficient d 1 . The formulas for the coefficients are given in [11, eqs. (4.12) and (4.13)]. We remark that the weight factor (2.1) respects center symmetry, similar to standard Yang-Mills theory and its discretizations. This is due to the fact that the weight factor is formulated in terms of the plaquette U p , which itself is a center-symmetric object. Center symmetry is of fundamental relevance for the order of the deconfinement transition, which we will study in section 5.

Simulation setup and parameter tuning
We would like to test the conjecture that the continuum limit of IPG reproduces continuum Yang-Mills theory for d > 2 when the continuum limit exists, i.e., if eq. (2.2) is fulfilled. To this end, we compare results obtained from IPG and WPG while taking the continuum limit. The first step along the way is to match the bare parameters such that the simulations in IPG and WPG are done at similar lattice spacings a. In this section we will discuss the matching for the test case of three-dimensional SU(2) gauge theory. This is the computationally cheapest non-trivial case to test the conjecture. From the theoretical point of view there is nothing special about this particular case, so that the results obtained here can be expected to be relevant also for other values of d > 2 and N c > 2.

Simulation setup
We consider pure SU(2) gauge theory discretized on a hypercubic lattice, for which the expectation value of an observable O is given by Here, Z normalizes 1 to unity, and ω[U ] is the weight factor, which is given in (2.1) for IPG and by for WPG with the Wilson plaquette action (2.4). In IPG a simulation point is characterized by two parameters, N b and α, while WPG has only one parameter, β.
For IPG, we are interested in the continuum limit for a fixed small value of N b , i.e., we approach the continuum limit along the blue arrows in figure 1. In particular, we will perform the tests for N b = 1 and 2. These two cases are convenient to test the bounds in eqs. (2.2) and (2.3) because both cases satisfy the first bound, which guarantees the existence of a continuum limit, while N b = 1 violates the second bound, which guarantees equivalence with YM theory in two dimensions. Following the arguments of section 2, we do not expect the latter bound to be relevant for d > 2.
While WPG can be simulated efficiently with the standard combination of heatbath [17] and overrelaxation updates [18], such algorithms are not available for IPG. We thus use a simple Metropolis algorithm [19] for the simulations, discussed in detail in appendix A.1.

Scale setting and parameter matching
We set the scale using the Sommer parameter r 0 [12]. This definition of the lattice scale relies on the physical properties of the static qq potential, which can be computed using Polyakov-loop correlation functions, for instance. The details of the extraction of the potential V (R) and the Sommer scale r 0 are discussed in appendix A.2. Since the change from WPG to IPG amounts to a change of the gauge action only, the operators relevant for the measurement of the potential, together with their spectral representation (cf. appendix A.2) remain unchanged.
We define equivalent lattice spacings by equivalent values of r 0 /a. This means that, for fixed value of N b , we match the bare parameters β and α so that r 0 /a WPG (β) = r 0 /a IPG (α). When tuned in such a way, observables are expected to differ by lattice artifacts only, as long as we are close enough to the continuum. Note that the resulting functional dependence β(α) is not unique. A similar matching could have been obtained using another observable, such as the string tension σ, for instance. The resulting matching relations will then differ by lattice artifacts.
The strategy for the matching is the following: We start by fitting the data for r 0 obtained from simulations of WPG to the expression 1 r 0 (β) =r 0 +r 1 β +r 2 β 2 . (3.3) Performing a simulation of IPG at a fixed value of α and computing the corresponding value of r 0 then gives us, via inversion of (3.3), the value of β to which this particular α should be matched. This procedure results in a number of pairs (α, β) for a given value of N b , which we fit to a suitable parameterization β(α). The only piece of information we include in this parameterization is the fact that β → ∞ should correspond to α → 1. We write the parameterization in the form of an asymptotic series, Perturbation theory supports the validity of such an asymptotic expansion around α = 1 and suggests n = 1, see eqs. (2.5) and (2.6). However, away from perturbation theory, it is not clear that the asymptotic expansion still provides a good parameterization of β(α). This has to be clarified by comparison with the numerical data, and we shall see below that n = 1 indeed provides the best description, in agreement with perturbation theory.

Numerical results for the matching
In order to compare the continuum approach of the two theories we have performed simulations at four β-values for WPG for which high-precision results for the quark-antiquark potential are available [20][21][22][23][24][25][26] Table 1. Simulation parameters and results for the measurements of r 0 in pure SU(2) gauge theory with Wilson action (WPG) and induced action (IPG) for N b = 1 and 2. Here, R gives the range of qq separations used in the analysis of Polyakov-loop correlation functions, t s is the temporal extent of the Lüscher-Weisz sublattices, n t is the number of sublattice updates, ∆ sw is the number of sweeps separating two sublattice measurements, N sw is the number of sweeps between two measurements, and is the size of the ball for the link proposal. For more details on the algorithms, e.g., the choice of ∆ sw and N sw , see appendix A. simulation parameters are given in table 1. All statistical error bars quoted in the following are Jackknife errors obtained with 100 bins. We have checked explicitly that the error bars and estimates for secondary quantities do not change significantly if we vary the binsize.
The methodology for the extraction of r 0 introduced in appendix A.2 relies on the assumption that IPG is a confining gauge theory for α < 1. This is not guaranteed, but the existence of a minimal coupling after which the theory is confining in the approach to the continuum limit is a necessary criterion for the approach to continuum Yang-Mills theory. The simulations have shown that this is the case for all couplings of table 1 so that we can extract r 0 and use it for scale setting. The results for r 0 are also listed in table 1. The first error is statistical, and the second error is the uncertainty associated with the interpolation. Note that we have kept the volume at L/r 0 7 in order to ensure that finite size effects for r 0 are small. That this is indeed the case can be seen from a comparison of the data presented in table 1 and the results for r 0 given in [24,25], where L/r 0 10. The different results are in very good agreement within the statistical accuracy.
We start the matching of the two theories by fitting the WPG data to the form given in eq. (3.3). The results are given in table 2. Note that we have added statistical and systematic uncertainties in quadrature when fitting the data, which likely overestimates the true uncertainties and can thus explain the rather small value of χ 2 /dof in the fit. Using these results as a definition for the behavior of r 0 with β, we then try to find a matching between β and α which leads to identical values of r 0 in the two theories. The guideline for the functional form of β(α) is eq. (3.4). First, it is necessary to determine the leading-order of the divergence of β in the limit α → 1. To this end the data for r 0 obtained in IPG are parameterized by eq. (3.3) and the parameters from 0.623( 4) -1.78(11) 3.6(7) Table 2. Results for the parameters of the matching fits (see text). If χ 2 /dof is not given, the parameterizations contain as many parameters as data points.
replaced by where b, b 0 , and n are free parameters. Note that this is not a fit, since there are as many free parameters as data points. The results for the parameters are listed in table 2. The important information is that in both cases, N b = 1 and N b = 2, we have n ≈ 1 to good accuracy, implying that we can expect the divergence to be a simple pole. This suggests that a suitable parameterization of β(α) is given by We test this parameterization by comparison to the data, using either all three terms (in which case there is no fit) or setting b 1 = 0 (in which case we have a fit with one degree of freedom). The resulting parameters are also listed in table 2. We see that the parameters with b 1 = 0 and b 1 = 0 are in good agreement for N b = 2 while there are some deviations for N b = 1. It appears that this is due to the large correction of the term including b 1 , which also leads to a large value of χ 2 /dof for the fit where b 1 = 0 with N b = 1. Even though a χ 2 /dof of 2.2 is not satisfactory, the fit for the N b = 2 case, in contrast, works reasonably well. Since in both cases the correction term including b 1 is not negligible (for N b = 2 signaled by χ 2 /dof > 2) we take the parameterization with b 1 = 0 for the matching between β and α, This relation will be updated in the next section, see eq. (4.2) below. Figure 2 shows the data for r 0 vs β, including also the data from the induced theory with β obtained from (3.7). The results from this type of matching can also be compared to the results from perturbation theory obtained in [11]. The comparison between the leading-order coefficient b −1 to the perturbative coefficient has already been done in section 4.3 of [11] and shows an excellent agreement between numerical data and the perturbative prediction, surprisingly even for N b = 1. The details of the comparison and the results for larger values of N b > 2 are summarized in appendix B.

The static qq potential
After matching the bare couplings of WPG and IPG we are now in a position to compare the results for other observables. In this section we will focus on the static qq potential, which is not only important to show that the theory is confining but is also related to an effective string theory for the QCD flux tube (for a recent review see [27]). The latter contains non-universal parameters, which can be used to distinguish between different microscopic theories.

Effective string theory and analysis strategy
A possible way to analyze the static qq potential is a comparison to the predictions from the associated effective low-energy theory, namely the effective string theory (EST) for the QCD flux tube, which governs the potential at intermediate and large qq separations R. The interested reader is referred to the reviews [27,28] for further reading about the foundations of the EST and a thorough list of references.
The main result from the EST that we will use in the analysis is the prediction for the R dependence of the static qq potential for large R [29][30][31][32][33], where σ is the string tension and ξ = 6 or ξ = 7, depending on whether the next term in the large-R expansion originates from another boundary term or a bulk term. The first term on the right-hand side is the light-cone (LC) spectrum [29], which is expected to appear due to the integrability of the leading-order S-matrix in the analysis using the thermodynamic Bethe ansatz [32,33]. The appearance of the full square-root formula is also in good agreement with numerical results for the potential (see [27] for a compilation of results). The parameterb 2 in eq. (4.1) is the leading-order boundary coefficient [30,31], which has been found to be non-universal [25,26]. In the spectrum, possible corrections to the standard EST energy levels, such as the rigidity term, first proposed by Polyakov [34], and corrections due to massive modes [32,33], have been left out. Our basic strategy is to try to reproduce, and to compare to, the high-precision results from [25], performing the same analysis steps. Since our aim is to perform a like-by-like comparison, rather than validating the string picture, we focus on the basic EST analysis, i.e., sections 3 and 4 in [25].

Simulation points and scale setting
For the comparison we have performed simulations at bare parameters α and volumes that are matched to the parameters in [25] using the matching relations (3.7). The new set of parameters is shown in table 3.
For the purpose of scale setting and to check the matching of eq. (3.7) we have computed the Sommer parameter using the same strategy as before. The results are given in table 4. We can use these results to update the scaling relation (3.7). For the parameterization with b 1 = 0 we obtain For a first comparison of observables in the approach to the continuum we can look at the expectation value of the plaquette. Even though trivial in the continuum limit, its dependence on the lattice spacing, i.e., the bare coupling, is non-trivial and can be computed in lattice perturbation theory. The three-loop result for SU (2)   given by [35] The numerical results for the plaquette are listed in table 4 and shown in figure 3. To make the small differences visible we do not display the raw data but rather the difference between the plaquette expectation values and the perturbative result. The plot shows that the WPG results are already very close to the perturbative result, i.e., lattice artifacts with respect to 3-loop perturbation theory for Wilson's gauge action are small. From the plot one might be led to the conclusion that lattice artifacts for IPG are larger. However, the coefficients of the perturbative result which we subtracted have been obtained from Wilson's plaquette action and are expected to be different for IPG in general, and for different values of N b in particular. Indeed, corrections to eq. (4.3) in WPG theory start at order g 8 0 , while for IPG theory corrections start at lower orders. Another option to set the scale is via the string tension σ, which governs the linear rise with R of the potential for R → ∞. As in [25] we extract σ with two different methods: (i) We fit the force to the form  motivated by the expansion of the EST potential to next-to-leading-order in 1/R.
(ii) We fit the potential to corresponding to the (full) leading-order EST prediction with an additive normalization constant V 0 .
This particular combination of ansätze is especially useful since corrections with respect to the full EST prediction, eq. (4.1), appear at different orders in the 1/R expansion. Consequently, we can determine the region where higher-order terms in the EST are negligible by comparing the results for σ from methods (i) and (ii). The basic strategy is to investigate the dependence of σ on the minimal value of R included in the fit, R min . In particular, in the region where the results from the two methods agree within errors and show a plateau, higher-order terms are expected to be negligible, and we can use any of the results for σ.
In figure 4 we show the results for the extraction of the string tension in IPG for N b = 1 (left) and N b = 2 (right). We see that for N b = 1 and α = 0.931 and 0.946 the results from method (i) already start to leave the plateau region when the results from method (ii) reach the plateau. Nevertheless, the plateau values agree within uncertainties, as the For comparison we also show the results for SU(2) WPG from [25] (red circles) and the associated continuum-extrapolated value (black circle).
results from method (i) agree within errors at the point where the results from method (ii) reach the plateau region. As the final results for the string tension, we will use the result from method (ii) obtained with the value of R min where the results of the two methods become consistent. In the cases without a common plateau we use the value of R min where the result from method (ii) starts to agree with the plateau from method (i). The final results are indicated by the red bands in figure 4. The results for σ are listed in table 4. The quoted uncertainties for σ include the systematic uncertainties for r 0 , which have been conservatively added in quadrature to the statistical uncertainties.
To compare the results for σ with the results for WPG with gauge group SU(2) from [25] we show the two sets of results in figure 5 for N b = 1 (left) and N b = 2 (right). For N b = 2 we observe good agreement between the results from IPG and WPG, while some slight differences are visible for N b = 1. The latter could either be fluctuations, remnants of uncontrolled systematic uncertainties connected to the extraction of σ, or simply due to the different lattice artifacts in the two theories. Eventually, we expect the results to agree in the continuum limit. To test this, we have performed a continuum extrapolation of the As in [25] we perform two fits, either with b σ,2 = 0 or with b σ,2 = 0. The resulting extrapolations are also shown in figure 5 together with the extrapolation for WPG from [25] (black circle). We see that for N b = 2 the continuum extrapolations are all in very good agreement with the WPG result, even though less precise, which can be expected since only three points are available for the extrapolation. For N b = 1 the extrapolation linear in a 2 overshoots the WPG result. However, the data also indicate the importance of higher-order terms. Including the a 4 term leads to an extrapolation which is fully consistent with the WPG result, albeit with large uncertainties. For comparisons we will use the results from the continuum extrapolation with b σ,2 = 0. This extrapolation has larger errors, but the central value agrees well with the continuum-extrapolated WPG result in both cases.

Results for the static potential
We will now look at the results for the static potential itself. The results for the potential, rescaled according to are shown in figure 6. As in [25] the results for the individual couplings have been rescaled using the string tension for this value of α, while the solid line is the rescaled version of the leading-order EST prediction, the LC spectrum, with the continuum limit of the string tension.
Results for the static potential in WPG were already presented in figure 3 of Ref. [25]. We do not reproduce that figure here. The comparison with figure 6 shows one of the main problems we are facing, namely the reduced accuracy of the present study, which is mainly due to the less efficient algorithm for IPG compared to the heat-bath algorithm used in the WPG simulations. This leads to a reduction in the range of R and reduces the precision in the extraction of σ. Concerning the results for the potential itself, the agreement with the WPG results in [25, fig. 3] is eminent for N b = 2 and also visible for N b = 1, although the lattice with α = 0.903 appears to be outside the scaling region. 3 As in WPG, the corrections to the LC potential in IPG are positive and tend to become stronger when approaching the continuum limit.
The next step in the analysis is to check whether the leading-order correction to the square-root formula in eq. (4.1) is indeed of order R −4 . To this end we fit the data to the form where η, m, σ and V 0 are fit parameters. If the predictions of the EST are correct, we will obtain m = 4. Otherwise we will find m < 4 if the square-root formula in eq. (4.8) is incorrect, or m > 4 if the corrections start at higher order. We show the results for m vs R min in figure 7. As in SU(2) WPG, cf. [25, fig. 4], we typically observe a plateau around 0.5 R/r 0 1.0. For N b = 1 the uncertainties are typically larger for larger R-values, so that the plateau does not last as long as for N b = 2. The plateau value is typically around m = 3.6. This slight discrepancy with expectations has also been found in WPG [25] and indicates a possible mixing with other correction terms. At finite N c the EST will receive corrections from virtual glueball exchange, for instance.

Extraction of the boundary coefficient
To make the comparison of the subleading properties of the potential more quantitative, we will now extract the boundary coefficientb 2 . As in [25], we include higher-order terms in the fit formula and fit the potential to the form Specifically, we perform the following fits [25]: A We use σ and V 0 from method (ii) as input and fit withb 2 , γ 0 and γ (2) 0 as free parameters and setb 2 = 0. The fits are performed for several values of R min , and we extract the final result from the second smallest R min which provides a χ 2 /dof < 1.5. The quality of the agreement with the data is then indicated by the value of R min (smaller values mean better agreement) in context with the number of higher-order terms included in the fit. Fit C, for instance, should allow for a smaller value of R min compared to fit B since the latter does not contain higher-order correction terms. To check for the systematic uncertainty associated with the fit interval we compare the result to the ones obtained with R min ± a. We list the results of the different fits for N b = 1 and 2 in tables 5 and 6, respectively. Comparing the results to those of [25, tab. 5], in particular for β = 5.0, 7.5 and 10.0, we see that both the possible fit ranges and the resulting parameters are very similar. The only exception is fit A, for which σ and V 0 have been taken over from section 4.2. Since the extraction of σ and V 0 was not as accurate as in [25] it is reasonable that this is the main reason for deviations in this particular fit. We can thus follow the discussion of [25] and conclude that fits B, C and D can be used in the following analysis. Fit E typically requires larger values for R min compared to fits C and D, even though it also includes two higher-order terms. Thus we conclude that the agreement with the data is worse for fit E than for the other fits. In any case, fit E only serves as a check that, given the data, b 2 does not vanish. Compared to [25], where fit E clearly showed less agreement with the data, direct conclusions are more difficult here since the IPG data are less precise than those of WPG.
As in [25], we determine the final results forb 2 on the individual lattice spacings via the weighted average over the results from fits B to D, where the weight is given  Table 5. Results of the fits for the extraction ofb 2 for N b = 1 IPG.   Table 7. Final results forb 2 in WPG [25] and IPG with N b = 1 and 2 for the individual couplings. The first error is the statistical uncertainty, the second the systematic one due to the unknown higher-order correction terms, estimated by computing the maximal deviations of the results from fits B to D, and the third is the systematic one associated with the choice of R min . ¹¼º¼ ¹¼º¼¿ ¹¼º¼¾ ¹¼º¼½ ¼ ¼º¼½ comparison in table 7 as well, we see that the results are similar in magnitude and have similar uncertainties. This is particularly true for N b = 2. For N b = 1 the systematic uncertainties are somewhat larger, but the overall agreement is still good. We show the results forb 2 in comparison to the results of [25] in figure 8. One can clearly see the similar behavior in the approach to the continuum and the good agreement between the results from the different simulations. In particular, the results are significantly different from the ones for SU(3) gauge theory, showing the discriminating power of the results.
Note that fit (1) is a parameterization of the results rather than a fit. In [25] a third fit has been performed, which included only the data for the finest lattice spacings. Such a fit does not make sense here due to the limited number of available couplings. To estimate the propagation of systematic uncertainties we follow the strategy of [25] and perform the fits (1) and (2) for the results from fits B to D and the fits with a minimal R-value of R min ± a individually. The final results have been extracted using a weighted average of the results from fits B to D, once more with the individual uncertainties as weights. As before, the individual systematic uncertainties are computed from the maximal deviations to the different fits. The curves for fit (2) with the main value of R min are shown in figure 9. The continuum results forb 2 from fits (1) and (2) are given in table 8. We use the results from fit (1) to estimate the systematic uncertainty associated with the continuum limit. The final results, which we take from fit (2), are thus given by    (1) and (2) (see text). The first error is the statistical uncertainty, the second the one associated with the unknown higher-order correction terms in the potential and the third the one due to the choice of R min .
Those can be compared to the final continuum estimate for SU(2) WPG from [25], In eqs. (4.11) and (4.12) the first error is purely statistical, the second is the systematic uncertainty due to the unknown higher-order terms in the potential, the third is the one associated with the particular choice for R min and the fourth is the systematic uncertainty due to the continuum extrapolation. We see that the results from eqs. (4.11) and (4.12) agree well within uncertainties. In addition, the sizes of the individual uncertainties are similar, except for the continuum extrapolation, which, however, is expected since in the IPG analysis ensembles at fewer lattice spacings are available. We briefly summarize the findings of this section. We have repeated the analysis of the potential in [25] for IPG and found that in every individual step the results agree extremely well with WPG, for both N b = 1 and 2. The whole analysis indicates that the fine structure of the potential in the continuum is indeed identical in IPG and WPG. Hence we can conclude that, at least for the potential, both theories lead to the same continuum limit (up to the accuracy of this study). This is also reflected in the significant difference to the results for SU(3) WPG, which shows the discriminating power of this comparison.

The finite-temperature phase transition
So far we have compared properties of IPG and WPG for observables at vanishing temperature and found good agreement. We will now show that the agreement prevails for thermodynamic observables. In particular, we will consider the deconfinement transition temperature T c and the ratio of critical exponents γ/ν. The latter can be regarded as a measure of the universality class of the transition (we expect a phase transition of 2nd order). The fundamental lattice observable that can be used to investigate the deconfinement transition is the absolute value of the Polyakov loop, |L| , the order parameter associated with the breaking of center symmetry. In particular, we choose a setup which is similar to that in [36] so that we can directly compare to this study.

Simulation parameters and results for the Polyakov loop
As before, we perform simulations in IPG using N b = 1 and 2. In addition, we also simulate in WPG for a direct comparison of the results. To test the approach to the continuum limit we use two different temporal extents, N t = 4 and 6, for which we vary the temperature T = 1/(aN t ) by varying the lattice spacing a via the lattice couplings α and β. For scale  Table 9. Simulation parameters of the finite-temperature runs in 3d SU(2) IPG and WPG theories. For each simulation point we have performed at least 100,000 measurements and increased the number of measurements to about 400,000 in the vicinity of T c , where the autocorrelations increase due to the approach of a second-order critical point. The results for |L| and χ L vs the temperature in units of the Sommer parameter are shown in figures 10 and 11 for N t = 4 and 6, respectively. We always show the two extremal cases of smallest and largest available volume. The plots show the remarkable agreement between the results of WPG and IPG with N b = 1 and 2. This already indicates the similarity between the corresponding phase transitions. In particular, the volume scaling is equivalent in the different cases so that the universality classes can be expected to be equivalent as well. In the following we will investigate this expectation more quantitatively.

The transition temperature
We define the critical temperature T c through the peak of the Polyakov-loop susceptibility. We determine it by fitting a Gaussian to the points in the vicinity of the peak. This definition assumes a Gaussian form of the susceptibility peak, which will, generically, not be the case. To account for the associated systematic uncertainty we use the distance between the two points surrounding the maximum as a conservative error estimate for T c .    Figure 11. Results for the absolute value of the Polyakov loop |L| (top) and its susceptibility χ L (bottom) at N t = 6 with volumes 32 2 (left) and 96 2 (right). The legend is the same for all plots. We have also checked that the systematic uncertainty associated with the number of points (symmetrically distributed around the maximum) entering the fit is much smaller than this estimate for the systematic uncertainty of the results. The results for T c in units of r 0 are listed in table 10 for the different values of N t and the different volumes. We have not listed the other fit parameters, such as the width of the Gaussian, since they are not relevant for the following analysis. As expected from the results of the previous section, the results for WPG and IPG with N b = 1 and N b = 2 agree well within uncertainties. For N t = 6 the uncertainties for WPG are larger due to the larger separation of simulation points. The results can be compared to the results of [36] in the V → ∞ limit. To convert to units in terms of r 0 we use the results for β c (∞) (given in [36, table 1]) and convert them to T c r 0 (∞) using the interpolation for r 0 (β) from section 3.3. The results are also listed in table 10, where the uncertainties include the uncertainties for β c (∞) and for r 0 /a. To be able to directly compare the results we need to perform a V → ∞ extrapolation of our data for T c . For N s → ∞, T c is expected to obey the scaling behavior where ν is the associated critical exponent. The Potts model associated with the SU(2) transition yields ν = 1 [37], i.e., a linear behavior of T c with 1/N s , which has been found to be in good agreement with the numerical data for T c [38]. Figure 12 is clearly in agreement with eq. (5.3), even though the uncertainties are too large to draw any final conclusions. Assuming that the three largest volumes are already in the scaling region, we perform a linear extrapolation to N s → ∞. The results are also listed in table 10. We show the results vs the inverse box length 1/N s (the spatial volume is given as N 2 s in table 10) in figure 12. The plot indicates good agreement between IPG, with both N b = 1 and 2, and WPG at finite and infinite volume. Note that our uncertainties are very conservative estimates for the systematic uncertainties associated with the fits for the extraction of T c and thus could be overestimated. These large uncertainties also propagate to T c (∞).  Figure 12. Results for the deconfinement transition temperatures T c r 0 vs the inverse box length 1/N s for N t = 4 (left) and N t = 6 (right). The data have been slightly displaced to improve the visibility of the different sets of points. The black data point at 1/N s = 0 is the result in the infinite-volume limit from the Bielefeld group [36], while the magenta and blue open symbols are the results for T c (∞) from IPG with N b = 1 and 2, respectively.

Order and universality class of the transition
Owing to the scaling of T c discussed above, we expect the transitions in IPG and WPG to be in the same universality class. Moreover, as already discussed in section 2, center symmetry is a good symmetry for both actions so that the deconfinement transition is accompanied by center-symmetry breaking. Another strong indication for the transitions being in the same universality class comes from the similar volume scaling of the susceptibility peaks of the Polyakov loop (cf. figures 10 and 11). We will now test whether these expectations are correct and both transitions are indeed in the 2d Ising universality class, which is known to be the case at least for 3d SU(2) WPG theory [36].
There are several types of analyses one can employ to determine the critical exponents which distinguish the different universality classes. Here we follow the strategy of [36] and use the χ 2 -method [39]. The starting point is the finite-size scaling formula for the susceptibility of the Polyakov loop expanded around the critical point, Here, c 0,1,2 are unknown coefficients, γ i is an unknown exponent, is the reduced temperature with respect to the critical temperature in the thermodynamic limit T c (∞), N s is the spatial extent of the lattice, and γ/ν is the desired ratio of critical exponents. Exactly at T c (∞), eq. (5.4) reduces to T r 0 Figure 13. Results for γ/ν from the fits discussed in the text, for 3d SU(2) gauge theory at N t = 4. The plot on the right displays the details of the T c (∞) region. The black dashed line is the prediction for γ/ν from a reduced model (RM) (cf. [36]) and the black data point, together with the orange band, the result from the Bielefeld group [36].
At large N s the second term (proportional to c 3 ) is a correction and can be neglected. We thus arrive at the scaling relation This scaling law can be tested by fitting the data at a given coupling to (5.7). Since at T = T c (∞) there are corrections of the form given in eq. (5.4), we will get the best fit when the coupling is equal (or very close) to the critical coupling. The critical temperature T c (∞) can thus be extracted by looking at χ 2 /dof, in principle.
The main problem of this analysis is the finite resolution of simulation points in the region around T c and the different accuracy of the results for the susceptibility. For WPG theory this problem can be overcome by the use of the multi-histogram method [40], which can be used for a well-controlled interpolation between simulation points. In addition, the method leads to enhanced and balanced statistics for all simulation points. In IPG theory this method cannot be applied since α, unlike β, does not appear as a simple prefactor in front of an observable (the average plaquette for WPG) in the action. For IPG this leads to the problem that χ 2 /dof fluctuates strongly and cannot be used as a conclusive indicator for T c as in [36]. Instead, we will compare to the results for γ/ν obtained for the individual temperatures by fitting to eq. (5.7) and making use of the results for T c (∞) from table 10. In this way we obtain a number of possible results for γ/ν in the region of T c (∞). To be conservative, we use the full spread of results as the uncertainty interval which encloses the final result for γ/ν and define the central value of our final result to be the midpoint of this interval. These results are also listed in table 10. Unfortunately, the uncertainties in IPG are too large to draw definite conclusions from the comparison. Nonetheless, the results are in agreement with those of WPG within errors.
The results of the analysis are shown in figures 13 and 14 for N t = 4 and 6, respectively. The plot on the right in the figures shows the details in the region around T c r 0 (∞). In the plots we do not show the results for γ/ν in IPG since the uncertainties are rather large. However, when we assume that a hypothetical high-precision result for T c (∞) in IPG would be similar to the result in WPG, the plot indicates that we would get a similar result for γ/ν, too. All in all we have compelling evidence that the transitions in WPG and IPG theory are in the same universality class. Moreover, cutoff effects are similar, so that we can expect this statement to be true for the continuum limit as well.

Conclusions
In this paper we have tested the conjecture of the equivalence of the continuum limit of induced pure gauge theory, formulated originally by Budczies and Zirnbauer in [10], and pure gauge theory with Wilson's gauge action [16]. To this end we have performed simulations with both discretizations in three dimensions with gauge group SU(2) at matched couplings to achieve similar lattice spacings. The matching via the Sommer scale was discussed in detail in section 3. It is found to be in good agreement with perturbation theory, both concerning its functional form and the numerical results for the matching coefficients (see also [11, sec. 4.3] and appendix B).
Using the matching, we have performed simulations at similar lattice spacings for a high-precision comparison of observables. In particular, we have looked at the leading and subleading properties of the static qq potential at intermediate and long distances in section 4. The leading properties are characterized by the string tension and the subleading ones by the non-universal boundary coefficientb 2 . Both observables show excellent agreement between WPG and IPG in the approach to the continuum for N b = 1 and N b = 2 and are similar already at finite lattice spacing. This indicates that the continuum potential is identical, at least at the current level of precision.
The thermodynamic properties of IPG theory have been investigated in section 5 for two temporal extents, N t = 4 and 6. Once more the agreement between the two theories for the behavior of the Polyakov loop and its susceptibility is remarkable for all volumes and over the full range of temperatures. Accordingly, the ratio of critical exponents γ/ν agrees very well with the known result for WPG from the Bielefeld group [36], indicating that the transition in both theories is in the same universality class. Furthermore, the transition temperature in the thermodynamic limit is also in agreement, which is another important crosscheck since T c is a non-universal quantity.
One of the problems of the Kazakov-Migdal model [8], an earlier model supposed to induce Yang-Mills theory, is the existence of a local center symmetry which constrains all Wilson lines to vanish. The present model for induced Yang-Mills theory does not have this problem. This is also reflected in the results for the Polyakov loop obtained in section 5. The results for the Polyakov loop fall in the two sectors singled out by the center symmetry of SU(2), i.e., they are located around |L| and −|L| on the real axis. This shows that the model also recovers the correct symmetry, whose breaking is associated with the transition in Yang-Mills theory. As mentioned already in section 2, this is also expected from the symmetries of the weight factor, eq. (2.1).
Concerning numerical efficiency, our simulations of IPG theory are up to a factor of 100 slower than the associated simulations in WPG theory, depending on the bare parameters in the simulations. We would like to stress that this inefficiency is an attribute of the choice of the simulation algorithm alone, as explained in detail in appendix A. We are confident that one can find an algorithm similar to the heat-bath algorithm for WPG [17], which should then lead to a similar performance. In particular, if one is interested in simulations including fermionic degrees of freedom, the Hybrid Monte Carlo (HMC) algorithm [41] is supposedly the algorithm of choice. We have tested the HMC algorithm for IPG theory in its bosonized version and found it to perform almost as well as the HMC for WPG theory in both SU(2) gauge theory for d = 3 and SU(3) in four dimensions. A general problem working in the bosonized version is the increase in autocorrelations, which can be enhanced up to an order of magnitude. However, this is compensated to some extent by the speed-up of the individual HMC steps. It is interesting to note that the HMC algorithm for IPG theory can be set up in such a way that the auxiliary boson fields are drawn from the exact distribution for the starting configuration. In this case no communication is needed to evolve the link variables in pure gauge theory, which results in a very efficient parallelization.
All in all, we find that the results from the two theories agree very well already away from the continuum limit and that the agreement improves as the continuum limit is approached. This is true already for N b = 1, which leads us to conclude that the bound in eq. (2.3) can indeed be relaxed for d > 2. Our results thus support the universality argument given in [10]. Therefore the modified BZ method, combined with a suitable strong-coupling approach, can be used to reformulate lattice gauge theories in a number of different ways. We leave these reformulations for future work.

A.1 Update algorithm
The weight for a link U µ (x) in the weight factor (2.1) in IPG theory is local so that we can use a local version of the Metropolis algorithm [19] for the simulations. More precisely, we propose a new link U µ (x) and accept it with probability P = min 1, Here, the product over p is taken over all plaquettes that include the link U µ (x), U p denotes the plaquette including the old link U µ (x), and U p the plaquette including the new link U µ (x). The crucial part affecting the efficiency of the algorithm is to find new links U µ (x) that lead to a large acceptance rate in the step defined by eq. (A.1). Since in our pilot study of IPG theory efficiency of the algorithm was not our primary concern, we simply took the new links U µ (x) to be random SU(2) matrices in an -surrounding of the old links U µ (x). In practice, we have generated a random element X = a x a T a ∈ su(2) (the Lie algebra of SU (2)), where the T a are the generators of SU(2) and the coefficients x a are taken from the interval [− , ], and constructed the new link via U µ (x) = exp(X)U µ (x).
To make sure that in each step a sufficient number of links is updated, we have tuned so that the overall acceptance rate is around 80%. To further decorrelate two measurements we have separated them by N sw sweeps, where N sw is chosen to be much larger than the integrated autocorrelation time in units of lattice sweeps.

A.2 Extraction of the static qq potential
To compare the two theories we mostly use quantities that are related to the potential between a static quark and antiquark. The cleanest way, in terms of excited-state contaminations, to extract the potential in numerical simulations is given by measuring the correlation function of two Polyakov loops, L( x), defined by (cf. eq. (5.1)) where n 0 denotes the temporal coordinate, N t is the number of lattice points in the temporal direction, and x is the vector including the spatial coordinates. The spectral representation of a correlation function of two Polaykov loops separated by a distance R is given by for R L/2, where T = aN t and L = aN s are the temporal and spatial extents of the lattice, b i denotes the overlap between the operator and the energy eigenstate, and E n (R) is the n-th energy level. Here the energies are ordered in ascending order, i.e., E 0 < E 1 < E 2 < . . . . In the limit T → ∞ the right-hand side of eq. (A.3) is dominated by the ground state with E 0 (R) = V (R). Excited states are suppressed exponentially with exp{−[E i (R) − E 0 (R)]T }. This means that excited states can be neglected for large values of T . In this case we can extract V (R) via The correlation functions in eq. (A.3) suffer from a well-known exponential decay of the signal-to-noise ratio with the area enclosed by the two loops. This renders the extraction of the potential difficult for large R. For simulations in WPG this problem can be overcome by the use of a multilevel algorithm introduced by Lüscher and Weisz [42]. The same algorithm can also be applied in IPG since the locality properties of the action are similar to those of the plaquette action. The details are discussed in appendix A.3.
A suitable observable to set the scale with the static potential is the Sommer parameter r 0 [12], which is defined implicitly by To extract r 0 /a we use the following four methods (see also [25]): (a) a numerical polynomial interpolation of R 2 F (R), (c) a parameterization of the form [12] F for the values of R corresponding to the four nearest neighbors of r 0 (motivated by the EST to LO), (d) the parameterization of (A.7) with f 2 = 0 for the two nearest neighbors of r 0 .
The final estimate for r 0 /a is obtained from method (d), while methods (a) to (c) are used to compute the systematic uncertainty associated with the interpolation.

A.3 Error reduction for large loops
A suitable algorithm to overcome the exponential decay of the signal-to-noise ratio for large loops is the multilevel algorithm proposed by Lüscher and Weisz [42]. The algorithm relies on a key property of the theory, the locality of the configuration weight, which ensures that sublattices, i.e., lattice domains separated by a time slice with fixed spatial links, are independent during local updates. IPG theory also has this form of locality since, as for the Wilson action, the weight for a particular link U µ (x) only depends on the plaquettes including this link. The error-reduction efficiency of the algorithm, however, depends on the properties of the transfer matrix since this is the object for which the uncertainty is decreased in the course of the sublattice updates. If we are close enough to the continuum and both theories indeed approach the same continuum limit, we can expect the transfer matrix to be similar in both theories (given by the continuum result plus lattice artifacts) so that the algorithm should lead to a comparable error reduction also in the case of IPG theory. A way to estimate the amount of error reduction for Polyakov-loop correlation functions, and thereby the optimal number of sublattice updates for loops of a particular size, has been proposed in [43]. Following this proposal we define the norm N (R) of a local two-link operator with link separation R on a particular sublattice as in [43, eq. (4)]. The two-link operator plays the role of a transfer matrix for Polyakov-loop correlators [42]. The decay of N (R) provides an estimate for the residual fluctuations of the transfer matrix. For a large number n t of sublattice updates we expect N (R) to fall off as 1/ √ n t .
To compare the error-reduction efficiency of the algorithm in the IPG and WPG theories we have performed dedicated simulations in the SU(2) theory on a 24 3 lattice at β = 5.0 and the roughly equivalent α = 0.65 in IPG with N b = 2. We have fixed the temporal extent t s of the sublattices to 2 throughout and measured two-link operators with link separations between 2 and 10. In figure 15 we show the results for link separation 10. The "optimal" value of n t is taken to be the point at which the expected asymptotic behavior N (R) ∼ 1/ √ n t sets in for the largest value of R considered. After this point no further exponential error reduction is achieved. The plot indicates that for WPG theory n t = 15000 to 20000 is sufficient. For IPG theory, N (R) × √ n t decreases more slowly if we perform measurements at each sublattice update. The main reason is that the sublattice configurations after one sublattice sweep of local Metropolis updates in IPG theory are more correlated than the configurations after one sublattice sweep of heat-bath updates for WPG theory. We therefore separate the measurements in IPG theory by ∆ sw sublattice sweeps (which saves simulation time since the measurement of the observable is more costly than a sweep of Metropolis updates). We see from figure 15 (left) that the decay of the norm becomes stronger when we increase ∆ sw . For the values shown in figure 15 (left), where = 0.07, ∆ sw = 50 − 100 appears to be a good choice. Moreover, the optimal value of ∆ sw depends on the value of , as shown in figure 15 (right). Since in the final simulations we have used a value of = 0.14, for which the integrated autocorrelation time is a factor of 3 smaller than for = 0.07, the optimal value of ∆ sw can be taken to be around 20 for α = 0.65. A similar tuning can be done for the other α-values in the simulation. Since the optimal values of and ∆ sw at a particular lattice spacing are a property of the algorithm, we can use the same setup for other values of N b as long as we have tuned α so that we simulate at similar lattice spacings. Note that the correctness of the algorithm does not depend on the particular values of and ∆ sw . The only effect of a suboptimal tuning of or ∆ sw is less error reduction.

B Comparison of the matching with perturbation theory
We want to compare the matching between the lattice couplings β and α discussed in section 3.2 to the perturbative results for the matching obtained in [11]. The perturbative result (after analytic continuation to the region where α → 1) as given in [11, eq. (4.15)], reformulated in terms of the lattice couplings β and α, is given by The N b -dependence of d 0 (N b ) and d 1 (N b ) can be computed in perturbation theory in the limit N b → ∞. For a direct comparison to perturbation theory it is convenient to replace the matching function from eq. (3.6) by We are particularly interested in the coefficient d 0 , for which the perturbative prediction is given in [11, eq. (4.78)].
To compare to perturbation theory we have used the parameterization (B.2) for a comparison with the data for N b = 1 and 2 from section 3.2 and performed additional simulations for N b = 3, 4 and 5, for which the simulation parameters are listed in table 11. We observe that this parameterization works well and leads to a similar description of the data as the parameterization used in section 3.2. The results for d 0 (N b )/N b , which have already been shown in [11, Figure 2], are listed in table 12. 4 For the comparison itself we refer to [11].
We can also compare directly to the matching results of section 3.2 by noting that the coefficients b −1 and d 0 are related by  Table 11. Simulation parameters and results for the measurements of r 0 in pure SU(2) IPG for N b = 3, 4 and 5. Here, R gives the range of qq separations used in the analysis of Polyakov-loop correlation functions, t s is the temporal extent of the Lüscher-Weisz sublattices, n t is the number of sublattice updates, ∆ sw is the number of sweeps separating two sublattice measurements, N sw is the number of sweeps between two measurements, and is the size of the ball for the link proposal. For more details on the algorithms, e.g., the choice of ∆ sw and N sw , see appendix A.
Using this relation to convert the results of eq.