Averaged transverse momentum correlations of hadrons in relativistic heavy-ion collisions

In relativistic heavy-ion collisions, the hot nuclear matter is created at early collision stage by intensively inelastic collisions of colliding nucleons (Heinz and Jacob, 2000; Arsene et al., 2005; Adcox et al., 2005; Back et al., 2005; Adams et al., 2005; Gyulassy and McLerran, 2005). Subsequently, the matter expands, cools and finally decomposes into hadrons scattering out. The evolution of hot nuclear matter is a complex process governed by non-perturbative QCD and is mainly modeled by hydrodynamic models (Kolb and Heinz, 2003) and transport models (Bass et al., 1998; Bass and Dumitru, 2000; Lin et al., 2005) at present. Hadrons produced from hot nuclear matter always have certain transverse momentum $p_{T}$, the component of momentum which is perpendicular to the beam direction. The $p_{T}$ distribution of hadrons carries lots of information on hot nuclear matter such as thermalization, transverse collective flow generated by system expansion in both partonic and hadronic stage, and is an important physical observable in experiments of relativistic heavy-ion collisions.

Rich experimental data for $p_{T}$ spectra of identified hadrons at mid-rapidity are successively reported in heavy-ion collisions at RHIC and LHC over the past decade (Adcox et al., 2002a, b; Abelev et al., 2006, 2009a, 2009b; Aamodt et al., 2011; Abelev et al., 2013a, 2015, b, 2014; Adam et al., 2020; Adamczyk et al., 2017). Based on these experimental data, lots of studies on properties of hadronic $p_{T}$ distribution are carried out, which greatly improve people’s understandings on property of the created hot nuclear matter and the mechanism of hadron production in relativistic heavy-ion collisions (Fries et al., 2003; Greco et al., 2003; Hwa and Yang, 2003; Huovinen and Ruuskanen, 2006; Chen and Ko, 2006; Gyulassy et al., 2004; Lokhtin and Snigirev, 2006; Zhang et al., 2004; Tang et al., 2009; Song et al., 2011; Cassing and Bratkovskaya, 2008; Karpenko et al., 2013). The average transverse momentum ($\left\langle p_{T}\right\rangle$) of hadrons is obtained by integrating over $p_{T}$ spectra of hadrons. It is dominated by property of hadronic $p_{T}$ spectra in the low $p_{T}$ range and therefore it reflects the property of soft hadrons and correspondingly that of hot nuclear matters.

In this paper, we study the property of $\left\langle p_{T}\right\rangle$ of identified hadrons produced in relativistic heavy-ion collisions. We compile experimental data for $\left\langle p_{T}\right\rangle$ of $\phi$, proton, $\Lambda$, $\Xi^{-}$ and $\Omega^{-}$ at mid-rapidity in Au+Au collisions at $\sqrt{s_{NN}}=$ 200, 39, 27, 19.6, 11.5, 7.7 GeV and in Pb+Pb collisions at $\sqrt{s_{NN}}=$ 2.76 TeV. We search the regularity in $\left\langle p_{T}\right\rangle$ data of these hadrons and, in particular, their dependence on hadron species and collision energy. We discuss what the underlying physics is responsible to the observed regularity. In particular, we study the effect of hadronization by an equal-velocity combination (EVC) mechanism of quarks and antiquarks (Song et al., 2017, 2021a) in explaining the experimental data of $\left\langle p_{T}\right\rangle$. We derive analytic expression for $\left\langle p_{T}\right\rangle$ of identified hadrons in EVC mechanism so as to give clear quark flavor dependence of hadronic $\left\langle p_{T}\right\rangle$ and provide intuitive explanations of experimental data.

The paper is organized as follows. In Sec. II, we briefly introduce our EVC model and derive $\left\langle p_{T}\right\rangle$ of identified hadrons for the simplified quark distributions at hadronization. In Sec. III, we show our finding for the systematic correlation among experimental data for $\left\langle p_{T}\right\rangle$ of hadrons in relativistic heavy-ion collisions and give an intuitive explanation using our EVC model. In Sec. IV, we show the regularity on $\left\langle p_{T}\right\rangle$ of hadrons in central heavy-ion collisions as the function of $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$. In Sec. V, we discuss the influence of resonance decays on $\left\langle p_{T}\right\rangle$ correlations of hadrons. Finally, the summary and discussion are given in Sec. VI.

In this section, we apply a particular quark combination model (Song et al., 2017; Gou et al., 2017) to describe the production of hadrons at hadronization and derive $\left\langle p_{T}\right\rangle$ of different hadrons. The model assumes the constituent quarks and antiquarks as the effective degrees of freedom for the final parton system created in collisions at hadronization stage. Based on constituent quark model of internal structure of hadrons at low energy scale, the model assumes that the equal velocity combination of these constituent quarks and antiquarks is main feature of hadron formation at hadronization. The quark masses are taken as the constituent masses so that the equal velocity combination of these constituent quarks and antiquarks can correctly construct the on-shell hadron. At mid-rapidity (i.e., taking $y=0$), $p_{T}$ distribution of a hadron ($dN/dp_{T}$) is the product of those of (anti-)quarks

	$\displaystyle f_{B_{i}}\left(p_{T}\right)$	$\displaystyle=\kappa_{B_{i}}f_{q_{1}}\left(x_{1}p_{T}\right)f_{q_{2}}\left(x_{2}p_{T}\right)f_{q_{3}}\left(x_{3}p_{T}\right),$		(1)
	$\displaystyle f_{M_{i}}\left(p_{T}\right)$	$\displaystyle=\kappa_{M_{i}}f_{q_{1}}\left(x_{1}p_{T}\right)f_{\bar{q}_{2}}\left(x_{2}p_{T}\right).$		(2)

Here, moment fraction $x_{i}=m_{i}/(m_{1}+m_{2}+m_{3})$ ($i=1,2,3)$ in baryon formula satisfies $x_{1}+x_{2}+x_{3}=1$ and $x_{i}=m_{i}/(m_{1}+m_{2})$ ($i=1,2)$ in meson formula satisfies $x_{1}+x_{2}=1$. $m_{i}$ is constituent mass of quark $q_{i}$ and we take $m_{u}=m_{d}=0.3$ GeV and $m_{s}=0.5$ GeV. Coefficients $\kappa_{B_{i}}$ and $\kappa_{M_{i}}$ are independent of $p_{T}$, see Refs. (Song et al., 2021a) for their detailed expressions, and therefore are not involved in derivation of averaged transverse momentum $\left\langle p_{T}\right\rangle$ in the following text.

The value of $\left\langle p_{T}\right\rangle$ is dominated by the $p_{T}$ spectrum of particles in the low $p_{T}$ range. Therefore, here we focus on $p_{T}$ spectra of quarks and antiquarks with low $p_{T}$. Unfortunately, quarks of low $p_{T}$ is governed by non-perturbative QCD dynamics and their distributions are difficult to calculate from first principle. Considering that experimental data for $p_{T}$ spectra of hadrons at mid-rapidity in the low $p_{T}$ range in relativistic heavy-ion collisions are generally well fitted by exponential function and/or Boltzmann distribution (Abelev et al., 2015, 2013a, 2009a, 2009b; Adamczyk et al., 2017; Adam et al., 2020), in this paper we take the following parameterization for quark $p_{T}$ spectra at mid-rapidity

$$f_{q_{i}}\left(p_{T}\right)=\mathcal{N}p_{T}^{k}\exp\left[-\frac{\sqrt{p_{T}^{2}+m_{i}^{2}}}{T_{i}}\right],$$

(3)

which is convenient to derive analytic results of hadronic $\left\langle p_{T}\right\rangle$. Here, $\mathcal{N}$ is coefficient to quantity the number of $q_{i}$, which is irrelevant to $\left\langle p_{T}\right\rangle$ calculations. $T_{i}$ is the slope parameter to quantify the exponential decrease of the spectrum. Exponent $k$ tunes the behavior of the spectrum at small $p_{T}$. In the case of two-dimensional Boltzmann distribution in the rest frame we have $k=1$, and in the one-dimensional case we have $k=0$. If we directly apply Eq. (3) to fit the experimental data of $p_{T}$ spectra of hadrons by Eq. 1, we should take $k\approx 1/3$ to properly describe baryon and $k\approx 1/2$ to properly describe meson (mainly $\phi$). In addition, the effect of strong collective radial flow should be included in the quark spectrum in the laboratory frame, which is dependent on collision energies in relativistic heavy-ion collisions. With these considerations, we take $k$ as an relatively-free parameter in the range $[0,1]$ in this study of hadronic $\left\langle p_{T}\right\rangle$ in relativistic heavy-ion collisions at RHIC and LHC energies .

Substituting Eq. (3) into Eqs. (1) and (2), we obtain

$$\left\langle p_{T}\right\rangle_{B}=\frac{\int f_{B}\left(p_{T}\right)p_{T}dp_{T}}{\int f_{B}\left(p_{T}\right)dp_{T}}=\frac{\int p_{T}^{3k+1}\exp\left[-\left(\frac{x_{1}}{T_{1}}+\frac{x_{2}}{T_{2}}+\frac{x_{3}}{T_{3}}\right)\sqrt{p_{T}^{2}+m_{B}^{2}}\right]dp_{T}}{\int p_{T}^{3k}\exp\left[-\left(\frac{x_{1}}{T_{1}}+\frac{x_{2}}{T_{2}}+\frac{x_{3}}{T_{3}}\right)\sqrt{p_{T}^{2}+m_{B}^{2}}\right]dp_{T}},$$

(4)

$$\left\langle p_{T}\right\rangle_{M}=\frac{\int f_{M}\left(p_{T}\right)p_{T}dp_{T}}{\int f_{M}\left(p_{T}\right)dp_{T}}=\frac{\int p_{T}^{2k+1}\exp\left[-\left(\frac{x_{1}}{T_{1}}+\frac{x_{2}}{T_{2}}\right)\sqrt{p_{T}^{2}+m_{M}^{2}}\right]dp_{T}}{\int p_{T}^{2k}\exp\left[-\left(\frac{x_{1}}{T_{1}}+\frac{x_{2}}{T_{2}}\right)\sqrt{p_{T}^{2}+m_{M}^{2}}\right]dp_{T}}$$

(5)

where $m_{B}=m_{q_{1}}+m_{q_{2}}+m_{q_{3}}$ and $m_{M}=m_{q_{1}}+m_{\bar{q}_{2}}$. We use the integral formula

$$\int_{0}^{\infty}p_{T}^{n}\exp\left[-a\sqrt{p_{T}^{2}+m^{2}}\right]dp_{T}=m^{n+1}\frac{2^{n/2}\Gamma(\frac{n+1}{2})}{\sqrt{\pi}}\frac{K_{n/2+1}(\alpha)}{\alpha^{n/2}}$$

(6)

where $\alpha=a\,m$, $\Gamma(z)$ is Gamma function and $K_{n}(z)$ is the modified Bessel function of the second kind. We obtain

$\displaystyle\left\langle p_{T}\right\rangle_{B}$

$\displaystyle=\left(m_{q_{1}}+m_{q_{2}}+m_{q_{3}}\right)\sqrt{\frac{2}{\alpha_{B}}}\frac{\Gamma(3k/2+1)}{\Gamma(3k/2+\frac{1}{2})}\frac{K_{3k/2+3/2}(\alpha_{B})}{K_{3k/2+1}(\alpha_{B})},$

(7)

$\displaystyle\left\langle p_{T}\right\rangle_{M}$

$\displaystyle=\left(m_{q_{1}}+m_{\bar{q}_{2}}\right)\sqrt{\frac{2}{\alpha_{M}}}\frac{\Gamma(k+1)}{\Gamma(k+\frac{1}{2})}\frac{K_{k+3/2}(\alpha_{M})}{K_{k+1}(\alpha_{M})}$

(8)

with

$$\alpha_{B}=\frac{m_{q_{1}}}{T_{1}}+\frac{m_{q_{2}}}{T_{2}}+\frac{m_{q_{3}}}{T_{3}}=\alpha_{q_{1}}+\alpha_{q_{2}}+\alpha_{q_{3}},$$

(9)

and

$$\alpha_{M}=\frac{m_{q_{1}}}{T_{1}}+\frac{m_{\bar{q}_{2}}}{T_{2}}=\alpha_{q_{1}}+\alpha_{\bar{q}_{2}}.$$

(10)

As shown by these expressions, $\left\langle p_{T}\right\rangle$ of different hadrons is correlated by the simple combination of slope parameter $\alpha_{q_{i}}$ of quarks at hadronization.

Bessel function $K_{\nu}(\alpha)$ and Gamma function $\Gamma(z)$ usually have complex expressions. Here, we present the numerical approximations for $\left\langle p_{T}\right\rangle_{B}$ and $\left\langle p_{T}\right\rangle_{M}$

	$\displaystyle\left\langle p_{T}\right\rangle_{B}$
	$\displaystyle\approx\left(m_{q_{1}}+m_{q_{2}}+m_{q_{3}}\right)\left(0.26+0.024k+\frac{0.96+2.99k}{\alpha_{B}}\right),$		(11)
	$\displaystyle\left\langle p_{T}\right\rangle_{M}\approx\left(m_{q_{1}}+m_{\bar{q}_{2}}\right)\left(0.25+0.03k+\frac{0.97+1.99k}{\alpha_{M}}\right)$		(12)

in order to see their dependence on $\alpha$ and $k$ in a numerically intuitive way. The relative errors of two approximations are less than about $3\%$ for the physical range of $\left\langle p_{T}\right\rangle_{B}$ and $\left\langle p_{T}\right\rangle_{M}$ in heavy-ion collisions at RHIC and LHC energies studied in this paper.

In this section, we study the correlation among $\left\langle p_{T}\right\rangle$ of different hadrons. In Fig. 1 (a), we present $\left\langle p_{T}\right\rangle$ of proton as horizontal axis and $\left\langle p_{T}\right\rangle$ of $\phi$ at correspondingly collision energy and centrality as the vertical axis to study the correlation between them. As we know, the mass of proton is close to that of $\phi$ but the quark flavor composition of proton ($uud$) is completely different from that of $\phi$ ($s\bar{s}$). Except a few datum points in central Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV, we see an relatively stable correlation between $\left\langle p_{T}\right\rangle$ of proton and that of $\phi$. In Fig. 1 (b) and (c), we show correlations among proton, $\Lambda$ and $\Xi^{-}$ in the manner of successive strangeness. In Fig. 1 (d), we show the correlation between $\phi$ and $\Xi^{-}$ which both have two strange (anti-)quarks. We see the systematic correlations among these hadrons. In Fig. 1 (e-h), we show the correlation among $\left\langle p_{T}\right\rangle$ data of $\phi$ and anti-baryons and also find systematic correlations among them.

As we know, the hot nuclear matters created in heavy-ion collisions at different collision energies and centralities have different size and geometry, different evolution times in partonic phase and in the subsequent hadronic re-scattering stage, etc. The correlations shown in Fig. 1 seem to assign these difference into a systematic manner and therefore may indicate some underlying physics which is universal in heavy-ion collision at both RHIC and LHC energies. We think that the universal hadronization mechanism may be a possible physical reason.

Therefore, we apply the EVC model in Sec. II to understand the above correlations in experimental data of hadronic $\left\langle p_{T}\right\rangle$. Hadrons in Fig. 1 are all made up of up, down, strange quarks and their antiquarks. In our model, exponent parameter $k$ and slope parameters $\alpha_{u}$, $\alpha_{d}$, $\alpha_{s}$, $\alpha_{\bar{u}}$, $\alpha_{\bar{d}}$ and $\alpha_{\bar{s}}$ are needed to fix in order to calculate $\left\langle p_{T}\right\rangle$ of hadrons according to Eqs. (7-10). Here, we can assume the isospin symmetry between up and down quarks $\alpha_{u}=\alpha_{d}$ at mid-rapidity in relativistic heavy-ion collisions. We also assume the charge conjugation symmetry for strange quarks and antiquarks $\alpha_{s}=\alpha_{\bar{s}}$, which is found to be a good approximation at mid-rapidity in relativistic heavy-ion collisions at RHIC and LHC energies (Song et al., 2021a). Finally, only three slope parameters $\alpha_{u}$, $\alpha_{\bar{u}}$ and $\alpha_{s}$ are left in addition to the exponent parameter $k$.

If we know the correlation between $\alpha_{u}$ and $\alpha_{s}$, we can calculate correlations among $\left\langle p_{T}\right\rangle$ of baryons and $\phi$, and therefore use them to explain experimental data in Fig. 1(a-d). Applying Eqs. (7) and (8) to proton and $\phi$, we have

	$\displaystyle\left\langle p_{T}\right\rangle_{p}$	$\displaystyle=3m_{u}\sqrt{\frac{2}{\alpha_{p}}}\frac{\Gamma(3k/2+1)}{\Gamma(3k/2+\frac{1}{2})}\frac{K_{3k/2+3/2}(\alpha_{p})}{K_{3k/2+1}(\alpha_{p})},$		(13)
	$\displaystyle\left\langle p_{T}\right\rangle_{\phi}$	$\displaystyle=2m_{s}\sqrt{\frac{2}{\alpha_{\phi}}}\frac{\Gamma(k+1)}{\Gamma(k+\frac{1}{2})}\frac{K_{k+3/2}(\alpha_{\phi})}{K_{k+1}(\alpha_{\phi})},$		(14)

where

	$\displaystyle\alpha_{p}$	$\displaystyle=2\alpha_{u}+\alpha_{d}=3\alpha_{u},$		(15)
	$\displaystyle\alpha_{\phi}$	$\displaystyle=\alpha_{s}+\alpha_{\bar{s}}=2\alpha_{s}$		(16)

according to Eqs. (9) and (10). By fitting experimental data of $\left\langle p_{T}\right\rangle$ of proton and $\phi$ shown in Fig. 1(a), we can reversely extract the correlation between $\alpha_{u}$ and $\alpha_{s}$, which can be parameterized as

$$\alpha_{u}(\alpha_{s})=c_{0}+c_{1}\alpha_{s}$$

(17)

where $c_{0}$ and $c_{1}$ are two coefficients. Because the extraction is also dependent on exponent parameter $k$, we list values of $c_{0}$ and $c_{1}$ at several different $k$ in Table. 1. The resulting fitted correlations between proton and $\phi$ with these different extractions are shown as lines with different types in Fig. 2 (a). In order to reduce to the bias in choice of $k$, we let these different extractions to represent the same correlation between $\left\langle p_{T}\right\rangle$ of proton and $\phi$, i..e, these lines are coincident with each other. By fitting experimental data of $\left\langle p_{T}\right\rangle$ of anti-proton and $\phi$, we also obtain the correlation between $\alpha_{\bar{u}}$ and $\alpha_{s}$ with the parameterized form

$$\alpha_{\bar{u}}(\alpha_{s})=d_{0}+d_{1}\alpha_{s}.$$

(18)

The values of coefficients $d_{0}$ and $d_{1}$ at several $k$ are shown in Table. 1 and the corresponding fit is shown in Fig. 2 (b). At low RHIC energies where $\alpha_{s}$ is large, $\alpha_{\bar{u}}$ is different from $\alpha_{u}$ to a certain extent, which is because the finite baryon density at low collision energies will cause the asymmetry between up/down quarks and their antiquarks.

With these two relationship, we can calculate correlations among $\left\langle p_{T}\right\rangle$ of various hadrons by Eqs. (7) and (8) with Eqs. (15), (16) and

	$\displaystyle\alpha_{\Lambda}$	$\displaystyle=2\alpha_{u}+\alpha_{s},$		(19)
	$\displaystyle\alpha_{\Xi}$	$\displaystyle=2\alpha_{s}+\alpha_{u},$		(20)
	$\displaystyle\alpha_{\Omega}$	$\displaystyle=3\alpha_{s}$		(21)

and the corresponding anti-baryons according to Eqs. (9) and (10). In Fig. 3, we present theoretical results for $\left\langle p_{T}\right\rangle$ correlations of $p\Lambda$, $\Lambda\Xi^{-}$, $\Xi^{-}\Omega$ and $\Xi^{-}\phi$ pairs, and these of anti-baryons, and we compare them with experimental data. Theoretical results at different $k$ with corresponding coefficients in Table. 1 are shown as lines with different types. These lines are different to a certain extent, which shows the theoretical uncertainties due to the selection of exponent parameter $k$. Overall, we see that the systematic feature of the correlations exhibited by $\left\langle p_{T}\right\rangle$ data of these hadrons can be described by our theoretical model.

This result is quite interesting. Here we only consider the effect of hadronization by EVC mechanism without any considerations on other dynamical ingredients such as system size and geometry, evolution time, hadronic re-scattering, etc. We run the event generators URQMD 3.4 Bass et al. (1998) and AMPT 2.26 (1.26) Lin et al. (2005) which practically include those dynamical processes and we do not find better description on systematic correlations in Fig. 1 when we use default parameter values of event generators. Therefore, our results in Fig. 3 indicate the important role of hadronization by EVC mechanism in describing $\left\langle p_{T}\right\rangle$ of those hadrons in relativistic heavy-ion collisions at both RHIC and LHC energies.

Experimental data for hadronic $\left\langle p_{T}\right\rangle$ shown in the form of Fig. 1 reveal correlations among $\left\langle p_{T}\right\rangle$ of different hadrons, which are mainly relevant to hadronization mechanism according to our studies in the previous section. In this section, we study another aspect of $\left\langle p_{T}\right\rangle$ of hadrons, i.e., their absolute values, and search some regularity underlying these experimental data in relativistic heavy-ion collisions at RHIC and LHC energies.

There are many physical ingredients that influence $\left\langle p_{T}\right\rangle$ of hadrons. Generally speaking, there are two main sources of generating the transverse momentum of hadrons. The first source is the intensive parton interactions at early collision stage which form the thermal bulk nuclear matter and generate primordial thermal or stochastic momentum of particles. Another source is the expansion of hot nuclear matter in both partonic phase (if exists) and hadronic phase, which generates the collective radial flow and therefore strengths the $\left\langle p_{T}\right\rangle$ of hadrons. The effects of two sources are both influenced by collision parameters such as collision energy, collision centrality and collision system. These collision parameters influence the size and geometry of the bulk nuclear matter, the intensity of soft parton/particle interactions, the time of system expansion and correspondingly the magnitude of collective radial flow. In view of these complex ingredients, it seems to be difficult to find a simple and perfect regularity for $\left\langle p_{T}\right\rangle$ of hadrons by directly analyzing the experimental data of $\left\langle p_{T}\right\rangle$ of hadrons in relativistic heavy-ion collisions.

Here, we try to take $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$ to quantify the excitation of hadronic $\left\langle p_{T}\right\rangle$. $dN_{ch}/dy$ is the rapidity density of charged particles at mid-rapidity. It can characterize the size of the created hot nuclear matter. In general case, i.e., as other conditions (such as collision energy and collision system) are not changed, the larger system means more intensive particle excitation (i.e., higher stochastic momentum or higher temperature) and more expansion (i.e., more radial flow). Experimental observation have shown that $\left\langle p_{T}\right\rangle$ of hadrons generally positively responds to the $dN_{ch}/dy$ at given collision energy and collision system (Olimov et al., 2021; Petrovici et al., 2017; Abelev et al., 2015, 2013a; Estienne, 2005; Abelev et al., 2009a, b; Adamczyk et al., 2017; Adam et al., 2020). Therefore, we take $dN_{ch}/dy$ as the main relevant ingredients parameterizing $\left\langle p_{T}\right\rangle$ of hadrons. $N_{part}$ is the number of participant nucleons calculated in Glauber model (Miller et al., 2007), which depends on collision energy and collision system and impact parameter. It can characterize the total amount of energy deposited in the collision region and therefore characterize the initial size and energy density of the created nuclear matter. The ratio $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$ quantifies the average number of charged particles produced by a pair of participant nucleons. It can roughly characterize the average number of charged particles produced by an “unit” effective energy deposited by the collision of a pair of nucleons. In general, the higher $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$ means more intensive particle excitation which needs more intensive parton interactions and also means more momentum generation. Therefore, we expect $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$ should positively correlate with $\left\langle p_{T}\right\rangle$ of hadrons.

The geometry property of hot nuclear matter, mainly controlled by impact parameter, also influence $\left\langle p_{T}\right\rangle$ of hadrons. In particular, in peripheral collisions where impact parameter is large, various-order anisotropic flows are generated and will influence the $\left\langle p_{T}\right\rangle$ to a certain extent by, for example, asymmetric distribution of $p_{x}$ and $p_{y}$. Therefore, in order to remove the freedom of impact parameter which is quite complex to parameterize, we only use experimental data of hadronic $\left\langle p_{T}\right\rangle$ in central collisions to search their possible regularity with respect to $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$.

In Fig. 4, we compile experimental data for $\left\langle p_{T}\right\rangle$ of $\phi$ and (anti-)baryons at mid-rapidity in central heavy-ion collisions at different collision energies. We see that these data of hadronic $\left\langle p_{T}\right\rangle$ exhibit a clear regularity when we plot them as the function $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$. Here data of $dN_{ch}/dy$ at mid-rapidity and $N_{part}$ are taken from Refs. (Abelev et al., 2013a; Aamodt et al., 2011; Abelev et al., 2009b; Adamczyk et al., 2017). In calculation of $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$, only experimental uncertainties of $dN_{ch}/dy$ are included.

According to behavior of experimental data in Fig. 4 and the approximated formula of hadronic $\left\langle p_{T}\right\rangle$ in Eqs. (11) and (12), we parameterize the $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$ dependence of slope parameter $\alpha_{q}$ of quarks as the following form

$\displaystyle\alpha_{q}$

$\displaystyle=\left[g_{q}+h_{q}\left(\frac{dN_{ch}/dy}{N_{part}/2}\right)^{2/3}\right]^{-1}.$

(22)

Coefficients $g$ and $h$ of $u$, $\bar{u}$ and $s$ quarks can be fixed by using Eqs.(7) and (8) to fit experimental data of proton, anti-proton and $\Omega^{-}+\bar{\Omega}^{+}$. Values of $g$ and $h$ at different $k$ are shown in Table. 2. These fittings to experimental data of proton, anti-proton and $\Omega^{-}+\bar{\Omega}^{+}$ are shown in Fig. 5 (a), (e) and (h) as lines of different types. In order to avoid the bias in selection of exponent parameter $k$, we let these different fitting groups generate the same $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$ dependence for $\left\langle p_{T}\right\rangle$ of proton, anti-proton and $\Omega^{-}+\bar{\Omega}^{+}$, that is, lines of different types are coincident with each other in Fig. 5 (a), (e) and (h). This treatment can enable us to study theoretical uncertainty in prediction of other hadrons. Results for $\left\langle p_{T}\right\rangle$ of other hadrons as the function of $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$ are shown in Fig. 5(b-d) and (f-g) and are compared with experimental data.

We see that theoretical results of hyperons at different $k$ in Fig. 5(b-c) and (f-g) are coincident with each other and they are in good agreement with experimental data. Results of $\phi$ in Fig. 5(d) are dependent on $k$ to a certain extent. In these results at different $k$, we see that the results of $\phi$ at $k=0.1,0.3$ are globally better than others. This feature is similar with that in correlations of hadronic $\left\langle p_{T}\right\rangle$ in Fig. 3.

In previous calculations, results of hadronic $\left\langle p_{T}\right\rangle$ are those for initially produced hadrons by hadronization and effects of resonance decay are not yet included. For proton, $\Lambda$ and $\Xi^{-}$, a certain fraction of these hadrons observed in experiments comes from decay of higher-mass resonances such as those from octet baryons $\Delta$, $\Sigma^{*}$ and $\Xi^{*}$, respectively. $\Omega^{-}$ and $\phi$ are generally expected to be less influenced. In this section, we study the influence of resonance decay on the correlations among $\left\langle p_{T}\right\rangle$ of different hadrons.

We apply the quark combination model developed in previous works (Song et al., 2017; Gou et al., 2017) to calculate the influence of resonance decay on $\left\langle p_{T}\right\rangle$ of hadrons. Following experimental corrections, results of $\Lambda$ and $\Xi^{-}$ do not include weak decay contributions but results of proton and anti-proton include them. We adopt the following strategy to quantify the effect of resonance decays. First, we use the model to calculate the $\left\langle p_{T}\right\rangle$ of final-state (anti-)proton and that of $\phi$ with the parameterized quark $p_{T}$ spectra in Eq. (3). We apply the model to fit the $\left\langle p_{T}\right\rangle$ correlation between experimental data of (anti-)protons and those of $\phi$ to obtain the correlation between $\alpha_{u}$ and $\alpha_{s}$ with the parameterization form Eq. (17) and that between $\alpha_{\bar{u}}$ and $\alpha_{s}$ with Eq. (18). The newly obtained coefficients $c_{0}$, $c_{1}$, $d_{0}$, $d_{1}$ are slightly different from those in Table 1 due to the effect of resonance decays. In the fitting process, we keep the same $p\phi$ and $\bar{p}\phi$ correlations shown as lines in Fig 2. Second, we calculate $\left\langle p_{T}\right\rangle$ correlations among other hadron pairs and compare them with results in Fig. 3 to study the effect of resonance decays.

In Fig. 6, we show $\left\langle p_{T}\right\rangle$ correlations among different final-state hadrons at given slope parameter $k=0.3$ and compare them, dashed lines, with results of directly-produced hadrons at same $k$, the dot-dashed lines, which are borrowed from Fig. 3. Experimental data (Abelev et al., 2013b, 2014, 2015, a; Estienne, 2005; Abelev et al., 2009a, b; Adamczyk et al., 2017; Adam et al., 2020) are also presented. We see that two sets of results are only slightly different, which indicates the weak influence of resonance decays on $\left\langle p_{T}\right\rangle$ correlations of these hadrons. Results at other values of $k$ are quite similar and therefore not presented.

In the similar way, we further study the effect of resonance decays on hadronic $\left\langle p_{T}\right\rangle$ as the function of $\left(dN_{ch}/dy\right)/\left(0.5N_{part}\right)$ in central heavy-ion collisions. In Fig. 7, the dashed lines denotes results for $\left\langle p_{T}\right\rangle$ of hadrons with resonance decays and dot-dashed lines denote results without resonance decays. Symbols are experimental data (Abelev et al., 2013b, 2014, 2015, a; Estienne, 2005; Abelev et al., 2009a, b; Adamczyk et al., 2017; Adam et al., 2020; Aamodt et al., 2011). Experimental data of protons (a), anti-protons (e) and $\Omega^{-}$ (h) are used to determine parameters of (anti-)quarks at hadronization, and we keep the same extent in reproducing these experimental data, i.e., dashed lines are coincident with the dot-dashed lines in three panels. In Fig. 7 (b-d) and (f-g), we see the difference between dashed lines and dot-dashed lines is small, which indicate the weak influence of resonance decays on $\left\langle p_{T}\right\rangle$ correlations of these hadrons.

In this paper, we have applied a quark combination model with equal-velocity combination approximation to study the averaged transverse momentum ($\left\langle p_{T}\right\rangle$) correlations of proton, $\Lambda$, $\Xi^{-}$, $\Omega^{-}$ and $\phi$ in relativistic heavy-ion collisions. We derived analytic formulas of hadronic $\left\langle p_{T}\right\rangle$ in the case of exponential form of quark $p_{T}$ spectra at hadronization, which can clarify correlations among $\left\langle p_{T}\right\rangle$ of identified hadrons based on the constituent quark structure of hadrons. We used these analytic formulas to explain the systematic correlations exhibited in $\left\langle p_{T}\right\rangle$ data of $p\Lambda$, $\Lambda\Xi^{-}$, $\Xi^{-}\Omega^{-}$ and $\Xi^{-}\phi$ pairs and those of anti-baryons. We discussed the regularity for $\left\langle p_{T}\right\rangle$ of these hadrons as the function of $(dN_{ch}/dy)/(N_{part}/2)$ at mid-rapidity in central heavy-ion collisions at both RHIC and LHC energies, and used our model to self-consistently explain $\left\langle p_{T}\right\rangle$ of these hadrons as the function of $(dN_{ch}/dy)/(N_{part}/2)$. In these studies, we use experimental data of (anti-)protons to fix the property of up/down (anti-)quarks and those of $\Omega$ or $\phi$ to fix that of strange quarks at hadronization. Then we predict correlations among $\left\langle p_{T}\right\rangle$ of other hadron pairs and compare with experimental data to test the theoretical consistency. Moreover, we studied the effects of resonance decays on $\left\langle p_{T}\right\rangle$ correlations of hadrons and find they are weak in comparison with hadronization.

Our studies shown that the $\left\langle p_{T}\right\rangle$ correlations among experimental data of proton, $\Lambda$, $\Xi^{-}$, $\Omega^{-}$ and $\phi$ in Au+Au collisions at RHIC energies and Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV can be self-consistently described by the equal-velocity combination mechanism of constituent quarks and antiquarks at hadronization. This indicates the important role of constituent quarks and antiquarks as the effective degrees of freedom of the hot nuclear matter at hadronization stage and their equal-velocity combination as the main feature of hadron formation in relativistic heavy-ion collisions. The current study is consistent with our previous works in studying the elliptic flow of these hadrons and the quark number scaling property of $p_{T}$ spectra of $\Omega^{-}$ and $\phi$ in relativistic heavy-ion collisions at RHIC and LHC energies using the same quark combination mechanism (Song et al., 2020, 2021b).

In relativistic heavy-ion collisions, re-scatterings of hadrons after hadronization will influence momentum of hadrons to a certain extent. For example, the signal of $\phi$ may be lost by the scattering of their decay daughters with the surrounding hadrons and $\phi$ may be also generated by the coalescence of two kaon. We will study this hadronic rescattering effect in the future works. In addition, we will also carry out a systematic study on $p_{T}$ spectra of identified hadrons at mid-rapidity in different centralities in Au+Au collisions at STAR BES energies. $p_{T}$ spectra of identified hadrons contain more dynamical information than their $\left\langle p_{T}\right\rangle$, which can be used to further test our quark combination model at low RHIC energies.

We thank Prof. X. L. Zhu for providing us the experimental data of hadronic $\left\langle p_{T}\right\rangle$ in Au+Au collisions at RHIC energies. This work is supported in part by Shandong Provincial Natural Science Foundation (ZR2019YQ06, ZR2019MA053), the National Natural Science Foundation of China under Grant No. 11975011, and Higher Educational Youth Innovation Science and Technology Program of Shandong Province (2019KJJ010).