Statistical Behavior of Lepton Pair Spectrum in the Drell-Yan Process and Signal from Quark-Gluon Plasma in High- Energy Collisions

We analyze the transversemomentum (pT) spectra of lepton pairs (l l) generated in the Drell-Yan process, as detected in proton-nucleus (pion-nucleus) and proton-(anti)proton collisions by ten collaborations over a center-of-mass energy ð ffiffiffiffiffiffiffi sNN p or ffiffi s p if in a simplified form) range from ~ 20GeV to above 10TeV. Three types of probability density functions (the convolution of two Lévy-Tsallis functions, the two-component Erlang distribution, and the convolution of two Hagedorn functions) are utilized to fit and analyze the pT spectra. The fit results are approximately in agreement with the collected experimental data. Consecutively, we obtained the variation law of related parameters as a function of ffiffi s p and invariant mass ðQÞ. In the fit procedure, a given Lévy-Tsallis (or Hagedorn) function can be regarded as the probability density function of transverse momenta contributed by a single quark (q) or anti-quark ( q). The Drell-Yan process is then described by the statistical method.


Introduction
There are more than one processes that can generate a pair of charged leptons (ℓ ℓ) in experiments of high-energy collisions. In 1970, Sidney Drell and Tung-Mow Yan firstly proposed ℓ ℓ production in a high-energy hadron scattering in a process we now call the "Drell-Yan" process [1]. In this process, a quark (q) in one hadron and an antiquark ( q) in another hadron are annihilated, and a virtual photon (γ * ) or Z boson is generated, which then decays into ℓ ℓ. This process is expressed as A + B ⟶ γ * /Z + X ⟶ ℓ + ℓ + X, where A and B are collision hadrons and X denotes other particles produced in the collisions. The Drell-Yan process has been extensively studied experimentally, theoretically, and phenomenologically.
Several phenomenological interpretations of experimental Drell-Yan data collected in the previous many years have been released by various groups, particularly in recent years where first extractions of quark transverse momentum distributions are becoming available from highly accurate theoretical descriptions of QCD perturbative ingredients. The factorization theorem for the Drell-Yan process allows to write the transverse momentum differential cross-section as a convolution of two transverse momentum-dependent (TMD) parton distribution functions (PDFs), which are, under certain conditions, very complicated. This complicated factorization involves soft factors that resum soft gluon radiation regularizing a certain class of divergences that arise in the theoretical formulae. The soft gluon resummation is especially important in the description of the quark-gluon plasma (QGP), where ℓ ℓ can be produced in a process similar to that of Drell-Yan but with different origins of quarks.
QGP is a new form of matter which is created in the central region of high-energy nucleus-nucleus collisions, where extreme density and high-temperature environment are developed. It has become one of the important areas of research in the field of nuclear and particle physics. The gradual maturity of QCD and gauge field theory provides a powerful explanation for this novel matter and phenomenon. In fact, QGP is particularly short-lived. In QGP, a quark q and anti-quark q can soon be annihilated into a virtual photon γ * or Z boson, which then decays to a pair of leptons ℓ ℓ. This happens in the QGP degeneration process in which most particles are produced. The yield, invariant mass, rapidity (y), and transverse momentum (p T ) distribution of ℓ ℓ depend on the momentum distribution of q q and gluons in QGP in the collision region. Therefore, the information of ℓ ℓ can be used to judge whether QGP is generated and further to study its thermodynamic status making the ℓ ℓ production one of the most important signals generated by QGP. Consequently, the study on ℓ ℓ becomes particularly critical.
From the above, it is clear that ℓ ℓ can be produced in high-energy hadronic/nucleonic collisions in two main ways: the Drell-Yan process and QGP degeneration process. To study the properties of QGP, we should remove the influence of the Drell-Yan process and vice-versa. Generally, we may use the same methodology to describe the two processes. At present, one mainly uses the statistical method to study the properties of QGP. Correspondingly, we may also use the statistical method to study the production of ℓ ℓ in the Drell-Yan process, especially because the factorization theorem is very hard to model. In short, the statistical description for the Drell-Yan process is necessary to better understand the properties of QGP.
The measurement of lepton-pair physical quantities (including energy, p T , and y) in experiments studying the Drell-Yan process provides lots of valuable information about the dynamic properties and evolution process of the produced particles. In particular, p T is Lorentz invariant in the beam direction and can be used to describe the particles' motion and system's evolution. There are different functions that can be used to describe the ℓ ℓp T spectra in statistics. For example, we can use the Lévy-Tsallis function [26][27][28][29][30], the (two-component) Erlang distribution [31][32][33], and the Hagedorn function [34,35] to fit the experimental data to obtain the analytical parameters of the p T spectrum. Since the Drell-Yan process is the result of the interactions of q and q, we can use the convolution of two functions to describe the p T spectra. The idea of convolution is concordant to the factorization theorem for the Drell-Yan process.
In this paper, we use three functions to fit and analyze the Drell-Yan ℓ ℓp T spectra obtained by ten collaborations from the experiments of high-energy proton-nucleus (pion-nucleus) and proton-(anti)proton collisions. These experimental studies provide a great resource for us to better understand the collision mechanism and dynamic characteristics of the mentioned process.

Formalism and Method
Naturally, the p T spectra of Drell-Yan ℓ ℓ depend on collision energy. For that reason, we should use different probability functions to study these spectra at different energies. Here, we briefly describe the three functions which will be used in this study. In the following, p t1 and p t2 are the transverse momenta of the two quarks, and p T is the transverse momentum of the two quark system, which equals the transverse momentum of the dilepton system at leading order.
2.1. The Lévy-Tsallis Function. The Boltzmann distribution is the most important probability density function in thermodynamic and statistical physics. We present the probability density function of p T as a simple Boltzmann distribution [36][37][38]: where N is the number of identical particles of mass m 0 produced in the collisions, C B is the normalization constant, and T B is the effective temperature of the collision system. The Boltzmann distribution is a special form of the Tsallis distribution, and the latter has a few alternative forms [26][27][28][29][30]. As one of the Tsallis distribution and its alternative forms, the Lévy-Tsallis function of the p T spectrum of hadrons [26][27][28][29][30] is used in this work. We have the following form to describe the transverse momentum (p t ) distribution of (anti-)quark: where N q is the normalization constant, T and n are the fitted parameters, and m q is the mass of q or q taking part in the reaction. In general, we use m u = m d = 0:3GeV/c 2 in the Drell-Yan process because q or q is from the participant hadrons. The same m u or m d is for sea quarks and those in baryons, where the sea quarks of higher mass are not considered in this work. In QGP, m u = 0:003 GeV/c 2 and m d = 0:007GeV/c 2 because the quarks are approximately bare [39]. It has been verified that the Tsallis distribution is just a special case of the Lévy distribution, but not the opposite [30].
2.2. The (Two-Component) Erlang Distribution. The Erlang distribution [31][32][33] is proposed to fit the p T spectra in the multisource thermal model [40]. Generally, a twocomponent Erlang distribution [31][32][33] is used to describe both the soft and hard processes. The contribution fractions of the two components are determined by fitting the experimental data. The numbers of parton sources participating in the soft and hard processes are represented by n S ≥ 2 and n H = 2, respectively. The contribution 2 Advances in High Energy Physics (p t ) of each parton source to p T of final-state particle is assumed to obey an exponential function: where hp t i represents the average p t contributed by the i-th source. Because hp t i is the same for different sources, the index i in hp t i i is omitted. The p T distribution contributed by n S (n H ) sources is the convolution of n S (n H ) exponential functions, which gives the Erlang distribution. Let k denote the contribution fraction of the first component (soft process). The two-component Erlang distribution is Fitting the data with the two-component Erlang distribution, we can get the changes of parameters hp t i S , hp t i H , and k.
We should discuss the values of n S and n H further. If n S > 2, the participant partons are expected to be q q and n S − 2 gluons in the soft or nonviolent annihilation process. Considering that the probability of multiparton participating together in the process is low, we have usually n S = 2 or 3 in this work. Generally, the larger the n S , the sharper the distribution peak. In many cases, n S = 3 means that q q and a gluon participate in the soft process. For all cases, n H = 2 (always true in this work) means that only q q participates in the violent annihilation in the hard process.
2.3. The Hagedorn Function. The Hagedorn function is an inverse power law [34,35] which is an empirical formula derived from perturbative QCD. Generally, this function can only describe the spectra at large p T , but not the entire p T interval. In the case of using the Hagedorn function in a wide range of p T , the probability density function of p T can be expressed as where A is the normalization constant, and p 1 and n 1 are the fitted parameters. The final-state particles with high momenta are mainly produced by the hard scattering process during the collisions. However, both the soft and hard processes contribute to the p T spectra. In some cases, the soft excitation process in the low p T range can also be described by the Hagedorn function. We try to use the Hagedorn function to describe the transverse momentum distribution of (anti-)quarks. That is, we may use p t instead of p T in Eq. (5) to obtain the transverse momentum distribution of (anti-)quarks, which is a new form of Eq. (5) and will be used in the following section.
We have tested the Hagedorn function with different revisions in which p T /p 1 in Eq. (5) is replaced by p 2 T /p 2 1 , A p T is replaced by Ap 2 T /m T , or p T /p 1 is replaced by p 2 T /p 2 1 , and Ap T is replaced by A, where p 1 and A vary in different revisions. The uses of the revised Hagedorn functions result in some overestimations in low (or high) p T region comparing to the Hagedorn function. Contrarily, these revisions result in some underestimations in high (or low) p T region due to the normalization. The revisions of the Hagedorn function are beyond the focus of the present work, and we shall not discuss them further.
2.4. The Convolution of Functions. The convolution of functions is an important operation process in functional analysis that can be used to describe the weighted superposition of input and system response (that is, two subfunctions). The Drell-Yan process is the result of the interactions of q and q in high-energy collisions, which means that we need the convolution of two functions to describe this process. Indeed, the above Eq. (2) or (5) can be used to describe the transverse momentum distribution, f 1 ðp t1 Þ, of a single (anti-)quark's contribution. The second (anti-)quark's contribution is where p T is still the transverse momentum of the ℓ ℓ system. So, the convolution of two probability density functions should be used to describe the p T spectrum of ℓ ℓ in the Drell-Yan process. We have the convolution of two Eq. (2) or (5) to be expressed as where It should be noted that the total transverse momentum before (of the two quarks system) and after (of the two leptons system) are equal. The assumption that the total transverse momentum is equal to the sum of the scalar transverse momenta of the two partons that is for the particular case in which the vectors p ! t1 and p ! t2 are parallel. Our recent works [41,42] show that this assumption is in agreement with many data. Naturally, we do not rule out other assumptions such as the particular case in which p ! t1 and p ! t2 are perpendicular and the general case which shows any azimuth for p ! t1 and p ! t2 . The particular case used in this work is more easier than other particular or general case in the fit to data. We are inclined to use the parallel case.
As a legitimate treatment, the convolution formula Eq. (6) can be used to fit the p T spectrum of ℓ ℓ in the Drell-Yan process, where f 1 ðp t1 Þ and f 1 ðp T − p t1 Þ are from empirical guess which is simpler than the factorization theorem based on perturbative QCD. On one hand, Eq. (6) can reflect the weighted contribution of the transverse momentum of each (anti-)quark to the p T spectrum in the process. On the other hand, Eq. (6) can also reflect the system in which two main participants take part in the interactions. Using the convolution to fit the data is a good choice for us, which allows us to 3 Advances in High Energy Physics more accurately understand the interaction process and mechanism between interacting partons, more completely describe the energy dependence and interdependence of the function parameters, and further better analyze the p T spectrum.

Results and Discussion
3.1. Comparison with Data. Figure 1 shows the p T spectra of ℓ ℓ [with different invariant masses (Q) or Feynman variables (x F )] produced by the Drell-Yan process in different collisions at different energies (with different integral luminosities if available in literature), where the concrete type of ℓ ℓ is also given in each panel. The symbols E and σ on the vertical axis denote the energy of ℓ ℓ and the cross-section of events, respectively. Among them, the data points presented in  Figure 1(f). Different collaborations have different intervals of Q and x F , while the detailed binning information is marked in the panels. In some cases, the range of rapidity y is not available due to other selection conditions such as the complex polar coverages and sensitivities of detector components or Feynman variable being used [45]. The total experimental uncertainties are cited from refs. [43][44][45] which include both the statistical and systematic uncertainties if both are available. The solid, dashed, and dotted curves in all panels are the results of our fittings with the convolution of two Lévy-Tsallis functions, the two-component Erlang distribution, and the convolution of two Hagedorn functions, respectively. The histograms in this and following figure correspond to QCD calculations which will be discussed later. We use the minimum-χ 2 to evaluate the goodness of the fits, where χ 2 = ∑ j ½ðData j − Fit j Þ 2 /Uncertainty 2 j and j are for the j-th data. We list the results of the fits (parameters), the χ 2 and the number of degrees of freedom (ndof) in Table 1. The numbers n S = 3 and n H = 2 which are not listed in the table to avoid trivialness. One can see that the three functions can approximately describe the p T spectra of ℓ ℓ produced by the Drell-Yan process in high-energy p-Cu and π − -W collisions. The two-component Erlang distribution describes better than the other two functions. Figure 2 shows the p T spectra of ℓ ℓ (with different invariant mass Q) generated by the Drell-Yan process in protonproton (p-p or pp) collisions and measured by three different collaborations. The data points in Figure 2(a) are from the experimental results measured by the R209 Collaboration [46]. The collision energy is ffiffi s p = 62GeV, and the Q range is 5 − 8GeV/c 2 . The data points in Figure 2(b) show the experimental results from the PHENIX Collaboration [47]. The ffiffi s p is 200 GeV, the Q range is 4:8 − 8:2GeV/c 2 , and the rapidity range is 1:2 < y < 2:2. The data points in Figure 2(c) are from the experimental results of the STAR Collaboration [24]. The ffiffi s p is 510 GeV, the Q range is 73 − 114GeV/c 2 , and the rapidity range is |y | <1. In some cases, the range of rapidity y is not available due to other selection conditions being used [46]. The curves of the three fits are also shown. and the extracted parameters are presented in Table 2. One can see that the three functions can approximately describe the p T spectra of ℓ ℓ produced by the Drell-Yan process in high-energy pp collisions.
Similar to Figure 2, Figure 3 shows the p T spectra of ℓ ℓ produced by the Drell-Yan process in proton-antiproton In some cases, the range of rapidity y is not available due to other selection conditions being used [48][49][50]. The curves of the three fits are also shown, and the extracted parameters are presented in Table 2. One can see that the three functions can approximately describe the p T spectra of ℓ ℓ produced by the Drell-Yan process in high-energy p p collisions.  Table 2. One can see that the three functions can approximately describe the p T spectra of ℓ ℓ produced by the Drell-Yan process in pp collisions at ultrahigh energies.
From the above comparisons, we see that some of the data sets are relatively poorly described by the fit. Notably, the FNAL-615 data, and the CDF data at 1.96 TeV, are not very well fitted. Specifically, the fit of the CDF data returns a poor χ 2 , unlike the other Tevatron experiments, as all LHC data sets have reasonably good χ 2 . We would like to point out that the relatively poor χ 2 of some data is understandable due to the fact that the fit function works well in most cases and does not work well in a few cases. In addition, the data at high transverse momentum has low statistics which causes large dispersion of the data from the fit. In other cases, the low statistics happen to be well fitted. To improve 4 Advances in High Energy Physics     5 Advances in High Energy Physics the relatively poor χ 2 of some fits, a better fit function could be used in the future; although, it is not feasible to find a fit function that works well in all cases. One could improve the quality of the fits by concentrating in the high-statistic regions or by rebinning the low-statistic data.

Tendency of Parameters.
In order to study the underlying law that governs the variation of the extracted parameters with ffiffi s p and Q, we preset the values of the fitted parameters as a function of the these two quantities.
x F = 0.5-0.6(×10 29 ) x F = 0.6-0.7(×10 28 )  5) and (6)], respectively, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations. 6 Advances in High Energy Physics     , there is a boundary between 60 and 500 GeV for hp t i S , between 500 and 1100 GeV for hp t i H , and between 200 and 500 GeV for k. Meanwhile, with increasing Q, hp t i S and hp t i H increase slowly and then quickly, and k decreases slowly and then quickly. There is a boundary at Q ≈ 20-60 GeV/c 2 .
In Figure 7, we show the parameters from the fits using the Hagedorn function and their dependence on collision energy  . The data points are quoted from the R209 (a) [46], PHENIX (b) [47], and STAR Collaborations (c) [24]. In panels (a) and (b), ℓ ℓ is μ + μ − . The solid, dashed, and dotted curves are our results of fitting the data points with the convolution of two Lévy-Tsallis functions [Eqs. (2) and (6)  Advances in High Energy Physics   Firstly, this is because the grouping of the quality in Figure 1(e) is different from the grouping of other data, and there are already many other groupings. Secondly, Figure 1(f) analyzes the p T spectra within different ranges of Feynman variables, which is different from others in terms of event sample. To avoid trivialness, we have not put the fitting results of Figures 1(e) and 1(f) in Figures 5-7; though, these results are also shown in Table 1. We find that, by d /dp

12
Advances in High Energy Physics analyzing these results, they do not contradict the trend of other results presented in Figures 5-7. The parameters T, hp t i S , hp t i H , and p 1 show monotonous increasing trend when ffiffi s p and Q increase. Naturally, the variation degrees are different. These increasing trends reflect that these parameters describe the violent degree of collisions between two (anti-)quarks in the Drell-Yan process. As the contribution fraction of the first component in the two-component Erlang distribution, k decreasing with increasing ffiffi s p and Q reflects naturally the increase of the contribution fraction of the second component. The parameter n increases firstly and then decreases, and the parameter n 1 decreases, with the increase of ffiffi s p and Q. This variation implies the change of interaction pattern and strength. A possible explanation is that the collision centrality between the two (anti-)quarks changes from periphery to center, or the violent degree of collisions increase, when ffiffi s p and Q increase. We should pay more attention on this variation in a future study.
In particular, the energy range considered in the present work is enough for the formation of QGP. It is possible that the p T spectra of ℓ ℓ in the Drell-Yan process are affected by QGP. The nonmonotonous change of n implies the maximum influence between Drell-Yan and QGP. According to the Tsallis statistics, n is opposite to the entropy index q because n = 1/ðq − 1Þ. The maximum n at~40 GeV implies that q is the closest to 1 at this energy. Meanwhile, at this energy, the system is the closest to the equilibrium. The soft-est point (~40 GeV) of equation of state from the excitation function of q is consistent with our very recent work [61] in which we stated that "the onset energy of the partial phase transition from hadron matter to QGP is 7.7 GeV and that of the whole phase transition is 39 GeV" with the uncertainty of 1-2 GeV.
We explain further the softest point here. In the energy range below~40 GeV, there is not enough volume for the interactions between Drell-Yan and QGP due to partial phase transition [61]. In the energy range above~40 GeV, there is not enough time for the interactions between Drell-Yan and QGP because the participants penetrate through each other quickly. Insufficient volume and time result in the system being not the closest to the equilibrium, though the system is still close to equilibrium at the softest point, Q ≈ 11-18 GeV/c 2 . At below (above) the softest point, we see naturally smaller (larger) Q, except for Q = 5-

13
Advances in High Energy Physics and n 1 show a trend of decrease or saturation. These obvious characteristics in parameters make a clear distinction between approaching to and far away from the onset energy of phase transition from hadron matter to QGP because QGP also affects the spectra of ℓ ℓ in the Drell-Yan process.    6), it reappears through the dependence of its parameters upon Q (and also the center-of-mass energy squared s). The dependence of the parameters on Q reflects the fact that the function f 1 ðp t1 Þ or f 1 ðp T − p t1 Þ must change its parameters with the hard scale. This result is consistent with the appropriate frameworks: collinear framework for large transverse momentum and TMD framework for small transverse momentum [2].
In addition, the formulae for the evolution of TMD parton distributions are rather complicated, and the resulting effect is difficult to parameterize with simple expressions. The present work is a preliminary attempt for the purpose of parameterization with simple expressions. Meanwhile, as a preliminary attempt, the present work is not accurate in some cases which show large χ 2 /ndof. In other cases, χ 2 /ndof is approximately 1 or not too large. In any case, we firmly believe that the underlying physics law is knowable and simple. A very complicated expression for the transverse momentum spectra of ℓ ℓ in the Drell-Yan process in highenergy collisions is not the final one. More work for simplifying the expressions is needed in the future.
After describing the spectra of ℓ ℓ in the Drell-Yan process, it is possible to subtract the contribution of the Drell-Yan process from the spectra of ℓ ℓ in the final state and leave behind only the contribution of QGP. Then, one may study the excitation function of related parameters from the spectra of ℓ ℓ contributed purely in QGP conditions and search for the softest points of equation of state from the excitation function. The energies corresponding to the softest points are expected to connect with the critical energy of phase transition from hadron matter to QGP. If both the spectra of ℓ ℓ in the Drell-Yan process and in QGP degeneration process are described by simple functions, one can study the phase transition more conveniently. We hope that the present work is significant in methodology for the study in the future.
We now discuss the results as they pertain to QCD calculations [23,24,46,55,[62][63][64]. To study more deeply, and in a visual way, the connection between our formalism and the standard QCD resummation, we present a direct comparison in Figures 1-4 as examples. As can be seen, the histograms presented in Figures 1-4 represent the results with different treatments based on QCD calculations, which will be discussed separately in the following.
The histograms in Figures 1(a)-1(d) are directly quoted from those in ref. [23] in which the next-to-next-to-leading order (NNLO) calculation is performed for nonperturbative structure of semi-inclusive deep-inelastic scattering and the Drell-Yan scattering at small transverse momentum within TMD factorization. The histograms in Figures 1(e) and 1(f) are directly quoted from those in ref. [62] in which the NNLO calculation is performed for the Drell-Yan processes within TMD factorization. In the figure, some histograms are renormalized to the data to give a better comparison.

Advances in High Energy Physics
The histogram in Figure 2(a) is transformed from the curve in ref. [46] to fit the style of other panels and to give a clear display. In the transformation from the curve to histogram, the areas under the histogram and curve in a given transverse momentum bin are kept being the same. In ref. [46], the calculation of QCD convoluted with the Gaussian function is used. The histogram in Figure 2(b) is directly quoted from that in ref. [63] in which the transverse momentum spectrum of low invariant mass Drell-Yan is produced at next-to-leading order (NLO) in the parton branching method. The histogram in Figure 2(c) is directly quoted from that in ref. [24] in which the TMD parton distributions are   17 Advances in High Energy Physics considered up to next-to-next-to-next-to-leading order logarithmic (N 3 LO or N 3 LL) from the Drell-Yan data.
The meanings of histograms in Figure 3 and Figure 4(a), Figure 4(b), Figure 4(c) with 46 ≤ Q < 66 and 116 ≤ Q < 150 GeV/c 2 , Figure 4(e), and Figure 4(f) are the same as Figures 1(a)-1(d), which will not be discussed anymore. The histogram in Figure 4(c) with 66 ≤ Q < 116 GeV/c 2 is directly quoted from that in ref. [64] in which the calculation is performed due to the Drell-Yan transverse momentum resummation of fiducial power corrections at N 3 LL. The histograms in Figure 4(d) are directly quoted from those in ref.
[55] in which the transverse momentum distributions of the Drell-Yan ℓ ℓ are calculated up to NNLO + N 3 LL.
From the above descriptions, one can see that the formalism of this paper is flexible in the fit to data. We may compare the fits with more QCD predictions, and it does not matter with or without the TMD PDFs. Except for several papers in which the Drell-Yan transverse momentum spectrum is predicted based on QCD factorization are referenced in the introduction [23,24] and the above discussion [46,55,[62][63][64], more works related to the PDFs or TMDs in QCD are available recently [65][66][67][68][69][70][71][72][73][74]. These QCD factorizations are complex in the calculation and show different forms of formalization from this paper but result in similar shapes for the dilepton spectra as observed in experiments. By adjusting the parameters, the formalism of this paper can flexibly fit the QCD predictions with or without the TMD PDFs. Compared to predictions from collinear PDFs, the formalism of this paper is closer to TMD PDFs due to the flexible parameter selection in the fit. The PDFs or TMDs in QCD and other QCD-based analyses reveal the dynamic aspect of the particle production process. The formalism used by us reflects the statistical behavior of the produced particles.

Summary and Conclusions
We have studied the transverse momentum spectra of lepton pairs generated by the Drell-Yan process in p-Cu, π − -W, and pp (p p) collisions over an energy range from~20 GeV to above 10 TeV. The low-energy data come from the E288, E605, R209, PHENIX, and STAR Collaborations. The highenergy data come from the CDF, D0, ATLAS, CMS, and LHCb Collaborations. The invariant mass range of the final-state particles produced in the collisions also has a large span of 4 < Q < 150GeV/c 2 . Three types of probability density functions are used to fit and analyze the collected experimental data. All the three functions are approximately in agreement with the experimental data and the QCD calculations. Some parameters are obtained.
In the convolution of two Lévy-Tsallis functions, as increasing ffiffi s p , there is a knee point for the trend of n at ffiffi s p ≈ 40-50 GeV. Meanwhile, there is a boundary at ffiffi s p ≈ 200-500 GeV above which T increases significantly.
With the increase of Q, there is a knee point for the trend of n at Q ≈ 14-15 GeV/c 2 . Meanwhile, there is a boundary at Q ≈ 20-60 GeV/c 2 above which T increases significantly. In the two-component Erlang distribution, there is a boundary at ffiffi s p ≈ 60-500 GeV for hp t i S , 500-1100 GeV for hp t i H , and 200-500 GeV for k, above which hp t i S , hp t i H , and 1 − k increase quickly. Meanwhile, there is a boundary at Q ≈ 20-60 GeV/c 2 above which hp t i S , hp t i H , and 1 − k increase quickly. In the convolution of two Hagedorn functions, p 1 increases obviously with increasing ffiffi s p and Q. There is a boundary for the trend of n 1 at ffiffi s p ≈ 40-200 GeV and Q ≈ 20-50 GeV/c 2 .
With increasing ffiffi s p and Q, the parameters T, hp t i S , hp t i H , 1 − k, and p 1 show monotonous increasing trend. These increasing trends reflect that these parameters describe the violent degree of collisions between quark and antiquark in the Drell-Yan process. The parameter n increases initially and then decreases, and the parameter n 1 decreases, with the increase of ffiffi s p and Q. This variation implies the change of interaction pattern and strength. A possible explanation is that the collision centrality between quark and antiquark changes from periphery to center, or the violent degree of collisions increases, when ffiffi s p and Q increase. The fit results in this paper are comparable to the QCD NNLO and N 3 LL results with TMD PDFs.
In conclusion, the large number of events collected by the investigated experiments allows us to study the statistical behavior of dilepton production in hadron and nuclei collisions. The lepton pairs can be produced through the Drell-Yan process and the QGP degeneration process. Both processes are predicted by QCD and are described well enough by our statistical thermodynamics models, especially in kinematic regions with sufficiently large numbers of events.

Data Availability
The data used to support the findings of this study are included within the article and are cited at relevant places within the text as references.

Ethical Approval
The authors declare that they are in compliance with ethical standards regarding the content of this paper.

Disclosure
The funding agencies have no role in the design of the study; in the collection, analysis, or interpretation of the data; in the writing of the manuscript; or in the decision to publish the results.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.