BI-TP 99/02

April 1999

Determination of Anisotropy Coefficients

for SU(3) Gauge Actions from the

Integral and Matching Methods

J. Engels, F. Karsch and
T. Scheideler^{1}^{1}1Present address: IBM, P.O. Box 265, 00101 Helsinki, Finland.

Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany

Department of Physics, P.O.Box 9, 00014 University of Helsinki, Finland

Abstract

We use two non-perturbative methods to obtain the anisotropy derivatives of the coupling constants (the anisotropy coefficients) of lattice gauge theory. These coefficients appear in the derivative formulae for the energy density and the pressure. We calculate them for the standard Wilson and two improved actions, the and the Square Symanzik action. Both methods lead for all investigated actions to compatible results which are clearly different from their known asymptotic perturbative limits. With increasing the limits are however approached in all cases. Our findings strongly support the equivalence of the integral and derivative methods for the calculation of energy density and pressure.

## 1 Introduction

For the investigation of the phase transition from ordinary hadronic matter to the quark gluon plasma and the approach of the plasma to its high temperature limit the calculation of the energy density and the pressure is of central importance. These observables are defined as derivatives of the partition function with respect to temperature and volume

(1.1) |

In the lattice formulation of QCD the temperature and the volume therefore have to be independent variables. On a lattice with a fixed number of points in each spatial direction, and points in the temporal direction, this is achieved by choosing corresponding independent lattice spacings and so that

(1.2) |

As a consequence two different coupling constants have to be introduced for the spatial and temporal contributions to the action. For lattice gauge theory we consider actions of the following form

(1.3) |

where is a generalized plaquette in the plane. Generalized plaquette expectation values are then defined by

(1.4) | |||||

(1.5) |

The couplings may be rewritten in terms of the coupling constants as

(1.6) |

Here is the anisotropy. One may also express the couplings with the bare anisotropy and a common coupling through [1]

(1.7) |

or equivalently

(1.8) |

The couplings depend on and , or, alternatively, on and . Following ref. [2] we adopt the latter two as independent variables. The energy density and pressure are then

(1.9) |

Obviously, the energy density contains now terms, which are proportional to derivatives with respect to the anisotropy , because

(1.10) | |||||

(1.11) |

After taking the derivatives one can evaluate the expressions on isotropic lattices. At we have

(1.12) |

and the relation for the function [3]

(1.13) |

Energy density and pressure are now conveniently combined in the two equations

(1.14) | |||||

(1.15) |

Here denotes a generalized plaquette on a symmetric () lattice. Its inclusion normalizes both and to zero at .

We define the anisotropy coefficients as follows

(1.16) |

In the weak coupling limit, they are related to the derivatives of the anisotropic parts of the couplings

(1.17) |

by

(1.18) |

The quantities were calculated in ref.[3] for the Wilson action perturbatively. For the difference one obtains

(1.19) |

and, from the renormalization group equation and eq. (1.13)

(1.20) |

These perturbatively calculated anisotropy coefficients were used in early finite temperature lattice calculations. In gauge theory, however, this led to unphysical behaviour of the pressure: it could become negative and at the critical point a gap in was observed [4, 5]. A second, rather important consequence of taking only the perturbative coefficients is the rapid approach of and to their respective ideal gas limits soon after the transition, though the limits themselves are independent of . The actual form of the approach is of course dependent on the residual interaction in the quark gluon plasma, and therefore of physical relevance.

In ref. [6] an alternative method for the calculation of energy density and pressure was proposed, which is now in common use. Here the pressure is obtained without the anisotropy coefficients from an integral over plaquettes. In order to be consistent, the integral method should be equivalent to the use of non-perturbative anisotropy coefficients, i.e. the corresponding values for should be essentially independent on finite size and cut-off effects. First results in ref.[7] showed however considerable cut-off effects for the Wilson action. To clarify the situation, we started [8] a non-perturbative calculation of the coefficients with the matching method, which we want to complete in this paper.

In the following section we discuss the integral method in some more detail as in [7] and [8] and apply it to the Wilson action. The matching method is described and applied to the same action in section 3. For comparison we determine in section 4 the anisotropy coefficients as well for two improved actions, the and the Square Symanzik action [9]. We close with a summary and the conclusions.

## 2 The integral method

For large homogeneous systems with isotropic interactions, the volume derivative in the pressure formula (1.1) may be replaced by a volume division. The pressure and the free energy density are then related by

(2.1) |

Instead of calculating ln directly, one may take advantage of the relation

(2.2) |

and integrate the measured expectation value over a range. For we obtain then

(2.3) |

As in eqs. (1.14) and (1.15), the inclusion of the term serves to normalize at to zero. Obviously, the last equation will lead to a continuous pressure, even if there is a gap in the plaquettes at the critical point. The energy density is found from

(2.4) |

by just adding . Here we have replaced the sum of the coefficients by the function. This function and/or have to be known non-perturbatively, not only to evaluate the last equation, but also to determine the temperature on fixed lattices at any given coupling. Parametrizations of the function for the Wilson action have been derived by several groups [7, 10, 11]. They are rather similar above , but below this value they differ from each other. Probably this is due to the transition from the weak coupling to the strong coupling regime.

### 2.1 The high temperature limit

The two methods for the calculation of and on the lattice suffer from cut-off effects, especially for small values. Though these cut-off effects are of similar size, there is a non-negligible difference. For the same method, the high temperature limits of and have the same cut-off dependencies, because in both cases the same formula for is used and

(2.5) |

4 | 1.4952 | 1.3778 | 1.0852 |
---|---|---|---|

6 | 1.1816 | 1.1323 | 1.0435 |

8 | 1.0867 | 1.0659 | 1.0195 |

However, as different methods are based on the evaluation of different operators, which have different cut-off dependencies, the high limit derived in different approaches leads also to distinct cut-off dependencies. In the derivative method the high temperature limit is derived from the lowest order of the expansion of the energy density ()

(2.6) | |||||

(2.7) |

The quantities are the corresponding expansion coefficients of the generalized plaquette expectation values

(2.8) |

For the Wilson action one obtains the following deviation of the energy density, , in the high limit from the continuum result, [12]

(2.9) |

If the integral method is used, the calculation of the high temperature limit starts from eq.(2.1) and the lowest order term of ln. Here the subtraction has to be done explicitly. With the Wilson action one arrives at

(2.10) |

where of course . The full expansions can be calculated numerically. In Table 1 we list them for the two methods in the case of the Wilson action and . We show also the ratio between the two high temperature limits. We stress that this difference, which expresses itself in the ratio , will play a crucial role in the discussion of the anisotropy coefficients presented in the following sections.

### 2.2 Determination of the anisotropy coefficients

As mentioned already in the introduction, the integral method implies the use of non-perturbative anisotropy coefficients. We can determine them by comparing the pressure formulae from the two different approaches

(2.11) |

and by exploiting the function relation, eq.(1.13). In the pressure equation (2.11) we have inserted the factor as a compensation for the differing cut-off effects. In principle, . We approximate by its value at (Table 1). From eqs.(1.14) and (1.15) we obtain

(2.12) |

and therefore

(2.13) | |||||

(2.14) |

In the following we discuss the behaviour of . The behaviour of is similar, because both coeffcients are related by the smooth function. It is clear, that the numerical determination of the coefficients becomes problematic below the transition point, because there both the pressure and the plaquette differences are small. Indeed, this is seen in Fig. 2, where we show for the Wilson action as obtained from eq. (2.13) on lattices using interpolations of the and data of ref. [7]. Apart from the close vicinity of the critical

values the results for different are consistent with each other. For comparison we show also the previous result for from ref. [8], where the difference in the cut-off effects of the two methods was not taken into account, that is was used. Obviously the influence of the factor is important. Moreover it leads to slightly lower results for , though, as before, there remains a substantial deviation from the asymptotic value from eqs. (1.19) and (1.20).

It is difficult to estimate all possible error sources in the calculation of , as for example that of the parametrization of the function. The major errors are certainly due to the errors in the plaquette expectation values. In Fig. 2 we present the results including these errors, but omitting the respective critical regions. We can summarize these data with a simple Padé fit of the form

(2.15) |

where .

A further source of errors is the result for . Whereas the integration itself contributes only a negligible error, the size of the pressure at the lower integration boundary is unknown. We have identified with the point where the integrand is disappearing and assumed, as usual, that there is zero. This may be not true. In order to assess the effect of a non-zero pressure contribution we have repeated our calculation with a test value at for . As a result the value of is diminished by 0.022 in the whole range considered. The same shift in is observed also for and 8. For comparison we have plotted in Fig. 2 a fit of the form (2.15) to the shifted data as well. Here we find .

Some insight into the question of the pressure size below the transition is gained from the behaviours of the plaquette differences around in and . In Fig. 3 we compare them to each other. The main difference between the two gauge groups is evidently arising from the fact, that the respective phase transitions are of different order. Also, there is no doubt, that in the pressure is non-zero well below .

As we have seen in the inclusion of the cut-off factor considerably lowers the result for . In the values which were previously [14] determined

for with the integral method and are well in accord with points measured in an independent way at and by Ejiri et al. [15]. One has to check then, whether this agreement is lost if the correct is taken into account. We do this by using the new high precision data of ref. [13], already shown in Fig. 3, the function of ref. [14] and the value from Table 1. The resulting is plotted in Fig. 4. We find only a small change in this range, the result is lowered at by 0.0067 and at by 0.0142 . The reason for this is, that up to the difference is still large and only at higher couplings the pressure term with the factor is dominating in eq. (2.13) . We expect therefore a stronger decrease at higher values.

## 3 The matching method

This method was first used by Burgers et al. [1] to measure the anisotropy as a function of the bare anisotropy on anisotropic lattices. With this information one may as well determine the anisotropy coefficients. Consider the couplings

(3.1) |

and their derivatives with respect to the anisotropy

(3.2) | |||||

(3.3) |

At we may apply eq. (1.12) to find the difference

(3.4) |

and, using eq. (1.13), we obtain

(3.5) |

A measurement of the function in the neighbourhood of will therefore enable us to calculate the anisotropy coefficients, once we know the function. In principle, the anisotropy may be determined with two measurements of a physical observable which depends on a distance, which can be chosen in a spatial or in the temporal direction. At the same physical distance the expectation values have to be the same and thereby fix the anisotropy. This is the idea of matching.

### 3.1 Matching of Wilson loop ratios

Suitable quantities for the matching process are obtained from Wilson loops of size (in lattice units). The Wilson loops are related to the heavy quark potential

(3.6) |

for . Here, . The potential differs from the continuum potential by a term , the self-energy of the heavy quarks, which is dependent on lattice spacing, but independent of [16]. The natural way to proceed then is to build ratios of Wilson loops, which depend only on for large . We use the following ratios

(3.8) |

On an anisotropic lattice there are two different types of loops : space-space () and space-time () Wilson loops. The corresponding ratios and are measuring the same potential , if the matching condition for the physical distance

(3.9) |

is fulfilled, apart from an independent factor , which is due to the dependence of on and

(3.10) |

At the factor is of course . We do not implement any smearing for our Wilson loops, because then space and time links would have to be smeared by the same amount in physical units. However, in order to enhance the accuray of the Wilson loop measurements we are using the fast link integration technique of de Forcrand and Roiesnel [17] wherever this is feasible. It cannot be done for loops with an and at loop edges. In applying this algorithm we are saving in addition half of the computer time by taking advantage of the fact that the necessary link contour integrals are real and can as well be obtained from a half contour in the complex plane.

For the Wilson action simulations were performed at four values: 5.4, 5.7, 6.3 and 7.2 using a lattice for , a lattice for and a lattice for . The Wilson loops were calculated after every fourth sweep through the lattice, where a sweep consisted of one heatbath update succeeded by four overrelaxation steps. The integrated autocorrelation times which one finds then are increasing with increasing from 0.5, 1.5, 3.0 to over 3.0 for . On the average we took 5000-7000 measurements for and 600 for . Generally, for the same number of measurements, the relative error of a Wilson loop increases with the size of its area. The increase itself is much steeper at lower values than at higher ones. Larger area Wilson loops are therefore only included in the matching procedure for the higher values. In Fig. 5 we show a typical example of this behaviour. The difference in accuracy between link integrated and plain loops is clearly seen in the plot.

After measuring the Wilson loops at a fixed value of and , we compute the ratios and . For we connect the timelike ratios with spline interpolations. In order to improve the interpolations we include as well the ratios with . The spatial ratios are then shifted in by a factor and in height

by a factor such that the sum of the squared deviation () of the shifted values from the respective interpolations becomes minimal. For the roles of and are interchanged. The fitting is done for all ratios with and at the same time, with suitably chosen minimal values and . In so far our procedure differs from the matching prescription of Klassen [18], who matches single ratios and looks for a possible plateau of shift values. In Figs. 7 and 7 we show matching examples at and for the ratios and . The matching in Fig. 7 is optimized for , in Fig. 7 for . Both lead to the same value for the anisotropy. For we find and for , that is is always very close to one. One may as well perform fits, as Klassen[18] does, with fixed . The corresponding is then considerably larger, the value for increases slightly. We shall come back to this point again.

As has been mentioned already, we have included only ratios in the matching process which are built from Wilson loops with at least minimal extensions and . This has been done for two reasons. Once, should anyhow not be too small(see eq. (3.6)) and second, Wilson loops with a single link on one side cannot be link integrated and are therefore less accurate. That is why we disregard ratios with and/or in the final analysis. The influence of the chosen minimal values on the matching result for is demonstrated in Fig. 8 for at and . We observe, that apart from the data all other measurements for are consistent with each other.

### 3.2 Results for from matching

In order to obtain the anisotropy coefficients we measure now the function in the neighbourhood of with the aim of determining the derivative . In principle one would choose rather close to one. However, the difference between spatial and temporal Wilson loops becomes then very small as well and the matching will be inexact. Fortunately, it turns out, that is linear in in a wide range around at all values. A successful strategy is therefore to measure the Wilson loops with high precision at not too many values to determine the slope of . In Fig. 9 we show as an example the function for the Wilson action at as obtained from the matching of with , together with a linear fit to the data points.

As in the last section, we use now the function of ref. [7] to deduce the anisotropy coefficients from the measured derivative . The corresponding

results for are shown in Fig. 10 as a function of the minimal area of the Wilson loops, which were included in the ratio-matching. The data for are somewhat problematic in several respects and can only serve as an indication for the true value. As demonstrated already in Fig. 5 too large Wilson loops cannot be included in the matching for , because their errors are too large. The second handicap is the unknown or ambiguous function in the strong coupling region. We have therefore just taken the same value for at and . At we find that the results are essentially independent on the minimal area,

5.4 | 5.7 | 6.3 | 7.2 | |
---|---|---|---|---|

0.413(23) | 0.315(09) | 0.261(08) | 0.227(15) | |

-0.374(23) | -0.276(09) | -0.202(08) | -0.159(15) | |

0.397(17) | 0.348(06) | 0.296(06) | 0.261(12) | |

-0.358(17) | -0.309(06) | -0.237(06) | -0.193(12) |

at the errors are larger because of the lower statistics of the measurements. The final results for and are given in Table 2. We obtained them by fitting the single values to a constant with the method. Their average errors were estimated from the error estimate of the constant times the square root of the number of fitpoints. The error estimate was rescaled whenever . This was only necessary for . In order to be able to compare our results to those of Klassen [18], we have repeated the matching analysis under the assumption . Though the single matching results seem to have smaller errors then, a constant fit as before yields always a per degree of freedom which is larger than 2. The inclusion of small Wilson loops in the matching seems to be more crucial here and also to lead to higher values for . This is observed in Fig. 11, where we compare the two methods at .

Finally, we have gathered in Fig. 12 all available results for the Wilson action : those which we obtained from the integral method are represented by the fit, eq.(2.15), and the fit showing the possible influence of ; our matching results for and ; the data of Ejiri et al.[15] and those of Klassen[18], both for our function. We observe consistent behaviour of the integral and matching results for , those for are still within the error bars (see Fig. 2). There is also full agreement with the data of Ejiri et al. . Only the results of Klassen are definitely higher and insofar incompatible with all other measurements of .

## 4 Anisotropy coefficients for improved actions

Up to now we have discussed in detail the calculation of the anisotropy coefficients for the Wilson action only. In the following we consider the more general family of actions defined by Wilson loops of size up to 4 in the ()-plane

(4.1) | |||||

Here, the parameters and are constrained by the equation

(4.2) |

which ensures the correct continuum limit. The action is improved, if additionally

(4.3) |

Still, one parameter, for example , is free. For we obtain the Symanzik action (), for the action (). Alternatively, one may require that the propagator can be diagonalized. This facilitates the analytic calculation of the anisotropy coefficients in the weak coupling limit. In that case we have the condition

(4.4) |

with the free parameter . One can combine both objectives for a special set of parameters ()

(4.5) |

The corresponding action was introduced by García Pérez et al. [9] under the name Square Symanzik action. In ref. [19] the asymptotic anisotropy coefficients of this action have actually been calculated. For one obtains

(4.6) |

It is therefore obvious to investigate the properties of this action further and in particular to determine its anisotropy coefficients non-perturbatively. For comparison we have as well considered the somewhat simpler action, which is improved, but its perturbative anisotropy coefficients are unknown.

In order to apply the integral method we need again the high temperature limits of energy density and pressure. For the Square Symanzik action we have derived the expansions corresponding to eqs. (2.9) and (2.10)

By chance, the ratio is even improved. Since we intend to make use of Monte Carlo simulations on lattices we list in Table 3 the respective numerical results for the complete expansions of the ratios. The numbers for the action have been taken from ref. [12].

Action | |||
---|---|---|---|

1.088 | 0.99 | 1.099 | |

Square Sym. | 0.957 | 0.91 | 1.052 |

### 4.1 Results for the action

The thermodynamics of the action has already been investigated in detail by simulations on and lattices in ref. [12]. We can therefore immediately apply the integral method to compute the anisotropy coefficients using the corresponding plaquette expectation values and the function. For =4 the critical coupling was found to be . The behaviour of is similar to that shown in Fig. 2 for the Wilson action. With increasing is slowly decreasing. In the neighbourhood of the critical point the integral method leads again to numerical difficulties.

4.4 | 5.0 | 5.9 | |
---|---|---|---|

0.226(08) | 0.157(16) | 0.115(16) | |

-0.190(08) | -0.096(16) | -0.067(16) | |

0.236(03) | 0.195(07) | 0.158(13) | |

-0.200(03) | -0.134(07) | -0.110(13) |

In Fig. 13 we compare the integral method results with matching results, which were obtained on anisotropic lattices with 2000-4000 Wilson loop measurements for and 5.0 and 600 for . We note, that due to the inclusion of Wilson loops in the action, less links can be integrated to increase the accuracy of the measurements. In Table 4 we present the values for and from the matching procedure. Again we observe, that the matching with fixed leads to somewhat higher values. We find agreement between the different methods, but like for the case of the Wilson action only after taking into account the cut-off correction factor from Table 3.

### 4.2 Results for the Square Symanzik action

In order to determine the function for this action we have performed simulations on lattices and measured the plaquettes and the string tension as in ref. [20]. The lattice result and the physical string tension are related by . The function is then derived from

(4.9) |

In Fig. 14 we show our measured values for together with a fit from which we have deduced the numerical values of the function. On a lattice we have then measured the generalized plaquettes and as well as the Polyakov loop and its susceptibility in the range . On the average we took 4000-8000, close to the critical point 18000 measurements. From the peak of the susceptibility we estimate the following critical coupling for

(4.10) |

As a side product of our calculations we obtain the ratio

(4.11) |

which is in agreement with results found for other improved actions [21]. Using the above mentioned data and from Table 3 we have again applied the integral method to find the anisotropy coefficient . It is shown in Fig. 15 together with the asymptotic value , eq. (4.6), and the results of our matching analysis. The latter are also listed in Table 5. They were obtained from Wilson loop measurements on lattices. Like for the action less link integrations than in the Wilson action case can be done. Comparing the different results in Fig. 15 we find again the same general dependence as for the other actions. Also, a sizeable difference between the different matching options for is found: the results show an even stronger influence of the Wilson loop area as before, the results agree nicely with the integral method data. We note that in the Square Symanzik case is about half as big as in the Wilson case.

4.0 | 4.5 | 4.75 | 5.0 | 5.5 | |
---|---|---|---|---|---|

2000 | 800-1000 | 700 | 500-1500 | 200 | |

0.149(11) | 0.143(17) | 0.132(20) | 0.119(18) | 0.109(21) | |

-0.102(11) | -0.065(17) | -0.042(20) | -0.031(18) | -0.033(21) | |

0.175(03) | 0.158(06) | 0.155(09) | 0.157(09) | 0.120(13) | |

-0.128(03) | -0.080(06) | -0.065(09) | -0.069(09) | -0.044(13) |

## 5 Summary and conclusions

One of the major objectives of our paper was the verification of the equivalence of the integral and derivative methods for the calculation of thermodynamic quantities like energy density and pressure. The major obstacle in the comparison of the methods is the cut-off dependence on finite lattices. In the high temperature limit the deviations of energy density and pressure from their respective continuum limits due to the cut-off effects can be calculated. They manifest themselves as a dependence on in the thermodynamic limit, and they are different for the two methods.

The anisotropy coefficients as such must be finite size independent. Their calculation allows then for a check on the internal consistency of the integral method and of its equivalence to the derivative method. Moreover, its non-perturbative behaviour is of importance for the approach of thermodynamic quantities to the continuum limit.

By taking into account the major cut-off effects through the correction factor we have shown for the Wilson action that one obtains similar results for =4, 6 and 8 which agree, inside error bars, also with the matching results and with the results of Ejiri et al.[15]. The data of Klassen[18] are however definitely higher. This may be connected to the difference in our matching methods, the accuracy of the data and the lattice sizes used. We found it essential for the application of the matching procedure, that link integration of the Wilson loops was carried out, wherever this was possible, and that ratios instead of Wilson loops directly were used for the matching. Thus we avoided perimeter and corner dependencies. Remaining effects of this kind are probably responsible for the difference of to matching fits. The fact that was linear near helped in the determination of the derivative .

We have repeated our analysis for two improved actions, the and the Square Symanzik action, and find similar behaviours of the corresponding anisotropy coefficients. For all investigated actions the non-perturbative decrease with increasing and approach their respective weak coupling limits from above; the results from the matching and integral methods are in all cases compatible with each other.

In addition we determined for the Square Symanzik action the deconfinement transition point for , the ratio and the function in terms of the string tension. As an analytic result we derived the thermodynamic high temperature limits for this action.

Acknowledgements

We thank Maria García Pérez for providing us with details of her work in refs. [9] and [19]. The help of Burkhard Sturm (alias B. Beinlich) in the calculation of the determinant of the inverse propagator and the high temperature limits is gratefully acknowledged. The work was partly supported by the Deutsche Forschungsgemeinschaft under grants Pe 340/3-3 and Ka 1198/4-1, and the European Commission under the TMR-Network ”Finite Temperature Phase Transitions in Particle Physics”, ERBFMRX-CT-97-0122.

## References

- [1] G. Burgers, F. Karsch, A. Nakamura and I.O. Stamatescu, Nucl. Phys. B 304 (1988) 587.
- [2] J. Engels, F. Karsch, I. Montvay and H. Satz, Nucl. Phys. B 205[FS5] (1982) 545.
- [3] F. Karsch, Nucl. Phys. B 205 (1982) 285.
- [4] B. Svetitsky and F. Fucito, Phys. Lett. 131 B (1983) 165.
- [5] Y. Deng, Nucl. Phys. B (Proc. Suppl.) 9 (1989) 334; F.R. Brown, N.H. Christ, Y. Deng, M. Gao and T. Woch, Phys. Rev. Lett. 61 (1988) 2058.
- [6] J. Engels, J. Fingberg, F. Karsch, D. Miller and M. Weber, Phys. Lett. 252 B (1990) 625.
- [7] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lütgemeier and B. Petersson, Nucl. Phys. B 469 (1996) 419.
- [8] J. Engels, F. Karsch and T. Scheideler, Nucl. Phys. B (Proc. Suppl.) 63A-C (1998) 427.
- [9] M. García Pérez, J. Snippe and P. van Baal, 2nd Workshop on Continuous Advances in QCD, Minneapolis, March 1996, hep-lat/9607007.
- [10] K. Akemi et al. (QCDTARO Collaboration), Phys. Rev. Lett. 71 (1993) 3063.
- [11] R.G. Edwards, U.M. Heller and T.R. Klassen, Nucl. Phys. B 517 (1998) 377.
- [12] B. Beinlich, F. Karsch and E. Laermann, Nucl. Phys. B 462 (1996) 415.
- [13] J. Engels and T. Scheideler, Nucl. Phys. B 539 (1999) 557.
- [14] J. Engels, F. Karsch and K. Redlich, Nucl. Phys. B 435 (1995) 295.
- [15] S. Ejiri, Y. Iwasaki and K. Kanaya, Phys. Rev. D 58 (1998) 094505.
- [16] J.D. Stack, Phys. Rev. D 27 (1983) 412.
- [17] Ph. de Forcrand and C. Roiesnel, Phys. Lett. 151 B (1985) 77.
- [18] T.R. Klassen, Nucl. Phys. B 533 (1998) 557.
- [19] M. García Pérez and P. van Baal, Phys. Lett. B 392 (1997) 163.
- [20] B. Beinlich, C. Legeland, M. Lütgemeier, A. Peikert and T. Scheideler, Nucl. Phys. B (Proc. Suppl.) 63A-C (1998) 260.
- [21] B. Beinlich, F. Karsch, E. Laermann, and A. Peikert, Eur. Phys. J. C 6 (1999) 133.