Interacting color strings as the origin of the liquid behavior of the quark-gluon plasma

We study the radial distribution function of the color sources (strings) formed in hadronic collisions and the requirements to obtain a liquid. As a repulsive interaction is needed, we incorporate a concentric core in the strings as well as the probability that a string allows core-core overlaps. We find systems where the difference between the gas-liquid and confined-deconfined phase transition temperatures is small. This explains the experimentally observed liquid behavior of the quark-gluon plasma above the confined-deconfined transition temperature.


I. INTRODUCTION
The heavy-ion Au-Au collisions at RHIC obtained a liquid of quarks and gluons with a shear viscosity over entropy density ratio lower than any other material ever known (quark-gluon plasma) [1-3]. This discovery was confirmed later in Pb-Pb collisions at LHC through the study of all of the harmonics of the azimuthal distributions and different correlations showing the existence of a collective motion of quarks and gluons [4][5][6]. Most of the properties seen in heavy-ion collisions have also been observed in pp, pA collisions at LHC [7], and dAu and 3 HeAu at RHIC [8].
Multiparticle production in pp, pA, and AA collisions is currently described in terms of color strings stretched between the partons of the projectile and target, which decay into new strings and subsequently into hadrons. The strings are extended in the longitudinal space between the partons of the colliding projectiles, which transforms it into a rapidity space describing the rapidity differences of the partons. This difference is given by the energy fraction of the projectiles carried by the partons. On the other hand, the strings also have a transverse extension. Thus, color strings may be viewed as small areas in the transverse plane of the collision filled with color field created by the colliding partons. Due to confinement, the strings have transverse circular areas with radius r 0 =σ/2 around 0.2-0.3 fm (σ is the diameter). With the growing energy or size of the colliding systems, the number of strings grows, and they start to overlap and interact. In the color string percolation model (CSPM) the interaction between strings occurs forming clusters when they overlap; the color field is given by the SU(3) * jhony.ramirezcancino@viep.com.mx † bodiazj@math.uc3m.es ‡ pajares@fpaxp1.usc.es color sum. Due to the random direction of the color field in the color space, the intensity of the resulting color field is only the square root of the individual strings color field [9][10][11][12]. In rapidity, the interaction of strings is taken into account by incorporating the energy-momentum conservation of the formed clusters. This gives rise to the particle production in the forward region outside the kinematical limits, the so-called cumulative effect, observed at RHIC energies [13][14][15]. However, in our study, we are only concerned with the transverse size. At one critical string density a spanning cluster is formed through the collision surface. For a uniform profile of the string distribution, this critical density value, in the thermodynamic limit, matches the percolation threshold of the classical two dimensional (2D) continuum percolation model, given by N/S=1.128/πr 2 0 , where N and S are the number of strings and the collision area respectively. This critical value can change slightly for a finite and not large N and other profiles [16][17][18]. This critical percolation density is associated with a temperature T = 160 MeV, which corresponds to the confineddeconfined phase transition of the quark and gluon matter [19].
On the other hand, the physical structure of the systems can be determined by analyzing the behavior of the radial distribution function, g(r) (also known as the pair correlation function). It describes how the average number of particles varies as a function of the radial distance from a point. It is widely used to characterize packing structures and contains information about long-range interparticle correlations and their organization [20].
Since models considering strings like fully penetrable disks correspond to the picture of the classical ideal gas, it is expected that they will have a flat radial distribution function. In Fig. 1, we show for the CSPM the flat behavior of g(r), regardless of the string density. This happens even when the CSPM can explain most of the experimental data on pp, pA, and AA collisions, including the azimuthal distributions of the produced particles, as well as the temperature dependence of the ratio between the shear and bulk viscosities over the entropy density [21,22]. This ideal gas behavior of the CSPM has been observed by some authors. For example, in Refs. [23,24] the authors analyzed the electrical and thermal conductivity of the quark-gluon plasma (QGP) using the CSPM and found that the system corresponds to an ideal gas. Also, in Ref. [18] the authors performed a finite-size analysis of the speed of sound and concluded that the CSPM corresponds to a mean field theory. On the other hand, there are other strings models [25][26][27] where the strings interact in a different way, for instance, by color rearrangements [26] or by shoving when they are close to each other [27]. The repulsion between strings has also been used to study the harmonics of the azimuthal distributions, obtaining a reasonable agreement with the data [28,29].
In this paper, we propose a hybrid core-shell model together with the traditional CSPM to incorporate the excluding or repulsive interaction between strings. We do this by providing the strings with a concentric region of exclusion (core region) of diameter λσ (0 ≤ λ ≤ 1). The rest of the string area is called the shell region. We also consider a probability q λ that a string allows overlap in its core with the core of another string. We refer to the strings as soft or hard if they admit overlap in their core region or not, respectively. Notice that the overlap condition applies only to the core-core interactions, while all core-shell and shell-shell overlaps are allowed.
In what follows, we call this modification the core-shellcolor string percolation model (CSCSPM). Notice that i) if (λ = 1, q λ = 0) the system corresponds to a hard-disks fluid [30], ii) if λ = 0 or q λ = 1 the system returns to the traditional 2D continuum percolation [31,32], iii) if q λ = 0 this model reproduces the continuum percolation of disks with hard cores [33,34], and iv) if λ = 1 it resembles the random sequential absorption model [35,36]. In this sense, our model is a generalization of them. For the CSCSPM, two phase transitions are observed as a function of (λ, q λ ): ideal gas to non-ideal to liquid. This paper aims to determine if some combinations of (λ, q λ ) allow the system to simultaneously exhibit the gas-liquid transition and the emergence of the spanning cluster. To this end, we determine: a) the gas-liquid phase transition temperature, which is calculated through the observation of the oscillation of the radial distribution function, and b) the critical temperature associated with the percolation threshold, which corresponds to the QGP formation temperature.
The plan of the paper is as follows. In Sec. II, we present how to calculate the temperature in percolationbased string models. In Sec. III, we provide the simulation and data analysis methods. In Sec. IV, we show the phase diagram of the CSCSPM in terms of (λ, q λ ), and determine the region in the λ-q λ plane where the transition and critical temperatures are close. Finally, Sec. V contains our conclusions and perspectives.

II. TEMPERATURE FOR PERCOLATION-BASED STRING MODELS
The interaction among strings gives rise to a reduction in multiplicity and an increase in the average transverse momentum. A cluster of n strings (remember that each string is considered as a disk) that occupies an area S n behaves as a single color source with a higher color field Q n corresponding to the vectorial sum of the color charges of each individual string Q 1 . As Q n = n Q 1 and the individual string colors may be oriented in an arbitrary manner respective to each other, the average Q 1i Q 1j is zero and Q 2 n = n Q 2 1 . Knowing the color charge, we can calculate the multiplicity µ n and the mean transverse momentum squared p 2 T n of the particles produced by a cluster, which are proportional to the color charge and color field respectively, as where µ 1 and p 2 T 1 are the multiplicity and transverse momentum squared of the particles produced from a single string with transverse area S 1 = πr 2 0 . For strings just touching each other, S n = nS 1 and therefore µ n = nµ 1 and p 2 T n = p 2 T 1 . On the other hand, when strings fully overlap, S n = S 1 and therefore µ n = √ nµ 1 and p 2 T n = √ n p 2 T 1 , so that the multiplicity is maximally suppressed and the transverse momentum is maximally enhanced. In the thermodynamic limit, the average value of nS 1 /S n for all of the clusters is [11] with η = N S 1 /S being the filling factor, where N and S are the number of strings and total surface area, respectively. The function F (η) is called the color reduction factor, which can be written in terms of the area covered by strings φ(η) as (3) In the thermodynamic limit, for fully penetrable disks [9,37], However, the explicit form of φ depends on the considered string model, e.g., if the system is finite or without periodic boundary conditions with a particular geometry [18] or, as in our case, if the system is composed of a combination of disks allowing overlap or not in the core region. Now, Eqs. (1a) and (1b) can be written as We must recall that the strings decay into new ones through color neutral q −q or qq −qq pair production. The Schwinger QED 2 string-breaking mechanism produces these pairs at the time τ ∼ 1fm/c, which subsequently hadronize to produce the observed hadrons. The Schwinger distribution for the produced particles is dN/dP 2 T ∼ exp(−p 2 T π/x 2 ), where the average value of the string tension (color field) is x 2 . The chromoelectric field in the string is not constant but fluctuates around its average value. Such fluctuations lead to a Gaussian distribution for the chromoelectric field, transforming the Schwinger distribution into the thermal distribution [19] dN dp 2 with x 2 = π p 2 T 1 /F (η). Equation (6) shows that the effective temperature is This temperature is related to the Hawking-Unruh effect [38,39]. In AA and pp collisions the thermalization can occur through the existence of an event horizon due to a rapid deceleration of the colliding nuclei [40]. Barrier penetration of the event horizon leads to a partial loss of information and is the reason for the stochastic thermalization of the q−q pairs. In string percolation, as clusters grow and cover most of the collision area, this local temperature can be regarded as the temperature of the total thermal distribution. In the 2D continuum percolation, for a uniform density profile, the critical filling factor at which the spanning cluster emerges is η c ∼1.128 [32]. This value can change slightly for a finite N and other profiles [16][17][18]. Introducing this value into Eq. (7), for a reasonable value of p 2 T 1 , around 200 MeV, we obtain a critical temperature T c around 160 MeV, which corresponds to the confined-deconfined phase transition of the quark and gluon matter [19].

III. SIMULATION METHOD AND DATA ANALYSIS
In this section, we discuss the computational implementation to calculate the radial distribution function g(r) which is the average number density at a radial distance from a tagged string relative to the ideal gas case (N/L 2 ). The plot of this function reveals the phase at which the system is [30]. In the case of an ideal gas, it corresponds to a flat function, indicating that the system is well distributed and that any disk may be at any place. Another situation corresponds to a diluted system composed of interacting disks. In this case, g(r) shows a global maximum around the distance λσ as a consequence of the (short-distance) repulsive interaction. Nevertheless, as r increases, g(r) becomes flat due to the low number density. Finally, if the number density increases, the distance between the strings becomes smaller, the particle interactions are more frequent, and the radial distribution function starts to oscillate. This behavior occurs because the particles settle around each other. From the perspective of a tagged particle, it looks to be surrounded by a first shell of particles at a distance equivalent to the repulsive interaction range. They are the nearest neighbors. Subsequently, if the number density is high enough, a second shell could take arise: the nextto-nearest neighbors. Thus, particle interactions give rise to a structured system, identified as a liquid (observed in x-ray scattering [41] and non-vibrating magnetic granular systems [42] experiments).
In the computational implementation, we use the random sequential addition algorithm [43] to add disk by disk, with the corresponding (λ, q λ ). The simulation process is stopped according to the observable under calculation.
To determine g(r) we add until 200 strings (if possible) distributed into a square box of size L = 8σ and randomly assigned as soft or hard according to q λ . The first added string is allocated in the geometrical center of the box and it is taken as a trial disk. Then, g(r) is calculated as g(r) = n(r)L 2 N (2π∆r(r + 0.5∆r) + π∆r 2 ) , where n(r) is the average number of strings at a distance between r and r + ∆r from the trial disk [33,34]. It is computed over 1 × 10 6 simulated systems, starting with N = 5 and increasing it in steps of ∆N = 5, and ∆r = 0.035. We classify the pair values (λ, q λ ) for each η according to g(r) and its derivative with respect to r, g ′ (r), which is calculated using the five-point stencil method (with a spacing between points of ∆r). The classification criteria are as follows (cf. Ref. [30]): 1. Ideal gas: if g(r) is a flat function.
2. Non-ideal gas: if g(r) has a global maximum around r 1 = λσ and g ′ (r) has a sustained negative trend for λ < r/σ < 2λ.
3. Liquid-like structure: if the system has already shown the non-ideal gas structure for some η, and if for a higher density, g(r) has a local minimum before 2λ, and g ′ (r) has a sustained positive trend for 2λ < r/σ < 2.5λ.
The general characteristic contained in criteria 1, 2, and 3 are based on the behavior of g(r) for classical fluids [30]. In particular, the details of the third point have been set by observing the oscillation behavior of g(r) for the hard-sphere case (λ = 1, q λ = 0). This is relevant since it corresponds directly with the classical fluid of hard spheres and the oscillations of the radial distribution functions are the indication that the liquid-like structure has taken place [30]. The first maximum is originated by the next nearest neighbors around the tagged particle, the first coordination shell, which for the case of a hard string is at least to a radial distance λσ. Meanwhile, the second maximum corresponds to the location of the next-to-nearest neighbors. We do this because for other values of (λ, q λ ) we do not expect to exactly reproduce the radial distribution function for a classical fluid since, in this work, g(r) is computed as the average over system simulations where the tagged string may be soft or hard. In Fig. 2 we show some examples of this classification for different pairs (λ, q λ ).
We analyze g(r) in a segment of length λσ, but as r takes discrete values, it is necessary to establish a condition on how many consecutive points we are going to check. We select five points; this implies that a transition from an ideal gas to a non-ideal gas is detected if λ ≥ 0.2. The classification begins with the pair (λ = 1, q λ = 0) (hard disk fluid limit) and we search if the system shows a liquid-like structure for some η. Taking into account that the systems should approach to the ideal gas case as λ decreases and q λ increases, we determine that if for some (λ 0 , q λ0 ) the system does not manifest the liquidlike structure for all η, then for (λ, q λ0 ) with λ < λ 0 or (λ 0 , q λ ) with q λ > q λ0 it will not show the liquid-like structure anymore. We use the same considerations for the transition from a non-ideal gas to an ideal gas.
We use the Mertens-Moore method [32] to determine the percolation threshold in the CSCSPM. Unlike the determination of g(r), now we add disks until the emergence of the spanning cluster and save the number n c of added disks. Then, using the information of 10 4 simulated systems, we calculate the probability, f L (n), that a spanning cluster exists when n disks have been added. The percolation probability P L (η) at η is computed (as in Ref. [32]) by convoluting the f L (n) probability with the Poisson distribution with average α = ηL 2 /πr 2 0 , for several η values around the maximum of the distribution of η c . Then, each data set is fitted to the function where η cL is the estimated percolation threshold of the system of size L, and ∆ L is the amplitude of the transition region [44,45]. To take into account the finite-size effects on η cL , we perform simulations with systems of size L=12, 16, 24, 32, 48, 64, and 96. Furthermore, we determine the percolation threshold in the thermodynamic limit, η c , through the finite-size scaling η c −η cL ∝ ∆ −1/ν L . Here, ν is the exponent associated with the scaling of the correlation length of the cluster size, which is calculated using ∆ L ∝ L −1/ν [45,46]. From the analysis of ∆ L (as a function of L), we find ν ∼ 4/3 for all considered pairs (λ, q λ ). This is in good agreement with the results of 2D percolating systems [46].
To measure the area covered by strings we draw a square grid with spacing L/20. Then, we count the cells whose center is inside of at least one disk. In this way, the area is approximated as A L = N L 2 /400, where N is the number of counted cells. We measure the area in the conditions: a) the density of the non-ideal gas to liquid-like transition, and b) the emergence of the spanning cluster. In particular, for situation b), it is possible to compute the area in the thermodynamic limit. In this case, we determine the mean area A L (n c ) that covers the n c strings. Then, the area covered by strings in the percolation threshold, φ(η c ), is calculated as the convolution of A L (n c ) with the Poisson distribution with average α c = η cL L 2 /πr 2 0 . The finite-size effects on φ L (η c ) satisfy φ(η c ) − φ L (η cL ) ∝ L −1/ν [this relation is no longer valid as (λ, q λ ) → (1, 0), this zone does not belong to the cases discussed above].

IV. RESULTS
To analyze the deviation between the gas-liquid transition temperature [T * c := T (η * c )] and the critical transition temperature [T c := T (η c )] for a given (λ, q λ ), we define where η * c is the minimum density at which the system shows a liquid-like structure, while η c is the percolation threshold. Notice that: a) τ is independent of p 2 T 1 [which may dependent on (λ, q λ ) and on the size of the system], b) it only depends on the critical filling factors, η c and η * c , and the corresponding covered area, and c) it can be interpreted as the relative deviation between T * c and T c .
As we have a discrete set in the pairs (λ, q λ ), we interpolate τ with cubic splines to determine the regions wherein it takes values less than τ 0 = 0.005, 0.02, 0.05, 0.1. These results allow us to determine the regions where T * c is bounded as Figure 3 contains the counter lines of τ values discussed above, together with the obtained phase diagram for the pairs (λ, q λ ) according to: i) an ideal gas if the system shows the ideal gas structure for all inspected η, ii) a non-ideal gas if the system only shows a transition from the ideal gas to the non-ideal gas case, and iii) a liquid-like structure, if the system presents both ideal gas to non-ideal and non-ideal to liquid-like transitions.  Let us make contact with the experimental data. The critical temperature at which the QGP is formed takes values between 150 and 170 MeV for zero chemical potential. In this case, for τ < 0.005 we have |T c − T * c | < 1 MeV, which means that both temperatures are very close. Then, our model predicts the liquid-like structure of the QGP (see Fig. 3), which is in agreement with the experimental observations that suggest the liquid behavior of the QGP. In Fig. 4 we show an example of the mentioned systems. . This corresponds to the pair (λ=0.7, q λ =0.25) inside the region where τ <0.005, with a filling factor near to the percolation threshold (η=0.861).

V. CONCLUSIONS
In summary, we have presented a percolation-based string model that incorporates a repulsive interaction. In this model, the systems can have three structures: ideal gas, non-ideal gas, and liquid-like. Our main result is the existence of systems that simultaneously show a gasliquid transition and the emergence of a spanning cluster. Then, our model describes the experimentally observed liquid behavior of the QGP. We must remark that the inclusion of only the core region to generate hard-core systems (q λ =0) is not enough to find the aforementioned systems. In our study, we used the color string percolation model, but the main results could be obtained in most of the string models if a repulsive interaction is incorporated properly.
In the simulations, we observed conglomerations of hard or soft strings that act as "droplets" and "bubbles," respectively. This is a consequence of a jamming-like effect produced by the hard strings. Then, there exists a particular filling factor from which can only be added soft strings in the system. This could be an explanation of the second rise of ε/T 4 at T ∼1.3-1.4 T c reported in Ref. [47]. In these kinds of systems, there would be an excess of soft strings at high η values, which could indicate a new transition from liquid to gas and subsequently to an ideal gas. In the CSPM, at this temperature, the mean distance d between the overlapping strings is smaller than the string radius (computed like in Ref. [48] for 2D systems, d/r 0 ∼0.8, 0.9). For models that consider strings with cores, like our model, this means that the color field of the strings penetrates the core region until λ ∼0.4, 0.45, starting to see undressed color sources.
Finally, several questions remain open, for instance: •Notice that we worked with τ , instead of calculating T c and T * c . For the temperatures, it is mandatory to determine F (η) and p 2 T 1 . With this information, it is possible to derive all of the phenomenology associated with the CSCSPM as a function of (λ, q λ ). Even when we do not expect large changes, it would be interesting if the results shifted toward the experimental data.
•Since the third harmonic of the azimuthal distribution, v 3 , is very sensitive to the fluctuations of the string locations, we expect significant effects in the implementation of this model. In particular, as the relative weight of the edge is larger in small systems like pp or pA collisions than in heavy-ions collisions, the effects on v 3 can be larger.
•Our model opens the possibility of incorporating other interactions between strings. For example, it is well known that the strings can interact through a Yukawatype potential, where the correlation length can be associated with the saturation scale R s = 1/Q s in the color glass condensate context.
•To simulate a more realistic parton distribution we can implement different profiles in the CSCSPM, such as the Gaussian or Woods-Saxon distributions.