#### Abstract

We analyze the transverse momentum () spectra of lepton pairs () 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 or if in a simplified form) range from GeV to above 10 TeV. 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 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 and invariant mass . 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 () or anti-quark (). The Drell-Yan process is then described by the statistical method.

#### 1. 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 () in one hadron and an antiquark () in another hadron are annihilated, and a virtual photon () or boson is generated, which then decays into . This process is expressed as , where and are collision hadrons and denotes other particles produced in the collisions. The Drell-Yan process has been extensively studied experimentally, theoretically, and phenomenologically.

The literature about the theoretical description of the Drell-Yan process within quantum chromodynamics (QCD) is well known and settled [2–9]. The framework for the description of the transverse momentum dynamics [sometimes, indicated as Collins-Soper-Sterman (CSS) formalism] is summarized in the well-known book by John Collins [2]. Some recent reviews on the subject of transverse momentum distributions in the Drell-Yan process can be found in refs. [4–6] in which many works were cited. More theoretical works at both small and large are available in literature [7–9]. At the same time, lots of phenomenological works were published in the past [10–15], recent [16–21], and very recent years [22–25].

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 and anti-quark can soon be annihilated into a virtual photon or 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 (), and transverse momentum () distribution of depend on the momentum distribution of 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, , and ) 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, 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 spectra in statistics. For example, we can use the Lévy-Tsallis function [26–30], the (two-component) Erlang distribution [31–33], and the Hagedorn function [34, 35] to fit the experimental data to obtain the analytical parameters of the spectrum. Since the Drell-Yan process is the result of the interactions of and , we can use the convolution of two functions to describe the 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 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.

#### 2. Formalism and Method

Naturally, the 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, and are the transverse momenta of the two quarks, and 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 as a simple Boltzmann distribution [36–38]: where is the number of identical particles of mass produced in the collisions, is the normalization constant, and 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–30]. As one of the Tsallis distribution and its alternative forms, the Lévy-Tsallis function of the spectrum of hadrons [26–30] is used in this work. We have the following form to describe the transverse momentum () distribution of (anti-)quark: where is the normalization constant, and are the fitted parameters, and is the mass of or taking part in the reaction. In general, we use in the Drell-Yan process because or is from the participant hadrons. The same or is for sea quarks and those in baryons, where the sea quarks of higher mass are not considered in this work. In QGP, and 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–33] is proposed to fit the spectra in the multisource thermal model [40]. Generally, a two-component Erlang distribution [31–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 and , respectively. The contribution () of each parton source to of final-state particle is assumed to obey an exponential function: where represents the average contributed by the -th source. Because is the same for different sources, the index in is omitted.

The distribution contributed by () sources is the convolution of () exponential functions, which gives the Erlang distribution. Let 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 , , and .

We should discuss the values of and further. If , the participant partons are expected to be and gluons in the soft or nonviolent annihilation process. Considering that the probability of multiparton participating together in the process is low, we have usually or 3 in this work. Generally, the larger the , the sharper the distribution peak. In many cases, means that and a gluon participate in the soft process. For all cases, (always true in this work) means that only 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 , but not the entire interval. In the case of using the Hagedorn function in a wide range of , the probability density function of can be expressed as where is the normalization constant, and and 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 spectra. In some cases, the soft excitation process in the low 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 instead of 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 in Eq. (5) is replaced by , is replaced by , or is replaced by , and is replaced by , where and vary in different revisions. The uses of the revised Hagedorn functions result in some overestimations in low (or high) region comparing to the Hagedorn function. Contrarily, these revisions result in some underestimations in high (or low) 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 and 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, , of a single (anti-)quark’s contribution. The second (anti-)quark’s contribution is , where is still the transverse momentum of the system. So, the convolution of two probability density functions should be used to describe the spectrum of in the Drell-Yan process. We have the convolution of two Eq. (2) or (5) to be expressed as where [] is shown as Eq. (2) if we use the Lévy-Tsallis function or Eq. (5) if we use the Hagedorn function.

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 and 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 and are perpendicular and the general case which shows any azimuth for and . 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 spectrum of in the Drell-Yan process, where and 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 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 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 spectrum.

#### 3. Results and Discussion

##### 3.1. Comparison with Data

Figure 1 shows the spectra of [with different invariant masses () or Feynman variables ()] 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 and on the vertical axis denote the energy of and the cross-section of events, respectively. Among them, the data points presented in Figures 1(a)–1(c) are quoted from the proton-copper (-Cu) collision experiments performed by the E288 Collaboration [43], and the collision energy per nucleon pair or if in a simplified form) is 19.4, 23.8, and 27.4 GeV, respectively. The data points shown in Figure 1(d) are the results of the -Cu collision experiment performed by the E605 Collaboration [44] at a collision energy of 38.8 GeV. For the E288 Collaboration, the invariant mass ranges from 4 to , while the corresponding invariant mass ranges of the E605 Collaboration are from 7 to . The experimental data points in Figures 1(e) and 1(f) are from negative pions () induced wolfram (W) (-W) collisions at 21.7 GeV performed by the FNAL-615 Collaboration [45]. The different symbols in Figure 1(e) represent invariant mass in the range of 4.05– with different scalings, where the units GeV/ are not shown in the panel due to crowding space. The Feynman variables range from 0 to 1 with a step of 0.1, as shown in Figure 1(f). Different collaborations have different intervals of and , while the detailed binning information is marked in the panels. In some cases, the range of rapidity 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–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- to evaluate the goodness of the fits, where and are for the -th data. We list the results of the fits (parameters), the and the number of degrees of freedom (ndof) in Table 1. The numbers and which are not listed in the table to avoid trivialness. One can see that the three functions can approximately describe the spectra of produced by the Drell-Yan process in high-energy -Cu and -W collisions. The two-component Erlang distribution describes better than the other two functions.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figure 2 shows the spectra of (with different invariant mass ) generated by the Drell-Yan process in proton-proton (- or ) 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 , and the range is . The data points in Figure 2(b) show the experimental results from the PHENIX Collaboration [47]. The is 200 GeV, the range is , and the rapidity range is . The data points in Figure 2(c) are from the experimental results of the STAR Collaboration [24]. The is 510 GeV, the range is , and the rapidity range is . In some cases, the range of rapidity 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 spectra of produced by the Drell-Yan process in high-energy collisions.

**(a)**

**(b)**

**(c)**

Similar to Figure 2, Figure 3 shows the spectra of produced by the Drell-Yan process in proton-antiproton (- or ) collisions with GeV/(Figures 3(a) and 3(b)),GeV/ (Figure 3(c)), and GeV/ (Figure 3(d)) at (Figures 3(a) and 3(c)) and (Figures 3(b) and 3(d)). The data points in Figures 3(a) and 3(b) are from the experiments of the CDF [48, 49] and D0 Collaborations [50–52], respectively. In some cases, the range of rapidity is not available due to other selection conditions being used [48–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 spectra of produced by the Drell-Yan process in high-energy collisions.

**(a)**

**(b)**

**(c)**

**(d)**

Figure 4 shows the spectra of produced by the Drell-Yan process in high-energy collisions, where the data used in the figure are all examples and there is no bias towards any specific CERN experiment. From Figures 4(a)–4(f), the event samples [with different conditions (, , and )] are shown in the panels. The data points are quoted from the experiments performed by the ATLAS (Figures 4(a)–4(d)) [53–55], CMS (Figure 4(e)) [56, 57], and LHCb Collaborations (Figure 4(f)) [58–60]. The extracted fit parameter values are listed in Table 2. One can see that the three functions can approximately describe the spectra of produced by the Drell-Yan process in collisions at ultrahigh energies.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

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 , unlike the other Tevatron experiments, as all LHC data sets have reasonably good . We would like to point out that the relatively poor 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 the relatively poor 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.

##### 3.2. Tendency of Parameters

In order to study the underlying law that governs the variation of the extracted parameters with and , we preset the values of the fitted parameters as a function of the these two quantities. Figures 5(a) and 5(b) show the trend of parameter , and Figures 5(c) and 5(d) show the trend of parameter , obtained by the fitting, using the Lévy-Tsallis function. Comparing Figures 5(a) and 5(c), we can analyze the dependence of parameters on . On average, it can be seen that as increases, the parameter increases quickly and then decreases slowly, and the parameter increases slowly and then significantly. There is a knee point for the trend of at –50 GeV. Meanwhile, there is a boundary at –500 GeV above which increases significantly. Similarly, we compare Figures 5(b) and 5(d) and analyze the variation of parameters with . It can be clearly seen that the parameter increases firstly and then decreases, and the parameter increases slowly and then significantly, with the increase of . There is a knee point for the trend of at –15 GeV/. Meanwhile, there is a boundary at –60 GeV/ above which increases significantly.

**(a)**

**(b)**

**(c)**

**(d)**

Figure 6 is similar to Figure 5, but shows the dependence of parameters (Figures 6(a) and 6(b)), (Figures 6(c) and 6(d)), and (Figures 6(e) and 6(f)) on (Figures 6(a), 6(c), and 6(e)) and (Figures 6(b), 6(d), and 6(f)) obtained from the two-component Erlang distribution. One can see that with increasing , and increase slowly and then quickly, and decreases slowly and then quickly. For , there is a boundary between 60 and 500 GeV for , between 500 and 1100 GeV for , and between 200 and 500 GeV for . Meanwhile, with increasing , and increase slowly and then quickly, and decreases slowly and then quickly. There is a boundary at –60 GeV/.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

In Figure 7, we show the parameters from the fits using the Hagedorn function and their dependence on collision energy and invariant mass. In Figures 7(a) and 7(b), one can see that the parameter has a slight tendency to increase with the increase of and . In Figures 7(c) and 7(d), one can see that the parameter decreases quickly and then slowly with the increase of ; also decreases slowly and then quickly with increase of . There is a boundary for the trend of at –200 GeV and –50 GeV/. We note that the variation of parameter is obtained from the the average values of the extracted parameters, for each or .

**(a)**

**(b)**

**(c)**

**(d)**

It should be pointed out that the values of the parameters in Figures 5–7 are all obtained by fitting the experimental data in Figures 1–4 using the convolution of two Lévy-Tsallis functions, the two-component Erlang distribution, and the convolution of two Hagedorn functions, where the values obtained from Figures 1(e) and 1(f) are not included. 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 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 analyzing these results, they do not contradict the trend of other results presented in Figures 5–7.

The parameters , , , and show monotonous increasing trend when and 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, decreasing with increasing and reflects naturally the increase of the contribution fraction of the second component. The parameter increases firstly and then decreases, and the parameter decreases, with the increase of and . 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 and 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 spectra of in the Drell-Yan process are affected by QGP. The nonmonotonous change of implies the maximum influence between Drell-Yan and QGP. According to the Tsallis statistics, is opposite to the entropy index because . The maximum at GeV implies that is the closest to 1 at this energy. Meanwhile, at this energy, the system is the closest to the equilibrium. The softest point ( GeV) of equation of state from the excitation function of 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 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 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, –18 GeV/. At below (above) the softest point, we see naturally smaller (larger) , except for –6 and 5–8 GeV/ which show below the softest point at 20–30 GeV and also show above the softest point at 60–200 GeV, which should be studied in the future.

At the energy below 200 GeV and the invariant mass below 20 GeV/, the system is close to the onset energy of phase transition. The parameters , , , , and show similar small and almost flat values, respectively, and and show different trends. At the energy of TeV and the invariant mass around 100 GeV/, the system is far away from the onset energy of phase transition. The parameters , , , , and show large values, and and 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.

##### 3.3. Further Discussion

In the exact factorization formula, the true TMD parton distributions depend on the hard scale of the process (for the Drell-Yan process, the invariant mass of ) and obey renormalization group equations. Although this scaling is not included in the function or of Eq. (6), it reappears through the dependence of its parameters upon (and also the center-of-mass energy squared ). The dependence of the parameters on reflects the fact that the function or 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 /ndof. In other cases, /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 high-energy 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–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.

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 considered up to next-to-next-to-next-to-leading order logarithmic (LO or LL) from the Drell-Yan data.

The meanings of histograms in Figure 3 and Figure 4(a), Figure 4(b), Figure 4(c) with and GeV/, 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 GeV/ 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 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 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–64], more works related to the PDFs or TMDs in QCD are available recently [65–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.

#### 4. Summary and Conclusions

We have studied the transverse momentum spectra of lepton pairs generated by the Drell-Yan process in -Cu, -W, and () collisions over an energy range from GeV to above 10 TeV. The low-energy data come from the E288, E605, R209, PHENIX, and STAR Collaborations. The high-energy 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 . 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 , there is a knee point for the trend of at –50 GeV. Meanwhile, there is a boundary at –500 GeV above which increases significantly. With the increase of , there is a knee point for the trend of at –15 GeV/. Meanwhile, there is a boundary at –60 GeV/ above which increases significantly. In the two-component Erlang distribution, there is a boundary at –500 GeV for , 500–1100 GeV for , and 200–500 GeV for , above which , , and increase quickly. Meanwhile, there is a boundary at –60 GeV/ above which , , and increase quickly. In the convolution of two Hagedorn functions, increases obviously with increasing and . There is a boundary for the trend of at –200 GeV and –50 GeV/.

With increasing and , the parameters , , , , and 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 increases initially and then decreases, and the parameter decreases, with the increase of and . 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 and increase. The fit results in this paper are comparable to the QCD NNLO and 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.

#### Acknowledgments

X.-H. Zhang and F.-H. Liu’s work was supported by the National Natural Science Foundation of China under Grant Nos. 12047571, 11575103, and 11947418, the Scientific and Technological Innovation Programs of Higher Education Institutions in Shanxi (STIP) under Grant No. 201802017, the Shanxi Provincial Natural Science Foundation under Grant No. 201901D111043, and the Fund for Shanxi “1331 Project” Key Subjects Construction.