Lattice gauge theory with fluctuating temperature

We study the possibility to implement the canonical Tsallis distribution for lattice field theory simulations. Formally, the application of the Tsallis distribution can be interpreted as introducing a fluctuating temperature. We give arguments for the approach and present our simulation method as well as our first numerical results in determining the equation of state for pure SU(2) lattice gauge fields.


Introduction
Lattice field theory, a systematic non-perturbative method has been successfully applied for studying the special features of strong interactions, such as quark confinement, the deconfining phase transition and, in general, the thermodynamics of strongly interacting matter [1][2][3][4][5].Usually, lattice simulations of Quantum Chromodynamics (QCD) are based on a conventional canonical ensemble approach where exponential distribution plays an important role.However, experimental particle spectra are not purely exponential.In heavy-ion collisions both exponential and power-law regions have been observed in the measured pion, kaon, etc. spectra [6][7][8][9][10][11][12][13][14].The experimentally measured hadron spectra may reflect statistical properties of the observed system.This suggests that the strongly interacting matter in extreme conditions may show non-conventional thermodynamical behaviour which manifests itself in non-conventional distributions.Such distribution with a power-law like tail is the so-called Tsallis distribution which may be derived from a non-extensive energy addition rule [15].This distribution has been widely used recently in many areas where statistical models have been applied.For instance, it leads naturally to the mass spectrum of hadrons when they are combined from massless partons with Tsallis-like energy distributions, and results in a Hagedorn-type limiting temperature [16].Fits using the Tsallis distribution to the observed particle spectra in heavy-ion reactions have also been remarkably successful [17,18].Non-conventional distributions are based in general on nonadditive composition rules [19], or on a non-conventional entropy formula (which replaces the Boltzmann entropy).Such an entropy formula is the Tsallis entropy [15], which is not an extensive quantity.One can find a monotonic function of it, however, which is proven to be extensive.Non-extensive thermodynamics can be viewed as an efa e-mail: tsbiro@rmki.kfki.hub e-mail: schram@phys.unideb.hufective theory for non-equilibrium phenomena and for systems where long-range interactions and long-range correlations are present.Applying it to describe the thermodinamical properties of the strongly interacting matter under extreme conditions is particularly interesting.
If E denotes the energy of a state of the system, the Tsallis distribution is given with the probability distribution (β is the inverse temperature).In the c → ∞ limit it restores the Gibbs factor: The quantity q = 1 + 1/c is the Tsallis index.It is easy to show that the following identity holds: In (3) one can recognize the Gamma distribution: The identity ( 3) is remarkable: it shows that averages with the Tsallis distribution can be evaluated as averaging over different β valued Gibbs expectation values.This is a specific example of the superstatistical approach [20], where the inverse temperature is not a simple quantity: it follows a Gamma distribution.
the lattice extension in time direction, a t is the timelike lattice spacing), the lattice fields satisfying periodic boundary conditions.In a usual lattice simulation the timelike and spacelike lattice spacings equal, and the volume of the system is V = N s 3 a s 3 (N s is the lattice size in space direction, a s is the spacelike lattice spacing).In our approach, so as to realize a fluctuating temperature, we use an anisotropic lattice where the fluctuation of the inverse temperature is simulated by a Gamma distributed anisotropy parameter, t = a t /a s .The mean value of t is 1, the finite width, 1/ √ c, is a parameter in this approach which can be taken from the observed particle spectra.In the following we apply the method for SU(2) gauge fields.The partition function in the superstatistical approach is given by where the SU(2) lattice action is The contribution of spacelike and timelike plaquettes, P s and P t , in a generic form is the following with U plaq being the ordered product of the SU(2) link variables around the elementary spacelike or timelike plaquette.Below we will also use a shorthand notation of the action, S (t, U) = a(U)t + b(U)/t = at + b/t, where a and b denotes the spacelike and timelike contribution to the action, respectively.The expectation value of an operator which includes a timelike link with power ν, Â = t ν A can be calculated as The evaluation of the above integral can be achieved in several ways.By integrating over the asymmetry parameter t we can derive an effective action: where The integration can be evaluated analytically giving Asymptotically, for any finite c the action S e f f → 2

√
ab as a and b becomes large.This assures that the above integral is well defined in the thermodynamical limit.On the other  2).Of course, this follows also from the fact that the Gamma distribution becomes a δ-function in the c → ∞ limit.In this report, instead of using the effective action, we calculate the expectation values via a direct numerical simulation.We outline the method in the next section.

Numerical Approach
In the direct numerical approach we start with the usual canonical Gibbs simulation, in our instance using a Metropolis update algorithm.Of course we have to modify the algorithm according to our previous discussion.Now we have an additional variable, the asymmetry, which has to be generated according to a Gamma distribution.Equation (3) assures us that if we cover the parameter space sufficiently then we may obtain the desired results.The characteristic values of the distribution are determined via the Tsallis parameter c.
We performed simulations for different c values: for c = 5.5, 13.5, 32.0 and for c = 1024.0(see Fig. 1).The lower values correspond to parameters derived from highenergy collision experiments, c = 1024 is meant to approximate the c → ∞ (Gibbs) limit, and c = 32 is chosen as an intermediate value to show the transition.For a chosen c we generate numerically random asymmetry values of t.Distributions of the realized asymmetry parameter values for different couplings at c = 5.5 can be seen in Fig. 2. The upper figure shows that for a smooth reconstruction of the Euler-Gamma distribution one should generate random values in the order of 10 5 (200000 is shown in our example).Our actual numerical simulations were performed for 1000 random values of the asymmetry parameter at each c value and at each coupling.It is obvious that about 1000 random values cannot follow precisely the required Gamma distributions.In fact, there are apparent fluctuations.However, in average the values cover the Gamma distributions quite well and Euler-Gamma fits to the 1000 generated numbers are rather good, reproducing the initial c parameters of the distributions.For example, in a simulation with c = 13.5, the fitted value of the shape parameter is 13.379 (instead of 13.5) and the rate parameter is 0.075243 (which, in our case should be compared to 1/c = 0, 074074).Nevertheless, longer simulation runs would be preferable, but these should be limited to reasonable lengths of computing time.
In general, once an asymmetry parameter is determined we can apply the standard Metropolis sweep to update the gauge fields.Then a new asymmetry parameter is thrown and so on.The ensemble averages are calculated with averaging the gauge configurations over the generated asymmetry values, just as in the canonical Monte-Carlo simula-Fig.3. The ratio at : b/t and the square of the asymmetry parameter t 2 (green crosses) as a function of the Monte-Carlo steps.The numerical average of t 2 and the theoretical expectation value is also shown.tions.Naturally, in this numerical approach several technical questions arise, partly connected to the sequential organization of the algorithm.For example one can generate a new asymmetry parameter t at each Monte-Carlo sweep over the lattice or perform several sweeps before choosing a new one.Going further, it is also possible to throw a new asymmetry value at each link update meaning local temperature fluctuations in space, or just to use the same t for the whole lattice, which would correspond to global thermal fluctuations.In the following we summarize our choices.
We perform simulations for the pure SU(2) gauge theory in 3+1 dimensions on 10 3 x2 and 10 4 lattices.We use the action given in equation (6).After heating up the initial configurations (performing 1000 Monte Carlo sweeps through the lattice) we throw a new asymmetry parameter for the whole lattice before each update sweeps to measure the observables.
In order to get more information about the reliability of the algorithm we investigated several quantities: the average spacelike and timelike plaquettes (containing the asymmetry factors t and 1/t, respectively), their differences, ratios etc. both as time series of Monte-Carlo sweeps and as functions of the coupling.Based on these studies we conclude that our algorithm is stable and able to provide the necessary quantities.As an illustration, in Fig. 3 we plotted the ratio of at and b/t compared to the t values in a simulation with 10000 sweeps on a 10 4 lattice.The average values of t 2 are also indicated, both the numerical one t 2 , and the theoretical average which is 1.17902.The ratio fluctuates approximately around the theoretical average, as it should, while the fluctuations are large.It is even more remarkable how large fluctuations the t parameter has in this, c = 5.5 case.One should note, however, that these fluctuations are inevitable for small c values.In these cases the Euler-Gamma distribution is wide, temperature values larger than double of the average or less than half of it give also important contributions to the expectation values.

Equation of State
To determine the equation of state for the superstatistical SU(2) Yang-Mills lattice gauge theory we calculate the energy density and pressure similar way as it is done in the conventional case.There also an asymmetry parameter is used to perform the volume and temperature differentiations independently.However, in our calculation the asymmetry parameter t = a t /a s cannot be set to unity after obtaining the derivatives.Rather, we have to write it in the expressions explicitely.
The pressure and energy density is obtained from the partition function (5) according to standard thermodynamical rules with β = 1/T inverse temperature [2,4]: It is well known that for an ideal gas p = e/3.Therefore it is particularly interesting to consider the interaction measure, From the above equations we obtain These formulae are extensions of the standard expressions used so far in lattice gauge theory in Gibbs-Boltzmann thermodynamics [22].In our calculations here, for the energy density (and pressure) we consider the first term only, so we can directly compare the values to the results of the pioneering work of [23].For the interaction measure, as the first terms cancel, we use the perturbative result in leading order for the derivative of the coupling.There is one more important issue regarding the above equations.Namely, (15) and ( 16) give significantly non-zero values even if we consider a system at zero temperature, i.e. if we use a symmetric, N 4 lattice at small inverse couplings (10 4 lattice in our case).This is a consequence of the fluctuating asymmetry parameter which causes that P s t P t /t even in the space-time symmetric, i.e.O(4) symmetric case.Naturally, this effect is remarkable for the physically interesting c values and disappears as c → ∞ (see Fig. 4).For this reason we use a simple renormalization method where we calculate the quantities on an asymmetric lattice (10 3 x2 lattice) first and subtract from them the corresponding quantities obtained on the symmetric, 10 4 lattice.We do the same for the interaction measure as well, which is given by the following expression: ∂ ∂a 4 g 2 tP s + P t /t asym − tP s + P t /t sym (17) Our numerical results for the energy density and for the interaction measure in the above approximation are shown in Fig. 5.One can see a slight shift of the curves to higher 4/g 2 values as c is getting smaller.Of course it is not easy to draw a precise conclusion about the consequences of that shift at this stage of the simulations.One needs at least more statistics and it would also be desirable to investigate the neglected terms in the expressions for the thermodynamical observables.Values at larger 4/g 2 should also be taken with caution: at inverse couplings, let say larger than about 3, the symmetric lattice does not represent any more a system at zero temperature.The lattice spacing at these couplings is so small that the lattice ceases to represent a real thermodynamical system.However, it is clear that even after the subtraction of the symmetric contributions there is an apparent effect due to the applied nonconventional distribution.
In the weak coupling limit, which means high temperature when N t is kept fixed, the energy density e should approach the value for the free gluon gas, e ∼ T 4 .This is the well known Stefan-Boltzmann relation.In conventional lattice calculations, for the energy density expressed in lattice units it means that ea 4 ∼ 1/N 4 t should be constant for fixed N t .Therefore it is usual (and interesting) to consider values of e/T 4 in lattice simulations.This raises immediately an important question in our non-conventional approach.How should one interpret and calculate this quantity when the temperature fluctuates?As we evaluate expectation values on the lattice, we might consider possibilities like e/ T 4 , e/ T 4 or e/T 4 .These expressions give significantly different results in our cases because the temperature is given via the asymmetry parameter t as T /t and the expectation values of different powers of t depend strongly on the Tsallis parameter c (see Fig. 6).Estimations for free massless Tsallis gas show that 05004-p.4According to these considerations we plot our numerical results for (e − 3p)/T 4 and for e/T 4 in Fig. 7.The t −3 c−1 values for different c parameters are shown in Table 1.Similar ratios (at least approximately) can be deduced from the data given in Fig. 7. Nevertheless, besides the lack of sufficient statistics (which we mentioned earlier), simulations on larger lattices are needed for exploring the weak coupling, high temperature limit and the Stefan-Boltzmann relation for the SU(2) Tsallis gas.
Alltogether, our first results seem promising.The superstatistical method we use to evaluate expectation values with Tsallis distribution on the lattice is stable and functioning.However, there are still some methodical questions we should address.One of them is the problem of thermalization: in the calculations presented here we have choosen a new random asymmetry parameter at each Monte-Carlo sweeps.Whether it is enough to thermalize the system or it is necessary to perform several sweeps at each asymmetry parameter, remains a question.Work on this and similar problems, as well as on improving the statistics of our calculations is in progress.

Fig. 1 .
Fig. 1.Gamma distributions for the values of the Tsallis parameter c which were used in our numerical simulations.

Fig. 2 .
Fig. 2. 200000 asymmetry parameters randomly generated via our algorithm compared to the Euler-Gamma distribution (upper figure) and the realized distributions of the t asymmetry in actual simulations at some values of the coupling (lower figure).The applied Tsallis parameter was c = 5.5.

Fig. 4 .
Fig. 4. The action difference, i.e. the first term of expression (15) on a symmetric, 10 4 lattice as a function of the coupling for different c values.For c = 1024 the values are practically zero as expected.

Fig. 5 .
Fig. 5.The energy density and the interaction measure in lattice units vs. coupling for the Tsallis parameters c = 5.5, 13.5, 32.0, and 1024.0.

Fig. 6 .
Fig. 6.The average values of the different powers of the asymmetry parameter at different couplings in our Monte-Carlo simulations for the Tsallis parameters c = 5.5, 13.5, 32.0, and 1024.0.