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

Authors: Xu-Hong Zhang; Fu-Hu Liu
Journal: Advances in High Energy Physics
 (2021)

---

## Abstract

We analyze the transverse momentum (pT) 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 sNN or s if in a simplified form) range from ~20 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 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 s 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.


---

## Body

## 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 (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.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 pT 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 quarkq 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 (pT) distribution of ℓℓ¯ depend on the momentum distribution of qq¯ 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,pT, 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, pT 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 ℓℓ¯pT 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 pT 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 pT 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ℓℓ¯pT 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, thepT 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, pt1 and pt2 are the transverse momenta of the two quarks, and pT 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 ofpT as a simple Boltzmann distribution [36–38]:
(1)fpTpT=1NdNdpT=CBpTexp−pT2+m02TB,where N is the number of identical particles of mass m0 produced in the collisions, CB is the normalization constant, and TB 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 pT spectrum of hadrons [26–30] is used in this work. We have the following form to describe the transverse momentum (pt) distribution of (anti-)quark:
(2)f1pt=Nqpt1+1nTpt2+mq2−mq−n,where Nq is the normalization constant, T and n are the fitted parameters, and mq is the mass of q or q¯ taking part in the reaction. In general, we use mu=md=0.3GeV/c2 in the Drell-Yan process because q or q¯ is from the participant hadrons. The same mu or md is for sea quarks and those in baryons, where the sea quarks of higher mass are not considered in this work. In QGP, mu=0.003GeV/c2 and md=0.007GeV/c2 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 pT 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 nS≥2 and nH=2, respectively. The contribution (pt) of each parton source to pT of final-state particle is assumed to obey an exponential function:
(3)fipt=1ptexp−ptpt,where pt represents the average pt contributed by the i-th source. Because pt is the same for different sources, the index i in pti is omitted.ThepT distribution contributed by nS (nH) sources is the convolution of nS (nH) 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
(4)fpT=kpTnS−1nS−1!ptSnSexp−pTptS+1−kpTptH2exp−pTptH.Fitting the data with the two-component Erlang distribution, we can get the changes of parametersptS, ptH, and k.We should discuss the values ofnS and nH further. If nS>2, the participant partons are expected to be qq¯ and nS−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 nS=2 or 3 in this work. Generally, the larger the nS, the sharper the distribution peak. In many cases, nS=3 means that qq¯ and a gluon participate in the soft process. For all cases, nH=2 (always true in this work) means that only qq¯ 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 pT, but not the entire pT interval. In the case of using the Hagedorn function in a wide range of pT, the probability density function of pT can be expressed as
(5)f1pT=ApT1+pTp1−n1,where A is the normalization constant, and p1 and n1 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 pT spectra. In some cases, the soft excitation process in the low pT 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 pt instead of pT 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 whichpT/p1 in Eq. (5) is replaced by pT2/p12, ApT is replaced by ApT2/mT, or pT/p1 is replaced by pT2/p12, and ApT is replaced by A, where p1 and A vary in different revisions. The uses of the revised Hagedorn functions result in some overestimations in low (or high) pT region comparing to the Hagedorn function. Contrarily, these revisions result in some underestimations in high (or low) pT 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 ofq 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, f1pt1, of a single (anti-)quark’s contribution. The second (anti-)quark’s contribution is f2pt2=f1pT−pt1, where pT is still the transverse momentum of the ℓℓ¯ system. So, the convolution of two probability density functions should be used to describe the pT spectrum of ℓℓ¯ in the Drell-Yan process. We have the convolution of two Eq. (2) or (5) to be expressed as
(6)fpT=∫0pTf1pt1f1pT−pt1dpt1,where f1pt1 [f1pT−pt1] 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 vectorsp→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 pT spectrum of ℓℓ¯ in the Drell-Yan process, where f1pt1 and f1pT−pt1 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 pT 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 pT spectrum.


## 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 ofpT as a simple Boltzmann distribution [36–38]:
(1)fpTpT=1NdNdpT=CBpTexp−pT2+m02TB,where N is the number of identical particles of mass m0 produced in the collisions, CB is the normalization constant, and TB 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 pT spectrum of hadrons [26–30] is used in this work. We have the following form to describe the transverse momentum (pt) distribution of (anti-)quark:
(2)f1pt=Nqpt1+1nTpt2+mq2−mq−n,where Nq is the normalization constant, T and n are the fitted parameters, and mq is the mass of q or q¯ taking part in the reaction. In general, we use mu=md=0.3GeV/c2 in the Drell-Yan process because q or q¯ is from the participant hadrons. The same mu or md is for sea quarks and those in baryons, where the sea quarks of higher mass are not considered in this work. In QGP, mu=0.003GeV/c2 and md=0.007GeV/c2 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 pT 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 nS≥2 and nH=2, respectively. The contribution (pt) of each parton source to pT of final-state particle is assumed to obey an exponential function:
(3)fipt=1ptexp−ptpt,where pt represents the average pt contributed by the i-th source. Because pt is the same for different sources, the index i in pti is omitted.ThepT distribution contributed by nS (nH) sources is the convolution of nS (nH) 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
(4)fpT=kpTnS−1nS−1!ptSnSexp−pTptS+1−kpTptH2exp−pTptH.Fitting the data with the two-component Erlang distribution, we can get the changes of parametersptS, ptH, and k.We should discuss the values ofnS and nH further. If nS>2, the participant partons are expected to be qq¯ and nS−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 nS=2 or 3 in this work. Generally, the larger the nS, the sharper the distribution peak. In many cases, nS=3 means that qq¯ and a gluon participate in the soft process. For all cases, nH=2 (always true in this work) means that only qq¯ 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 pT, but not the entire pT interval. In the case of using the Hagedorn function in a wide range of pT, the probability density function of pT can be expressed as
(5)f1pT=ApT1+pTp1−n1,where A is the normalization constant, and p1 and n1 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 pT spectra. In some cases, the soft excitation process in the low pT 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 pt instead of pT 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 whichpT/p1 in Eq. (5) is replaced by pT2/p12, ApT is replaced by ApT2/mT, or pT/p1 is replaced by pT2/p12, and ApT is replaced by A, where p1 and A vary in different revisions. The uses of the revised Hagedorn functions result in some overestimations in low (or high) pT region comparing to the Hagedorn function. Contrarily, these revisions result in some underestimations in high (or low) pT 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 ofq 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, f1pt1, of a single (anti-)quark’s contribution. The second (anti-)quark’s contribution is f2pt2=f1pT−pt1, where pT is still the transverse momentum of the ℓℓ¯ system. So, the convolution of two probability density functions should be used to describe the pT spectrum of ℓℓ¯ in the Drell-Yan process. We have the convolution of two Eq. (2) or (5) to be expressed as
(6)fpT=∫0pTf1pt1f1pT−pt1dpt1,where f1pt1 [f1pT−pt1] 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 vectorsp→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 pT spectrum of ℓℓ¯ in the Drell-Yan process, where f1pt1 and f1pT−pt1 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 pT 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 pT spectrum.


## 3. Results and Discussion

### 3.1. Comparison with Data

Figure1 shows the pT spectra of ℓℓ¯ [with different invariant masses (Q) or Feynman variables (xF)] 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 Figures 1(a)–1(c) are quoted from the proton-copper (p-Cu) collision experiments performed by the E288 Collaboration [43], and the collision energy per nucleon pair sNN or s 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 p-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 14GeV/c2, while the corresponding invariant mass ranges of the E605 Collaboration are from 7 to 18GeV/c2. 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 Q in the range of 4.05–13.05GeV/c2 with different scalings, where the units GeV/c2 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 Q and xF, 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–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=∑jDataj−Fitj2/Uncertaintyj2 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 nS=3 and nH=2 which are not listed in the table to avoid trivialness. One can see that the three functions can approximately describe the pT 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 1
Transverse momentum spectra ofℓℓ¯ (with different invariant masses Q or Feynman variables xF) produced by the Drell-Yan process in different collisions at different energies. The data points in (a)–(c) are quoted from the E288 Collaboration [43] and obtained in p-Cu collisions at s=19.4 GeV (4≤Q<9GeV/c2), 23.8 GeV (4≤Q<13GeV/c2), and 27.4 GeV (5≤Q<13GeV/c2), respectively. The data points in (d) come from p-Cu collisions at s=38.8 GeV performed by the E605 Collaboration [44]. The data points in (e) (4.05≤Q<13.05GeV/c2) and (f) (0.0≤xF<1.0) come from π−-W collisions at s=21.7 GeV measured by the FNAL-615 Collaboration [45], where the units GeV/c2 are not shown in (e) due to crowding space. In panel (d), ℓℓ¯ is e+e−, and in other panels, ℓℓ¯ 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)], the two-component Erlang distribution [Eq. (4)], and the convolution of two Hagedorn functions [Eqs. (5) and (6)], respectively, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations.
(a)(b)(c)(d)(e)(f)Table 1
Values of parametersn and T in the Lévy-Tsallis function (solid curves), ptS, ptH, and k in the two-component Erlang distribution (dashed curves), as well as p1 and n1 in the Hagedorn function (dotted curves) in Figure 1. The χ2 is for each fit. The ndof is for the fits of the two-component Erlang distribution. For the fits of Lévy-Tsallis function and Hagedorn function, the ndof needs one more. For Figure 1(f), the selection is not the invariant mass Q, but the Feynman variable xF. nS=3 and nH=2 for Figure 1, which are not listed in the table.
Figures (GeV)Q (GeV/c2)nT (GeV)χ2ptS (GeV/c)ptH (GeV/c)kχ2p1 (GeV/c)n1χ2ndofFigure1(a)19.44-53.910±0.0890.161±0.0051220.294±0.0030.611±0.0060.961±0.003521.028±0.0066.994±0.022164105-63.969±0.0890.166±0.005820.303±0.0040.608±0.0060.974±0.003241.031±0.0066.992±0.02299116-74.008±0.0900.169±0.005340.327±0.0060.588±0.0050.978±0.003101.035±0.0066.987±0.02243107-84.139±0.0930.171±0.00590.331±0.0080.511±0.0040.973±0.00311.046±0.0066.944±0.0211288-94.155±0.0980.176±0.00560.348±0.0080.494±0.0040.961±0.00331.049±0.0066.921±0.02087Figure1(b)23.84-53.571±0.0870.181±0.0061550.318±0.0040.642±0.0060.964±0.003131.002±0.0056.264±0.016157115-63.892±0.0890.182±0.006960.339±0.0050.598±0.0050.967±0.003221.003±0.0046.264±0.01599116-73.909±0.0920.195±0.007670.364±0.0070.579±0.0050.961±0.003151.005±0.0046.230±0.01574117-83.988±0.0890.196±0.007380.379±0.0080.543±0.0050.970±0.003141.010±0.0046.194±0.01340108-94.135±0.1020.197±0.007420.380±0.0080.494±0.0040.974±0.003311.011±0.0046.193±0.01442911-124.410±0.1140.241±0.01030.427±0.0110.433±0.0040.982±0.00331.025±0.0056.021±0.01246Figure1(c)27.45-63.264±0.0600.187±0.0061030.367±0.0060.773±0.0070.947±0.003351.030±0.0056.164±0.010159116-73.760±0.0840.202±0.0081840.386±0.0080.711±0.0060.953±0.003881.041±0.0056.144±0.009180117-83.653±0.0880.213±0.0091440.403±0.0080.688±0.0060.953±0.003381.067±0.0056.136±0.010165118-93.905±0.0920.217±0.0091090.409±0.0080.638±0.0060.953±0.003451.082±0.0066.125±0.0101221111-124.247±0.1060.230±0.012210.442±0.0100.602±0.0050.931±0.003101.081±0.0056.093±0.009271112-134.466±0.1130.253±0.014150.464±0.0110.564±0.0050.944±0.003111.086±0.0066.089±0.00916103-144.986±0.1330.256±0.01440.468±0.0110.552±0.0050.944±0.00331.084±0.0066.092±0.00948Figure1(d)38.87-83.211±0.0640.254±0.013190.446±0.0100.654±0.0060.983±0.00361.257±0.0096.191±0.0104278-94.371±0.0800.262±0.013470.469±0.0120.602±0.0060.992±0.003111.180±0.0086.201±0.011521011.5-13.54.694±0.1180.267±0.013320.473±0.0130.529±0.0050.995±0.00371.194±0.0076.190±0.010441113.5-18.04.949±0.1310.278±0.015940.480±0.0130.434±0.0070.986±0.003341.197±0.0076.171±0.0129810Figure1(e)21.74.05-4.503.221±0.0850.155±0.0054120.397±0.0080.473±0.0130.958±0.0031070.572±0.0062.175±0.012563134.50-4.953.073±0.0730.163±0.0052270.429±0.0090.314±0.0050.968±0.003940.505±0.0052.696±0.010311144.95-5.403.663±0.0890.176±0.0061220.411±0.0080.481±0.0140.966±0.003880.512±0.0052.357±0.014180135.40-5.853.348±0.0880.165±0.005470.425±0.0090.632±0.0160.955±0.003570.595±0.0052.803±0.01279115.85-6.753.841±0.1120.186±0.008440.428±0.0090.358±0.0060.977±0.003550.551±0.0052.540±0.01468146.75-7.653.453±0.0880.188±0.008320.442±0.0100.499±0.0150.971±0.003480.563±0.0052.746±0.01570117.65-9.003.965±0.0900.202±0.008200.433±0.0100.423±0.0090.997±0.003280.552±0.0052.323±0.01433129.00-10.353.321±0.0880.167±0.005150.379±0.0070.482±0.0140.983±0.003100.506±0.0052.904±0.018201110.35-11.703.013±0.0730.172±0.005170.372±0.0070.219±0.0030.984±0.003130.561±0.0052.767±0.01620911.70-13.052.707±0.0450.153±0.00530.354±0.0061.491±0.0570.979±0.00330.521±0.0053.312±0.02537Figure1(f)21.7xF (1/c)0.00-0.102.702±0.0450.148±0.005320.408±0.0080.645±0.0160.974±0.003150.189±0.0052.035±0.00847110.10-0.202.873±0.0480.113±0.004690.374±0.0070.503±0.0150.976±0.003250.188±0.0052.806±0.01292150.20-0.303.978±0.0920.187±0.006900.386±0.0080.487±0.0140.979±0.003460.190±0.0052.865±0.012134130.30-0.402.938±0.0580.136±0.0051200.393±0.0080.580±0.0150.963±0.003620.184±0.0052.468±0.011159130.40-0.503.046±0.0610.152±0.0051050.393±0.0080.512±0.0140.975±0.003570.185±0.0053.581±0.018190130.50-0.602.820±0.0520.158±0.0051060.399±0.0080.491±0.0140.959±0.003540.293±0.0103.302±0.016162130.60-0.702.721±0.0450.180±0.0061410.383±0.0080.403±0.0120.969±0.003520.298±0.0102.841±0.012198130.70-0.802.938±0.0480.158±0.0051190.376±0.0070.427±0.0130.956±0.003550.202±0.0083.102±0.016166130.80-0.903.022±0.0630.126±0.004740.284±0.0040.394±0.0110.989±0.003460.367±0.0123.673±0.018101110.90-1.003.143±0.0680.095±0.012250.195±0.0020.281±0.0080.985±0.00360.386±0.0153.868±0.022177Figure2 shows the pT spectra of ℓℓ¯ (with different invariant mass Q) generated by the Drell-Yan process in proton-proton (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 s=62GeV, and the Q range is 5−8GeV/c2. The data points in Figure 2(b) show the experimental results from the PHENIX Collaboration [47]. The s is 200 GeV, the Q range is 4.8−8.2GeV/c2, 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 s is 510 GeV, the Q range is 73−114GeV/c2, 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 pT spectra of ℓℓ¯ produced by the Drell-Yan process in high-energy pp collisions.Figure 2
Transverse momentum spectra ofℓℓ¯ (with different invariant masses Q) produced by the Drell-Yan process in pp collisions at s=62 (a), 200 (b), and 510 GeV (c). 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)], the two-component Erlang distribution [Eq. (4)], and the convolution of two Hagedorn functions [Eqs. (5) and (6)], respectively, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations.
(a)(b)(c)Table 2
Values of parametersn and T in the Lévy-Tsallis function (solid curves), ptS, ptH, and k in the two-component Erlang distribution (dashed curves), as well as p1 and n1 in the Hagedorn function (dotted curves) in Figures 2–4. The χ2 is for each fit. The ndof is for the fits of the two-component Erlang distribution. For the fits of the Lévy-Tsallis function and Hagedorn function, the ndof needs one more. For Figures 4(a) and 4(b), different y ranges are included. nS=3 and nH=2 for Figures 2 and 4, as well as nS=2 and nH=2 for Figure 3, which are not listed in the table.
Figures (GeV)Q (GeV/c2)nT (GeV)χ2ptS (GeV/c)ptH (GeV/c)kχ2p1 (GeV/c)n1χ2ndofFigure2(a)62.05-83.887±0.0800.136±0.006210.287±0.0060.913±0.0750.989±0.003261.001±0.0056.829±0.018189Figure2(b)2004.8-8.23.187±0.0680.258±0.015180.617±0.0170.277±0.0100.977±0.003181.111±0.0055.089±0.016158Figure2(c)51073-1143.162±0.0650.792±0.01841.924±0.0621.103±0.1100.801±0.003111.077±0.0123.137±0.01943Figure3(a)180066-1162.618±0.0280.820±0.016252.631±0.1107.092±0.3200.511±0.002241.088±0.0122.916±0.0153733Figure3(b)196066-1162.661±0.0180.840±0.0182702.690±0.1288.120±0.4030.524±0.0021341.222±0.0133.030±0.01647653Figure3(c)180075-1052.588±0.0210.779±0.016162.721±0.1018.151±0.4000.539 ± 0.002161.040±0.0122.934±0.0153014Figure3(d)196070-1102.745±0.0170.904±0.017142.830±0.1108.580±0.4500.539±0.002131.359±0.0153.106±0.017269Figure4(a)700066-1160≤∣y∣<12.415±0.0160.996±0.01633.521±0.20110.710±0.5010.450±0.00231.241±0.0132.780±0.01511101≤∣y∣<22.396±0.0161.003±0.01623.520±0.20110.793±0.5010.439±0.00231.231±0.0122.760±0.0157102≤∣y∣<2.42.336±0.0120.980±0.01613.482±0.20011.742±0.5070.434±0.00211.121±0.0122.679±0.014110Figure4(b)800066-1160.0≤∣y∣<0.42.385±0.0231.023±0.01823.520±0.20110.571±0.5050.425±0.00221.232±0.0122.738±0.015670.4≤∣y∣<0.82.405±0.0231.031±0.02023.549±0.20310.713±0.5090.432±0.00221.266±0.0132.761±0.016670.8≤∣y∣<1.22.411±0.0311.031±0.01923.561±0.21011.760±0.5150.438±0.00221.267±0.0132.767±0.016671.2≤∣y∣<1.62.384±0.0231.031±0.02023.663±0.21311.106±0.5240.438±0.00211.267±0.0132.752±0.016671.6≤∣y∣<2.02.298±0.0270.990±0.02013.700±0.22011.817±0.5280.432±0.00211.178±0.0122.686±0.015372.0≤∣y∣<2.42.259±0.0220.953±0.02013.525±0.20011.482±0.5200.419±0.00211.122±0.0102.658±0.01517Figure4(c)46-662.918±0.0110.987±0.021242.890±0.1758.152±0.4120.555±0.002181.164±0.0102.787±0.01650766-1162.402±0.0221.038±0.02043.493±0.20910.540±0.5050.415±0.00241.290±0.0122.769±0.0151212116-1502.183±0.0111.008±0.011163.960±0.22112.841±0.6010.410±0.002101.371±0.0122.684±0.014407Figure4(d)1300066-1162.177±0.0120.962±0.018113.862±0.22312.570±0.6030.428±0.00261.035±0.0102.550±0.0142210Figure4(e)700060-1202.487±0.0201.071±0.020143.563±0.20510.390±0.5010.433±0.00251.373±0.0172.831±0.01124680002.341±0.0241.001±0.01973.441±0.20310.851±0.5050.409±0.002191.210±0.0152.719±0.014116Figure4(f)700060-1202.483±0.0110.936±0.017193.114±0.1839.211±0.4320.436±0.002321.214±0.0112.846±0.01536880003.438±0.0100.941±0.017153.121±0.1839.263±0.4310.433±0.002301.211±0.0112.810±0.015478130003.394±0.0450.996±0.020172.590±0.1294.653±0.2090.444±0.002361.170±0.0102.859±0.015168Similar to Figure2, Figure 3 shows the pT spectra of ℓℓ¯ produced by the Drell-Yan process in proton-antiproton (p-p¯ or pp¯) collisions with 66≤Q<116 GeV/c2(Figures 3(a) and 3(b)),75≤Q<105GeV/c2 (Figure 3(c)), and 70≤Q<110GeV/c2 (Figure 3(d)) ats=1.8TeV (Figures 3(a) and 3(c)) and s=1.96TeV (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 y 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 pT spectra of ℓℓ¯ produced by the Drell-Yan process in high-energy pp¯ collisions.Figure 3
Same as Figure2, but showing the pT spectra of ℓℓ¯ (with different invariant masses Q) produced by the Drell-Yan process in pp¯ collisions at s=1.8 (a, c) and 1.96 TeV (b, d). The data points are quoted from the CDF (a, b) [48, 49] and D0 Collaborations (c, d) [50–52]. In all panels, ℓℓ¯ is e+e−, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations.
(a)(b)(c)(d)Figure4 shows the pT spectra of ℓℓ¯ produced by the Drell-Yan process in high-energy pp 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 (s, Q, and y)] 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 pT spectra of ℓℓ¯ produced by the Drell-Yan process in pp collisions at ultrahigh energies.Figure 4
Same as Figure2, but showing the pT spectra of ℓℓ¯ [with different conditions (s, Q, and y)] produced by the Drell-Yan process in pp collisions at the LHC energies. The data points are quoted from the experiments performed by the ATLAS (a)–(d) [53–55], CMS (e) [56, 57], and LHCb Collaborations (f) [58–60]. The flavor of ℓℓ¯ is shown in the panels, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations.
(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χ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 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.


### 3.2. Tendency of Parameters

In order to study the underlying law that governs the variation of the extracted parameters withs and Q, 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 n, and Figures 5(c) and 5(d) show the trend of parameter T, 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 s. On average, it can be seen that as s increases, the parameter n increases quickly and then decreases slowly, and the parameter T increases slowly and then significantly. There is a knee point for the trend of n at s≈40–50 GeV. Meanwhile, there is a boundary at s≈200–500 GeV above which T increases significantly. Similarly, we compare Figures 5(b) and 5(d) and analyze the variation of parameters with Q. It can be clearly seen that the parameter n increases firstly and then decreases, and the parameter T increases slowly and then significantly, with the increase of Q. There is a knee point for the trend of n at Q≈14–15 GeV/c2. Meanwhile, there is a boundary at Q≈20–60 GeV/c2 above which T increases significantly.Figure 5
The variation of parametersn (a, b) and T (c, d) in the Lévy-Tsallis function with energy s (a, c) and invariant mass Q (b, d). The parameter values are taken from Figures 1–4 and recorded in Tables 1 and 2. Since the invariant mass grouping in Figure 1(e) is different from others, and the varying quantity in Figure 1(f) is the Feynman variable, Figure 5 does not include the parameters from Figures 1(e) and 1(f) to avoid trivialness.
(a)(b)(c)(d)Figure6 is similar to Figure 5, but shows the dependence of parameters ptS (Figures 6(a) and 6(b)), ptH (Figures 6(c) and 6(d)), and k (Figures 6(e) and 6(f)) on s(Figures 6(a), 6(c), and 6(e)) andQ (Figures 6(b), 6(d), and 6(f)) obtained from the two-component Erlang distribution. One can see that with increasing s, ptS and ptH increase slowly and then quickly, and k decreases slowly and then quickly. For s, there is a boundary between 60 and 500 GeV for ptS, between 500 and 1100 GeV for ptH, and between 200 and 500 GeV for k. Meanwhile, with increasing Q, ptS and ptH increase slowly and then quickly, and k decreases slowly and then quickly. There is a boundary at Q≈20–60 GeV/c2.Figure 6
Similar to Figure5, but showing the variation of parameters ptS (a, b), ptH (c, d), and k (e, f) in the two-component Erlang distribution with s (a, c, e) and Q (b, d, f).
(a)(b)(c)(d)(e)(f)In Figure7, 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 p1 has a slight tendency to increase with the increase of s and Q. In Figures 7(c) and 7(d), one can see that the parameter n1 decreases quickly and then slowly with the increase of s; n1 also decreases slowly and then quickly with increase of Q. There is a boundary for the trend of n1 at s≈40–200 GeV and Q≈20–50 GeV/c2. We note that the variation of parameter is obtained from the the average values of the extracted parameters, for each s or Q.Figure 7
Similar to Figure5, but showing the parameters p1 (a, b) and n1 (c, d) in the Hagedorn function with (a) (c) s and (b) (d) Q.
(a)(b)(c)(d)It should be pointed out that the values of the parameters in Figures5–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 pT 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 parametersT, ptS, ptH, and p1 show monotonous increasing trend when s 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 s 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 n1 decreases, with the increase of s 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 s 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 thepT 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 softest 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 ~40GeV, 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/c2. At below (above) the softest point, we see naturally smaller (larger) Q, except for Q=5–6 and 5–8 GeV/c2 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/c2, the system is close to the onset energy of phase transition. The parameters T, ptS, ptH, 1−k, and p1 show similar small and almost flat values, respectively, and n and n1 show different trends. At the energy of TeV and the invariant mass around 100 GeV/c2, the system is far away from the onset energy of phase transition. The parameters T, ptS, ptH, 1−k, and p1 show large values, and n and n1 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 massQ of ℓℓ¯) and obey renormalization group equations. Although this scaling is not included in the function f1pt1 or f1pT−pt1 of Eq. (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 f1pt1 or f1pT−pt1 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 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 Figures1(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 Figure2(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 (N3LO or N3LL) from the Drell-Yan data.The meanings of histograms in Figure3 and Figure 4(a), Figure 4(b), Figure 4(c) with 46≤Q<66 and 116≤Q<150 GeV/c2, 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/c2 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 N3LL. 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+N3LL.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.


## 3.1. Comparison with Data

Figure1 shows the pT spectra of ℓℓ¯ [with different invariant masses (Q) or Feynman variables (xF)] 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 Figures 1(a)–1(c) are quoted from the proton-copper (p-Cu) collision experiments performed by the E288 Collaboration [43], and the collision energy per nucleon pair sNN or s 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 p-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 14GeV/c2, while the corresponding invariant mass ranges of the E605 Collaboration are from 7 to 18GeV/c2. 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 Q in the range of 4.05–13.05GeV/c2 with different scalings, where the units GeV/c2 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 Q and xF, 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–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=∑jDataj−Fitj2/Uncertaintyj2 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 nS=3 and nH=2 which are not listed in the table to avoid trivialness. One can see that the three functions can approximately describe the pT 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 1
Transverse momentum spectra ofℓℓ¯ (with different invariant masses Q or Feynman variables xF) produced by the Drell-Yan process in different collisions at different energies. The data points in (a)–(c) are quoted from the E288 Collaboration [43] and obtained in p-Cu collisions at s=19.4 GeV (4≤Q<9GeV/c2), 23.8 GeV (4≤Q<13GeV/c2), and 27.4 GeV (5≤Q<13GeV/c2), respectively. The data points in (d) come from p-Cu collisions at s=38.8 GeV performed by the E605 Collaboration [44]. The data points in (e) (4.05≤Q<13.05GeV/c2) and (f) (0.0≤xF<1.0) come from π−-W collisions at s=21.7 GeV measured by the FNAL-615 Collaboration [45], where the units GeV/c2 are not shown in (e) due to crowding space. In panel (d), ℓℓ¯ is e+e−, and in other panels, ℓℓ¯ 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)], the two-component Erlang distribution [Eq. (4)], and the convolution of two Hagedorn functions [Eqs. (5) and (6)], respectively, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations.
(a)(b)(c)(d)(e)(f)Table 1
Values of parametersn and T in the Lévy-Tsallis function (solid curves), ptS, ptH, and k in the two-component Erlang distribution (dashed curves), as well as p1 and n1 in the Hagedorn function (dotted curves) in Figure 1. The χ2 is for each fit. The ndof is for the fits of the two-component Erlang distribution. For the fits of Lévy-Tsallis function and Hagedorn function, the ndof needs one more. For Figure 1(f), the selection is not the invariant mass Q, but the Feynman variable xF. nS=3 and nH=2 for Figure 1, which are not listed in the table.
Figures (GeV)Q (GeV/c2)nT (GeV)χ2ptS (GeV/c)ptH (GeV/c)kχ2p1 (GeV/c)n1χ2ndofFigure1(a)19.44-53.910±0.0890.161±0.0051220.294±0.0030.611±0.0060.961±0.003521.028±0.0066.994±0.022164105-63.969±0.0890.166±0.005820.303±0.0040.608±0.0060.974±0.003241.031±0.0066.992±0.02299116-74.008±0.0900.169±0.005340.327±0.0060.588±0.0050.978±0.003101.035±0.0066.987±0.02243107-84.139±0.0930.171±0.00590.331±0.0080.511±0.0040.973±0.00311.046±0.0066.944±0.0211288-94.155±0.0980.176±0.00560.348±0.0080.494±0.0040.961±0.00331.049±0.0066.921±0.02087Figure1(b)23.84-53.571±0.0870.181±0.0061550.318±0.0040.642±0.0060.964±0.003131.002±0.0056.264±0.016157115-63.892±0.0890.182±0.006960.339±0.0050.598±0.0050.967±0.003221.003±0.0046.264±0.01599116-73.909±0.0920.195±0.007670.364±0.0070.579±0.0050.961±0.003151.005±0.0046.230±0.01574117-83.988±0.0890.196±0.007380.379±0.0080.543±0.0050.970±0.003141.010±0.0046.194±0.01340108-94.135±0.1020.197±0.007420.380±0.0080.494±0.0040.974±0.003311.011±0.0046.193±0.01442911-124.410±0.1140.241±0.01030.427±0.0110.433±0.0040.982±0.00331.025±0.0056.021±0.01246Figure1(c)27.45-63.264±0.0600.187±0.0061030.367±0.0060.773±0.0070.947±0.003351.030±0.0056.164±0.010159116-73.760±0.0840.202±0.0081840.386±0.0080.711±0.0060.953±0.003881.041±0.0056.144±0.009180117-83.653±0.0880.213±0.0091440.403±0.0080.688±0.0060.953±0.003381.067±0.0056.136±0.010165118-93.905±0.0920.217±0.0091090.409±0.0080.638±0.0060.953±0.003451.082±0.0066.125±0.0101221111-124.247±0.1060.230±0.012210.442±0.0100.602±0.0050.931±0.003101.081±0.0056.093±0.009271112-134.466±0.1130.253±0.014150.464±0.0110.564±0.0050.944±0.003111.086±0.0066.089±0.00916103-144.986±0.1330.256±0.01440.468±0.0110.552±0.0050.944±0.00331.084±0.0066.092±0.00948Figure1(d)38.87-83.211±0.0640.254±0.013190.446±0.0100.654±0.0060.983±0.00361.257±0.0096.191±0.0104278-94.371±0.0800.262±0.013470.469±0.0120.602±0.0060.992±0.003111.180±0.0086.201±0.011521011.5-13.54.694±0.1180.267±0.013320.473±0.0130.529±0.0050.995±0.00371.194±0.0076.190±0.010441113.5-18.04.949±0.1310.278±0.015940.480±0.0130.434±0.0070.986±0.003341.197±0.0076.171±0.0129810Figure1(e)21.74.05-4.503.221±0.0850.155±0.0054120.397±0.0080.473±0.0130.958±0.0031070.572±0.0062.175±0.012563134.50-4.953.073±0.0730.163±0.0052270.429±0.0090.314±0.0050.968±0.003940.505±0.0052.696±0.010311144.95-5.403.663±0.0890.176±0.0061220.411±0.0080.481±0.0140.966±0.003880.512±0.0052.357±0.014180135.40-5.853.348±0.0880.165±0.005470.425±0.0090.632±0.0160.955±0.003570.595±0.0052.803±0.01279115.85-6.753.841±0.1120.186±0.008440.428±0.0090.358±0.0060.977±0.003550.551±0.0052.540±0.01468146.75-7.653.453±0.0880.188±0.008320.442±0.0100.499±0.0150.971±0.003480.563±0.0052.746±0.01570117.65-9.003.965±0.0900.202±0.008200.433±0.0100.423±0.0090.997±0.003280.552±0.0052.323±0.01433129.00-10.353.321±0.0880.167±0.005150.379±0.0070.482±0.0140.983±0.003100.506±0.0052.904±0.018201110.35-11.703.013±0.0730.172±0.005170.372±0.0070.219±0.0030.984±0.003130.561±0.0052.767±0.01620911.70-13.052.707±0.0450.153±0.00530.354±0.0061.491±0.0570.979±0.00330.521±0.0053.312±0.02537Figure1(f)21.7xF (1/c)0.00-0.102.702±0.0450.148±0.005320.408±0.0080.645±0.0160.974±0.003150.189±0.0052.035±0.00847110.10-0.202.873±0.0480.113±0.004690.374±0.0070.503±0.0150.976±0.003250.188±0.0052.806±0.01292150.20-0.303.978±0.0920.187±0.006900.386±0.0080.487±0.0140.979±0.003460.190±0.0052.865±0.012134130.30-0.402.938±0.0580.136±0.0051200.393±0.0080.580±0.0150.963±0.003620.184±0.0052.468±0.011159130.40-0.503.046±0.0610.152±0.0051050.393±0.0080.512±0.0140.975±0.003570.185±0.0053.581±0.018190130.50-0.602.820±0.0520.158±0.0051060.399±0.0080.491±0.0140.959±0.003540.293±0.0103.302±0.016162130.60-0.702.721±0.0450.180±0.0061410.383±0.0080.403±0.0120.969±0.003520.298±0.0102.841±0.012198130.70-0.802.938±0.0480.158±0.0051190.376±0.0070.427±0.0130.956±0.003550.202±0.0083.102±0.016166130.80-0.903.022±0.0630.126±0.004740.284±0.0040.394±0.0110.989±0.003460.367±0.0123.673±0.018101110.90-1.003.143±0.0680.095±0.012250.195±0.0020.281±0.0080.985±0.00360.386±0.0153.868±0.022177Figure2 shows the pT spectra of ℓℓ¯ (with different invariant mass Q) generated by the Drell-Yan process in proton-proton (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 s=62GeV, and the Q range is 5−8GeV/c2. The data points in Figure 2(b) show the experimental results from the PHENIX Collaboration [47]. The s is 200 GeV, the Q range is 4.8−8.2GeV/c2, 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 s is 510 GeV, the Q range is 73−114GeV/c2, 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 pT spectra of ℓℓ¯ produced by the Drell-Yan process in high-energy pp collisions.Figure 2
Transverse momentum spectra ofℓℓ¯ (with different invariant masses Q) produced by the Drell-Yan process in pp collisions at s=62 (a), 200 (b), and 510 GeV (c). 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)], the two-component Erlang distribution [Eq. (4)], and the convolution of two Hagedorn functions [Eqs. (5) and (6)], respectively, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations.
(a)(b)(c)Table 2
Values of parametersn and T in the Lévy-Tsallis function (solid curves), ptS, ptH, and k in the two-component Erlang distribution (dashed curves), as well as p1 and n1 in the Hagedorn function (dotted curves) in Figures 2–4. The χ2 is for each fit. The ndof is for the fits of the two-component Erlang distribution. For the fits of the Lévy-Tsallis function and Hagedorn function, the ndof needs one more. For Figures 4(a) and 4(b), different y ranges are included. nS=3 and nH=2 for Figures 2 and 4, as well as nS=2 and nH=2 for Figure 3, which are not listed in the table.
Figures (GeV)Q (GeV/c2)nT (GeV)χ2ptS (GeV/c)ptH (GeV/c)kχ2p1 (GeV/c)n1χ2ndofFigure2(a)62.05-83.887±0.0800.136±0.006210.287±0.0060.913±0.0750.989±0.003261.001±0.0056.829±0.018189Figure2(b)2004.8-8.23.187±0.0680.258±0.015180.617±0.0170.277±0.0100.977±0.003181.111±0.0055.089±0.016158Figure2(c)51073-1143.162±0.0650.792±0.01841.924±0.0621.103±0.1100.801±0.003111.077±0.0123.137±0.01943Figure3(a)180066-1162.618±0.0280.820±0.016252.631±0.1107.092±0.3200.511±0.002241.088±0.0122.916±0.0153733Figure3(b)196066-1162.661±0.0180.840±0.0182702.690±0.1288.120±0.4030.524±0.0021341.222±0.0133.030±0.01647653Figure3(c)180075-1052.588±0.0210.779±0.016162.721±0.1018.151±0.4000.539 ± 0.002161.040±0.0122.934±0.0153014Figure3(d)196070-1102.745±0.0170.904±0.017142.830±0.1108.580±0.4500.539±0.002131.359±0.0153.106±0.017269Figure4(a)700066-1160≤∣y∣<12.415±0.0160.996±0.01633.521±0.20110.710±0.5010.450±0.00231.241±0.0132.780±0.01511101≤∣y∣<22.396±0.0161.003±0.01623.520±0.20110.793±0.5010.439±0.00231.231±0.0122.760±0.0157102≤∣y∣<2.42.336±0.0120.980±0.01613.482±0.20011.742±0.5070.434±0.00211.121±0.0122.679±0.014110Figure4(b)800066-1160.0≤∣y∣<0.42.385±0.0231.023±0.01823.520±0.20110.571±0.5050.425±0.00221.232±0.0122.738±0.015670.4≤∣y∣<0.82.405±0.0231.031±0.02023.549±0.20310.713±0.5090.432±0.00221.266±0.0132.761±0.016670.8≤∣y∣<1.22.411±0.0311.031±0.01923.561±0.21011.760±0.5150.438±0.00221.267±0.0132.767±0.016671.2≤∣y∣<1.62.384±0.0231.031±0.02023.663±0.21311.106±0.5240.438±0.00211.267±0.0132.752±0.016671.6≤∣y∣<2.02.298±0.0270.990±0.02013.700±0.22011.817±0.5280.432±0.00211.178±0.0122.686±0.015372.0≤∣y∣<2.42.259±0.0220.953±0.02013.525±0.20011.482±0.5200.419±0.00211.122±0.0102.658±0.01517Figure4(c)46-662.918±0.0110.987±0.021242.890±0.1758.152±0.4120.555±0.002181.164±0.0102.787±0.01650766-1162.402±0.0221.038±0.02043.493±0.20910.540±0.5050.415±0.00241.290±0.0122.769±0.0151212116-1502.183±0.0111.008±0.011163.960±0.22112.841±0.6010.410±0.002101.371±0.0122.684±0.014407Figure4(d)1300066-1162.177±0.0120.962±0.018113.862±0.22312.570±0.6030.428±0.00261.035±0.0102.550±0.0142210Figure4(e)700060-1202.487±0.0201.071±0.020143.563±0.20510.390±0.5010.433±0.00251.373±0.0172.831±0.01124680002.341±0.0241.001±0.01973.441±0.20310.851±0.5050.409±0.002191.210±0.0152.719±0.014116Figure4(f)700060-1202.483±0.0110.936±0.017193.114±0.1839.211±0.4320.436±0.002321.214±0.0112.846±0.01536880003.438±0.0100.941±0.017153.121±0.1839.263±0.4310.433±0.002301.211±0.0112.810±0.015478130003.394±0.0450.996±0.020172.590±0.1294.653±0.2090.444±0.002361.170±0.0102.859±0.015168Similar to Figure2, Figure 3 shows the pT spectra of ℓℓ¯ produced by the Drell-Yan process in proton-antiproton (p-p¯ or pp¯) collisions with 66≤Q<116 GeV/c2(Figures 3(a) and 3(b)),75≤Q<105GeV/c2 (Figure 3(c)), and 70≤Q<110GeV/c2 (Figure 3(d)) ats=1.8TeV (Figures 3(a) and 3(c)) and s=1.96TeV (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 y 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 pT spectra of ℓℓ¯ produced by the Drell-Yan process in high-energy pp¯ collisions.Figure 3
Same as Figure2, but showing the pT spectra of ℓℓ¯ (with different invariant masses Q) produced by the Drell-Yan process in pp¯ collisions at s=1.8 (a, c) and 1.96 TeV (b, d). The data points are quoted from the CDF (a, b) [48, 49] and D0 Collaborations (c, d) [50–52]. In all panels, ℓℓ¯ is e+e−, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations.
(a)(b)(c)(d)Figure4 shows the pT spectra of ℓℓ¯ produced by the Drell-Yan process in high-energy pp 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 (s, Q, and y)] 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 pT spectra of ℓℓ¯ produced by the Drell-Yan process in pp collisions at ultrahigh energies.Figure 4
Same as Figure2, but showing the pT spectra of ℓℓ¯ [with different conditions (s, Q, and y)] produced by the Drell-Yan process in pp collisions at the LHC energies. The data points are quoted from the experiments performed by the ATLAS (a)–(d) [53–55], CMS (e) [56, 57], and LHCb Collaborations (f) [58–60]. The flavor of ℓℓ¯ is shown in the panels, see text in subsection 3.3 for the meanings of histograms which are quoted from QCD calculations.
(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χ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 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.


## 3.2. Tendency of Parameters

In order to study the underlying law that governs the variation of the extracted parameters withs and Q, 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 n, and Figures 5(c) and 5(d) show the trend of parameter T, 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 s. On average, it can be seen that as s increases, the parameter n increases quickly and then decreases slowly, and the parameter T increases slowly and then significantly. There is a knee point for the trend of n at s≈40–50 GeV. Meanwhile, there is a boundary at s≈200–500 GeV above which T increases significantly. Similarly, we compare Figures 5(b) and 5(d) and analyze the variation of parameters with Q. It can be clearly seen that the parameter n increases firstly and then decreases, and the parameter T increases slowly and then significantly, with the increase of Q. There is a knee point for the trend of n at Q≈14–15 GeV/c2. Meanwhile, there is a boundary at Q≈20–60 GeV/c2 above which T increases significantly.Figure 5
The variation of parametersn (a, b) and T (c, d) in the Lévy-Tsallis function with energy s (a, c) and invariant mass Q (b, d). The parameter values are taken from Figures 1–4 and recorded in Tables 1 and 2. Since the invariant mass grouping in Figure 1(e) is different from others, and the varying quantity in Figure 1(f) is the Feynman variable, Figure 5 does not include the parameters from Figures 1(e) and 1(f) to avoid trivialness.
(a)(b)(c)(d)Figure6 is similar to Figure 5, but shows the dependence of parameters ptS (Figures 6(a) and 6(b)), ptH (Figures 6(c) and 6(d)), and k (Figures 6(e) and 6(f)) on s(Figures 6(a), 6(c), and 6(e)) andQ (Figures 6(b), 6(d), and 6(f)) obtained from the two-component Erlang distribution. One can see that with increasing s, ptS and ptH increase slowly and then quickly, and k decreases slowly and then quickly. For s, there is a boundary between 60 and 500 GeV for ptS, between 500 and 1100 GeV for ptH, and between 200 and 500 GeV for k. Meanwhile, with increasing Q, ptS and ptH increase slowly and then quickly, and k decreases slowly and then quickly. There is a boundary at Q≈20–60 GeV/c2.Figure 6
Similar to Figure5, but showing the variation of parameters ptS (a, b), ptH (c, d), and k (e, f) in the two-component Erlang distribution with s (a, c, e) and Q (b, d, f).
(a)(b)(c)(d)(e)(f)In Figure7, 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 p1 has a slight tendency to increase with the increase of s and Q. In Figures 7(c) and 7(d), one can see that the parameter n1 decreases quickly and then slowly with the increase of s; n1 also decreases slowly and then quickly with increase of Q. There is a boundary for the trend of n1 at s≈40–200 GeV and Q≈20–50 GeV/c2. We note that the variation of parameter is obtained from the the average values of the extracted parameters, for each s or Q.Figure 7
Similar to Figure5, but showing the parameters p1 (a, b) and n1 (c, d) in the Hagedorn function with (a) (c) s and (b) (d) Q.
(a)(b)(c)(d)It should be pointed out that the values of the parameters in Figures5–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 pT 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 parametersT, ptS, ptH, and p1 show monotonous increasing trend when s 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 s 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 n1 decreases, with the increase of s 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 s 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 thepT 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 softest 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 ~40GeV, 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/c2. At below (above) the softest point, we see naturally smaller (larger) Q, except for Q=5–6 and 5–8 GeV/c2 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/c2, the system is close to the onset energy of phase transition. The parameters T, ptS, ptH, 1−k, and p1 show similar small and almost flat values, respectively, and n and n1 show different trends. At the energy of TeV and the invariant mass around 100 GeV/c2, the system is far away from the onset energy of phase transition. The parameters T, ptS, ptH, 1−k, and p1 show large values, and n and n1 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 massQ of ℓℓ¯) and obey renormalization group equations. Although this scaling is not included in the function f1pt1 or f1pT−pt1 of Eq. (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 f1pt1 or f1pT−pt1 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 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 Figures1(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 Figure2(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 (N3LO or N3LL) from the Drell-Yan data.The meanings of histograms in Figure3 and Figure 4(a), Figure 4(b), Figure 4(c) with 46≤Q<66 and 116≤Q<150 GeV/c2, 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/c2 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 N3LL. 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+N3LL.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 inp-Cu, π−-W, and pp (pp¯) 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 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 4<Q<150GeV/c2. 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 increasings, there is a knee point for the trend of n at s≈40–50 GeV. Meanwhile, there is a boundary at s≈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/c2. Meanwhile, there is a boundary at Q≈20–60 GeV/c2 above which T increases significantly. In the two-component Erlang distribution, there is a boundary at s≈60–500 GeV for ptS, 500–1100 GeV for ptH, and 200–500 GeV for k, above which ptS, ptH, and 1−k increase quickly. Meanwhile, there is a boundary at Q≈20–60 GeV/c2 above which ptS, ptH, and 1−k increase quickly. In the convolution of two Hagedorn functions, p1 increases obviously with increasing s and Q. There is a boundary for the trend of n1 at s≈40–200 GeV and Q≈20–50 GeV/c2.With increasings and Q, the parameters T, ptS, ptH, 1−k, and p1 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 n1 decreases, with the increase of s 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 s and Q increase. The fit results in this paper are comparable to the QCD NNLO and N3LL 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.



---

*Source: 9548737-2021-06-01.xml*