Studying the quark-gluon plasma (QGP) is nowadays one of the most important goals in high-energy physics. In recent years, experiments at the Relativistic Heavy-Ion Collider (RHIC) at Brookhaven National Laboratory  and the Large Hadron Collider (LHC) at CERN have provided a wealth of data. The most striking result, obtained from analysis of these experimental data   , is that the deconfined quark-gluon matter behaves as an almost perfect fluid rather than a perfect gas, as it could be expected from the asymptotic freedom.
The most fundamental way to compute the properties of strongly interacting matter is provided by the lattice quantum chromodynamics    . However interpretation of these very complicated numerical computations requires the application of various quantum chromodynamics (QCD) motivated, albeit schematic, models simulating various aspects of the full theory and allowing for a deeper physical understanding. Moreover, such models are needed in cases when the lattice QCD fails, e.g. at large quark chemical potentials and out of equilibrium. It is, therefore, crucial to devise reliable and manageable theoretical tools for a quantitative description of non-Abelian QGP both in and out of equilibrium. Kinetic theory for the QGP can be formulated in two ways, namely: the color degrees of freedom are treated quantum mechanically and the distribution function of plasma constituents is a matrix in color space; in the second approach, on the other hand, the color may be considered as a continuous classical variable  . Here we use the latter approach that describes a particle carrying a classical color charge interacting with the chromodynamic field. The used approach is based on a quasiparticle picture and is motivated by the expectation that the main features of non-Abelian plasmas can be understood in simple semi-classical terms without the difficulties inherent to a full quantum field-theoretical analysis.
In this work, we extend the previous classical nonrelativistic simulations  to take into account quantum and spin effects. This is done in the frame of quantum Monte Carlo simulations where we rewrite the Wigner function of this system in the form of color path integrals. For the integration over color variables we have developed a procedure of sampling the color quasiparticle variables according to the SU(3) group Haar measure with the quadratic and cubic Casimir conditions. The developed approach self-consistently takes into account the Fermi (Bose) statistics of quarks (gluons). The main goal of this work is to calculate the quark and gluon pair distribution functions and momentum distribution functions and to treat influence of the strong interparticle interaction on these quantities.
2. Basics of the Model
The basic assumptions of the model are similar to those of Ref.   :
1) Quasiparticles masses (m) are of order or higher than the mean kinetic energy per particle. This assumption is based on the analysis of QCD lattice data    . For instance, at zero net-baryon density it amounts to , where T is a temperature.
2) In view of the first assumption, interparticle interaction is dominated by a color-electric Coulomb potential. Magnetic effects are neglected as subleading ones.
3) Relying on the fact that the color representations are large, the color operators are substituted by their average values, i.e. by Wong’s classical color vectors [eight-dimensional (8D) in SU(3)] with the quadratic and cubic Casimir conditions  .
4) We consider the 3-flavor quark model. For the sake of simplicity we assume the masses of “up”, “down” and “strange” quarks to be equal. As for the gluon quasiparticles, we allow their mass to be different (heavier) from that of quarks.
Thus, this model requires the following quantities as a function of temperature (T) and quark chemical potential ( ) as an input:
a) quasiparticle masses, for quarks and gluons , and
b) the coupling constant , or .
Input quantities should be deduced from lattice QCD data or from an appropriate model simulating these data. Applicability of such approach was discussed in Refs.   in detail. Our approach differs from that of Ref.   by a quantum treatment to quasiparticles instead of the classical one.
We consider a multi-component QGP consisting of color quasiparticles: gluons, quarks and antiquarks. The Hamiltonian of this system is with the kinetic and color Coulomb interaction parts
Here i and j summations run over quark and gluon quasiparticles, are their chemical potentials, , , and are total numbers of quarks and antiquarks of all flavours (up, down and strange), respectively, 3D vectors are quasiparticle dimensionless spatial coordinates, is coupling constant, the denote the Wong’s quasiparticle color variable (8D-vector in the group ), denote scalar product of color vectors. Nonrelativistic approximation for potential energy is used, while for kinetic energy we still keep relativistic form as the quasiparticle masses are not negligible as compared with temperature. The eigenvalue equation of this Hamiltonian is usually called the “spinless Salpeter equation”.
The grand canonical ensemble with given temperature, net-quark-number ( ) and strange ( ) chemical potentials, and fixed volume V are fully described by the grand partition function
Here . In Equation (2) we explicitly wrote sum over different quark flavors (u,d,s). Below the sum over quark degrees of freedom is understood in the same way. Usual choice of the strange chemical potential is (nonstrange matter), such that the total factor in front of is zero. Therefore, below we omit from the list of variables.
The partition function in canonical ensemble and the related thermodynamic properties of many particle systems of particles are defined by the diagonal matrix elements of the density-matrix operator
where x, and Q denote the multi-dimensional vectors related to spatial, spin and color degrees of freedom of quasiparticles with related flavor indexes respectively. The summation and spacial ( ) and color ( ) integrations run over all individual degrees of freedom of the particles, denotes integration over SU(3) group Haar measure   .
3. The Wigner Function for Canonical Ensemble
The Wigner function of the multiparticle system in canonical ensemble is defined as the Fourier transform of the off-diagonal matrix element of density matrix in coordinate representation   :
To avoid discussed in literature problems with definition of the relativistic Wigner function     in this paper we take in Equation (4) the non relativistic limit for kinetic energy operator (1).
Further we are going to obtain a new representation of the Wigner functions in the path integral form in the color phase space, which allows the numerical Monte Carlo simulations of the strongly coupled quantum systems of particles in canonical ensemble  -  . Average value of arbitrary quantum operator can be written as the Weyl’s symbol averaged over the color phase space ( ) with the Wigner function :
where the Weyl’s symbol of operator is:
Weyl’s symbols for usual operators like etc. can be easily calculated directly from definition (6).
4. Path Integral Representation of the Matrix Elements of the Density Matrix Operator
1For the sake of notation convenience, we ascribe superscript (0) to the original variables.
The exact density matrix of interacting quantum systems can be constructed using a path integral approach  based on the operator identity , where the r.h.s. contains identical factors with , which allows us to rewrite1 the integral in Equation (3) as follows
where , , spin gives rise to the spin part of the density matrix (S) with exchange effects accounted for by the permutation operators , and acting on the quark, antiquark and gluon degrees of freedom and the spin projections . The sum runs over all permutations. In Equation (7) and are permutation parity, while
is the off-diagonal element of the density matrix. Since the color charge is treated classically, we keep only diagonal terms in color degrees of freedom. Accordingly each quasiparticle is represented by a set of coordinates (“beads”) and an 8-dimensional color vector in the group. Thus, all “beads” of each quasiparticle are characterized by the same spin projection, flavor and color charge. Notice that masses and coupling constant in each are the same as those for the original quasiparticles, i.e. these are still defined by the actual temperature T.
The main advantage of decomposition (7) is that it allows us to use perturbation theory to obtain approximation for density matrices , which is applicable due to smallness of artificially introduced factor . Each factor should be calculated with the accuracy of order of with , as in this case the error of the whole product in the limit of large M will be equal to zero. In the limit can be approximated by a product of two-particle density matrices . Generalizing the electrodynamic plasma results  to the quark-gluon plasma case, we write approximate
In Equation (9) the effective total color interaction energy
is described in terms of the off-diagonal elements of the effective potential approximated by the diagonal ones by means of
Here the diagonal two-particle effective quantum Kelbg potential is
with , , Other quantities in Equation (9) are defined as follows: with
being a thermal wavelength of an a type quasiparticle ( ). The antisymmetrization and symmetrization takes into account quantum statistics and results in appearing permanent for gluons and determinants for quarks/antiquarks.
Functions ( are defined by modified Bessel functions. Gluon matrix elements are , while quark and antiquark matrix elements depend additionally on spin variables and flavor index of the particle, which can take values “up”, “down” and “strange”, and are the Kronecker’s deltas. These functions allow to exclude the Pauli blocking for particles with different spins, flavors and colors. Here arguments of modified Bessel functions are .
The coordinates of the quasiparticle ‘‘beads” , ( ) are expressed in terms of and vectors between neighboring beads of an i particle, defined as , while are vector variable of integration in Equation (7).
In the limit of functions describe the new relativistic measure of developed color path integrals. This measure is created by relativistic operator of kinetic energy . Let us note that in the limit of large particle mass the relativistic measure coincide with the Gaussian one used in Feynman and Wiener path integrals.
According to Equation (4) the antisymmetrized Wigner function can be written as the Fourier transform of the off-diagonal matrix element of the density matrix:
where is constant and . The antisymmetrization for quarks and symmetrization for gluons takes into account quantum statistics. Here we have replaced variables of integration for any given permutation by relation
In Equation (12) E is the unit matrix, while the matrix presenting permutation is equal to unit matrix with appropriately transposed columns.
5. Harmonic and Linear Approximation for the Wigner Function
The expression (12) for the Wigner function contains complex-valued function and so is inconvenient for Monte Carlo simulations, usually making use real valued functions. To overcome this difficulty we have to get explicit analytical real valued expression for the Wigner function. This can be done for free particles ( ) and linear or harmonic potentials . For general this integral cannot be analytically calculated but proper approximation of the Wigner function for any potential can be obtained from the Taylor expansion up to the first or second order in the variables :
Here . The second term means scalar product of the vector related to combination with the multidimensional gradient of pseudopotential, while the third term means quadratic form with matrix of the second derivatives.
This approximation for the Wigner function takes the form of the gaussian integral and can be calculated analytically. Here for simplicity let us consider expressions related to the linear approximation accounting for the linear term in this expansion. It is straightforward to account for the quadratic term. In the linear approximation, the Wigner function can be written in the form:
To test the developed approximations calculations of thermodynamic values and the ground state wave functions for quantum particle in 1D and 3D potential field, which strongly differs from harmonic one have been carried out in   . The used approximation gives practically exact results even for potentials, which have no matter with harmonic potential.
In degenerate system, average distance between fermions is less than the thermal wavelength and virtual trajectories in path integrals (15) are strongly entangled. This is the reason why permutations cannot significantly affect the potential energy in (15) in comparison with the case of the identical permutation. So we can replace permutation in the potential energy in (15) by the identical permutation. Now all permutations in (15) are acting only on variables x and and can be taken out of the path integral.
As was done in  for electromagnetic plasma it is enough for our purpose to take into account the pair permutations. These permutations cannot significantly affect the potential energy in (15) in comparison with the case of the identical permutation, as virtual trajectories in path integrals (15) are strongly entangled. So we can replace permutation in the potential energy in (15) by the identical permutation    . In this approximation permutations in (15) are acting only on variables x and and can be taken out of the path integral. So the Wigner function is determined by the path integral over all closed trajectories and can be presented in the form in (15) in comparison with the case of the identical permutation, as.
So the Wigner function is determined by the path integral over all closed trajectories ( ) and can be presented in the form
Here in the delta function of color vectors we have introduced a fit parameter z and and are the Kronecker’s deltas depending on spin and flavor indexes of quasiparticle. The later can take values “up”, “down” and “strange”. These functions allow to realize the exchange interaction for particles with the same spins, flavors and colors.
To regularize integration over momenta and to avoid difficulties arising at Monte Carlo simulations due to the presence of delta-function in expression (16) let us consider the positive Husimi distributions being a coarse-graining Wigner function with a Gaussian smoothing  . So the Wigner function has to be averaged over arbitrary “small” phase space cell. Let us introduce the function
with the equivalent cell area in the color phase space
For small and making use of the mean value theorem for integrals the Wigner function can be written in the form:
where the final expression for the phase space pair pseudopotentials accounting for the quantum statistical effects look like:
Here , and . To extent the region of applicability
of obtained phase space pair pseudopotential the and can be considered as a fit function small in comparison with unity. Our test calculations  have shown that the best fit for can be written in the form , while and were of order 0.1. Note that the expression (17) explicitly contains term related the classical Maxwell distribution however modified by terms accounting for influence of interaction on the momentum distribution function.
6. Simulation of QGP
The main idea of the Monte Carlo simulations of QGP consists in constructing a Markov chain of different quasiparticle states in the color configuration or color phase spaces  . The computational procedure comprises two stages. At the first stage, a dominant, i.e., maximal, -term in the sum of (2) is determined by calculations in the grand canonical ensemble. This term is indeed dominant in the thermodynamic limit of the box volume ( ). In the grand canonical ensemble, the quasiparticle numbers in the simulation box are varied, i.e., the consecutive states of the Markov chain can differ from each other by numbers of quarks, antiquarks, and gluons. Transitions between these states are the first type of Markovian elementary step. At the second type of elementary step, the coordinates of a randomly chosen quasiparticle were changed. The color variables are changed according to the SU(3) group Haar measure at the third type of Markovian elementary step. The Markov chain is generated until a full convergence of calculated values is achieved. This allows one to determine the average numbers of quarks, antiquarks, and gluons in the box at fixed temperature. Here, only densities of each species, i.e., the ratios of these average numbers with respect to the box volume, have physical meaning. Usually, after several million elementary steps, the average numbers of these quasiparticles become stable, and, for example, at the zero baryon chemical potential, the average number of quarks practically equals the average number of antiquarks. This equality can be considered as an inherent test of the consistency of the calculations.
At the second stage, the fixed number of quarks, antiquarks, and gluons is chosen to be equal to the obtained at the first stage average values of quasiparticles and calculations are performed in the canonical ensemble. However now according to (17) it is necessary do to integration in the color phase space and beside the second and third types of elementary step described above it is necessary to introduce new type of Markovian elementary step changing the momentum of quasiparticle  .
The input parameters of the model should be deduced from the lattice QCD data or can be taken as HTL values of and . In the present simulations we take only a possible set of parameters  :
where is the number of quark flavors which can be excited, for SU(3) group, and is the QCD running coupling constant squared, generally depending on T and all . All masses depend on combinations
and rather than on two independent variables T and .
It is also reasonable to assume that is a function of this single variable . This choice is done because like gluons is related to the whole system rather than one specific quark flavor. Then we can use the same “one-loop analytic coupling”
and substitute Q by to use this coupling in our simulations.
7. The QGP Mometum Distribution Functions
Here we present preliminary results of our QGP simulations. Figure 1 shows dependences of the quark and gluon densities on temperature at barion chemical potential equal to zero ( ) obtained according to the Equation (2) at the first stage of Monte Carlo simulation in grand canonical ensemble  .
Figure 1. (Color online) The quark and gluon densities versus temperature at barion chemical potential equal to zero.
Making use of the quark and gluon densities (see on Figure 1) obtained in grand canonical ensemble we present the Monte Carlo calculations in canonical ensemble of the QGP pair distribution and momentum distribution functions of quark, antiquark and gluons averaged over spin, color and flavours of quasiparticles.
First of all for physical analysis of interparticle interaction let us consider spatial arrangement of the quasiparticles in the QGP by studying a pair distribution function (PDF) defined as
where and are types of the particles ( or g). The PDF gives a probability density to find a pair of particles of types a and b at a certain distance from each other. The PDF depends only on the difference of coordinates because of the translational invariance of the system. In a non-interacting classical system, , whereas interactions and quantum statistics result in a redistribution of the particles. At temperatures , 350, and 525 MeV the PDF averaged over the quasiparticle spin, colors and flavors are shown in Figure 2.
The quasiparticle interaction is dominated by attraction at short distances. Indeed, the QGP lowers its total energy by minimizing the color Coulomb interaction energy via a spontaneous “anti-ferromagnetic”-like ordering of color vectors, i.e. the color vectors of nearest neighbor quasiparticles become anti-parallel. This short-distance attraction is stronger for gluon-gluon and gluon-(anti)quark pairs than for (anti)quark-(anti)quark ones because of the corresponding difference in values of quadratic Casimir invariants, which determine the maximal values of the effective color charge products in color Kelbg (Coulomb) potentials. For gluon-gluon pairs , for gluon-(anti)quark pairs , and for (anti)quark-(anti)quark pairs . Stronger gg attraction additionally enhances correlation of the gluon-gluon pairs at short distances. At the same time
Figure 2. (Color online) The Monte Carlo results for averaged over color, flavor and spin variables pair distribution functions for quarks, antiqurks and gluons of the non ideal (solid lines 1 - 6) and ideal (dashed lines 1 - 6) QGP plasma. Here the dimensionless average distance between quasiparticles is Wigner-Seitz radius ( , n is the density of all quasiparticles, ). For top, middle and bottom panels temperatures are , 350, and 525 MeV respectively. Small oscillations indicate the statistical error of Monte Carlo calculations.
the short-distance attraction is the only reason of the gluon-(anti)quark short-distance correlation.
The short-distance correlation implies formation of the bound states of the gluon clusters (glueballs), if the product related to gluons has maximum (not shown here)  . Glueballs and strongly correlated pairs of the other quasiparticles are uniformly distributed in space. The last conclusion comes from the fact that the gg, qq and PDF’s at distance larger than 0.7 fm are practically equal to unity. Possible existence of the medium-modified bound states was actively discussed some time ago, e.g., in  and later in  based on results from lattice QCD calculations of spectral functions   .
Obtained at the first time the Monte Carlo results for strongly coupled QGP momentum distribution functions have been calculated at the mentioned before second stage of simulation. From physical point of view the momentum distribution function gives a probability density for particle of type a to have momentum p:
The non ideal classical systems of particles due to the commutativity of the kinetic and potential energy operators have Maxwell distribution function (MD) proportional to even at strong coupling. Here at Figure 3 and Figure 4 we present results for both ideal QGP plasmas and results for
Figure 3. (Color online) The Monte Carlo results for averaged over color, flavor and spin variables momentum distribution functions for quarks and gluons of the non ideal (solid lines) and ideal (dashed lines) QGP plasma. The dotted line presents the Maxwell distributions. The momentum distribution functions for quarks and antiquarks practically coincide with each other (as example, the solid and dash-dot-dot lines 1). Here , is the thermal wave length and is the Plank’s constant. For top, middle and bottom panels temperatures are , 350, and 525 MeV respectively. All momentum distribution functions are normalized to unity.
Figure 4. (Color online) The momentum distribution functions for quarks and gluons in logarithmic scale of the non ideal and ideal QGP plasmas (solid lines 1 and 2 respectively) for . The dotted line 3 presents the Maxwell distributions. The momentum distribution functions for quarks and antiquarks practically coincide with each other (as example, the solid and dash-dot-dot lines 1). High energy quantum “tail” (the solid line 4) is approximated by sum of the Maxwell distributions and product of and the Maxwell distributions with effective temperature that exceeds the temperature of medium  . So in approximation the two fit constants are made use of. For top, middle and bottom panels temperatures are , 350, and 525 MeV respectively.
strongly coupled QGP plasmas. On these figures the dotted lines present the standard Maxwell distributions, while the dashed lines present Monte Carlo calculations for ideal QGP plasma.
The momentum distribution functions for quarks and antiquarks practically coincide with each other (as example, the solid and dash-dot-dot lines 1).
Quantum effects can affect the shape of kinetic energy distribution function. Quantum ideal systems of particles due to the quantum statistics have Fermi or Bose momentum distribution functions. In addition interaction of a quantum particle with its surroundings restricts the volume of configuration space available for particles, which, can also affect the shape of momentum distribution function due to the uncertainty relation, i.e., in a rise in the fraction of particles with higher momenta  .
Solid lines on Figure 3 and Figure 4 present the Monte Carlo results supporting these phenomena for non ideal QGP plasma. Difference in momentum distribution functions for both ideal and non ideal QGP plasma is decreasing with increasing temperature and lowering coupling plasma parameter. As before the main physical reason responsible for difference in behaviour of momentum distribution function of quarks and gluons is that the quadratic Casimir values responsible for interparticle interaction is significantly larger for gluons in comparison to quarks.
The peculiarities in asymptotic region of the quark momentum distribution function relates in arising out “quantum tails” due to the uncertainty relation as we have mentioned above. High energy quantum “tail” (the solid line 4) is approximated by sum of the Maxwell distributions and product of and the Maxwell distributions with effective temperature that exceeds the temperature of medium  . So in this approximation the two fit constants are made use of.
To verify the relevance of all above discussed trends, a more refined color-, flavor-, spin-resolving analysis of the distribution functions is necessary. This work is presently in progress.
In this paper, we have derived the new color path integral representation of Wigner function for the QGP plasma quasiparticle model for canonical ensemble. We have obtained explicit expression of the Wigner function in linear approximation resembling the Maxwell-Boltzmann distribution on momentum variables, but with quantum corrections. This approximation contains also the oscillatory multiplier describing quantum interference. The Monte Carlo calculations of the quark and gluon densities, pair distribution function and the momentum distribution function for strongly coupled QGP plasma in thermal equilibrium at barion chemical potential equal to zero have been done. Comparison with classical Maxwell-Boltzmann shows the significant influence of the interparticle interaction on the high energy asymptotics of the momentum distribution functions resulting in appearance of quantum “tails”. As a continuation of this work we are going to do detailed calculations of the spacial pair distribution functions, color correlation functions, investigate screening properties of QGP and compare Debye mass with available lattice QCD calculations.
We acknowledge stimulating discussions with A.N. Starostin, D. Blaschke and Yu.B. Ivanov. Authors acknowledge the Presidium RAS program number 6 “New approaches to the creation and study of state of matter under extreme conditions” for financial support.
 Shen, C. and Heinz, U. (2012) Collision Energy Dependence of Viscous Hydrodynamic Flow in Relativistic Heavy-Ion Collisions. Physical Review, 85, Article ID: 054902.
 Borsanyi, S., Endrodi, G., Fodor, Z., Jakovac, A., Katz, S.D., Krieg, S., Ratti, C. and Szabo, K.K. (2010) The QCD Equation of State with Dynamical Quarks. Journal of High Energy Physics, 2010, 77.
 Gelman, B.A., Shuryak, E.V. and Zahed, I. (2006) Classical Strongly Coupled Quark-Gluon Plasma. II. Screening and Equation of State. Physical Review C, 74, Article ID: 044908.
 Petreczky, P., Karsch, F., Laermann, E., Stickan, S. and Wetzorke, I. (2002) Temporal Quark and Gluon Propagators: Measuring the Quasiparticle Masses. Nuclear Physics B—Proceedings Supplements, 109, 513-515.
 Karsch, F. and Kitazawa, M. (2009) Quark Propagator at Finite Temperature and Finite Momentum in Quenched Lattice QCD. Physical Review D, 80, Article ID: 056001.
 Filinov, V.S., Ivanov, Yu.B., Bonitz, M., Fortov, V.E. and Levashov, P.R. (2013) Color Path-Integral Montecarlo Simulations of Quark-Gluon Plasma: Thermodynamic and Transport Properties. Physical Review C, 87, Article ID: 035207.
 Larkin, A.S. and Filinov, V.S. (2014) Wigner’s Pseudo-Particle Relativistic Dynamics in External Potential Field. Physics Letters A, 378, 1876-1882.
 Larkin, A.S., Filinov, V.S. and Fortov, V.E. (2016) Path Integral Representation of the Wigner Function in Canonical Ensemble. Contributions to Plasma Physics, 56, 187-196.
 Larkin, A.S. and Filinov, V.S. (2017) Phase Space Path Integral Representation for Wigner Function. Journal of Applied Mathematics and Physics, 5, 392-411.
 Zamalin, V.M. and Norman, G.E. (1973) The Monte-Carlo Method in Feynman's Formulation of Quantum Statistics. USSR Computational Mathematics and Mathematical Physics, 13, 169.
 Larkin, A.S., Filinov, V.S. and Fortov, V.E. (2018) Peculiarities of Momentum Distribution Functions of Strongly Correlated Charged Fermions. Journal of Physics A: Mathematical and Theoretical, 51, Article ID: 035002.
 Datta, S., Karsch, F., Petreczky, P. and Wetzorke, I. (2004) Behavior of Charmonium Systems after Deconfinement. Physical Review D, 69, Article ID: 094507.
 Starostin, A.N., Grjaznov, V.K. and Petrushevich, Yu.V. (2017) Development of the Theory of Momentum Distribution of Particles with Regard to Quantum Phenomena Page. Journal of Experimental and Theoretical Physics, 125, 940-947.