# 1 Introduction

## 1 Introduction

Lattice gauge fixing is unavoidable to compute non perturbatively the propagators of the fundamental fields appearing in the QCD Lagrangian . It is also necessary in non-perturbative renormalization schemes  which use gauge dependent matrix elements to renormalize composite operators, and it becomes a fundamental technical ingredient in the so called non gauge invariant quantizations of chiral theories .   
Up to now, the Landau gauge is the only covariant gauge for which an efficient numerical algorithm is available.

A few years ago, an unconventional method to implement a generic covariant gauge on the lattice was proposed . In principle this procedure takes into account the contributions from each gauge orbit with appropriate weight. With this method a gauge dependence of the gluon propagator, in particular at zero momentum, had been found . On the other hand the implementation of this method for physical lattices can be numerically demanding.

Here we follow a more conservative procedure to fix a generic covariant gauge on the lattice that stems from ref. . It is based, as in the Landau case, on the minimization of a functional $`H_A[G]`$ chosen in such a way that its absolute minima correspond to a gauge transformation $`G`$ satisfying the appropriate gauge condition (see section 2). In the continuum it is easy to show  that this procedure is equivalent to the Faddeev-Popov quantization for covariant gauges, i.e. the two procedures lead to the same matrix elements for a generic gauge dependent operator. As usual the Faddeev-Popov factor can be written as a Gaussian integral of local Grassman variables, the resulting effective action is invariant under the BRST transformations and the correlation functions of the operators satisfy the appropriate Slavnov-Taylor identities. On the compact lattice the presence of Gribov copies is an obstruction to define partition functions with exact BRST invariance as it has been shown by H. Neuberger a few years ago . The Neuberger’s argument concerns a fundamental theoretical point whose relevance for the measurement of gauge dependent operators on the lattice must be clarified. Nevertheless, we believe that new numerical gauge fixing algorithms represent a concrete tool to improve the understanding of the lattice gauge properties.

In the following we propose a new numerical gauge fixing procedure which generalizes the usual Landau gauge fixing to generic covariant gauge. The numerical implementation of our procedure on the lattice uses a straight-forward generalization of the standard Landau algorithm, i.e. a steepest descent iterative algorithm minimizes a discretization of the new functional $`H_A[G]`$. In order to have an efficient and feasible gauge fixing procedure, the discretized functional $`H_U[G]`$ must be chosen carefully. In fact a naïve discretization of this functional, due to its complicated structure, would lead to an algorithm which either does not converge or takes too much computer time .

The main purpose of this paper is to show that a simple discretization $`H_U[G]`$ exists and leads to a minimization algorithm which converges with a good efficiency. The application of this algorithm to the non-perturbative evaluation of gauge dependent Green functions like the quark and gluon propagator will clarify the gauge dependence of the fitted parameters.

In order to show the feasibility of our method, we have computed the gluon propagator at different values of the gauge parameter at small volumes, finding a sensitive gauge dependence.

The plan of the paper is as follows. In section 2 we review the method proposed to fix a generic covariant gauge on the lattice. In section 3 we describe a very simple and efficient discretization of the gauge fixing functional. In section 4 we give the details of our numerical simulations and report our main numerical results.

## 2 Covariant Gauge Fixing

In this section we set the notation and briefly formulate the covariant gauge fixing method in the continuum. In the following we will neglect the problem of Gribov copies, assuming that a gauge section in the space of gauge fields intersects all gauge orbits once and only once. In the Landau gauge the expectation value of a gauge dependent operator is given by

$$𝒪=\delta A_\mu \delta \eta \delta \overline{\eta }𝒪e^{S(A)S_{ghost}(\eta ,\overline{\eta },A)}\delta (_\mu A_\mu ),$$

and the gauge fixing condition can be enforced non-perturbatively by minimizing Gribov’s functional 

$$F_A[G]A^G^2=\text{Tr}\left(A_\mu ^GA_\mu ^G\right)d^4x.$$
(1)   
The Landau gauge is readily extended to a general covariant gauge-fixing condition of the form

$$_\mu A_\mu ^G(x)=\mathrm{\Lambda }(x),$$
(2)   
where $`\mathrm{\Lambda }(x)=\lambda ^a(x)\frac{T^a}{2}`$ belongs to the Lie algebra of the group and $`\text{Tr}(T^aT^b)=2\delta ^{ab}`$. Since gauge-invariant quantities are not sensitive to changes of gauge condition, it is possible to average over $`\mathrm{\Lambda }(x)`$ with a Gaussian weight

$$𝒪=\delta \mathrm{\Lambda }e^{\frac{1}{\alpha }{\scriptscriptstyle d^4xTr[\mathrm{\Lambda }^2]}}\delta A_\mu \delta \eta \delta \overline{\eta }𝒪e^{S(A)S_{ghost}(\eta ,\overline{\eta },A)}\delta (_\mu A_\mu \mathrm{\Lambda }),$$
(3)   
obtaining the standard formula

$$𝒪=\delta A_\mu \delta \eta \delta \overline{\eta }𝒪e^{S(A)S_{ghost}(\eta ,\overline{\eta },A)}e^{\frac{1}{\alpha }{\scriptscriptstyle d^4xTr[(_\mu A_\mu )^2]}}.$$

This formula is adopted as a definition for the expectation value of gauge dependent operators. The invariance under BRST transformations of the effective action in eq. (2) leads to the Slavnov-Taylor identities among different correlation functions. In particular the longitudinal part of the gluon propagator has to be equal to the free one, i.e. 

$$\frac{2}{N^21}\text{Tr}_\mu A_\mu (x)_\nu A_\nu (y)\alpha \delta (xy),$$
(4)   
where $`N`$ is the number of colors. Following the usual technique, the gauge-fixing condition (2) is obtained non perturbatively by minimizing the new functional 

$$H_A[G]d^4x\text{Tr}\left[(_\mu A_\mu ^G\mathrm{\Lambda })(_\nu A_\nu ^G\mathrm{\Lambda })\right],$$
(5)   
which obviously reaches its absolute minima ($`H_A[G]=0`$) when eq. (2) is satisfied. Therefore in this case the Gribov copies of the equation (2) are associated with different absolute minima of eq. (5). All the stationary points of $`H_A[G]`$ correspond to the following gauge condition

$$D_\nu _\nu (_\mu A_\mu ^G\mathrm{\Lambda })=0,$$
(6)   
where $`D_\nu `$ is the covariant derivative. Hence, due to the complexity of the new functional, its minimization could find “spurious” solutions which correspond to zero modes of the operator $`D_\nu _\nu `$ and do not satisfy the gauge condition in eq. (2). Of course the numerical minimization of the discretized version of eq. (5) can reach relative minima (spurious solutions) with $`H_A[G]0`$ which are not distinguishable from the absolute minima. Hence this could simulate the effect of an enlarged set of numerical Gribov copies. Within the limitations of this preliminary study we do not find any practical difference between the use of the new functional with respect to the standard Landau one (see results given in Fig. 2 below).

On the lattice, the expectation value of a gauge dependent operator $`𝒪`$ in a generic covariant gauge is

$$𝒪=\frac{1}{Z}𝑑\mathrm{\Lambda }e^{\frac{1}{\alpha }_xTr[\mathrm{\Lambda }^2]}𝑑U𝒪(U^{G_\alpha })e^{\beta S(U)},$$
(7)   
which is the straight-forward discretization of Eq. (3), where $`\mathrm{\Lambda }`$ is dimensionless on the lattice. $`S(U)`$ is the Wilson lattice gauge invariant action and $`G_\alpha `$ is the gauge transformation that minimizes the discretized version of the functional (5). On the lattice, the correct adjustment to the measure is included in eq. (7) by evaluating the operator over the gauge rotated links. Therefore it is not necessary to introduce ghost fields but it is mandatory to fulfill the gauge fixing condition numerically in order to get $`G_\alpha `$.

## 3 The Driven Discretization

In the Wilson discretization of gauge theories, the fundamental fields are the links $`U_\mu `$ which act as parallel transporters of the theory. Hence the lattice fields $`A_\mu `$ are derived quantities which tend to the continuum gluon field as the lattice spacing vanishes. As a consequence on the lattice it is possible to choose different definitions of $`A_\mu `$ formally equal up $`O(a)`$ terms. In quantum field theory this ambiguity is well understood because any pair of operators, differing from each other by irrelevant terms, will tend to the same continuum operator. This feature, checked in perturbation theory, has been verified numerically at the non-perturbative level in ref. , where it has been shown that different definitions of the gluon field give rise to Green’s functions proportional to each other, guaranteeing the uniqueness of the continuum gluon field.

The freedom to choose the lattice definition of $`A_\mu `$ can be used to build discretized functionals which lead to efficient gauge-fixing algorithms. In the standard Landau gauge fixing, for example, the discretization of the functional (1) is given by

$$F_U[G]=\frac{1}{VTa^2g^2}Tr\underset{x,\mu }{}\left[U_\mu ^G(x)+U_\mu ^G(x)2I\right],$$
(8)   
where $`V`$ is the 3-dimensional volume and $`T`$ the time size of the lattice. This formula corresponds, only up to $`O(a)`$ terms, to the naïve discretization of eq. (1) that is obtained from the standard lattice definition of the gluon field

$$A_\mu (x)\left[\frac{U_\mu (x)U_\mu ^{}(x)}{2iag}\right]_{Traceless}.$$
(9)   
On the other hand $`F_U[G]`$ has the important property that it depends only linearly on $`G(\overline{x})`$, when the iterative algorithm visits the lattice point $`\overline{x}`$. This feature would be spoiled if $`F_U[G]`$ were defined assuming literally the naïve discretization (9) in eq. (1).

In order to study the convergence of the algorithm, two quantities are usually monitored as a function of the number of iteration steps: $`F_U[G]`$ itself and

$$\theta _F=\frac{1}{VT}\underset{x}{}Tr[\mathrm{\Delta }_F\mathrm{\Delta }_F^{}],$$
(10)   
where

$`\mathrm{\Delta }_F(x)=\left[X_F(x)X_F^{}(x)\right]_{Traceless}{\displaystyle \frac{\delta F_U[G]}{\delta ϵ}},`$ (11)   
being $`G=e^{iϵ^a\frac{T^a}{2}}`$ and

$`X_F(x)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}\left(U_\mu (x)+U_\mu ^{}(x\mu )\right).`$ (12)   
$`\mathrm{\Delta }_F(x)`$ is proportional to the first derivative of $`F_U[G]`$ and reaches zero as the functional is extremized. Therefore the quality of the gauge fixing is determined by the parameter $`\theta _F`$ which corresponds to the continuum quantity $`d^4x\text{Tr}(_\mu A_\mu )^2`$.

A naïve discretization of $`H_A[G]`$ will generate a quadratic dependence on $`G`$, which could prevent the convergence of the algorithm. This obstacle has been overcome by taking advantage of the freedom to choose the gluon field definition. In fact, as in the Landau case, it has been possible to find a discretization of $`H_A[G]`$ (”driven discretization”) that depends linearly on $`G(\overline{x})`$ and corresponds, up to $`O(a)`$ terms, to the continuum limit, eq. (5). This aim can be reached by choosing each different term of $`H_U[G]`$ in order to guarantee the local linear dependence on $`G(\overline{x})`$, instead of being the algebraic consequence of a particular $`A_\mu `$ definition. We propose the following compact form of $`H_U[G]`$:

$$H_U[G]=\frac{1}{VTa^4g^2}Tr\underset{x}{}J^G(x)J^G(x),$$
(13)   
where

$`J(x)`$ $`=`$ $`N(x)ig\mathrm{\Lambda }(x),`$   
$`N(x)`$ $`=`$ $`8I+{\displaystyle \underset{\nu }{}}\left(U_\nu ^{}(x\nu )+U_\nu (x)\right).`$ (14)   
It is easy to see that locally $`H_U[G]`$ transforms linearly in $`G(\overline{x})`$ and its continuum limit is the functional (5). $`H_U[G]`$ is positive semidefinite and, unlike the Landau case, it is not invariant under global gauge transformations.   
The functional $`H_U[G]`$ can be minimized using the same numerical technique adopted in the Landau case. In order to study the convergence of the algorithm, two quantities can be monitored as a function of the number of iteration steps: $`H_U[G]`$ itself and

$$\theta _H=\frac{1}{VT}\underset{x}{}Tr[\mathrm{\Delta }_H\mathrm{\Delta }_H^{}],$$
(15)   
where

$`\mathrm{\Delta }_H(x)`$ $`=`$ $`\left[X_H(x)X_H^{}(x)\right]_{Traceless}{\displaystyle \frac{\delta H_U[G]}{\delta ϵ}}`$ (16)   
and

$`X_H(x)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}\left(U_\mu (x)J(x+\mu )+U_\mu ^{}(x\mu )J(x\mu )\right)`$ (17)   
$``$ $`8J(x)72I+igN(x)\mathrm{\Lambda }(x).`$   
$`\mathrm{\Delta }_H`$ is the driven discretization of the eq. (6); it is proportional to the first derivative of $`H_U[G]`$ and, analogously to the continuum, it is invariant under the transformations $`\mathrm{\Lambda }(x)\mathrm{\Lambda }(x)+C`$, where $`C`$ is a constant matrix belonging to the $`SU(3)`$ algebra. During the minimization process $`\theta _H`$ decreases to zero and $`H_U[G]`$ becomes constant. The quality of the convergence is measured by the final value of $`\theta _H`$.

## 4 Numerical Simulations and Results

We have generated 50 $`SU(3)`$ thermalized link configurations using the Wilson action with periodic boundary conditions at $`\beta =6.0`$ for $`8^4`$ and $`8^3\times 16`$ volumes. Following the prescription contained in eq. (7), for each Monte Carlo configuration and for each lattice site we have extracted a matrix $`\mathrm{\Lambda }(x)`$ according to a Gaussian distribution at a fixed $`\alpha `$ value. The gauge-fixing code implements an iterative overrelaxed minimization algorithm for $`F_U[G]`$ and $`H_U[G]`$. We have monitored the quantities $`F_U`$ and $`\theta _F`$ for the standard Landau gauge-fixing algorithm and $`H_U[G]`$ and $`\theta _H`$ for the new one after every lattice sweep. We have found that the value of the overrelaxing parameter $`\omega `$ adopted for the Landau algorithm $`\omega =1.72`$ is a good choice also for the new gauge fixing.

In Fig. 1 we report the values of $`\theta _H`$ for different $`\alpha `$ values as a function of the gauge fixing sweeps for a typical thermalized configuration at $`\beta =6.0`$ and $`VT=8^4`$.

Each sweep of the new algorithm takes $`15\%`$ more computer time with respect to the Landau case. For $`\alpha =0`$, which in this new procedure corresponds to the Landau gauge and is obtained by taking $`\mathrm{\Lambda }=0`$, the number of sweeps necessary to fix the gauge is $`1.2`$ times that of the standard algorithm. In the case of $`\alpha 0`$, the number of sweeps to minimize $`H_U[G]`$ with a given precision increases when $`\alpha `$ decreases.

Once the configuration have been rotated in a given gauge, we have tested our procedure computing the following two point correlation functions

$`𝒜_0𝒜_0(t)`$ $``$ $`{\displaystyle \frac{1}{V^2}}{\displaystyle \underset{𝐱,𝐲}{}}\text{Tr}A_0(𝐱,t)A_0(𝐲,0),`$ (18)   
$`𝒜_i𝒜_i(t)`$ $``$ $`{\displaystyle \frac{1}{3V^2}}{\displaystyle \underset{i}{}}{\displaystyle \underset{𝐱,𝐲}{}}\text{Tr}A_i(𝐱,t)A_i(𝐲,0),`$ (19)   
$`𝒜𝒜(t)`$ $``$ $`{\displaystyle \frac{1}{4V}}{\displaystyle \underset{\mu ,\nu }{}}{\displaystyle \underset{𝐱}{}}\text{Tr}_\mu A_\mu (𝐱,t)_\nu A_\nu (𝐱,0),`$ (20)   
where $`\mu `$ and $`\nu `$ run from 1 to 4, $`i`$ from 1 to 3 and the trace is over the color indices. In eqs. (18)-(20) the definition (9) is adopted for the gluon field on the lattice and $`_\mu `$ indicates the usual backward derivative. The correlators in eqs. (18) and (19) are relevant to the investigation of the QCD gluon sector, while we use the correlation in eq. (20) to check the Slavnov-Taylor identity for the longitudinal component of the gluon propagator. The statistical errors for the correlation functions have been estimated by the jacknife method. A thorough study of these operators will be presented in a forthcoming paper .

In Fig. 2 the Green functions computed with the standard Landau gauge fixing $`𝒜_i𝒜_i_F`$ and $`𝒜_i𝒜_i_H`$ evaluated with the new algorithm with $`\alpha =0`$, are reported.

The remarkable agreement shows that the new algorithm reproduces very well the results of the standard one within an overall scale factor close to 1. Moreover, the possible spurious solutions of the new gauge condition apparently do not affect this operator at least at $`\alpha =0`$.

Although in this paper we did not perform a systematic study of the Gribov copies for this gauge fixing procedure, a rough preliminary search shows, for $`\alpha =0`$, the same pattern of copies of the Landau case.

In Fig. 3 the data for the correlation (18) are reported. For $`\alpha =0`$ they show the well known flat behaviour forced by the Landau gauge.

For $`\alpha 0`$ the data show an almost flat curve with an enhancement at $`t=1`$ probably due to the effect of contact terms. In this case the flatness stems from the fact that $`d^3x_0A_00`$ because the average value of $`\mathrm{\Lambda }`$ on each time slice is negligible.

In Fig. 4 our results for the correlations $`𝒜_i𝒜_i_H`$ for different $`\alpha `$ values as a function of the time slice $`t`$ are reported for a set of $`50`$ SU(3) configurations at $`\beta =6.0`$ and $`VT=8^4`$.

Although in the case of small volume and a small number of configurations, the gauge dependence of the gluon propagator is clearly shown. For increasing $`\alpha `$ values, the time dependence of the gluon propagator becomes flatter since it approaches a limit where the gauge fixing effects disappear. It is also interesting to note that the $`\alpha `$ dependence of the gluon propagator shown in Fig. 4 does not seem to be re-absorbed by an overall scaling factor. This plot shows the feasibility of this procedure to study the gauge dependence of physically interesting correlators.

We complete the numerical checks of the gauge-fixing we are proposing, showing the numerical behaviour of a particular Slavnov-Taylor identity. This identity is a direct consequence of our way to implement the covariant gauge fixing and it must be understood as a numerical check of the procedure. This check measures the correlator of the operator $`𝒜`$ with itself, and if this operator has been well gauge-fixed, it is related to the autocorrelator of the $`\mathrm{\Lambda }`$’s and then it must be a delta function. In fact, for a finite volume $`VT`$ with periodic boundary conditions, the naïve Slavnov-Taylor identity gives for the correlator of eq. (20) the following relation:

$$𝒜𝒜(t)\alpha (\delta _{t,0}\frac{1}{VT})+𝒪(a).$$
(21)   
Eq. (21) shows that the shape of the correlation function (20) is expected to be proportional to a $`\delta `$-function up to the square of the $`_\mu A_\mu (x)`$ renormalization constant, apart from volume and $`𝒪(a)`$ effects. In order to verify eq. (21), in Fig. 5 we plot the correlation functions (20) versus the time $`t`$. We believe that the good agreement between these data and eq. (21) is an indication of the absence of possible spurious solutions induced by our method. The small signal at $`t=2,8`$ is likely due to the residual correlation induced by the enhancement in the point $`t=1`$.

## 5 Conclusions

In this paper we have described an efficient method to fix non-perturbatively a generic covariant gauge on the lattice. This procedure is equivalent to the Faddeev-Popov quantization for covariant gauges. Therefore it can be used to compare the lattice renormalized correlation functions with the same quantities computed in continuum perturbation theory.   
We have tested the algorithm on different SU(3) lattices for different values of the gauge parameter $`\alpha `$. Our numerical data show that the statistical fluctuations of measured quantities is comparable to the standard Landau gauge fixing algorithm.   
We have computed the correlation functions relevant for the investigations of the gluon propagator at a few $`\alpha `$ values and we have found a sensitive dependence on the gauge parameter.   
In collaboration with the Boston University Center for Computational Science, we are applying this method to the study of the gluon propagator on physical lattices.

## Acknowledgments

We warmly thank Massimo Testa for many fruitful discussions on this subject.

L.G. thanks Claudio Rebbi for interesting discussions on this subject.

We thank the Center for Computational Science of Boston University where part of this computation has been done.

L. G. was supported in part under DOE grant DE-FG02-91ER40676.

S. P. has been partially supported by the INFN-MIT “Bruno Rossi” Exchange Program while part of this paper was written. He thanks all the members of the MIT Center for Theoretical Physics for their warm hospitality.