Particle yields from numerical simulations

In this paper we use numerical simulations to calculate the particle yields. We demonstrate that in the model of local particle creation the deviation from the pure exponential distribution is natural even in equilibrium, and an approximate Tsallis-Pareto-like distribution function can be well fitted to the calculated yields, in accordance with the experimental observations. We present numerical simulations in classical $\Phi^4$ model as well as in the SU(3) quantum Yang-Mills theory to clarify this issue.


I. INTRODUCTION
In collider experiments the observed hadron yields are surprisingly far from the Boltzmann distribution expected from black body radiation of the hot plasma. It is true for all yields coming from hadron collisions, starting from p-p collisions in ALICE, CMS, STAR and PHENIX collaborations, respectively [1][2][3][4]. Part of these particle spectra can be explained by perturbative QCD calculations [5], but for a complete description of all of the fits a Tsallis-Pareto [6] like Ansatz is necessary. These fits have just one additional parameter compared to a Boltzmannian, yet they work very nicely in experimental fits [2,[7][8][9] although in certain cases the soft and hard physics has to be treated separately [10].
There are several interpretation of these results. A QCDbased generalized Ansatz [11][12][13] results in distribution functions that are similar to Tsallis distribution. More natural explanation is to assume some collective behaviour that leads to these type of distributions. The reason for the deviation from the Boltzmann distribution can be the finite volume [14], fluctuating temperature [15,16]. In fact the leading deviation from Boltzmann distribution is of Tsallis form [17]. It is also possible that the event-by-event distribitions are Boltzmannian, but the hadron multiplicities fluctuate according to negative binomial distribution, and only in the cumulative yields do we see Tsallis distribution [10].
In this paper we suggest another possible mechanism to observe non-Boltzmann distribution. It is possible, namely, that the energy levels of the system are Boltzmann distributed, still the observed particle yields follow a different distribution. To understand this possibility, we have to rethink the mechanism, how particles emerge from a strongly interacting plasma.
In weakly interacting gases, like the photon gas that interacts only with the wall of the cavity, the particle states are (almost) the same as the energy eigenstates: they have definite momentum and energy, and they extend to the whole cavity. If a photon escapes from the system it carries information about the distribution of the energy eigenstates, and so, correspond- * Electronic address: homor.marietta.m@gmail.com † Electronic address: jakovac@caesar.elte.hu ingly, we obtain a photon yield that is distributed according to the Bose-Einstein distribution.
In strongly interacting plasmas, however, the situation changes. Particle states are no longer energy eigenstates, in fact they consist of a lot of energy eigenstates with a Lorentzian envelope: they are quasiparticles. Moreover, they are basically local objects, they do not extend to the whole plasma. The phenomenon of jet quenching [18] clearly indicates that high energy particles are created locally, in a volume of at most of order one fm. As a consequence, if they escape the plasma, they do not carry information about the occupation of the global energy levels, but about the local energy density. The fluctuation of the energy in a small volume, however, is different than in a large subsystem, it need not (and actually do not) follow Boltzmann distribution.
According to this picture to assess the particle yields coming from a strongly interacting plasma we have to measure the distribution of the local energy density. This is a measurable quantity even in numerical simulations, making possible to give predictions on the observable yields which are otherwise very hardly accessible quantities.
In this paper we have considered two models basically to demonstrate the method: a classical Φ 4 model and a quantum SU(3) Yang-Mills gauge theory. As it turns out, in both cases the distribution of the local energy density stabilizes relatively fast, well before the actual thermal equilibration, and a Tsallis distribution is an excellent fit to them. In the process of approaching equilibrium, the parameters of the Tsallis distribution change. In classical scalar model case it is possible to follow real time evolution, and so the variation of the temperature and Tsallis parameter in real time. In the followings we discuss the histogram method to determine the local energy distribution function, consider the classical quartic model and the quantum Yang-Mills model, and finally we close the paper with a Conclusion section.

II. LOCAL ENERGY DENSITY DISTRIBUTION
As stated in the introduction, our aim is to determine the local energy density distribution. In this section, we recall some basic definitions of probability theory.
First of all, let X be a stochastic variable. The indicator of X being in the [x, x + ∆x] interval is I [x,x+∆x] (X) = Θ(X− arXiv:1709.00230v1 [hep-ph] 1 Sep 2017 x)Θ(x + dx − X), where Θ is the Heaviside function. The expectation value of this indicator equals the probability of X being in the interval: In statistical approach we take the expectation value above some configuration space, then we assume that X is a function(al) of the configurations X(A). If the distribution of the configurations A is given, let us denote it by f (A), then the expectation value of a general quantity R(A) can be computed as In canonical ensemble f (A) = exp(−βH(A))/Z. In practice, however, we generate a lot of configurations according to the distribution function f (A), either by solving the equation of motion (classical theory), by following a Markov-process (Monte Carlo simulations) or considering small subsystems of a configuration in a large volume. In any case we have A 1 , A 2 , . . . A n configurations and we take the expectation value by summing above them: Therefore the expectation value of the indicator is proportional to the number of configurations n i , where the X(A i ) quantity have value between x and x + ∆x. Therefore To determine the complete distribution function, therefore, we divide the possible outputs of X to bins, each of them is ∆x wide. Then we scan over all available configurations, and each configuration contributes to the bin that contains X(A i ). This provides a histogram that is exactly the desired P(X ∈ [x, x+ ∆x]) distribution function. The limit ∆x → 0 provides the probability density of the variable X: If we write (1) into the definition (5) of f (x), then we get a Dirac-δ approximation and we can write rather formally: In this work, the quantity in question -in other words the stochastic variable -is the local energy density x where x is an arbitrary space-time coordinate. With this (6) becomes: It is important to note that the local energy density, or the energy in a small volume need not to follow Boltzmann distribution, even if we use canonical ensemble to generate configurations. The reason is that in a small volume the energy contains a considerable amount of surface energy, too. The surface energy depends both on the state of the singled out volume and on its environment, and so we cannot associate it to any of them. Only in largish volumes which are small compared to the complete system, but large enough to neglect the surface energy terms, can we deduce that the probability density of measuring a given energy value is Boltzmanndistributed.

A. Local energy distribution in free systems
It is worth to think about the form of the local energy distribution when it is Boltzmann-like. So let us assume that p( ) = N e −c with some constants. We have two constraints: these fix the constants to be Therefore we do not expect e −β form, only if x = T . In simple systems we in fact obtain such a form. Most simply, in a system built up from local independent systems we have for all configurations σ In this case because at sites j = i the corresponding Z j factors drop out.
We also have the same results in case of free systems, even when the energy is the sum of the momentum states, while the local energy density is localized in real space. To prove this statement we first rewrite the energy density distribution as where we assumed x > 0 for all configurations. Then we expand the exponential To avoid UV divergences we renormalize the above expression taking the normal ordered product, so we calculate The local energy density can be defined in a number of ways, each definitions differ from each other in total divergences. We will choose a simple representation, where the local energy density can be written as function of the creationannihilation operators a p and a † p in d dimensions as It is simple to see that Now we can compute the expectation value in question at x = 0: where E n = n| : H : |n . Since the same state stands in the left and right hand side of the expectation value, we must have the same number of creation and annihilation operators for each momenta. We will omit the possibility that more than two operators have the same momenta, since these contributions are suppressed by factors of V d the volume of the d dimensional space. This means that we have to make pairs (Wick theorem). It is easy to see that all pairings give the same contribution, so finally we have (18) Substituting back this result we see that the expectation value of the -times local energy density is proportional to the expectation value of the local energy density to the th power: Therefore and so the inverse Fourier transform yields: which means that in the free systems the local energy density is indeed Boltzmann-distributed. We see from this calculation that the validity of the Boltzmann distribution depends on very sensitive details, for example that the expectation value of powers of the local energy density is proportional to powers of the expectation value of the local energy density (cf. eq.19). In a general theory it will not be true anymore, resulting that e −iλ x is not a simple pole and then p( ) is no longer exponential. We can not determine the actual form, but based on very general arguments [17] we expect that if the deviation is small, then it must be a Tsallis-Pareto distribution.
In the following sections, we consider the real time simulation of the classical Φ 4 theory in 3 dimensions and perform a standard Monte Carlo simulation with heat-bath algorithm for the Euclidean SU(3) gauge theory.

III. A TOY MODEL: CLASSICAL Φ 4 THEORY
Our first toy model is the well-known classical Φ 4 theory. One of the advantages of classical theories is that we can perform real-time simulations by successively solving the canonical equations and we can calculate physical quantities that are hardly accessible in other methods. Classical theories are used to approach the full theory in a lot of contexts [19][20][21][22][23][24][25][26][27][28].
The discretized version has the Hamiltonian [28] where U denotes the discretization mesh (in our case a cubic lattice with N sites in all directions, N = 40, 50), and x is the local energy density In this expression we have to use the discretized gradient where a is the discretization spacing, and e i is the unit vector pointing to the ith direction. The corresponding equations of motion reaḋ is the discretized Laplacian. The continuous equation of motion preserves energy, but in the time discretized version the energy conservation depends on the algorithm. We used leapfrog and Runge-Kutta methods; for further discussion cf. [28].
We note that the system can be rescaled as t → t/a, Φ → √ λ aΦ, Π → √ λ a 2 Π, then we have the same equations of motion with a = 1, m → am and λ = 1. This means that the value of λ does not modify the classical dynamics, it can be compensated by the normalization of the fields. The energy density rescales as → λa 4 . In the simulations we have used the a = 1 unit, but we have kept the value of λ to test the numerical effects.
After thermalization we can use the thermodynamical notions. The temperature (T ) of the system is defined as [28]: where N 3 is the number of lattice sites and Π k is the Fouriertransformed momentum field. We use this formula to check whether the system reached thermal equilibrium by verifying that |Π k | 2 is independent of k (equipartition). It turned out, that we can distinguish two time scales, as higher modes thermalise much faster than low ones. After 10000 time steps, the system can be considered fully thermalised. We have solved the classical EoM on 40 3 and 50 3 lattices. As initial conditions we have chosen Φ x = 0 for every point, and we have assumed pointwise independent distributions for the canonical momenta. In one case the momentum distribution was a uniform distribution in the [0, 1] range, in the other case we had a 1/ cosh(x) distribution. In the simulation we have chosen dt = 0.1 time step in lattice spacing units.
After starting the simulation, the distribution of the local energy density very quickly stabilizes. Already after the 17th time step the histograms reach the characteristic form which remained true in all later times. We can see these distributions in figures 1 and 2. Each data point in the histogram is the average of 50 runs with the same initial condition and the errorbars represent the standard error of the mean (SEM). It is apparent that the distribution deviates from the Boltzmannian exponential form, but a Tsallis-Pareto distribution proved to be a very good fit Note that for q → 1 it gives back the Boltzmann distribution.
The actual values, q = 0.974 and q = 1.094 respectively, are very close to the Boltzmannian case, and so it can be revealed only by a thorough analysis with at least 10 6 independent data points.
In the 17th time step the system is very far from equilibrium, but the Tsallis-Pareto distribution of the energy density remained true, with time dependent Tsallis parameter q(t). This function is plotted in Figures 3 and 4 for two different time intervals, starting with different total energy (corersponding to different temperatures after thermalization and different lattice sizes). We can observe that even in that cases when q started from a value smaller than 1, finally in all runs it reached a value that is consistently larger than 1. It is interesting that this value seems to be independent on the temperature as well as on the lattice sizes we studied. The Tsallis-parameter takes its equilibrium value already in the pre-thermalised state (i.e. when only higher modes are thermalised) within error. The actual value of the equilibrium Tsallis parameter is q = 1.024 is in the order of the experimental values.   At this stage it is not clear, whether the deviation from the Boltzmann distribution is a property of the energy density only, or the microscopic distribution function f (A) over the configurations (cf. (2)) is also non-Boltzmannian. To this end we have carefully studied the distribution of Π 2 x for different energy values and different lattice sizes with large statistics. If the thermal ensemble has a distribution function e −βH , then the x = Π 2 x /2 values for a given x must follow e −βx expontential distribution, i.e. we must have q = 1 in the Tsallis fit. The result for the averaged data   is q momentum = 0.999 ± 0.001 for all the cases we studied. This means that the distribution over the configurations is in fact the standard Boltzmann distribution, and only the energy density has a Tsallis-like distribution. Learning this fact we could proceed to the quantum field theory case, where only the equilibrium can be studied with Monte Carlo methods, but the distribution of the energy density is expected to be Tsallis-like even in this case.

IV. EUCLIDEAN SU(3) PURE GAUGE THEORY
Our second model is the quantum SU(3) Yang-Mills theory. The action of the model in Euclidean formalism is the following: where F µν (x) = −igF a µν (x)T a is the gluon field strength tensor, T a are the generators of the Lie-algebra and g is the coupling constant.
We use the Wilson-action for the lattice formulation of the theory: where U (x+µ, x) is the parallel transporter from lattice coordinate x to x + µ and the plaquette variable corresponding to x is the product of four parallel transporters along the closed It is well-known that the Wilson-action corresponds to the continuum theory if one chooses β = 2N/g 2 (in our case N = 3) and the connection between the parallel transporters and the gauge fields is U (x, µ) = e −aAµ(x) . Continuum limit is reached as β → ∞.
The local energy density is now given by the plaquette energy: We determine the histogram of the local energy density in the same way as we have done in the classical theory, picking out independent configurations from the thermal ensemble. An advantage of the histogram method is that renormalization can be explicitly traced in the distributions. In case of the local energy density, being a composite operator, we expect a multiplicative renormalization as well as an eventual mixing with the unit operator (additive renormalization). In the histograms the two types of renormalization show up as a dilatation and a position shift. Neither of these effects modify the power of the high-energy tail: therefore the Tsallis parameter is not renormalized.
We use Monte Carlo simulation with the well-known heatbath algorithm to determine the distribution of with zero energy initial condition (i.e. p = 0 for all plaquettes) at lattices N t × N 3 s . For the fits of the resulting distributions we have taken into account the reduction of the phase space at low energy densities. As a result, we modify our previous fit function (26) by a factor of x n where n ≈ 3: The error of the histogram is assumed to be a Gaussian (meaning √ N standard deviation for a bin which contains N points). However, in more realistic models it would be advised to perform a maximum likelihood parameter estimation (in the context of Tsallis distributions cf. [29]).

A. Plaquette energy histogram
A typical plaquette energy histogram is presented in Fig. 5. In the figure the best Tsallis (30) and Boltzmannian g(x) = Ax N e −Bx fits can also be seen. It is evident that the Tsallis fit is better than the Boltzmannian one, similarly as in the case of classical Φ 4 theory.  The evolution of the Tsallis parameter (with its fit error) is presented in Fig. 6 as a function of the Monte Carlo time. We covered a wide range of β, starting from β = 6 to as high values as β = 20. Note, that the plaquette energy has an upper bound (due to the properties of the trace of SU(3) matrices), and this severely distorts the histogram below β = 6.
It is interesting that the shape of the distribution shows up very early, well before thermalization: already after the first MC sweep we find a roughly Tsallis-like distribution, although the Tsallis fit converges well after about τ = 10 (first vertical line in Fig. 6). Thermalization time is τ ≈ 30 (second vertical line in Fig. 6). We have calculated the equilibrium Tsallis parameter from configurations after τ = 100 (third vertical line) for each simulation.

C. Tsallis q for various Nt
To connect the numerical observations to physical units, we have to fix the scale. We have chosen the Sommer-scale [30] and interpolated the β and lattice constant a relation based on the data from [31][32][33]. The dimensionless temperature is as follows: r 0 T = 1 Nt r0 a , where r 0 ≈ 0.5 fm is the Sommerscale parameter.
The temperature dependence of the q parameter is shown in Fig. 7. Five different N t value is considered. To check the thermodynamic limit, we repeated the simulation for N s = 60, N s = 50 and N s = 40.

D. Continuum limit
At T = 0, collecting all the measured values at different lattice spacing, we can determine the continuum limit of the Tsallis parameter. Our result is presented in Fig. 8. We used second order polinomial to fit the numeric data and aquired  The particles emerging from a strongly interacting plasma are created locally and so they carry information about the local energy density. The distribution of this quantity is in general different from the canonical energy level distribution, except for the free (or very weakly interacting) theories. Therefore in the particle yields coming from a strongly interacting plasma we should not expect Boltzmann distribution. For the actual expectation we have to measure the distribution of the local energy density.
In this work we have determined the local energy density distribution with histogram method in the classical Φ 4 and in the quantum SU(3) Yang-Mills theory. In both cases the energy level distribution is Boltzmannian, but we have found that the Boltzmann distribution does not fit well to the local energy distribution in either case. However, the Tsallis distribution is a good fit, similarly to experimental data. The corresponding Tsallis parameter differs significantly from 1. Thermodynamic limit analysis is performed in both cases and we carried out the continuum limit analysis as well for SU(3) gauge theory. We remark that the renormalization of the local energy density does not affect the power of the power law tail, i.e. the Tsallis parameter q receives no renormalization correction.
We have found that in case of the classical Φ 4 theory, the Tsallis parameter q = 1.024 ± 0.001 -this is in the order of the experimental values obtained from heavy ion collisions. Regarding the SU(3) gauge theory, q = 0.9835 ± 0.0005. Interestingly it is smaller than 1, although |1 − q| is in the same order of magnitude in both models.
These results encourage us to proceed to our main goal, namely to perform similar analysis for QCD and to compare the results with experimental data. This work was supported by the Hungarian Research Fund (OTKA) under contract No. K104292.