Noise reduction algorithm for Glueball correlators

We present an error reduction method for obtaining glueball correlators from monte carlo simulations of SU(3) lattice gauge theory. We explore the scalar and tensor channels at three different lattice spacings. Using this method we can follow glueball correlators to temporal separations even up to 1 fermi. We estimate the improvement over the naive method and compare our results with existing computations.


Introduction
Glueballs are colourless states in Quantum Chromodynamics (QCD) made entirely out of gluons. QCD predicts the existence of glueballs. However, no glueball has yet been discovered unambiguously even though there are several candidate glueball resonances, such as f 0 (1370), f 0 (1500), f 0 (1710), f J (2220) etc. [1]. One reason is that glueball states can mix with mesons in the same J P C channel and so it is very difficult to unambiguously extract glueball masses experimentally. It remains nevertheless a very exciting proposition and for a recent review on the status of glueball searches we refer the reader to ref. [2].
Glueball masses can be computed in lattice quantum chromodynamics and a lot of effort has been directed towards this computation. However there is still no consensus regarding the mass spectrum. It is a difficult computation in lattice QCD with dynamical fermions due to the high masses of the glueballs (> 1 GeV) and their mixing with mesonic operators in the same symmetry channels. In recent times computation of glueball masses in lattice QCD with dynamical fermions have been attempted in references [3][4][5].
Glueball masses are often computed in pure Yang-Mills theory. Advantages are that there is no mixing with mesonic operators and the glueballs are stable as they cannot decay. Thus it is much easier to extract the glueball masses from monte-carlo simulations of pure Yang-Mills theory than lattice QCD with dynamical fermions. Nevertheless, even in this theory, glueball correlation functions are dominated by statistical noise at large temporal separations and contribution from excited states at short separations. Hence it becomes very hard to extract single masses reliably from the correlators. Fitting of glueball correlators with reasonable accuracy becomes very difficult and one is forced to compute the "effective mass" which is the logarithm of the ratio of the values of the correlation function between successive time slices. If the effective mass does not vary over a significant temporal range then one assumes that the effective mass is the same as the globally fitted mass.
To remove the effect of excited states, conventional methods involve computing correlation matrices with matrix elements between a large set of interpolating operators constructed from smeared or fuzzed links [6] in the relevant symmetry channel 1 . The ground state is obtained by diagonalizing the correlation matrix in each channel [8,9]. As it is difficult to follow the correlator signal to large physical distances, even using the above techniques, one often uses asymmetric lattices with a significantly smaller temporal spacing compared to the spatial lattice spacing with the expectation to observe a flat behaviour of the effective masses [10,11] over several time slices.
A different approach is to use noise reduction algorithms. Such algorithms have been used in the past for computing the glueball spectrum for U(1), SU(2) and SU(3) lattice gauge theories [12][13][14][15]. Barring the U(1) case, it has so far not been possible to follow the correlators to large enough distances (significantly beyond 0.5 fermi) to extract masses from them using global fits.
In this article we follow the latter approach. We restrict ourselves to pure Yang-Mills theory with gauge group SU(3) and employ only the standard operators in each J P C channel (scalar and tensor) but try to follow the correlator to large temporal separations using a new noise reduction algorithm. Since the contamination due to excited states fall off exponentially, we expect that correlators at distances beyond 0.5 fermi to be dominated by the ground state.
In section 2 we describe the algorithm. Section 3 is devoted to our results on the correlators and masses for the scalar and tensor channels. In section 4 we discuss the improvement obtained over existing conventional methods. Finally in section 5 we draw our conclusions and outline directions for future studies.

Algorithm
We compute glueball correlators using Mote-Carlo simulations of SU(3) lattice gauge theory with the Wilson gauge action at three different lattice spacings for both the scalar and the tensor channels. For updating the links we use the usual Cabibo-Marinari heat-bath for SU(3) and use three over-relaxation steps for every heat-bath step. We set the scale on the lattice through the Sommer parameter r 0 [16]. Our simulation parameters are given in table 1. The Sommer parameter for our lattices have been computed in [17] and we use those values.
The noise reduction scheme we implement follows the philosophy of the multilevel algorithm. The multilevel algorithm was introduced in [18] as an exponential noise reduction technique for measuring polyakov loop correlators in lattice gauge theories with a local action. However the principle is general and can be applied to other observables as well. In addition to polyakov loop correlators, it has been used to measure observables such as the polyakov loop [19], wilson loop [20], components of the energy-momentum tensor [21] as well as the glueball mass spectrum [12][13][14][15].  Table 1: Simulation parameters for all the lattices. Lattices A,B and C were used for the scalar channel while D,E and F were for the tensor channel.
The main principle of the multilevel algorithm is to compute expectation values in a nested manner. Intermediate values are first constructed by averaging over sub-lattices with boundaries and then the full expectation values are obtained by averging over the intermediate values with different boundaries. We refer the reader to ref. [14] for details.
For our implementation, we slice the lattice along the temporal direction by fixing the spatial links and compute the intermediate expectation values of the glueball operators by performing several sub-lattice updates. Individual correlators are created using products of the averaged operators at different time slices. The scheme is depicted in figure 1 and the extents of the sub-lattices and number of sub-lattice updates along with other simulation parameters are shown in table 1.
The glueball operators between which we compute our correlation function (source and sink) are extended wilson loops denoted by P ab where a, b go over the three spatial directions x, y, z. The operators are projected to zero momentum states as usual. We denote the temporal separation between the source and sink operator by ∆t. The sizes of the loops used for the different lattices are given in table 1. Correlation functions between large loops have the advantage that they have much less contamination from the higher excited states compared to those between elementary plaquettes. Such an approach was reported in [22]. There, however, single exponential fits to the correlators were not possible as the data was too noisy. Nevertheless it was observed that glueballs seemed to have the largest overlap with loops of spatial extent 0.5 fermi in each direction. We therefore choose loops of roughly the extent r 0 × r 0 to construct our glueball operators.
Our first noise reduction step was to use a semi-analytic multihit on the SU(3) links [23] with which the wilson loops were constructed and in addition use sub-lattice updates to obtain the expectation values of the glueball operators with very little noise. The choice of the number of sub-lattice updates "iupd" is an important parameter of the algorithm. For the tensor channel, the rule of the thumb we follow is that the operator expectation value over the sub-lattice updates should be the same order as the square root of the correlator at a large value of ∆t. For the scalar channel the same holds but for the connected parts. We compare the overall noise reduction of our algorithm with the naive method (where operators are constructed from elementary plaquettes and only full lattice updates using  heat-bath and over-relaxation is used) in section 4. The multilevel algorithm is very efficient for calculating quantities with very small expectation values. While the operators in the tensor channel viz. E 1 = Re (P xz − P yz ) and E 2 = Re (P xz + P yz − 2P xy ) have zero expectation values and are therefore ideal for direct evaluation using the multilevel scheme, the scalar operator A = Re (P xy + P xz + P yz ) has a non-zero expectation value which has to be subtracted to obtain the connected correlator. For the scalar channel, an alternative is to evaluate the derivative of the glueball correlator directly as that does not need a subtraction. This was carried out in [14] for the U(1) case. In our current calculations we found that evaluating the derivative is more profitable at stronger couplings and we did that for our lattice A. However for the other lattices we found it was better to do the subtraction.
We observed one more phenomenon which is particular to this algorithm. For the smaller values of ∆t where most contributions come from slices which are within the same sub-lattice, there are strong effects due to the short temporal extent of the sub-lattice itself. In such cases we were forced to take into account only correlators between those time slices which lay in different sub-lattices. We found this effect to be significant only in the tensor channel (probably because of the larger value of "iupd" in those cases).

Results -masses
In this section we describe our fitting procedures and the masses we obtain. All the correlators except for the scalar channel at β = 5.7, where we compute the derivative directly and not the correlator, are fitted to the form where m is the glueball mass and T is the full temporal extent of the lattice. Since the correlator is symmetric about T /2, as usual, we fold the data about T /2 and use only one half of the temporal range for the fits. For the scalar channel at β = 5.7, following [14], we use the fit form The fit range was decided on the following two criteria: (i) the range should extend to as large a value of ∆t as possible and (ii) the fit to the form in eq. 3.1 should have f stands for degrees of freedom) to be around one. The fit range for all the different channels and couplings are indicated in table 2. For the tensor channel, we compute the masses using both E 1 and E 2 . In table 2 we report masses only from the operator E 1 , but we checked that masses computed from E 2 agreed within error bars with those from E 1 in each case.
In addition to masses from global fits, we also compute the effective masses from the correlators as To estimate the error on the effective masses we take C(∆t) to be a jackknife bin and we compute masses for 20 such bins. The error on the effective mass is the jackknife error computed from the spread of the masses from the different bins.
In figures 3 and 4 we plot the correlators along with the effective masses for each channel and coupling. Even though the fits were done on the folded data, in the figures we plot the fitted correlator on the full range. It can be clearly seen, especially in the tensor channel that the correlators have contamination form the excited states for the smaller values of ∆t. The same thing is seen for the effective masses. The masses fall at first and then stabilize to a plateau. In our figures we have plotted the mass range obtained from fitting the correlator along with the effective masses to check the extents to which the two results agree. Except for the scalar channel at β = 5.8, we find good agreement. For this case we need a longer lattice in the temporal direction to see the single mass plateau. By the time the excited states die out ∆t is too close to T /2 to obtain the effective mass through equation (3.3). A numerical solution to C(∆t+1) yields the fitted mass though.
We cross-check our data by comparing them against known results. Unfortunately we do not have continuum results yet so we cannot use existing data from asymmetric lattices. We therefore compare our masses with results in [12,13,15]. In [12] the scalar effective mass is computed to be 0.929(49) at β = 5.7 and in [15] the mass is computed from the ratio of partition functions also at the same β to be 0.935(42). This compares quite well with our global fit value of 0.924(33) at β = 5.7. In the tensor channel we could not unfortunately find any data for comparison at exactly the same β value. In [13] the effective mass at β = 6.0 has been computed between time slices 2 and 3 to be 1.078 (16). At β = 6.07 between the same time slices we obtain the effective mass to be 1.061(3). Our global mass for this β value is 0.487 (18) but that is reached only beyond a ∆t of 7.  Figure 2: The derivative of the scalar correlator at β = 5.7. Adopting the notation in [14] δ and δ * represent the forward and backward derivatives at the source and sink time slices.

Results -algorithmic gains
Let us start by discussing the advantage of the current algorithm over the naive method.
To get an idea of the advantage this algorithm brings, we did a few runs for the same computer time using both the methods. Since it is not yet clear how the algorithm behaves as either ∆t or β changes we report our experience for different values of ∆t and β. For the tensor channel at β = 5.8 we carried out runs for 200 hours. The multilevel algorithm had an error of 3% at ∆t = 3 which is just below r 0 (see table 1), while the naive algorithm had an error of 81%. It would be interesting to compare the performance at a value of ∆t between r 0 and 2r 0 . So we choose points around 1.5r 0 (in this case ∆t = 6). Even after 200 hrs of runtime we do not have a signal at that distance for the naive method. So to estimate the % error we multiply the naive correlator at the largest value of ∆t where we have a signal (i.e. ∆t = 3) by corr multilevel (∆t = 6)/corr multilevel (∆t = 3). Doing this we get the % error to be 850% for the naive method while it is 29% for the multilevel   scheme. Thus for the tensor channel at β = 5.8 we estimate that the error reduction algorithm produces an error which is between 27 times smaller than the naive method at both ∆t = 3 and 6. Since the error ∝ time 2 we estimate the new method is more efficient by at least a factor of 729 or so.
At β = 5.95, the runs were for about 100 hours. There at ∆t = 3 the multilevel algorithm produced an error of about 4% while the naive algorithm had an error of 75%. Doing a similar estimate as β = 5.8 we estimate that at 1.5r 0 (∆t = 8) the errors are 150% for the multilevel algorithm while it is about 3000% for the naive method. Thus we estimate that at this β value, the gains in terms of % error is about a factor 20. At β = 6.07, we did not get a signal for the naive algorithm for any ∆t other than ∆t = 1 even after about 300 hours of runs. Thus we see that the gain has very little dependence on ∆t but does depend on β.
For the scalar channel at β = 5.8 the runs were carried out for about 1000 mins. In this case we have a signal at 1.5r 0 for both methods and we get an error of 13% for the multilevel scheme while it is about 70% for the naive method. Thus the gain in terms of % error is about 5.5 or in terms of time about 30. At β = 5.95 in the scalar channel, again we do not have a signal at 1.5r 0 using the naive method and are forced to use the same method as in the tensor channel to estimate the errors. At ∆t = 3 we obtain the errors to be 2% and 37% for the multilevel and the naive methods respectively while at 1.5r 0 they are 29% and 500% (estimated) respectively. Thus the ratio of errors is about 18 or gain in terms of time 324. At β = 5.7 we got a modest gain of about 2.5 in terms of error or about 6 in terms of time.
A compilation of our results is presented in table 3.

Conclusions
Extraction of glueball masses from correlators is a difficult problem in lattice QCD due to a very low signal to noise ratio at large Euclidean times. In this article we present a new method, based on the multilevel scheme, to enhance the signal to noise ratio in glueball correators. We observe that this error reduction technique works quite well at least in pure gauge theories. For a given computational cost, the improvement in the signal to noise ratio is several times to even a couple of orders of magnitude. We are able to follow the correlator to temporal separations of about 1 fermi and can perform global fits to the correlators between 0.5 and 1 fermi. Our effective masses also show a plateau in the same range obtained from the global fits. The current simulations are on relatively coarse lattices and we see a significant dependence of the masses on the lattice spacing. We therefore cannot yet extrapolate to the continuum. It is of course of interest to reach the continuum limit and we are continuing our runs at finer and larger lattices and will report our results in subsequent publications..