Fractals, non-extensive statistics and QCD

In this work we analyse how scaling properties of Yang-Mills field theory manifest as self-similarity of truncated n-point functions by scale evolution. The presence of such structures, which actually behaves as fractals, allow for recurrent non-perturbative calculation of any vertex. Some general properties are indeed independent of the perturbative order, what simplifies the non-perturbative calculations. We show that for sufficiently high perturbative orders a statistical approach can be used, the non extensive statistics is obtained, and the Tsallis index, $q$, is deduced in terms of the field theory parameters. The results are applied to QCD in the one-loop approximation, where $q$ can be calculated, resulting in a good agreement with the value obtained experimentally. We discuss how this approach allows to understand some intriguing experimental findings in high energy collisions, as the behavior of multiplicity against collision energy, long-tail distributions and the fractal dimension observed in intermittency analysis.


Introduction
Yang-Mills field (YMF) theory is a prototype theory to describe three among the four known natural interactions, namely, strong, weak and electromagnetic, with only gravitational force left aside. The importance of such cathegory of field theory is summarized in the Satandard Model, which is a single framework to describe all three interactions in a unified formalism. Scale invariance is a fundamental aspect of the Yang-Mills field theory, playing an important role in the renormalization of the theory after the ultraviolet divergences are subtracted [1][2][3]. The theory was extremely succesful in QED and Weak interactions, providing accurate and precise descriptions for the observed phenomena. QCD has proved to be more challenging, since the calculational methods applied to the other two forces is not appropriate to provide accurate results for strong interacting systems. In this scenario. QCD has been tested in high energy elementary particle

Formalism
The simplest scale free non-abelian gauge field theory has Lagrangian density including bosons and fermions given by where and where ψ and A represent, respectively, the fermion and the vector fields, with f abc being the structure constants of the group and T a the matrices of the group generators in the fermion representation. The UV regularized vertex functions are related to the renormalized vertex functions with renormalized parameters,m andḡ, as Γ(p, m, g) = λ −D Γ(p,m,ḡ) .
This property is mathematically described by the renormalization group equation, which introduce the beta-function, which allowed to show that QCD is asymptotically free [9][10][11][12]. Such equation is known as Callan-Symanzik equation [13][14][15], and is given by where M is the scale parameter, and the β-functions are defined as D = D o + d, with D o being the natural dimension of the phase-space. In general, d is not an integer, therefore the scaling dimension, D, may be fractionary. The parameter γ is the anomalous dimension given by a combination of the scaling dimensions of the fields ψ and A.
In the irreducible one-particle representation, self-energy is taken into account by using the effective parton mass,m, what allows a reduction of the calculational complexity, since only proper vertexes must be considered. As schematically shown in Fig. 2, the scale invariance means that, after proper scaling, the loop in a higher order graph is identical to a loop in lower orders [9][10][11][12]. The scale invariance of the vertex function is a direct consequence of the Callan-Symanzik equation, and it is of fundamental importance in what follows. When effective charges and masses are used, the line of the respective field in Feynman diagram represent an effective particle or state, since parts of the diagram representing self-energy contribution are omitted because they are already taken into account by the renormalization. Such diagrams are called irreducible graphs. We will next refer to the physical system represented by irreducible graphs as effective parton. In irreducible graphs, the only allowed vertexes, called proper vertexes, are related to the creation of an effective parton. Despite the great reduction in complexity when irreducible graphs are used, still complex graphs are unavoidable in performing perturbative calculations, specially for QCD, since vertex functions may include several orders in perturbative approximation. We will show that some additional simplification may be achieved under some circumstances that are relevant in hadron structure and in multiparticle production. Preliminary results of the approach developed here were discussed in Ref. [16].

Statistical description of the partonic state
In this work the quantity of interest is the trace of the evolution operator where Ψ o is an eigenvector of the Hamiltonian H corresponding to the initial state of the fields. H governs the effective parton evolution in the presence of the strong interaction. From now on we refer to proper vertexes as interactions.
The time evolution of an initial partonic state is The state Ψ can be written as with the basis states |Ψ n corresponding to a fixed number, n, of interactions in the vertex function. In the perturbative method each proper vertex gives rise to a term in the Dyson series, and at any time, t, the partonic state is given by where g represents the interaction and t n > t n−1 > · · · > t 1 > t o , and {n} runs over all possible terms. with n interaction vertexes. Observe that Ψ n |Ψ n = δ n n . The symbol {n} indicates that the summation is performed over all possible configuration of fields with n interactions.
The number of particles in the state |Ψ n is not directly related to n, since high order contributions to the N particles states can be important, but certainly N ≤ M (n) := n(Ñ −1)+1, whereÑ is the number of particles created or annihilated at each interaction. In Yang-Mills field theoryÑ = 2. 1 We can introduce states of well defined number of effective partons, |ψ N , so that The states |ψ N are autovectors of H o with fixed number of particles, N . Of course Ψ n |ψ N = 0 whenever N > M (n), and Ψ N |Ψ N = δ N N . Since the number of partons is fixed and they do not interact but by contact interaction, the states |ψ N can be understood as the states of an ideal gas of N partons. Therefore where m i and p i are the mass and momentum of the i partonic state, and γ i represents all relevant quantum numbers necessary to completely characterize the partonic state. S is the symmetrization operator acting over fermions and bosons. Since the mass of the effective partons varies continuously, S gives a negligible modification of the single parton states, so mass and momentum of each parton can vary independently, as far as the total energy is conserved. Notice that the states with N partons can be obtained in several ways, since with ψ N |Ψ n = 0 for M (n) ≥ N . For sufficiently high number of interactions, n, there are so many ways to obtain a N particle system that the possibility to get a particular configuration becomes insensitive to the initial state. This situation is similar to that of complex systems where a statistical approach is possible. Therefore we assume that the probability to get a particular state can be calculated on the sole basis of the number of configurations corresponding to that state, and supposed that all configurations have the same probability to be observed. This assumption is not really new, since it is in fact the same statistical assumption made in LQCD where the gluon field configuration is obtained by statistical approach following a Monte Carlo method. The probability to find a state where at least one parton has mass between m o and m o + dm o , and momentum coordinates between p oi and p oi + dp oi is given by matrix elements like so we have the probability For reasons that will become apparent later, we will suppose that the state |Ψ(t) has energy between E and E + dE that follows a probability distribution P (E). The first bracket in Eq. (14) is related to the probability to have n interactions up to the instant t, and depends on the intensity of the interaction, G, so Therefore, the bracket in Eq. (14) must be associated to the probability The second bracket depends on the relative number of possibilities to get the configuration with N particles after n interactions, so The last bracket in the expression above is calculated statistically, see appendix for technical details. The result is Notice that the density of states with one particle with energy p 0 j can be obtained by integrating on the vector momentum coordinates, resulting in Observe that if we write E = (4N − 5)µ, and take the limits N → ∞ and µ → 0 keeping constant N µ = kT , we obtain the exponential factor usually found in the thermodynamical limit. Here, however, we cannot take such limit because N < M (n). At the rest mass frame the total energy of the system is equal to its mass, i.e., E = M , and setting p 0 j = ε j , then Eq. (18) can be written as with and From the expression for A(N ) we see that the number of states increases with N 4 , so those configurations with large number of particles are favoured. The maximum number possible is M (n) = n(Ñ − 1)+1, so the probability to get a configuration with N particles will be, for sufficiently large n, approximately (N/(nÑ − n)) 4 , so using Eqs. (15) and (16) we obtain Observe that when N is sufficiently large and ε/M sufficiently smaller than unit, we have therefore we can set and finally obtain 3 Self-similarity and fractal structure The result obtained in Eq. (26) shows that for an ideal gas with finite number of particles, the probability depends on a power-law function of the ratio ε j /M . While this is valid for an ideal gas, where particles have no internal structure, the same result cannot be directly applied to the gas of effective partons because effective partons always have internal structure related to self-energy contributions coming from interaction with the vacuum.
We will now investigate how the scaling properties determined by the renormalization group equation can be used to obtain the result for effective partons. The scaling properties represented by the Callan-Simanzyk Equation demands that effective partons at any scale are similar to any other parton after the proper rescaling is performed. In other words, if partonic properties are expressed in terms of scale free variables they must be described by the same functions of those variables. The energy distribution of a parton, as given by Eq. (26), depends on the ratio χ = ε j /E. Now let us consider that the system with energy E in which the parton with energy ε j is one among N constituents, is itself a parton inside a larger system with energy M. Then the ration E/M is represented by the same variable, χ, that describes the ration ε j /E. In addition, we can write P (E) in Eq. (26) as P (E/M), and this probability density must follow that same function as P (ε j /E), that is where the sign ∼ indicates that the three functions above are indeed the same function of the scale free variable χ. The same conclusion is valid for any parton with energy ε which is a component of another parton with energy Λ, and to let it clear in the following calculations we write χ = ε/Λ. The reasoning just used here is the same that applies to fractal structures, so what we are doing is to introduce the mathematical tools common to fractals studies in the analysis of Yang-Mills fields. The self-similarity among the partons implies that the probability that the parent parton with mass E inside a larger system with mass M is similar to the probability given in Eq. (26), so the dependence on N that appears in the exponent must be changed to a parameter α, which remains to be determined. This parameter represents to the total number of degrees of freedom of the fractal, playing the same role of the exponent 4N − 5 in the case of the ideal gas. We can write and the function being integrated in Eq. (26) contains the term Notice that the parameters α and ν describe the complexity of the interaction involved in the gauge field interaction, and in this sense they measure the sensibility of the effective parton to its internal degrees of freedom. The self-similarity implies, from relation (29), that since the same probability governs any parton distribution.Therefore The parameter ν represents the fraction of total number of degrees of freedom of the state |ψ N that is involved in each interaction. Observe that we expect ν ≤ 1, therefore α is positive.
Observe that q is related to the resolution parameter ν and to the number of particles, N . But this two parameters are not independent, since N increases as the resolution increases. In fact, we can write where λ is a reduced scale, independent of the number of degrees of freedom relevant to the system. Since variable χ must follow a universal distribution for any parton that is independent of the level in the fractal structure, and since the smallest parton is the one obtained in 1-loop approximation, where N =Ñ , and therefore is constant, then also 1/α = λ/Λ is independent of the position the system occupy in the fractal structure, and so α is constant. We can introduce a parameter q such that with q > 1 been constant for any parton in the fractal system. Then we obtain This result shows that the distribution of parton energy created by a system governed by Yang-Mills fields depends only on the ratio between the parton energy, ε, and the energy scale per degree of freedom λ. Furthermore, it shows that the energy distribution follows the q-exponential function commonly found in Tsallis non extensive statistics. This result shows that the distribution of parton energy created by a system governed by Yang-Mills fields depends only on the ratio between the parton energy, ε, and the energy scale per degree of freedom λ. Furthermore, it shows that the energy distribution follows the q-exponential function commonly found in Tsallis non extensive statistics.
Similar results have been obtained through a different approach using the concept of thermofractals, introduced in Ref. [17] and studied in details in Ref. [18]. There it is shown that the fractal structure leads to the non extensive statistics, and it is discussed the relations between thermofractals and Hagedorn's self-consistent thermodynamics developed to study high energy collisions [19,20], and that was extended to non extensive statistics [21].
The results obtained here mean that an effective parton with energy Λ can be represented as a system of effective partons at an energy scale λ, which can be represented by states of the type |γ 1 , ε 1 , . . . , γ N , ε N with ε 1 , . . . , ε N ∼ λ. A particular state where one particle has energy between ε o and ε o + dε o , denoted by |γ o , ε o = |γ o , ε o , . . . , and the matrix element becomes This is obtained by increasing one order in the perturbative calculation, that is, this results from a vertex function Z to which one loop is added and an external line is added to the new loop. This means that the q-exponential plays the role of an effective coupling constant in the vertex function, that is with

Discussion and applications
As an application of the formalism proposed here, let us consider a vertex in two different orders, as depicted in Fig. 2. The vertex function is The next order in perturbative approximation is given by the vertex with one additional loop, which results in a vertex function whereḡ is the effective coupling given bȳ with Observe that the effect of the effective coupling as defined by Eq. (37) is to increase the contribution of low energy partons and decrease the contribution of high energy partons. At each vertex, the creation of a parton with energy much higher than that expected for partons at the fractal level as determined by the scale λ is strongly suppressed. The same is valid for the parton annihilation: partons with high energy are more likely to remain in the system with small probability to interact with the less energetic and therefore less massive partons in the medium. Depending on its color content, the heavy parton can reach the surface of the medium and scape. In this way, a complete description of the behaviour of those heavy partons must include the possibility of color exchange and coalescence, and some works in this direction, using power-law distributions, have already been done [22]. Keeping these aspects in mind, we continue here with a purely cinematic description of the partonic evolution of QCD. Due to energy-momentum conservation constraint, ε 5 + ε 6 ∼ λ, so the coupling strongly favors configurations where ε 5 ∼ ε 6 ∼ λ/2. In this way the effective coupling controls how the energy transferred to an effective parton is distributed among its internal degrees of freedom.
Notice that Γ M and Γ o are similar, differing only by the presence of the non interacting field (γ 4 p 4 ) and by the different energies of the fields at each vertex. The scaling properties of Yang-Mills fields allow us to relate Γ M to Γ o by an appropriate scale, λ. This scale must take into account that this is a 1-loop contribution, so λ/λ o =Ñ −1 , and the coupling constant and fields in Γ M must be modified accordingly to [1,2] where m is the mass of the parton at scale λ o andm is its mass at scale λ. Here φ i refers to the field of parton i, and Z i are multiplicative factors arising from the vertex renormalization.
The renormalization properties imposes that [1, 2] For one loop approximation we have λ = λ o + dλ, where dλ < 0, and the integration above turns into from where it follows that The field (γ 4 p 4 ) in Γ M does not interact at this vertex, so it does not need to be scaled, so it will be omitted in the following calculation for the sake of clarity. The scaling behavior of the Γ(M ) function is obtained from dimensional analysis, and it results to be Γ M (λ) = (λ/λ o ) 4 , after energy-momentum conservation is taken into account, so with d = 4 in the present case. From these considerations and from Callan-Symanzik equation, it results that The contribution to the scaling transformation fromḡ appears through g(λ), while the other two vertex are at the inital scale, that is, g = g(λ o ). Also, in order to compare with QCD results, we study the behavior of g(λ) at λ = λ o /µ. We have the logarithmic derivative The beta function can be easily calculated and we get and we emphasize thatÑ = 2. Scaling properties of QCD have been extensively studied in the 1-loop approximation. The beta-function for QCD is [31,32] where and therefore relating the entropic index, q, from Tsallis statistics to fundamental parameters of the field theory. Quantitatively, the parameters c 1 and c 2 are related to the number of colors and flavors by c 1 = N c and c 2 = N f /2. Using N c = N f /2 = 3 we get which leads to q = 1.14. From experimental data analysis [33][34][35] we have q = 1.14 ± 0.01, showing a good agreement between theory and experiments.  [31,32], and the dash-dotted line represents β as calculated in 4-loop approximation [36,37].

The beta function and the effective coupling
The result obtained here shows that in the 1-loop approximation we recover the behavior expected for QCD. Let us now consider the asymptotic limit, when µ → ∞. This limit is equivalent to say that the 1-loop approximation is made at higher orders of logarithmic approximation, i.e., when λ o → λ o λ o . In this case also λ → λ , such that the ratio λ /λ o ∼ λ/λ o , i.e., we are doing the same 1-loop approximation, but to a process that occurs at higher order of logarithmic approximation. At this point, we have to take into account that the fractal structure considered here, with all effective parton showing an internal structure, is broken since quarks are elementary particles. This implies the existence of an inferior limit to ε = ε o . When doing so, we have a minimum value for the scale, λ c = (q − 1)ε o , then the scale is redefined to λ /µ = λ o /µ − λ c . It results that g(µ ) → 0 and also the logarithmic derivative , µ∂g(µ)/∂µ → 0, as shown in Fig. 3. These results confirm that the effective coupling used here leads to the asymptotic freedom of quarks obtained in QCD, since g → 0 as µ → ∞. The beta function also approaches zero asymptoticaly, but for low values of the scale parameter µ it also goes to zero, presenting a peak around µ/λ o = 40. For a direct comparison with the QCD running coupling, we show in Fig. 3 the comparison with QCD calculations in 1-loop and in 4-loop calculations. We observe that in the asymptotic region, when g → 0, the behavior of the effective coupling proposed in the present work is similar to the QCD behavior, diverging from it for larger coupling values, when perturbative QCD is not expected to work.
It is interesting to note that the scaling invariance can be used in two ways: one can fix a value for λ, based on experimental resolution or some physical change in the system that makes the scale invariance invalid or useless, and then the several fractal layers are obtained by considering systems of large and large energies corresponding to different values of Λ; the other way is to fix the total energy and let λ vary, meaning that the system is observed with better resolution, and therefore deeper layers of the system are investigated, i.e., the better the resolution, the lower the λ. The first way finds application in high energy collisions, and the second, in studies of hadron structure.
The results obtained here have shown that a system with fractal structure, similar to the thermofractals introduced in Ref [17] and studied in [18], can be understood as a natural consequence of the scale invariance of gauge field theories. There is is shown that the fractal structure leads to the non extensive statistics, and it is discussed the relations between thermofractals and Hagedorn's self-consistent thermodynamics developed to study high energy collisions [19,20], and that was extended to non extensive statistics [21]. This fractal structure has been already used to investigate properties of hadrons [23], phasetransition in hot hadronic matter [24] and neutron stars [25]. The power-law distribution of energy and momentum, which is a direct consequence of the fractal structure, was used to describe p T distributions from high energy collisions experiments [34,35,38] and to describe hadron mass spectrum [26,34]. The results obtained here give a stronger basis for the interpretation of those experimental and phenomenological studies in terms of non extensive statistics.

Multiplicity
It is possible to understand, from the considerations made here, that the fractal structure of YMF is the basis for investigations of hadron properties [23], phase-transition in hot hadronic matter [24], neutron stars [25] and cosmic-ray [30]. These phenomenological approaches are, in fact, implementations of scaling symmetries observed in Yang-Mills fields. The fractal structure also allows the understanding of the self-similarity [27][28][29] and scaling properties observed in high energy experimental data. In fact, these findings are direct consequences of the scaling properties of YMF, as discussed here. Moreover, the fact that the entropic index, q, is obtained from well-known field-theoretical parameters, the results we have obtained allows a new interpretation of Tsallis statistics in terms of fractal structure in the same lines it was obtained in thermofractals approach [26].
The fractal structure presents at least one fractal dimension, and the Haussdorf dimension is a characteristic dimension that can be calculated by using the box-counting technique [6], where the dimension D is related to the number of boxes, N , necessary to completely cover all possible values for the measured quantity and D t is the topological dimension. At some scale r these quantities are related by [6] N r −D ∝ r −Dt .
In our case, D t = 1 since we are dealing with system energy as a measure. The procedure to obtain the Haussdorf dimension is similar to that followed in Ref. [17]. The average energy of the partons at the scale λ, already introduced as the energy scale per degree of freedom, is The ratio between the average energy of the components and the parent system energy is related to the level of the fractal structure relative to the scale λ by The number of boxes with length λ necessary to completely cover the possible range of energies in which the fractal components can be found is N =Ñ n . Then it follows from Eq. (57) that Eq. (57) also shows that, in terms of the scales, the energy of the system varies as E ∼ r −1 . Let us now write the dependence of the parton energies at scale λ as ε ∼ r −D , then therefore From Eq.-(58) it follows that From Eqs. (55) e (56), and using E = λ r /(q − 1), we get Using the value q − 1 = 1.14 it follows that D = 0.69. This result is in good agreement with the fractal dimension observed in intermittency analyses of high energy experimental data [39,40]. These analyses allows to access fractal dimensions by studying the behavior of cummulants of the measured distributions [41][42][43][44][45][46][47][48], and the systematics show that for pp collisions there is a good agreement between the value obtained from the theory with those resulting from experimental data analyses. The fractal dimension gives the behavior of the parton energy with the energy scale, r, that is, while the total energy goes as E ∝ r −1 , the partons observed at scale λ have energies that depend on the scale as ε ∝ r −D . A more direct way to access the fractal dimension is the particle multiplicity. In fact, being M the particle multiplicity, we have From the dimensional behavior obtained above, we get For the case of hadrons, as we have seen, q = 1.14 and D = 0.69, so we obtain M ∝ E 0.31 , which is in excellent agreement with the result obtained for pp collision at high energy [49], which gives, for a power-law fit, an exponent corresponding to 1 − D = 0.302.

Conclusions
In the present work the scaling properties f Yang-Mills fields are analyzed under the light of the concept of fractals. It is shown that those scaling properties lead to the formation of a fractal structure, and usig a thermodyanamical hypotheis similar to that used in LQCD we obtain Tsallis distributions that have been associated to the long-tail distributions observed in high energy data. The entropic index, q, which in Tsallis statistics is in most cases a parameter to be determined experimentally, in the present case can bedetermined completely in terms of the fundamental parameters of the field theory. The result is used to obtain, by a recurrence formula that reflects the self-similar features of the fractal, the effective coupling which is expressed in terms of a scale dependent formula. Applying this expression for the case of QCD in the 1-loop approximation, the entropic index is calculated, for the first time, from the numbers of colors and flavors. The result is shown to be in good agreement with the value obtained for q by fitting Tsallis distributions to data. From the analysis of the fractal structure we obtain the fractal dimension associated to Yang-Mills fields, which is determined completely in terms of the field theory parameters. The fractal dimension is calculated for QCD and the result is in good agreement with the value obtained from intermittency analysis of high energy distributions.
Finally, it is shown that the fractal dimension allows us to determine the behaviour of the particle multiplicity against the collision energy, and we obtain a result that is in good agreement with that observed in high energy collisions.
In summary, we have shown in the present work that renormalizable field theories lead to fractal structures, which can be studied, from a thermodynamical point of view, with Tsallis statistics. A recursive method allows to perform non perturbative calculations to describe the particles structure governed by the gauge theory. In the case of multiparticle production, the calculations lead to a thermodynamical description where non extensive statistics must be used. The results obtained here give a solid basis from QCD to the use of non extensive self-consistent thermodynamics to describe properties of strong interacting systems and to the use of thermofractal structure to describe hadrons.
where p µ i = (p 0 i , p i ) is the fourth momentum of particle i, the goal is to compute the density of states ρ(p), with the normalization condition In the following we won't assume a fixed value for the mass m i of particle i, where m 2 i = p µ p µ , so that p 0 i and p i are variables that may change independently each other. Let us compute the phase space volume of region H ≤ E, i.e.
where Θ(x) is the step function. We can express the integral in the form Note that for particle i one has p i 2 = (p 0 i ) 2 −m 2 i . Because m i is not fixed to a particular value, the limit of the integration in the d 3 p i integral is 0 ≤ | p i | ≤ p 0 i , and this has been expressed with the corresponding step function in Eq. (68). This integral can be easily computed to give d 3 p i Θ(p 0 i − | p i |) = 4 3 π(p 0 i ) 3 , so that Eq. (68) turns out to be From dimensional analysis one can see that the result of this integral should be of the form where c N are some coefficients to be determined, and the factor π N has been explicitly extracted for convenience. It is possible to obtain the coefficients c N by mathematical induction. In the case N = 1, one can easily check that so that c 1 = 1/3. Assuming that the expression of Ω N (E) is known, let us compute Ω N +1 (E). It writes where in the last equality we have performed the integral only in the momentum of particles i = 1, · · · , N . The integral in d 4 p N +1 can be easily performed following similar steps as above, i.e.
This result should be identified with Ω N +1 (E) = c N +1 π N +1 E 4(N +1) , and from it one obtains a relation between c N +1 and c N , which can be used iteratively to obtain In the last equality we have used that c 1 = 1/3. Finally, the result of Eq. (67) is Let us introduce a constant distribution in the microcanonical ensemble that should be normalized to 1 as indicated in Eq. (66). Then, one finds that so that C = 1 Ω N (E + ∆E) − Ω N (E) .
The probability distribution for particle j to have momentum p µ j is obtained by integrating the joint distrubution ρ(p 1 , · · · , p N ) over all the variables except p j . Then, one has f (p j )d 4 p j = d 4 p j d 4(N −1) p ρ(p) Θ(E ≤ H ≤ E + ∆E) where the prime in N i=1 means that the term i = j should be excluded in the summation. Using the result above for ρ(p), this can be expressed in the form Finally, by using the explicit expression of Ω N (E) given by Eq. (76) and taking the limit ∆E → 0, one arrives at the final result f (p j )d 4 p j = d 4 p j 1 8π Γ(4N ) Γ(4(N − 1))