{'fraction_non_alphanumeric': 0.04017559548004227, 'fraction_numerical': 0.013283472888383059, 'mean_word_length': 4.41664464993395, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Testing Proposals for the Yang-Mills Vacuum Wavefunctional by Measurement of the Vacuum

We review a method, suggested many years ago, to numerically measure the relative amplitudes of the true Yang-Mills vacuum wavefunctional in a finite set of lattice-regulated field configurations. The technique is applied in 2+1 dimensions to sets of abelian plane wave configurations of varying amplitude and wavelength, and sets of non-abelian constant configurations. The results are compared to the predictions of several proposed versions of the Yang-Mills vacuum wavefunctional that have appeared in the literature. These include (i) a suggestion in temporal gauge due to Greensite and Olejnik; (ii) the"new variables"wavefunction put forward by Karabali, Kim, and Nair; (iii) a hybrid proposal combining features of the temporal gauge and new variables wavefunctionals; and (iv) Coulomb gauge wavefunctionals developed by Reinhardt and co-workers, and by Szczepaniak and co-workers. We find that wavefunctionals which simplify to a"dimensional reduction"form at large scales, i.e. which have the form of a probability distribution for two-dimensional lattice gauge theory, when evaluated on long-wavelength configurations, have the optimal agreement with the data.


I. INTRODUCTION
Most of the key non-perturbative properties of non-abelian gauge theories, such as the static quark potential, the chiral condensate, and the topological charge density, are actually properties of the vacuum of the quantized theory. In the Hamiltonian formulation, the vacuum state is the ground state wavefunctional of the Hamiltonian operator, and all of the excited states of the theory, i.e. the mesons, baryons, and, in a pure gauge theory, the glueballs, are simply small excitations on top of that underlying ground state. For this reason, knowledge of the Hamiltonian ground state wavefunctional could be essential in understanding the infrared properties of a non-abelian gauge theory.
Proposals for the ground state of pure Yang-Mills theory go back over thirty years [1,2]. However, with only a few exceptions [3][4][5][6][7], very little work was done in this area after those initial efforts. In recent years, however, there has been a modest revival of interest in this area, and a number of plausible suggestions for the vacuum state have been advanced. These proposals will be described, along with their motivations, in the next section. Briefly, there are suggestions which have been put forward in temporal gauge [8], in Coulomb gauge [9][10][11][12][13] and, in 2+1 dimensions, in terms of gauge-invariant "new variables" [14]. Since these suggestions differ in various ways, it would be interesting to know which (if any) is the true vacuum state, or at least a reasonable approximation to the true vacuum state.
In this article we will apply an old method [15][16][17] for measuring, via lattice Monte Carlo simulations, the relative magnitudes of the true Yang-Mills wavefunctional in any given set of lattice gauge field configurations. The evaluations will be carried out for two types of lattice configurations: non-abelian constant gauge fields of varying amplitudes, which are constant in space but noncommutative [U i ,U j ] = 0, and abelian plane waves of various amplitudes and wavelengths, which are abelian in the sense that [U i ,U j ] = 0. The results are compared to the corresponding values obtained in each of the proposed vacuum wavefunctionals. The method can be applied in any number of space-time dimensions, but here we will work exclusively in 2+1 dimensions, since the new variables proposal [14] is formulated only in that case.
In section II below we will introduce and motivate each of the wavefunctionals to be tested. Section III reviews the method for measuring the true vacuum wavefunctional, and section IV compares the results obtained by this method with the predictions of each of the proposed ground states. Our conclusions are in section V, and some numerical details are found in the appendix.

II. VACUUM STATE PROPOSALS
The Yang-Mills Hamiltonian operator takes on its simplest form in temporal gauge, namely in the continuum theory in D + 1 dimensions, and on the lattice, where the sums are over links l and spatial plaquettes p, respectively. Physical states in temporal gauge must obey the Gauss law constraint D ab k E b k Ψ = 0, or more explicitly which implies that physical states must be invariant under infinitesimal gauge transformations. The Gauss law constraint in temporal gauge is a mixed blessing in the search for an approximate ground state. On the one hand, gauge invariance can be seen as an aid in selecting a good ansatz for the vacuum state. On the other hand, by severely limiting the choice, certain states which are perfectly acceptable in Coulomb gauge, and which may be much more amenable to an analytical treatment, must be discarded in temporal gauge. A very important relation, for our purposes, is the equality of the vacuum wavefunctionals in temporal and Coulomb gauge (see, e.g., ref. [18]), when evaluated on gauge fields satisfying the Coulomb gauge condition ∇ · A = 0, and which also lie in the first Gribov region. Since our numerical method, to be described in the next section, will generate the relative amplitudes of vacuum wavefunctionals in temporal gauge, in any finite set of gauge field configurations, we will be able to check proposals in Coulomb gauge by ensuring that the given set satisfies the Coulomb gauge condition, and lies within the first Gribov horizon. The ground state wavefunctional is known in two limits: the free-field g 2 = 0 limit, and also at strong lattice couplings g 2 ≫ 1. In the free-field limit, in either Coulomb or temporal gauge, while in the strong-coupling limit, in SU(N) gauge theory, it has been shown that [19] Ψ 0 [U] = N exp N g 4 (N − 1) ∑ P TrU(P) + c.c. , (6) to leading order in 1/g 2 . It was suggested long ago in ref. [1], by one of the present authors, that the Yang-Mills vacuum wavefunctional in 3+1 dimensions might have the form when evaluated on sufficiently long-wavelength, slowly varying field configurations. This wavefunctional has the property of dimensional reduction: If we write then R [A] has the form of the Euclidean Yang-Mills action in one lower dimension (three dimensions, in this case). It is clear that the strong-coupling vacuum state (6) does, in fact, have this property. The dimensional reduction vacuum (7) in 3+1 dimensions is confining, i.e.
if and only if Yang-Mills theory in three Euclidean dimensions has that property, where U(C) is a Wilson loop holonomy around the planar, spacelike loop C. Of course we have good reasons to believe that Yang-Mills theory is confining in three Euclidean dimensions. It was noted by Halpern [2] that a dimensional-reduction vacuum state in 2+1 dimensions must be confining, since Yang-Mills theory in two Euclidean dimensions is known to confine. Dimensional reduction was also suggested somewhat later, on rather different grounds, by Ambjorn, Olesen, and Peterson [20,21]. These authors were the first to make the connection between dimensional reduction and the property that has come to be known [22] as Casimir scaling. Strong evidence for Casimir scaling at intermediate distance scales was found in [23].
On the other hand, the dimensional reduction wavefunctional cannot be correct as it stands, because the short-distance structure is completely wrong. For example, equal-time twopoint correlators in D + 1 dimensions, at short distances, cannot be identical to short-distance two-point correlators in D Euclidean dimensions; the singularity structure in the approach to zero separation would be wrong. In general one would expect that the vacuum state evaluated on short wavelength configurations would agree with the perturbative ground state, whose zeroth order approximation is given by (5).
There are other reasons, apart from short-distance singularity structure, that dimensional reduction cannot be exact even for infrared physics. Dimensional reduction from 2+1 to two Euclidean dimensions would imply a non-vanishing string tension, and perfect Casimir scaling, for any color group representation. This cannot be right in 2+1 dimensions, because of color screening. 1 As argued in ref. [8], it is quite plausible that color screening is achieved by small corrections to the dimensional reduction form.
Another argument against exact dimensional reduction from 3+1 to three Euclidean dimensions was raised in refs. [26,27], which pointed out that this reduction would imply a match between the equal-time Coulomb gauge gluon propagator in 3 + 1 dimensions, and the Landau gauge propagator in D = 3 Euclidean dimensions. It was shown in the same references that these propagators actually do agree quite well in a range of low and intermediate momenta around 1 GeV (a range which is relevant for phenomenology), but the equivalence cannot hold in the far infrared.
For all of these reasons, a purely dimensional reduction vacuum wavefunctional is clearly inadequate. Corrections are essential, and what is really required is an approximation to the vacuum state which holds at all distance scales. There are now a number of proposals, which may or may not obtain the dimensional reduction form in some limit, but which do claim to approximate the ground state at all length scales. These we will briefly review.

A. Temporal gauge
It was suggested in ref. [8] that the Yang-Mills ground state wavefunctional, in D = 2 + 1 dimensions and in temporal gauge, is approximated by 2 where B a = F a 12 , D 2 is the covariant Laplacian, λ 0 is the lowest eigenvalue of −D 2 , and m 2 is a parameter which vanishes as g → 0. The motivation was to find the simplest possible gauge-invariant expression which would agree with the freefield (5) and dimensional reduction (7) wavefunctionals in the appropriate limits. In support of this conjecture, it was found that Ψ GO 1. solves the Yang-Mills Schrödinger equation in the strong-field, zero-mode limit; 2. confines if the mass parameter m > 0, and that m > 0 seems to be energetically preferred; 3. produces results for the mass gap, the Coulomb gauge ghost propagator, and the color Coulomb potential, which are in rather good agreement with results derived from standard lattice Monte Carlo simulations.
The subtraction of λ 0 is essential, and was introduced because −D 2 has a positive semi-definite spectrum, and in general the lowest eigenvalue tends to infinity for typical vacuum configurations in the continuum limit. This fact is obvious perturbatively, and is confirmed numerically. Without the subtraction (and this was the form originally suggested by Samuel [6]), the kernel joining B a (x) and B b (y) in (10) effectively vanishes in the continuum limit, and the corresponding string tension would be infinite. In contrast, the spectrum of −D 2 − λ 0 is well-behaved, and not far from that of the free-field Laplacian operator −∇ 2 [8].
If one drops all components of the vector potential apart from the zero mode (analogous to the "minisuperspace" approximation in quantum cosmology), then the Lagrangian and the Hamiltonian operators are simply 2 A factor of g has been absorbed into the definition of the gauge field, so that A k has units of inverse length. This accounts for the overall factor of 1/g 2 in the exponent of the wavefunction.
where V is the volume of 2-space, and the cross-product and dot-product are defined with respect to SU(2) color indices. Solving for the ground state is a problem in quantum mechanics, rather than quantum field theory, and to leading order in 1/V the solution is Now in the region of parameter space where the zero mode is much larger than all other modes, the covariant Laplacian is approximated by and m 2 is negligible. It is then found, after some algebra, that the proposed wavefunctional (10) reduces to the zero-mode solution (12). Dimensional reduction follows by expanding the B-field in eigenmodes φ a n of −D 2 . Then the part of the wavefunctional that depends only on the low-lying modes, with eigenvalues λ n − λ 0 ≪ m 2 has the form of the dimensional reduction wavefunctional (7), with µ = 1/m. If we assume that the asymptotic string tension is due to the low-lying modes, then calculation of the string tension is simply an exercise in twodimensional Yang-Mills theory, and the result is If we turn this around, and write m = 16σ /(3g 2 ), then we have a complete proposal for the vacuum wavefunctional, although the string tension must be supplied as an input. A method for obtaining equal time expectation values by numerical simulation, with a suitable lattice regularization, was also introduced in [8], and applied to calculate the mass gap. The Coulomb gauge ghost propagator and color Coulomb potential were derived via numerical simulation of Ψ 2 GO in [28], by the method of generating thermalized lattice configurations from the Ψ 2 GO distribution, and then transforming these configurations to Coulomb gauge. The results, as already mentioned, were in very good agreement with those obtained from standard lattice Monte Carlo simulations. For details, we refer the reader to the cited references.

B. New variables
While the temporal gauge ground state can be credited with some numerical success, it remains an educated guess, and requires the string tension as an input. A more ambitious program in 2+1 dimensions, which aims to calculate both the Yang-Mills vacuum state and the string tension analytically, was initiated by Karabali, Kim, and Nair [14], and has been further developed by Karabali and Nair in a series of papers, cf. [29] and references therein.
The starting point in the Karabali, Kim, Nair (KKN) approach is temporal A 0 = 0 gauge, and the remaining two components of the A-field are combined into a complex field A = (A 1 + iA 2 )/2, related to a matrix-valued field M via where z = x 1 − ix 2 , and z = x 1 + ix 2 are the usual holomorphic variables in the complex plane. The matrix-valued field M takes values in the group SL(2,C), and transforms covariantly, M → GM, under a gauge transformation G. This field can be used to define gauge-invariant field variables where C A is the quadratic Casimir in the adjoint representation. In terms of these gauge invariant variables, the Hamiltonian becomes where T is derived from the E 2 term in the standard Hamiltonian with and (∂ ≡ ∂ z ) and also Inner products are evaluated with respect to the integration measure where the second line is the new variables state converted back to usual variables. KKN assume that the dimensional reduction form is obtained for long-wavelength configurations by simply dropping −∇ 2 in the kernel, i.e. (26) and then the string tension for a spacelike Wilson loop is obtained from solving Yang-Mills theory in two Euclidean dimensions, with the result Very remarkably, this value is within a few percent of the value found by Bringoltz and Teper [30] in lattice Monte Carlo simulations of the 2+1 dimensional theory, after careful extrapolation to the continuum limit. 3

C. A hybrid wavefunctional
The problem with Ψ KKN is that, in terms of new variables, it is not holomorphic invariant, and in terms of the usual variables (second line of (25)) it is not gauge invariant. Therefore Ψ KKN , as it stands, is not a physical state. Of course, KKN do not claim that Ψ KKN [J ] in eq. (25) is exact, and presumably gauge and holomorphic invariance requires consideration of contributions to R[J ] involving higher powers of J . However, ignorance of the gauge/holomorphic-invariant wavefunctional calls into question the assumed dimensional reduction form (26), which was required for the successful prediction of the string tension. For example, suppose we assume that higher powers of J in the expansion of R[J ] would have, as its main effect, the conversion of the ordinary Laplacian into a covariant Laplacian; i.e. in the usual variables In that case, for configurations which are non-abelian ([A x , A y ] = 0) in general, dropping −D 2 is invalid even for configurations which vary very slowly compared to the length scale 1/g 2 , and indeed is invalid even for configurations which have no spatial variation whatever. As we have remarked above, in connection with Ψ GO , the covariant operator −D 2 has a positive semi-definite spectrum, and for typical lattice configurations the lowest eigenvalue diverges in the continuum limit. In that case, rather than replacing −D 2 by zero to obtain the dimensional reduction result, one should replace it by infinity! This is obviously nonsense.
Assuming that the KKN wavefunctional applies to abelian configurations ([A x , A y ] = 0), the corresponding vacuum state for more general configurations is still a mystery; one can only guess what the gauge and holomorphic invariant completion of Ψ KKN might be. But the gauge-invariant completion is essential, if one is going to invoke dimensional reduction to compute the string tension. At this stage there are an infinite number of possibilities, and the validity of the KKN prediction for the string tension depends on which of these possibilities is the correct one. One possible approach is to retain Ψ KKN for abelian configurations, and ask for the simplest gauge-invariant generalization which would lead to the dimensional reduction form (26). Then it is natural to merge features of Ψ GO and Ψ KKN into a conjectured "hybrid" form for the ground-state wavefunctional which we will include in our numerical tests below.
An alternative approach has been followed by Leigh, Minic, and Yelnikov (LMY) [31], who begin with the ansatz where L = −∆/m 2 , and ∆ is the holomorphic-covariant Laplacian. They then derive and solve a differential equation for K(L), where L is treated as a number, rather than an operator, and by solving this equation they arrive at where J 1,2 are Bessel functions. By construction, the LMY proposal is a physical state. If the infrared limit means L → 0, then K → 1, and Ψ 0 has the dimensional reduction form (26), leading to the same prediction for the string tension. Leigh et al. also obtain predictions for the glueball mass spectrum in 2+1 dimensions, which appear to be in good agreement with standard lattice Monte Carlo results. The reservation in this case is that the LMY approach assumes a certain operator identity (eq. (56) of ref. [31]) whose validity, in our opinion, is questionable. It would nevertheless be interesting to test Ψ LMY numerically, but unfortunately it is not clear to us that the method we will use in this article could be easily applied to the LMY proposal.

D. Coulomb gauge
In Coulomb gauge, after resolving Gauss' law, eq. (3), one obtains the Yang-Mills Hamiltonian [32] in terms of the transverse components of the gluon field, ∇ · A A A = 0, is the canonical momentum (electric field) operator and is the color charge of the gluons and is the so-called Coulomb kernel. The gauge fixed Hamiltonian eq. (32) is highly non-local due to the Coulomb kernel, eq. (35), and due to the FP determinant, eq. (33). In addition, the latter occurs also in the functional integration measure of the scalar product of Coulomb gauge wavefunctionals Any normalizable state, expressed as a functional of the transverse gauge field, is a physical state in Coulomb gauge. This means in particular that a wavefunctional which is Gaussian in the gauge field may be a viable proposal for the ground state. Unlike the GO and KKN/hybrid proposals, such a state cannot have the dimensional reduction property in general, since that property calls for a wavefunctional which, on large scales, is Gaussian in the field strengths rather than the gauge fields. On the other hand, also unlike the other proposals, the Gaussian wavefunctional is tractable analytically. Efforts in this direction were spearheaded by Szczepaniak and Swanson [9,33]. They used a Coulomb gauge ground state wavefunctional of the form The proposal was further developed in ref. [10], where the contribution from the Faddeev-Popov determinant was included at one-loop order. The field-independent function ω(k) was determined from a gap equation obtained by minimizing the energy expectation value. The gap equation depends on the so-called ghost dressing function d(k), which is defined in terms of the expectation value of the inverse Faddeev-Popov operator 4 and the Coulomb form factor, f (k), defined by In terms of d(k) and f (k) the expectation value of the Coulomb kernel in eq. (35), which determines the Coulomb potential V , is given by (41) Finally, inclusion of the Faddeev-Popov determinant at oneloop order introduces dependence on the function 5 which is related to the expectation value of J . In ref. [10] χ(k) (there denoted by F(k)) was derived in context of the gap equation, while the explicit representation of J in terms of χ(k) was derived by Reinhardt and Feuchter in ref. [12] (cf. eq. (47) below). The set of coupled Schwinger-Dyson equations for χ(k), d(k), f (k) and ω(k) is UV divergent and requires renormalization. In the variational approach this is achieved by adding relevant and marginal counter-terms to the Hamiltonian and, if needed, renormalizing the functional measure. The latter was obtained in [10] and reads where I χ (k) is given by the right hand side of eq. (42). In [10] the renormalization program was, however, not fully implemented. In particular a Hamiltonian counter-term proportional to AΠ, which defines the c 1 renormalization constant (cf. eq. (52) below), was omitted and thus only an approximate low-energy solution could be obtained. It was found, however to be qualitatively consistent with the results of [9] that used the J = 1 (χ(k) = 0) approximation. This hints that within the one-loop variational approach, contributions from the FP operator may be accounted for by the gaussian wavefunctional itself, with an appropriate choice of the gaussian parameter ω(k). Such a possibility was rigorously demonstrated by Reinhardt and Feuchter [12] (cf. eq. (46) below and the discussion that follows). Inspired by the wavefunctional of a spinless particle in an s-state of a spherical potential Feuchter and Reinhardt in [11] suggested to use the ansatz which has a number of technical advantages: The factor of J [A] in the integration measure (eq. (36)) cancels against J [A] −1 from the square of the wavefunction and thus drops out from the calculation of equal-time vacuum expectation values. As a consequence Wick's theorem can be applied directly, and in particular ω(k) appearing in eq. (44) is found to be directly related to the static gluon propagator In ref. [12] Reinhardt and Feuchter considered a general wavefunctional of the type (46) In the one loop approximation they showed that the Faddeev-Popov determinant, eq. (33), can be represented as where χ(k), thereafter referred to as the curvature, is given by which, to the order of approximation considered, after renormalization, coincides with the one given in eq. (43). Combining eq. (46) and eq. (47) leads to and establishes equivalence, at a one-loop level, between the ansatz of the Indiana group eq. (37), which corresponds to α = 0, and that of the Tüebingen group eq. (44), corresponding to α = 1/2. 6 However, using equivalent variational ansätze did not lead to the same results for the correlation functions, d(k), f (k), χ(k), ω(k). This is because the approaches of the Indiana and Tübingen groups differ in i) the approximation scheme used to evaluate the expectation value of the Hamiltonian and ii) the renormalization scheme. While the Tübingen group fully includes the Faddeev-Popov determinant to the order considered, the Indiana group set J = 1 throughout ref. [9] and neglected J in the Coulomb term in the numerical calculations of ref. [10]. (In the analytic calculation of ref. [10] J was, however, fully included.) Also, while the Indiana group considers the one-loop corrections to the Coulomb form factor f (k), the Tübingen group employs the d(k) = 1 approximation in the equation for f (k).
Ref. [10], in which the renormalization program was not fully implemented, missed a Hamiltonian counter-term proportional to AΠ, which defines the c 1 renormalization constant (cf. eq. (52) below). The existence of this term was realized by Feuchter and Reinhardt [11], who carried out the complete renormalization program. The c 1 counter-term missed in [10] plays an important role in determining the IR properties of the wavefunctional, as realized by Reinhardt and Epple [35], and will be crucial for the investigations given in the present paper. Therefore throughout this paper we will use the fully renormalized approach of the Tübingen group [11,35].
For later convenience we define where ω(k) corresponds to the wave functional in eq. (44), and write the wave functional of eq. (44) in the form The fully renormalized gap equation for ω, which ultimately determines ω, reads [11,35] with and V (k) given by eq. (41). The gap equation, together with eq. (43) and the Schwinger-Dyson equations for the ghost form factor, and Coulomb form factor, form a closed set of coupled integral equations for χ, d, f and ω. In the gap equation (52), c 1 and c 2 are (finite) renormalization constants. For the critical solution, where one imposes the horizon condition for the ghost dressing function, eq. (38), both ω(k) and χ(k) are infrared divergent, which implies that the transverse gluon propagator vanishes at k → 0, while [35] So even when enforcing the horizon condition, the quantity c 1 = ω(0) is undetermined and may be taken to be either infrared finite or zero. However, a perimeter law of the 't Hooft loop requires c 1 = 0 and this value is also favoured by the variational principle [35]. Furthermore, for c 1 = 0, in the IR limit k → 0, the wavefunctional eq. (51) becomes independent of the gluon zero mode which agrees with the behavior of the exact vacuum wavefunctional in 1 + 1 dimensions [36], and corresponds to the so-called ghost loop dominance in higher dimensions [37]. But although there is strong evidence to favor c 1 = ω(0) = 0, our numerical studies in Section IV B will also look at the case of a non-zero, but small, value for ω(0). The renormalization parameter c 2 , on the other hand, has no influence on the IR or UV behavior of the solutions of the gap equation (52). Only the mid momentum regime of ω(k) is weakly dependent on c 2 [11]. Since we are mainly interested in the IR properties we will put c 2 = 0 throughout this paper. The set of coupled integral equations can be solved analytically in the IR (for the critical solution) using the power law ansätze [11,38] while the full numerical solutions of the above equations were given, for D = 3 space dimensions, in [11,39,40]. For D = 2, the numerical solution was presented in ref. [13] and it will be used in Section IV B for comparison with lattice simulations.
One criticism that can be leveled at the Coulomb gauge proposal is that it is not clear how it could ever lead to an area law falloff for spatial Wilson loops. In order to address this issue, a modified version of a Gaussian ansatz, which incorporates monopole configurations, has been proposed by Matevosyan and Szczepaniak [41]. Furthermore, recently [42] Campagnari and Reinhardt have developed a method which allows to use non-Gaussian wavefunctionals in the variational approach. Specifically, a wavefunctional containing vertices with up to four gluon fields was considered. Tests of these modified versions are, however, deferred to future investigations.

III. THE MEASUREMENT METHOD
We begin with the identity where, in the infinite volume limit, Ψ 0 is the ground state of the operator H, defined via the transfer matrix T = exp[−Ha t ], with a t the lattice spacing in the time direction. In the continuous time limit, H is the Hamiltonian of the lattice gauge theory. Now consider a finite set of lattice configurations k (x), m = 1, 2, ..., M} at a fixed time, and define This is the partition function of a statistical system in which the lattice configurations at time t = 0 are restricted to the set U . The rescaled wavefunctional has the interpretation as the probability P n that, in this statistical system, a lattice configuration on the t = 0 time-slice is equal to the n-th configuration U (n) i (x) ∈ U in the given set. The probability P n can be computed numerically by a modified lattice Monte Carlo simulation. All links at t = 0 are updated in the usual way, which for the SU(2) gauge group with the Wilson action is a simple heat bath. On the t = 0 plane, however, one of the M configurations from the set U is selected at random, and then accepted or rejected by the Metropolis algorithm. Let N n represent the total number of times, in a given simulation, that the n-th configuration in the set is selected by the Metropolis algorithm, with N tot the total number of updates of the t = 0 plane. Then Since Ψ 0 [U (n) ] is simply a constant rescaling of Ψ 0 [U (n) ], it follows that the relative amplitudes of the vacuum wavefunctional Ψ 0 in the set U are given by Now suppose we have some theoretical proposal for the Yang-Mills vacuum wavefunctional If the proposal is correct, i.e. Ψ theory = Ψ 0 , and we make a plot of then the data points should fall on a straight line, with slope equal to one. The method just described was introduced and applied in refs. [15][16][17]. In that early work, however, the simulations were carried out on small lattices and relatively small values of β = 4/g 2 , while comparison to theory was limited to simple wavefunctionals, resembling (6), inspired by the strongcoupling expansion. It is now possible for us to greatly improve on these previous studies.
In this investigation we will consider sets of three different different types of configurations: • Abelian plane waves with fixed wavelength L and variable amplitude where m = 1, 2, ..., m max with L the lattice extension and α, γ some constants.

IV. RESULTS
Since the measurement method in the previous section relies on a lattice regularization, we must apply this regulator to the vacuum wavefunctionals under study. Let us begin with Ψ GO . The proposal is that where R 0 = − log(N 2 ), and in the continuum 7 The factor of 20 in the definition of a (m) is an arbitrary scaling of the parameters, which could of course be absorbed into α,γ.
In the special case of abelian plane waves with A a 1 (x) = A 1 (x)δ a3 , A a 2 (x) = 0, we have the simpler expression The engineering dimension of the kernel, in 2+1 spacetime dimensions, is also inverse length. We now latticize the theory and absorb dimensions into a lattice spacing a, with where ∂ L is the lattice finite difference operator, and all of the other subscript L quantities are dimensionless. All factors of a cancel in R [A], and the result is

A. The GO and KKN wavefunctionals for abelian plane waves
Now we specialize to the lattice abelian plane wave configurations listed in the previous section (lattice sites are x = (n 1 , n 2 ))

Substituting these configurations into R[A], the result is
with and where k and m are the momentum and the mass parameters in physical units, i.e. k 2 = k 2 /a 2 , m = m L /a. The same regularization applied to the KKN wavefunctional yields, for the abelian plane wave configurations, with The theoretical values for ω(k 2 ) are to be compared against the data obtained from the numerical simulation. For a given lattice coupling β E of the Wilson action, at a given lattice size L corresponding to a value of k 2 given in eq. (72), we obtain from the numerical simulation described in the previous section the values Then ω MC ( k 2 ) is obtained from a best linear fit of to the data points {r n }. Figure 1 shows a typical plot of r n vs. 2(α + γn) at β E = 9 and L = 24; ω MC ( k 2 ) is the slope of the line (best linear fit) shown. The values for α, γ used at each β E and L are listed in Table III of the Appendix. The theoretical expressions for ω GO (k 2 ) and ω KKN (k 2 ) involve two dimensionful parameters m and g 2 . Once these parameters are chosen, the results can be compared with the data obtained for ω MC ( k 2 ) on any lattice, providing the dimensionless squared momentum k 2 on the lattice is converted into physical units k 2 = k 2 /a 2 using the lattice spacing a. For a choice of lattice coupling β E , the lattice spacing in physical units is given by where σ L = σ L (β E ) is the D = 3 dimensional string tension in lattice units, and σ is the string tension in physical units. On grounds of tradition, we make the arbitrary choice σ = (440 MeV) 2 . Figure 2 is a plot of ω MC (k 2 ), extracted from a best fit of the data to eq. (78). Each data point is obtained at a particular β E = 6, 9 or 12 on a given lattice of extension L, with L = 16, 24, 32, 40 or 48, and the wavelength of the plane wave on each lattice is the largest wavelength λ = L available. This plot also displays the two theoretical curves with the parameters g 2 and m obtained, for each curve, from a best fit to the data points. Observe that in this range of momentum, the difference between the two fitting functions is essentially negligible, and in fact only becomes noticeable for k 2 > 4 GeV 2 .
Cumulative data for ω MC vs. p 2 in physical units, on lattices of extensions L = 16, 24, 32, 40, 48, and Euclidean lattice couplings β E = 6, 9, 12. The curves labeled "GO fit" and "KKN fit" (there are actually two curves, difficult to distinguish from one another), are the theoretical values for ω GO (p 2 ), and ω KKN (p 2 ), using the parameters of m and g 2 in Table I. The line labeled "Coulomb gauge" is obtained from the ansatz for the Coulomb gauge vacuum wavefunctional Ψ CG [A] (eq. 51) as described in Section IV B.
With the parameters obtained from the fit, we can use dimensional reduction (naively, in the KKN case, as explained in section II C) to compute the string tension, and compare it with our input value of (440 MeV) 2 . Dimensional reduction gives The parameters g 2 , m from the best fit, and √ σ from obtained dimensional reduction, in the GO and KKN cases are shown in Table I The product of m and g 2 , in either the GO or KKN approach, determines the string tension σ in either approach. The dimensionless ratio g 2 /m is an output of the KKN approach, where it is predicted to be π. If m and g 2 are determined from a best fit to the data, then the actual ratio is g 2 /m = 2.95. It is not clear, at this stage, whether this small discrepancy is significant, or should just be attributed to deviations from the continuum scaling due to a finite lattice spacing.

B. Tests of the Coulomb gauge wavefunctional
To test the wavefunctional eq. (51), we first have to transfer it to the lattice. We begin by rescaling the gauge field A i → A i /g so that a prefactor g −2 appears in the exponent of eq. (51), and A i (x) has engineering dimension of a mass. With these conventions, the Fourier transformed kernel ω(k) and curvature χ(k) also have dimensions of mass.
Next we latticize as in eq. (70) and rescale the gauge field again to obtain the dimensionless field 8 A c k (x) ≡ a A c k (ax). For Coulomb gauge fixed connections, it is, in principle, important to use the so-called midpoint rule when extracting the gauge fields from the lattice links U k : As compared to simpler prescriptions such as eq. (72), we have two modifications: 1. The shift in the argument on the lhs ensures that the resulting lattice connection is exactly lattice transversal if the link fields are, After Fourier transformation, the shift leads to a phase factor in the connection which affects general observables but happens to drop out in the (quadratic) exponent R[A] tested here.
2. The η-correction in eq. (82) comes from the SU(2) algebra for parallel transporters over a finite distance a, It is only relevant for very strong fields far from the continuum limit. (In our numerical studies, the correction never exceeded 5%.)

After Fourier transformation
where k i = (2π/L)ℓ i (with −L/2 ≤ ℓ i < L/2), a simple calculation leads to the lattice version of the CG wavefunctional, (84) Notice that the dimensionless momentum argument in the numerical continuum solution of the gap equation is k/g 2 , so that its lattice counterpart becomes To complete the lattice transcription, we only have to find an expression for the function From the leading terms of order O(β −2 ), we find b = 16σ /g 4 and therefore This completes the lattice transformation of the Coulomb gauge wavefunctional. Let us first look at the non-Abelian constant configurations (65). The corresponding lattice connection has the special colour structure A c i ∼ δ c i , but is otherwise constant in space, i.e. Fourier transformation projects out the zero frequency contribution, where the approximation in the second line comes from discarding the η-correction in eq. (82).
From eq. (56), the quantity ω(0) is given by the (finite) renormalization constant c 1 and, as already mentioned in sect. II D, the energetically preferred value is c 1 = 0, which is also required for a perimeter law in the 't Hooft loop [35]. Obviously, with this choice of renormalization constant ω(0) = c 1 = 0 the Coulomb gauge wavefunctional cannot account for the constant non-Abelian gauge field configurations. Whether this failure is important remains to be seen. At least it does not necessarily imply that the Coulomb gauge wavefunctional is a bad approximation to the true vacuum wavefunctional since constant configurations form a set of measure zero in field space. One could give up the preferred value c 1 = 0 and choose ω(0) = c 1 as a fitting parameter, cf. fig.3. This gives reasonable agreement with the lattice data for one set of constant non-Abelian configurations but does not cure the general problem. From the results presented in Sec. IV C below, it will become clear that constant non-Abelian gauge fields can only be accounted for if we include quartic terms ∼ (A A A × A A A) 2 in the exponent of the wavefunctional. The use of such non-Gaussian wavefunctionals in the variational principle has recently become feasible [42], but the solution for the wavefunctional has not yet been determined explicitly up to quartic terms in the exponent.
For these reasons, we will use the energetically favored value ω(0) = c 1 = 0 in the following. We will now show that the Coulomb gauge wavefunctional does a good job for Abelian plane waves of the type eq. (64). In this case we have carried out simulations at β = 6 on a fixed lattice volume of extension L = 24, and varied the amplitude of the plane waves, at given wavelength L/M, according to (90) This can be simplified considerably, if the η-correction in the definition of the connection, eq. (82), is discarded. Then the sums in eq. (90) can be performed explicitly and we obtain a much simpler expression where c M = 2 for the highest frequency M = L/2 and c M = 1 otherwise for L even (L = 24 in this case). From eq. (91), it is obvious that the plane wave configuration tests the kernel ω = ω/g 2 − χ/g 2 exactly at the lattice momentum p M which corresponds to the frequency of the plane wave. Figure 4 shows the result of the numerical evaluation of eqs. (90), (91) against the lattice MC data for Abelian plane wave configurations of varying wavenumber and amplitude. As can be clearly seen, the individual plane waves with fixed wavenumbers M and varying amplitude fall on a straight line, but the slope of that line differs from unity. (We have chosen the solution ω(k) of the variational problem with the preferred renormalisation constant c 1 = 0.) Morever, the slopes of the lines vary slightly with M, i.e. effectively with the momentum picked by the plane wave: For the smallest momentum M = 1, we find a slope of 1.19, which decreases down to 1.02 for M = 2, and then increases again up to 1.52 for the largest momentum M = 12 representable on a L = 24 lattice. If we relax the condition on the renormalisation constant c 1 and take it as a free parameter, we observe that the spread in the slope between the various wave numbers is increased, which is another hint that the choice c 1 = 0 should be preferred. Since the plane waves test the kernel ω(k) at varying momenta, we can use a fit to the MC data as explained in the previous section to find a numerical estimate ω MC (k). In the Coulomb gauge wavefunctional, this quantity corresponds to ω(k) = g −2 ω(k) − χ(k) . After rescaling to physical units (see eq. (86) and below), the result is plotted along with the values obtained by numerical simulation, ω MC (k), in fig. 2. It is evident that the variational solution for ω(k) fits the MC data very well, at least in the infrared region for momenta up to k ≈ 1.3 GeV. For larger momenta, ω(k) starts to deviate and becomes slightly larger than the numerical estimate, but at most by a few percent within the phenomenologically relevant mid-momentum regime. (For very large momenta not plotted here, ω(k) ∼ k is exact by asymptotic freedom.)

C. Non-abelian constant configurations: fixed amplitude, variable "non-abelianicity"
For general non-abelian configurations we have, in a lattice regularization, where with U(P x ) a product of links around a plaquette, starting with a link at site x. The lattice covariant Laplacian, in the adjoint representation, is given by In terms of the parameters g 2 , m in the GO row of Table I, we use β = 4/(g 2 a) and m L = ma, where a is the lattice spacing. For comparison with the Monte Carlo data generated at the lattice coupling β E of the Wilson action, we determine a from eq. (79). It is important to note that while we expect β /β E → 1 in the continuum limit, this ratio need not be exactly equal to one at any finite β E . In the same way, the latticized "hybrid" wavefunctional is with β , m L determined using the parameters g 2 , m in the KKN row of Table I, and the lattice spacing from eq. (79). We will consider first the configurations of eq. (66), with fixed amplitude and variable "non-abelianicity" specified by the θ parameter. If the amplitude is chosen small enough, then −D 2 − λ 0 is negligible compared to m 2 , and the kernel reduces to for the GO wavefunctional, and for the hybrid. This is the dimensional reduction limit, and in either case, for the configurations (66), For the Coulomb gauge wavefunctional, however, R[U] ∝ A 2 1 + A 2 2 , and hence, since the amplitudes of A 1 and A 2 are fixed in the set (66), independent of the angle θ n . If ω(0) = 0, which seems optimal for agreement with the plane wave data, then R CG would also be independent of the amplitude of the gauge fields. However, it is important to recall that the Coulomb gauge wavefunctional should not be evaluated outside the first Gribov horizon. So even if ω(0) = 0, the restriction to the Gribov region amounts to a cutoff in the amplitude of non-abelian constant configurations. The Monte Carlo simulation was carried out on a 12 3 lattice at β E = 6, with the t = 0 configurations chosen from U (n) with α = 0.193, and θ n = (n − 1)π/38. By explicitly calculating numerically the lowest lying eigenvalues of the lattice Faddeev-Popov operator, we have checked that these lattice configurations are all inside the first Gribov horizon.
In Fig. 5 it can be seen that the logarithm of the wavefunctional is indeed proportional to sin 2 (θ ), as one would expect from the GO and hybrid wavefunctionals in the dimensional reduction limit. The data does not seem to be compatible, however, with the θ -independence (99) of the CG wavefunctional (51).
We recall that if Ψ[U] = exp[− 1 2 R(U)] is the true vacuum state, then the data points for − log(N n /N T ) vs. R[U n ] should fall on a straight line, with unit slope. Plotting the data for − log(N n /N T ) against R GO [U n ], as in Fig. 6, we find the slope obtained from a linear fit through the data is indeed close to unity. In the GO case the slope is 1.02(6); a similar analysis for the hybrid wavefunctional results in a slope of 1.12 (7). Some numerical details concerning the simulations are found in the Appendix.

D. Non-abelian constant configurations: variable amplitude, maximal "non-abelianicity"
We now consider the non-abelian constant configurations of maximal "non-abelianicity,", i.e. θ = π/2, which are the configurations of eq. (65), with index m running from 1 to 20. All Monte Carlo calculations were carried out on lattices of volume 32 3 at β E = 6, 9, 12, and the corresponding values of β , m L at each β E are given in Table II, where the values for the hybrid wavefunctional are taken to be the KKN values, since the hybrid reduces to the KKN form on abelian configurations. The test of the GO and hybrid wavefunctionals is to see whether or not the data points for − log[N n /N tot ], when plotted against R[U (n) ], fall on a straight line whose slope is close to unity.  An example of the − log[N n /N tot ] vs. R GO [U (n) ] data at β E = 6 is shown in Fig. 7, for the choice α = 2, γ = 0.15. Although the data is nicely fit by a straight line which has a slope close to unity, this fact must be interpreted with caution because, since the number N n falls off exponentially with R GO [U (n) ], the range of R must necessarily be kept small; typically ∆R ≈ 4 − 5. This could mean that the tendency of the data to lie on a straight line is misleading, and perhaps we are simply looking at the tangent of a curve. It is therefore necessary to extract the slope of the straight line over small intervals centered around points over a wide range of R. The question  is whether those slopes are constant, in which case the linearity hypothesis is verified, or whether they vary significantly as R increases. This is the motivation to calculate − log[N n /N tot ] in sets of twenty configurations, using different values of the parameters (α, γ) in each set. The parameters we have used are shown in Table IV of the Appendix. Figure 8 is a plot of the slope vs. R at β E = 6, 9, 12, where the value of R at each data point is the midpoint of the range in which the slope was computed. Things are not perfect; there is some slight variation in the slope with R, there is a little variation with β , and the values of the slope are not exactly one (they seem to be closer to 1.1 at the large R values). On the other hand, we have made no claim that the GO wavefunctional is exact, nor is asymptotic scaling exact at these lattice couplings. The point is that scaling is not bad, and the slopes are fairly close to unity over a large range of R, using g 2 , m values that were extracted from fits to a completely different type of lattice configuration (i.e. abelian plane waves).
Results for the hybrid wavefunctional turn out to be quite close to those of the GO wavefunctional. The values for β E = 12, for both types of wavefunctionals, are shown in Fig. 9, with similar agreement at the two other β E values.

E. The ghost propagator and the Coulomb potential
Because of the equality (4) of the vacuum wavefunctionals in temporal and Coulomb gauges, when evaluated on transverse (∇ · A = 0) gauge fields, equal-time expectation values in Coulomb gauge can be derived from and we may use for Ψ 0 either of the temporal gauge proposals, Ψ GO , Ψ hybrid , or the Coulomb gauge proposal Ψ CG to calculate such objects as the ghost propagator aa xy |x−y|=R (102) and the color Coulomb potential 9 In eq. (101) there is an implicit restriction of the integration domain to the Gribov region. For the Coulomb gauge wavefunctional Ψ CG [A] the ghost propagator and the Coulomb potential are presented in [13].
In an ordinary Monte Carlo (MC) simulation, Coulomb gauge expectation values are obtained by first generating lattice configurations with the usual probability distribution exp[−S]/Z, where S is the standard lattice action, transforming those configurations to Coulomb gauge, and evaluating the observable Q in the ensemble of transformed configurations. In principle the same strategy applies to evaluating the right hand side of (101) numerically; the problem in that case is to generate configurations with the probability distribution Ψ 2 [U], and this problem was solved, for the Ψ GO proposal, in ref. [8]. The simulation method developed in [8] is also applicable (although it has not been applied until now) to the hybrid proposal. The lattice ghost propagator and Coulomb potential were calculated numerically from Ψ GO , and compared to the corresponding results in ordinary lattice Monte Carlo, in ref. [28]. In that work, however, the authors chose β = β E and m L = 4β σ L /3. In the present article the philosophy has changed somewhat. We have two parameters with dimensions of mass, g 2 and m, and a scale set (arbitrarily) by taking √ σ = 440 MeV. Then g 2 , m are chosen to give a best fit to the abelian plane wave data in Fig. 2. To compare wavefunctional results with standard Monte Carlo results we determine the lattice spacing a, at each β E , from σ L /σ , and then β = 4/(g 2 a) and m L = ma 9 More precisely, for color charges in some representation r, the Coulombic potential energy is obtained by multiplying V c (R) by the quadratic Casimir C r , and dividing by the dimension of the adjoint representation. are the corresponding dimensionless parameters to use in the latticized wavefunctional Ψ GO or Ψ hybrid . With the new procedure we have β = β E , and the obvious question is whether this fact will tend to destroy the agreement that was found previously, in [28], between ghost propagators and Coulomb potentials derived from simulation of Ψ 2 GO , and the corresponding quantities found in ordinary lattice Monte Carlo simulations. We would also like to calculate the Coulomb gauge ghost propagator and Coulomb potential for the hybrid wavefunctional proposal. Figure 10 shows the equal-times ghost propagator G(R) computed in a standard Monte Carlo simulation on a 32 3 lattice at β E = 9. On the same plot we see the corresponding results obtained by generating lattices with probability distribution Ψ 2 GO and Ψ 2 hybrid by the methods of [8], transforming to Coulomb gauge, and evaluating the ghost propagator, in each case using the appropriate values of β , m L corresponding to β E = 9. It can be seen that the agreement between Monte Carlo, GO, and hybrid results is almost perfect.
The agreement for the Coulomb potential V c (R) is not as good. In Fig. 11 we display the data from MC, GO, and hybrid simulations, again at β E = 9, with a cut in the data, discarding configurations with |V (0)| greater than some bound equal to 5, 10, 50, 300. If we restrict the data set to configurations with |V (0)| < 5, then the agreement between MC, GO, and hybrid results is again almost perfect. Roughly half of all configurations meet this criterion. The agreement is still fairly good for |V (0)| < 10, which accounts for about 80% of all configurations. However, as the cut is gradually removed, the Coulomb potential derived from GO and hybrid simulations, while roughly linear in R, deviates quantitatively from the MC result. But how can there be such a noticeable deviation when the ghost propagators agree so accurately, without any cuts at all? The explanation probably has to do with a discrepancy in the tail of the probability distribution. If two probability dis-  11. Data for the Coulomb potential at β E = 9 and L = 32, derived from MC, GO and hybrid simulations, with a cut on the data, discarding configurations for which |V 0 | is greater than 5, 10, 50, and 300, respectively. tributions agree in their lower moments, but disagree in higher moments, then it means that the two distributions agree pretty well where the probability is substantial, but disagree in the tail of the distributions. That is what seems to be going on here.
What was found already in ref. [28] is that the Coulomb potential is quite sensitive to a comparatively small number of "exceptional" configurations, in which the lowest eigenvalue of the Faddeev-Popov operator −∇·D is far below the average value for the lowest eigenvalue. The reason that such exceptional configurations are relevant for the Coulomb potential, but not the ghost propagator, is presumably because the ghost propagator involves only one factor of the inverse F-P operator, while the Coulomb potential involves two factors. Because the inverse F-P operator becomes singular as the lowest eigenvalue λ 0 approaches zero, higher powers of the inverse F-P operator (such as the Coulomb potential) will be more sensitive to infrequent configurations with exceptionally low values of λ 0 than lower powers (such as the ghost propagator). The probability distribution of infrequent configurations is, of course, governed by the tail of the probability distribution. So our interpretation of the ghost and Coulomb propagator results is that Ψ 2 GO and Ψ 2 hybrid agree quite closely with each other, and with the probability distribution of the true Yang-Mills vacuum wavefunctional Ψ 2 0 , in the "bulk" of the distribution. The Coulomb potential data suggests, however, there is some small disagreement in the tail of the distribution.
In general, our results for the Coulomb gauge ghost prop-agator and Coulomb potential with the new fitting procedure for β , m agree quite closely with our previous results (based on setting β = β E ) reported in ref. [28] (for a quantitative comparison, cf. [44]). The GO and hybrid results are, once again, virtually indistinguishable. Since both choices of parameters, and the GO and hybrid wavefunctionals, have about the same dimensional reduction limit, our results suggest that the quantities we have computed, at the couplings we have employed, are mainly sensitive to that limit.

V. CONCLUSIONS
We have compared several suggestions for the Yang-Mills vacuum wavefunctional to the true Yang-Mills vacuum wavefunctional in 2+1 dimensions, whose exact form is unknown, but whose relative magnitudes in any set of lattice configurations can be obtained numerically. Three types of lattice configurations were studied: abelian plane wave configurations, non-abelian constant configurations of fixed amplitude but varying "non-abelianicity," and non-abelian constant configurations of maximal abelianicity but of differing wavelengths and varying amplitudes. For purposes of comparison, the physical scale was set by taking the string tension to be √ σ = 440 MeV.
For abelian plane waves, up to the shortest wavelength corresponding to p 2 = 2.5 GeV 2 that we have investigated, the GO and Karabali-Kim-Nair proposals are almost indistinguishable, and both agree very well with the values obtained for the true vacuum wavefunctional, evaluated on these configurations. The Coulomb gauge wavefunctional can also fit the plane wave data with an appropriate choice of parameters, providing in particular that the renormalization constant c 1 in eq. (41) is set equal to zero. Both the GO and KKN wavefunctionals reduce to the dimensional reduction form exp[−µ F 2 ] at long wavelengths, and it seems likely that this is also true for the Coulomb gauge proposal, in this special case of abelian configurations, for the choice of renormalization constant c 1 = 0.
For non-abelian configurations, we have suggested a gaugeinvariant wavefunctional which reduces to the KKN proposal for abelian configurations, and incorporates the covariant Laplacian and eigenvalue subtraction of the GO proposal, which we have termed the "hybrid" wavefunctional. Both the GO and hybrid wavefunctionals have the dimensional reduction form when restricted to configurations which, when expanded in eigenstates of the covariant Laplacian, contain only low-lying eigenmodes. Once again, the GO and hybrid wavefunctionals are almost indistinguishable when evaluated on non-abelian constant configurations, and this is probably because they have almost the same dimensional reduction limit. We find that the GO and hybrid wavefunctionals are in good agreement with the true vacuum wavefunctional for non-abelian constant configurations, as well as for abelian plane waves. The Coulomb gauge wavefunctional, however, which does not have the dimensional reduction property for non-abelian lattices, does not seem compatible with the data for non-abelian constant configurations, particularly the data with variable non-abelianicity.
The Coulomb gauge wavefunctional has been used to compute Coulomb gauge ghost and gluon propagators, with results in 2+1 dimensions, reported in [13], indicating a Coulomb potential rising almost (but not quite) linearly. We have also computed these quantities by direct simulation of the GO and hybrid wavefunctionals. The GO and hybrid results agree with one another, and almost perfectly with the lattice Monte Carlo results for the ghost propagator. The GO and hybrid wavefunctionals also lead to an apparently linear Coulomb potential and agree very closely with each other. On the other hand there is some difference in the GO and hybrid Coulomb potentials in comparison to the lattice Monte Carlo results, and this can be attributed to a difference associated with exceptional configurations with unusually small values of the lowest Faddeev-Popov eigenvalue. Thus the GO and hybrid wavefunctionals would seem to agree with the true Yang-Mills vacuum wavefunctional for the bulk of the probability distribution, but there would appear to be a small disagreement in the tail of the distribution.
The main effort in this article has been to calculate the relative magnitudes of the true vacuum wavefunctional on particular sets of lattice configurations; namely, abelian plane waves and non-abelian constant configurations, and to compare those results with a number of proposals for the vacuum state. We have found that the lattice data for the abelian plane waves have been nicely reproduced by all proposals consid-ered, while good agreement with the data for non-abelian constant configurations appears to require wavefunctionals with the property of dimensional reduction.
where X is a matrix, and the coefficients are given by [45]  with an implicit summation over lattice sites x and color indices a. All the matrix operations, including the determination of λ 0 , can be carried out numerically using sparse matrix techniques, which results in a considerable savings in computation time, often by an order of magnitude or more in our calculations. We have checked the accuracy of the Zolotarev approximation by evaluating R numerically, in several cases, without this approximation, and have found the results with and without the approximation to differ only at the third significant digit. This is sufficient for our purposes. In the case of R hybrid the formula (A.2) is not directly applicable, and the numerical evaluation was carried out without the help of the Zolotarev approximation.
In the Monte Carlo simulations, we set up eight runs each time with the same parameters, but different seeds for the random number generator. Each run is itself a number of independent jobs, which we refer to as "cycles", whose results for − log(N n /N T ) are averaged together at the end of the run. At the beginning of each cycle the links are all set to the identity matrix, except for the spacelike links on the t = 0 plane, which are set to the first (n = 1) configuration out of the set of {U i (x,t = 0)} of non-abelian constant configurations. The lattice at t = 0 then thermalizes for 5000 sweeps with the n = 1 configuration at t = 0 held fixed. All timelike links are fixed to the unit matrix, except for the timelike links at t = L/2, which are updated in the usual way. After thermalization we carry out another 30000 sweeps, with the configuration at t = 0 updated only once every 40 sweeps. On reaching the t = 0 plane every 40th sweep, we carry out 20 Metropolis "hits"; i.e. the Metropolis algorithm is used to update the t = 0 plane, and at each hit the plane is changed to a new configuration (or not, depending on the result of the algorithm), and the appropriate configuration counter N n is incremented. At the end of each cycle the value for − log(N n /N T ) for each configuration n is recorded. At the upper end (higher n) it is usually the case that N n = 0 on one or more cycles; all data from these higher n configurations are deemed statistically unreliable, and discarded. The number of cycles (used for eight runs at the same set of parameters) varied from a minimum of 17 to a maximum of 70, but was mostly around 30. The result for the slope of − log(N n /N T ) vs. R[U (n) ] was obtained from the best fit to the data in each run, and the results from the eight independent runs were used to estimate the error.
Finally we record, in Tables III and IV, the values of α, γ used, in eqs. (72) and (65), to generate sets of abelian plane waves and non-abelian constant configurations with varying amplitudes. The aim, in choosing parameters, was to keep the variation of r n = − log(N n /N tot ) in a relatively small range ∆r n ≈ 4 (otherwise, because of the exponential falloff, there would be few or no data points at the larger values of n). In the case of non-abelian constant configurations, we choose different α values so as to sample the slope of − log(N n /N tot ) vs. R[U] in a small interval of ∆R, centered around a wide range of values of R, as explained in subsection IV D.