Using the Tsallis distribution for hadron spectra in $pp$ collisions: Pions and quarkonia at $\sqrt{s}$ = 5--13000 GeV

TL;DR

This paper introduces a thermal Tsallis-based model to accurately describe hadron spectra in proton-proton collisions across a wide energy range, providing predictions for future experiments.

Contribution

It presents a novel combination of the Tsallis distribution and blast-wave model to fit and predict hadron spectra in high-energy $pp$ collisions.

Findings

Successfully describes experimental pion and quarkonia spectra from 5 GeV to LHC energies.

Provides parametrizations for energy dependence of model parameters.

Enables predictions for particle yields at new collision energies.

Abstract

A thermal model, based on the Tsallis distribution and blast-wave model, is proposed to compute hadron double-differential spectra $d^{2}N/dp_{T}dy$ in $pp$ (also high energy $p\overline{p}$) collisions. It successfully describes the available experimental data on pion and quarkonia ($\phi$, $J/\psi$, $\psi(2S)$, $\Upsilon$ family) production at energies from $\sqrt{s} = 5$ GeV to the LHC ones. Simple parametrizations for the $\sqrt{s}$ dependence of the model parameters are provided allowing predictions for the yields of these particles at new collision energies. The model can be used also for the pion Bose-Einstein correlation studies.

Tables

Table 1. Table 1: Parameters of Eqs. ( 18 ) − - ( 21 ) for π ± , 0 superscript 𝜋 plus-or-minus 0 \pi^{\pm,0} , ϕ italic-ϕ \phi , J / ψ J 𝜓 {\rm J}/\psi , ψ ​ ( 2 ​ S ) 𝜓 2 S \psi({\rm 2S}) and Υ ​ ( 1 ​ S ) Υ 1 S \rm\Upsilon(1S) , obtained from the combined fits of data measured at different energies. The χ 2 superscript 𝜒 2 \chi^{2} and N ​ D ​ F 𝑁 𝐷 𝐹 NDF correspond to the fits when all the parameters, except V ~ ~ 𝑉 \widetilde{V} , are fixed to their central values. Additional parameters for higher Υ Υ \rm\Upsilon states and non-prompt J / ψ J 𝜓 {\rm J}/\psi and ψ ​ ( 2 ​ S ) 𝜓 2 S \psi({\rm 2S}) production are given in Sec.  V .

	$\pi$	$\varphi$	$J/\psi$	$\psi \left(2S\right)$	${\rm Y}\left(1S\right)$
$\tilde{V}$ (${\mathrm{GeV}}^{-3}$)	5030.9 $\pm$7.2	561.2 $\pm$5.2	107.4 $\pm$0.4	19.2 $\pm$0.1	0.241 $\pm$0.001
$e_1$ (GeV)	12.5	12.5	7.8	7.8	30.0
$p_1$	3.5	3.5	76.1	128.1	0
$p_2$	2.3	3.1	0	0	0
$p_3$	135.6	166.3	56.3	56.3	3.0
$p_4$	0	0	29.2	29.2	8.6
$p_5$	46.7	46.7	87.9	94.0	$-$4.3
$e_2$ (GeV)	-	8786	13900	35500	1$\cdot {10}^9$
$e_3$ (GeV)	-	225	63.1	63.1	16000
$p_6$	-	0.047	0.072	0.060	0.058
$p_7$	-	2.30	2.01	0.50	3.50
${\chi }^2$	8716	235	5985	1567	5773
$NDF$	2293	223	1707	974	622

Equations

$$E\frac{d^{3}N}{d^{3}p}=\frac{g}{(2\pi {)}^3}{p}_{\nu }{\int }_{{\Sigma }_f}{d}^3{\Sigma }^{\nu }f(X),$$

$$\begin{array}{rl}& N=\frac{g}{(2\pi {)}^3}{V}_f\int d^{3}pf(\frac{E-\mu }{T}),\\ & V_f={\int }_{{\Sigma }_f}{d}^3{\Sigma }^{\nu }{u}_{\nu }(x).\end{array}$$

$$u^{\nu }={\gamma }_r(\cosh \eta ,v_r\cos \varphi ,v_r\sin \varphi ,\sinh \eta ),$$

$$p_{\nu }{u}^{\nu }={\gamma }_r[m_{\mathrm{T}}\cosh (y-\eta )-v_r{p}_{\mathrm{T}}\cos ({\varphi }_p-\varphi )].$$

$$d^3{\Sigma }^{\nu }={\tau }_f(\cosh \eta ,0,0,\sinh \eta )d\eta d^{2}r.$$

$$R(\eta )=R_0\sqrt{1-{\eta }^2/{\eta }_{max}^2}$$

$$v_r(r)=v_s\cdot (r/R_0{)}^2$$

$$⟨v_r⟩=\frac{1}{V_f}{\int }_{{\Sigma }_f}{d}^3{\Sigma }^{\nu }{u}_{\nu }(x)v_r(r).$$

$$V_f$$

$$⟨v_r⟩=\frac{3}{2}\frac{V_0}{V_f}{\int }_0^{1}dx(1-3x^2)\arcsin [v_s(1-x^2)].$$

$$f(X)={\left[[1+(q-1)X{]}^{\frac{1}{q-1}}-\xi \right]}^{-q},$$

$$E\frac{d^{3}N}{d^{3}p}=\frac{d^{3}N}{d^2{p}_{\mathrm{T}}dy}=\frac{g{\tau }_f}{(2\pi {)}^3}\sum _{k=0}^{\infty }{\xi }^k\left(\genfrac{}{}{0ex}{}{q-1+k}{k}\right)\times$$

$${\int }_{-{\eta }_{max}}^{{\eta }_{max}}d\eta {\int }_0^{R(\eta )}rdr{\int }_0^{2\pi }d\varphi \frac{m_{\mathrm{T}}\cosh (y-\eta )}{[1+(q-1)X{]}^{\frac{q+k}{q-1}}}.$$

$$\frac{d^{2}N}{p_{\mathrm{T}}d{p}_{\mathrm{T}}dy}=g\frac{3V_0}{8{\pi }^2}\sum _{k=0}^{\infty }{\xi }^k\left(\genfrac{}{}{0ex}{}{q-1+k}{k}\right){\int }_{-{\eta }_{max}}^{{\eta }_{max}}\frac{d\eta }{{\eta }_{max}}$$

$$\times {\int }_0^{R(\eta )}\frac{rdr}{R_0^2}\frac{m_{\mathrm{T}}\cosh (y-\eta )a^{\frac{q+k}{q-1}}}{[1+\frac{{\gamma }_r{m}_{\mathrm{T}}\cosh (y-\eta )-\mu }{T/(q-1)}{]}^{\frac{q+k}{q-1}}}{P}_{\frac{k+1}{q-1}}(a),$$

$$b=\frac{{\gamma }_r{v}_r{p}_{\mathrm{T}}}{T/(q-1)+{\gamma }_r{m}_{\mathrm{T}}\cosh (y-\eta )-\mu }$$

$$\underset{q\to 1}{\lim }{a}^{\frac{k+1}{q-1}}{P}_{\frac{k+1}{q-1}}(a)=I_0(\frac{k+1}{T}{\gamma }_r{v}_r{p}_{\mathrm{T}})$$

$${\sigma }_{in}=(26.2+0.1717{\ln }^2\frac{s}{3.521}+\frac{53.2}{s^{0.40}}-\frac{27.0}{s^{0.48}}-\lambda \frac{33.8}{s^{0.545}})\text{mbox}mb.$$

$$T=T_{\infty }(1+\frac{1.33\sqrt{x}-0.21}{1+x^2}),$$

$${\eta }_{max}$$

$$v_s$$

$$V_0=\tilde{V}{\eta }_{max}{y}_b(T_{\infty }/T{)}^{2.06},$$

$$\tilde{V}[1+\frac{E-\mu }{T/(q-1)}{]}^{\frac{q}{1-q}}={\tilde{V}}_{\pi }[1+\frac{E-{\mu }^{\prime }}{(T-T\delta )/(q-1)}{]}^{\frac{q}{1-q}},$$

$$n=q/(q-1),$$

$$n=\frac{n_{\infty }}{1+p_{1}x}+\frac{p_3/\ln x}{\ln x-p_2}+p_4{x}^{0.37}-p_{5}x,$$

$$\mu =p_6(\ln \frac{e_2}{\sqrt{s}}-p_7\frac{1-\sqrt{s}/e_2}{1+\sqrt{s}/e_3})m,$$

$${\mu }_{\pi }=\frac{2\mathrm{GeV}}{\sqrt{s}}(1+0.6Q_{\pi })e^{-\sqrt{s}/e_0}{m}_{\pi },$$

$$\tilde{V}={\tilde{V}}_{NP}(1+\frac{c_1}{1+(c_2/p_{\mathrm{T}}{)}^4}),$$

Full text

I Introduction

Theoretical description of the hadron transverse momentum ($p_{\rm T}$) and rapidity ($y$) spectra produced in proton-proton ($pp$), proton-nucleus ($pA$) and nucleus-nucleus ($AA$) collisions is one of the important tasks of high-energy physics. Since its realization in the QCD is still not fully satisfactory (e.g., due to the parton hadronization complicated processes, especially at low $p_{\rm T}$), alternative phenomenological methods are also in use. For instance, the thermal models of the stationary fireball (hadronic gas) with conventional Boltzmann-Gibbs distribution (BGD) are widely used to explain the hadronic abundances and $p_{\rm T}$-spectra at low $p_{\rm T}$ (see, e.g., BRS ; Becat ; Ths ; Shr ). At high $p_{\rm T}$ the exponential BGD is not adequate since the spectra have a power-law form. The thermal models with expanding (also called flowing) fireball, like the blast-wave model (BWM) BW1 ; BW2 , are included in hadron generators Therm ; Amel ; Drag . They assume the physics scenario that the initial collision creates a thermalized quark-gluon fireball, which expands, cools, hadronizes and goes through the chemical freeze-out and finally the kinetic freeze-out, when it decays into the free-streaming hadrons. Hadron spectra are computed usually by the Cooper-Frye formula CF and flow-boosted BGD. The longitudinal flow helps us to explain the $y$-spectra (see, e.g., BW1 ; BW2 ; BSW ), while the radial (or transverse) flow flattens the $p_{\rm T}$-spectra and improves the data description up to $p_{\rm T}$ values of several GeV/$c$ (see, e.g., BW1 ; BW2 ; Star0 ; Ryb ; Ghosh ; Chatt ; Alice1 ).

Recent years the thermal models employing the Tsallis distribution (TD) Tsallis1988 have become very popular, especially after the LHC operation Bed ; Beck ; Wilk1 ; Wilk2 ; Biro1 ; Tang ; Conroy ; Cley1 ; Cley2 ; Azmi ; Wlod1 ; Wong1 ; Wong2 ; Wong3 ; Wilk3 ; Wilk4 ; Urm1 ; Wilk5 ; Li ; Depp1 ; Depp2 ; Depp3 ; Zheng ; Cley3 ; Urm2 ; Star1 ; Phen1 ; Phen2 ; CMS1 ; Atlas1 ; Alice2 111There are hundreds of papers that develop and/or use such models. We cite only some of them which include further references.. Its ability to describe the charged hadron $p_{\rm T}$-spectra in a large $p_{\rm T}$ range 0–200 GeV/$c$  Wong1 ; Wong2 ; Wong3 is very impressive. TD is a generalization of the BGD. Besides the temperature $T$ and chemical potential $\mu$ it has an additional parameter $q$ and reduces to BGD in the limit $q\rightarrow 1$. TD can be considered as a result of averaging of the temperature fluctuations in the BGD, where $q-1$ characterizes the strength of these fluctuations Wilk1 (for other interpretations, see Urm1 ; Wilk5 ). The relation of TD with the QCD hard-scattering formulas is discussed in Wong2 ; Wong3 . Thanks to parameter $q$, TD provides a smooth transformation of the $p_{\rm T}$-spectrum shape from the nearly exponential form at low $p_{\rm T}$, similar to BGD, to the power-law form at high $p_{\rm T}$, which is the usual domain of the perturbative QCD. Most of the TD-based models consider a stationary fireball and are devoted to the fits of hadron $p_{\rm T}$-spectra in different collisions. Papers Tang ; Cley3 use a flowing fireball of the BWM and study the radial flow effect on the $p_{\rm T}$-spectra. In Li ; Depp3 , the $y$-spectra of charged particles are also considered in the two-fireball models with a longitudinal flow. All these studies give different values for $T$ and $q$. Parameter $q$ increases slowly with the collision energy $\sqrt{s}$ and varies in the range 1–1.2, depending on the hadron and collision types. Some theoretical arguments give the upper limit $q=11/9$ Beck .

In the present paper we propose a new thermal model based on the TD and BWM with a flowing fireball. It utilizes thermodynamically consistent version of the TD Conroy ; Cley1 and differs from similar models by a suitable choice of the BWM ingredients (see Sec. II), allowing us to describe the shape and normalization of the hadronic $p_{\rm T}$ and $y$ spectra, measured in $pp$ collisions at energies from $\sqrt{s}=5$ GeV to the highest LHC one of 13 TeV and in $p\overline{p}$ collisions at $\sqrt{s}>500$ GeV. In our model, unlike others which also use TD, the kinetic freeze-out temperature is the same for all hadron species. Here we consider only pions and quarkonia since the pion data are the most abundant (in terms of statistics and $\sqrt{s}$ values) and the quarkonia (${\rm J}/\psi$, $\rm\Upsilon(1S)$, …) data cover large intervals of $p_{\rm T}$ and $y$, which are important for our fits to better fix the model parameters. We provide simple parametrizations for the $\sqrt{s}$-dependence of the model parameters allowing us to predict the pion and quarkonia yields $d^{2}N/dp_{\rm T}dy$ in $pp$ collisions at new energies. Other particles as well as $pA$ and $AA$ collisions will be considered elsewhere.

The paper is organized as follows: Sec. II gives details of the model. In Sec. III, we discuss the model parameters and fit procedure. Sections IV and V are devoted to the description of pion and quarkonia data, respectively. In the last section, our concluding remarks are given.

II Model description

In thermal models the single-particle invariant yield is usually defined by the Cooper-Frye integral over the kinetic freeze-out space-time hypersurface $\Sigma_{f}$ CF

[TABLE]

where $X=(p_{\nu}u^{\nu}-\mu)/T$. Here the integrand $f$ is the freeze-out distribution of particle four-momentum $p=(E,\vec{p})$ and four-coordinate $x=(t,\vec{x})$ with temperature $T$ and $x$-dependent collective flow four-velocity $u$, $u_{\nu}u^{\nu}=1$, $\mu$ is the particle chemical potential and $g=2J+1$ is its spin degeneracy factor. Generally, $T$ and $\mu$ may also depend on $x$, but in order to keep our model as simple as possible, we assume that they are constant on the $\Sigma_{f}$. Then, the invariant volume $V_{f}$, which is called the fireball effective volume of particle production and includes the flow effects, could be factored out due to the Lorentz invariance in the expression for the particle total integrated yield AkkSin1 ; Soll ; Heinz ; Cley4 ; Bron ; AkkSin2

[TABLE]

We further assume, according to the BWM BW1 ; BW2 , that the fireball flow and geometry are azimuthally symmetric and boost invariant along the longitudinal ($z$) direction, as expected at high-energy $pp$ (also $p\overline{p}$ and central $AA$) collisions. Now, instead of the Cartesian coordinates, it is convenient to introduce the radial vector $\vec{r}=(r\cos\phi,r\sin\phi)$ and the Bjorken longitudinal proper time $\tau=\sqrt{t^{2}-z^{2}}$ and space-time rapidity $\eta=\frac{1}{2}\ln\frac{t+z}{t-z}$. Then, the flow four-velocity could be written as BW1

[TABLE]

where $\gamma_{r}=1/\sqrt{1-v_{r}^{2}}$ and $v_{r}$ is the radial flow velocity. Expressing the particle four-momentum via the $y$ and $p_{\rm T}$ ($m_{\rm T}=\sqrt{m^{2}+p_{\rm T}^{2}}$ is transverse mass), $p^{\nu}=(m_{\rm T}\cosh y,~{}p_{\rm T}\cos\phi_{p},~{}p_{\rm T}\sin\phi_{p},~{}m_{\rm T}\sinh y)$, we get

[TABLE]

The hypersurface $\Sigma_{f}$ in the BWM is defined by the condition that the freeze-out happens at a constant value of the proper time: $\tau=\tau_{f}=const$. In this case the hypersurface element four-vector has a simple form BW2 ; Amel

[TABLE]

We fix the $\Sigma_{f}$ geometry as follows: in the longitudinal direction, it is limited in the interval $-\eta_{max}<\eta<\eta_{max}$, where a maximum longitudinal flow rapidity $\eta_{max}$ is required by the finite total energy (this breaks the exact longitudinal boost invariance). In the radial direction the upper boundary of $r$ is given by radius $R(\eta)$ that depends on $\eta$. This dependency plays a major role in our model for the proper description of the hadron rapidity spectra. We have tried different forms for it and found that the following simple one (see, e.g., BW2 )

[TABLE]

is very successful. Since $R_{0}$ is the radius at $\eta=0$, the fireball gets thinner with increase of $|\eta|$.

Now, we need to define the radial flow velocity $v_{r}$. Usually one assumes that it equals zero at $r=0$ and grows with $r$ according to a power-law dependence BW1 . We have found that the simple quadratic dependency

[TABLE]

allows us to correctly describe the hadron $p_{\rm T}$ spectra. Here $v_{s}=v_{r}(R_{0})$ is the surface velocity. A useful quantity is the mean value of $v_{r}(r)$, which can be defined as

[TABLE]

According to the Eqs. (2)$-$(7) one has

[TABLE]

where $V_{0}=4/3\,\pi R_{0}^{2}\tau_{f}\eta_{max}$. Performing similar calculations with Eq. (8) we obtain

[TABLE]

Fig. 1 shows that the ratios $V_{f}/V_{0}$ and $\langle v_{r}\rangle/(0.4v_{s})$ are equal unity at $v_{s}=0$ and grow with the $v_{s}$.

Thus, we defined the BWM ingredients of our model. Now we specify the function $f$ in Eq. (1) by choosing the thermodynamically consistent TD Conroy ; Cley1 (in contrast to the TD version, defined by Eq. (11) with the external power index $-1$ instead of $-q$)

[TABLE]

where $\xi$ equals 1 or $-1$ to account for the quantum statistics of bosons or fermions, respectively. This quantum correction matters only for pions due to their small mass. Expanding the right-hand side of Eq. (11) into the binomial series and substituting it in Eq. (1), one gets

[TABLE]

Using Eq. (4) and performing integrations over $\phi$ and $\phi_{p}$ (second integration gives $2\pi$), we obtain

[TABLE]

where $a=1/\sqrt{1-b^{2}}$  ,

[TABLE]

and $P_{\nu}(a)$ is the Legendre function of the first kind Erdelyi1 . Taking into account the relation ($I_{0}(x)$ is the modified Bessel function)

[TABLE]

one can easily verify that in the limit $q\rightarrow 1$ Eq. (13) reproduces usual BWM formulas BW2 based on the BGD. Eq. (13) (with Eqs. (6) and (7)) is the main formula of our model. We have checked that the series in this formula is convergent if $\mu<m$ (like for similar series in the thermal models with BGD Shr ). This condition is fulfilled according to Eqs. (21) and (22). Higher terms of the series are important only for pions (mostly for $\pi^{+}$ which has larger $\mu$) at low values of $\sqrt{s},\,|y|$ and $p_{\rm T}$. For example, for the case of $\sqrt{s}=30.6$ GeV, $y=0$ and $p_{\rm T}\sim 0$, considered in Sec. IV, first three terms of the series give together about 97% of the $\pi^{+}$ yield. At lower energies, more terms of the series should be used for the accurate computation of pion yields. For heavier hadrons, one can safely use $\xi=0$.

III Parameters and fit procedure

Here, we utilize Eq. (13) for fitting the hadron spectra in $pp$ and $p\overline{p}$ collisions. We follow two aims. First is to show that our model with a possibly minimum number of parameters is able to describe well the available data on $p_{\rm T}$ and $y$ spectra for different particles and energies $\sqrt{s}$. The second aim is to systematize the fit results for different $\sqrt{s}$ and provide simple parametrizations for the $\sqrt{s}$-dependence of model parameters, permitting predictions for the future experiments.

To fit the data given in terms of the cross section $\sigma$, we convert it to the invariant yield $N$ via the relation $\sigma=N\sigma_{in}$, where $\sigma_{in}$ is the $pp$ or $p\overline{p}$ inelastic cross section at the energy $\sqrt{s}$/GeV ($\lambda$ equals 1 or $-$1, respectively, see the L2 model of Table B1 in sigInel )

[TABLE]

As in other applications of TD for inclusive pions, we do not calculate explicitly the feed-down contribution from the resonance decays, assuming that directly produced pions and secondary ones have the same spectral shapes. Secondary pions are expected to dominate at low $p_{\rm T}$ (see, e.g., Star0 ).

Eq. (13) has six independent parameters: $T$, $q$, $\mu$, $V_{0}$, $\eta_{max}$ and $v_{s}$. Generally, they can depend on the $\sqrt{s}$ and hadron type. We assume that the kinetic freeze-out temperature T is the same for all hadron species (the chemical freeze-out temperature may rise with the hadron mass). Since the neutral pions and quarkonia ($\phi$, ${\rm J}/\psi$, $\psi({\rm 2S})$, $\rm\Upsilon$ family) do not have conserved quantum numbers, their chemical potential $\mu$ must equal zero in the chemical equilibrium BRS . We have verified that at $\sqrt{s}>50$ GeV the neutral and charged pion data can be successfully fitted with $\mu=0$, while this is not true for heavier hadrons. The nonzero $\mu$ can be interpreted as a measure of the non-equilibrium for the given particle. A similar fact is well known in the non-equilibrium thermal models based on BGD, where one introduces so-called phase-space occupancy $\gamma$, related to the chemical potential as $\mu=T\ln\gamma$ Shr . To ensure the same yield for the pion three charge states at high energies, as follows from the data, we assume that all the model parameters, except $\mu$, are the same for these states. Moreover, we will use for them a common averaged mass $m_{\pi}=(2m_{\pi^{\pm}}+m_{\pi^{0}})/3$.

Using the above-mentioned assumptions, we have done $\chi^{2}$ fits (in the ROOT framework ROOT ) of the existing data on pion and quarkonia $p_{\rm T}$ spectra for different values of $y$ and $\sqrt{s}$. We started with the pion fits and have observed that parameter $T$ increases with energy at low energies up to about $\sqrt{s}=10$ GeV. Then, it decreases and becomes practically constant at $\sqrt{s}>500$ GeV. This behavior can be parametrized as (see Fig. 2)

[TABLE]

where $x=\sqrt{s}/(16$ GeV) and $T_{\infty}=78$ MeV is the temperature at $\sqrt{s}\to\infty$. Similar energy dependence was observed for the kinetic freeze-out temperature in $AA$ collisions using thermal models with the BGD (see, e.g., Fig. 11a in Chatt ).

We then utilized Eq. (15) in the fits of all hadrons. The fit results for $\eta_{max}$ and $v_{s}$ are parametrized as

[TABLE]

where $e_{0}=45$ GeV, $m$ is the mass of given hadron, $y_{m}=\ln(\sqrt{s}/m)$ its maximum rapidity, $m_{p}$ the proton mass and $y_{b}=\ln(\sqrt{s}/m_{p})$ the beam rapidity in high energy $pp$ or $p\overline{p}$ collisions. In our model, the $y$-spectrum width is proportional to $\eta_{max}$ and grows logarithmically with $\sqrt{s}$. Besides, the larger $m$ is, the smaller $\eta_{max}$ is and hence the narrower the $y$-spectrum is. Parameter $v_{s}$ changes in the range 0–0.78, increases with $\sqrt{s}$ and decreases with increasing $m$. Eqs. (10) and (17) show that radial flow velocity for pions is significant even at $\sqrt{s}\sim 5$ GeV (while it vanishes at $\sqrt{s}\leq 540$ GeV in Tang ). Note that the last terms in Eqs. (16) and (17) are important only for pions at low energies.

The remaining fit parameters also demonstrate properties common for different hadrons. The volume parameter can be expressed as

[TABLE]

where $\widetilde{V}$ is $\sqrt{s}$-independent but strongly decreases with increase of the hadron mass (see Table 1)222 In principle, it is possible to redefine the Eq. (13) parameters and obtain for heavier hadrons the same $\widetilde{V}$ as for pions using the following identity transformation of the TD (see also Cley2 ):

$\widetilde{V}[1+\frac{E-\mu}{T/(q-1)}]^{\frac{q}{1-q}}=\widetilde{V}_{\pi}[1+\frac{E-\mu^{\prime}}{(T-T\delta)/(q-1)}]^{\frac{q}{1-q}}\,,$

where $\mu^{\prime}=\mu-T\delta/(q-1)$ and parameter $\delta=1-(\widetilde{V}/\widetilde{V}_{\pi})^{1-1/q}$ grows with $\widetilde{V}_{\pi}/\widetilde{V}$. As mentioned in Sec. I, the $q-1$ characterizes the temperature fluctuations around the mean value $T$. According to Wilk2 , the quantity $T\delta$ can be interpreted as a measure of the energy transfer, caused by these fluctuations, from the fireball region where the particle is produced to the surrounding regions. Note that $\delta\sim(q-1)$ at $q\to 1$, as expected in Wilk2 . . Note that $V_{0}\sim{\ln}^{2}{s}$ at high $\sqrt{s}$. The normalization constant $\widetilde{V}$ for inclusive pions, given in Table 1, includes the contribution of the resonance decays and hence is expected to be larger than the one for the directly produced pions.

The fitted values of $q$ grow with $\sqrt{s}$ and are different for different hadron species. However they vary in the very small interval 1–11/9, as noted in Sec. I. Therefore, it is more convenient to use the parameter $n$ instead

[TABLE]

which controls the large-$p_{\rm T}$ behavior of Eq. (13). Then, $n>n_{\infty}=$11/2 Beck and $n\to\infty$ at $q\to 1$. the resulting fitted values of $n$ can be parametrized for different hadron species by the formula, valid at $x=e_{1}/\sqrt{s}<1$,

[TABLE]

where parameters $e_{1}$ and $p_{1}$–$p_{5}$ are listed in Table 1 and $e_{1}$

is the energy when $n$ becomes infinity. So, at $\sqrt{s}\leq e_{1}$, the TD reduces to BGD. The limiting value $n_{\infty}$ provides that Eq. (13) at $\sqrt{s}\to\infty$ has the same large-$p_{\rm T}$ behavior as the jet production in the lowest-order perturbative QCD Wong3 . Fig. 3 shows the energy dependence of $n$ for $\pi$, ${\rm J}/\psi$ and $\rm\Upsilon(1S)$. The corresponding curves for $\phi$, $\psi({\rm 2S})$ and higher $\rm\Upsilon$ states are similar to the ones for $\pi$, ${\rm J}/\psi$ and $\rm\Upsilon(1S)$, respectively.

The obtained values of $\mu$ are always smaller $m$. They are proportional to $m$ and vanish with increasing $\sqrt{s}$. We parametrize $\mu$ for different quarkonium species by the formula, valid at $\sqrt{s}\leq e_{2}$,

[TABLE]

where parameters $e_{2}$, $e_{3}$, $p_{6}$ and $p_{7}$ are given in Table 1 and $\mu=0$ at $\sqrt{s}>e_{2}$ ($e_{2}$ is larger for heavier particles). For pions we use

[TABLE]

where $Q_{\pi}=0$ for $p\overline{p}$ collisions while for $pp$ collisions $Q_{\pi}$ equals the pion charge, to account for the difference of $\pi^{+}$, $\pi^{0}$ and $\pi^{-}$ yields in low energy $pp$ collisions, related to the charge-conservation effects. $\mu_{\pi}$ vanishes with increasing energy, in agreement with the fact that these yields almost coincide at $\sqrt{s}\geq 62.4$ GeV Phen1 ; Phen2 .

In Sections  IV and V we will discuss in more detail the results of combined fits of pion and quarkonia data using Eqs. (13)$-$(22). The parameter values as well as the $\chi^{2}$ and $NDF$ of the fits for each hadron type are given in Table 1. Additional parameters for ${\rm J}/\psi$ and $\psi({\rm 2S})$, produced via bottom hadron decays, and for higher $\rm\Upsilon$ states will be considered in Sec. V. Note that rather large ratios $\chi^{2}/NDF$ are due to the large amount of data included in the fits, which use $\sqrt{s}$-dependent parametrizations for the model parameters. Since the quality and normalization of different measurements for given hadron do not always agree well with each other, the combined fit gives larger $\chi^{2}/NDF$ than the individual fits for each measurement. To get not-too-large $\chi^{2}/NDF$, we have excluded some data samples from the combined fits.

IV neutral and charged pions

Here, we present the results of the combined fit of

$\pi^{0}$ Phen1 ; Kourk ; Phen3 ; Phen4 ; Star2 ; Star3 ; Phen5 ; UA2_1 ; Alice3 ; Alice4 and $\pi^{\pm}$ Phen2 ; Alice2 ; BHM ; Beier76 ; NA61 ; Alper ; Star4 ; UA2_2 ; Alice5 ; Alice6 inclusive production $p_{\rm T}$-spectra measured for different values of $y$ at energies $\sqrt{s}$ from 30.6 GeV Kourk to 7 TeV Alice4 for $\pi^{0}$ and from 4.93 GeV BHM to 7 TeV Alice2 for $\pi^{\pm}$. The used charged pion data are mostly for $\pi^{-}$. High-energy data at $\sqrt{s}=$ 0.54, 2.76 and 7 TeV are for the averaged ($\pi^{-}+\pi^{+}$)/2 production. From Kourk , we included in the fit only the data obtained with so-called retracted geometry, and from Alper , we included only the data measured at $\sqrt{s}=$ 30.6, 44.6 and 52.8 GeV which cover larger intervals of $p_{\rm T}$ and $y$. We did not include in the fit the $\pi^{+}$ $p_{\rm T}$-spectra from BHM , which give too large $\chi^{2}$; however, our model describes well the corresponding $p_{\rm T}$-integrated data (see Fig. 6). Since the charged pion measurement at $\sqrt{s}=200$ GeV by the STAR Collaboration Star4 is for the non-single diffractive (NSD) yield, it was converted to an inclusive cross section using $\sigma_{NSD}=30$ mb Star4 . Figures 4 and 5 (left) show examples of the fits of pion $p_{\rm T}$-spectra for mid-rapidity and different $\sqrt{s}$ values while Fig. 5 (right) shows fits for different values of rapidity at $\sqrt{s}=30.6$ GeV Alper . To demonstrate the quality of the fits, the data points have been divided by the corresponding values of the fit function, and the ratios are also plotted. Generally, the quality is always good. Only the $\pi^{0}$ data Star2 ; Phen5 show a large excess at $p_{\rm T}>10$ GeV/$c$. In Fig. 5 (left), the dashed lines represent the fit functions at $\sqrt{s}=$ 6.3 and 900 GeV for $v_{s}=0$ to illustrate the importance of the radial flow in our model.

Note that we did not include in the combined fit the $\pi^{\pm}$ high-statistics data NA49 measured in $pp$ collisions at $\sqrt{s}=17.3$ GeV since the published results do not quote the dominant systematic uncertainties. But we included the $\pi^{-}$ measurement at the same energy NA61 , and both data sets agree well, as shown in NA61 . We also did not use for the combined fit the charged pion data measured at $\sqrt{s}=200$ GeV and high rapidities $y=$ 2.95 and 3.3 Brahms1 but have checked that our model describes them well at $p_{\rm T}<2.5$ GeV/$c$. At higher $p_{\rm T}$, corresponding to pion energies larger than 25 GeV, the model overestimates the data. So, our model is not valid for such high rapidities at $\sqrt{s}=200$ GeV, related to diffractive processes, which is generally expected for a thermal model. Note also that the $\pi^{\pm}$ data Alper provide a large contribution into the $\chi^{2}$ and $NDF$ values shown in Table 1. A combined fit without these data gives $\chi^{2}/NDF=5959/1734$.

Fig. 6 presents an example of our predictions, based on Eqs. (13)$-$(22), for the pion $p_{\rm T}$-integrated yields in $pp$ collisions at mid-rapidity and varying $\sqrt{s}$. It shows a good agreement with the available data.

V Quarkonia ($\phi$, ${\rm J}/\psi$, $\psi({\rm 2S})$, $\rm\Upsilon(1S)$–$\rm\Upsilon(3S)$)

V.1 $\phi$ meson

The following results are for the combined fit of $\phi$ meson inclusive production data measured in $pp$ collisions at $\sqrt{s}=$ 17.3 GeV NA49phi , 200 GeV Phen1 ; Phen6 , 900 GeV Alice09phi , 2.76 TeV Alice3phi , 7 TeV LHCb7phi ; Alice7phiee ; Alice7phimm and in $p\overline{p}$ collisions at $\sqrt{s}=$1.96 TeV CDFphi . The $p_{\rm T}$-spectrum from E735phi is not included in the fit since its normalization is about six times lower the one in CDFphi at similar energy. It appeared that the fitted values of the $\phi$ meson $n$ parameter at different $\sqrt{s}$ are close to the pion ones. So, in the parametrization Eq. (20) for $\phi$, we have fixed some of the parameter values to the ones for pion (see Table 1). Examples of the $p_{\rm T}$-spectra fits are shown in Fig. 7 for mid-rapidity and different $\sqrt{s}$ values (left) and for different values of $y$ at $\sqrt{s}=7$ TeV (right). As an example of our predictions, Fig. 8 presents the $\phi$ meson $p_{\rm T}$-integrated cross section in $pp$ (also $p\overline{p}$) collisions versus $\sqrt{s}$ at mid-rapidity and

at forward rapidity of $y=3.25$ (for dimuon decay channel measurements of LHCb LHCb7phi and ALICE Alice3phi ; Alice7phimm ). At $y=3.25$ two values for the $p_{\rm T}$-integration lower limit are considered: 0 and 1 GeV/$c$. Comparison of calculations with the available data shows a reasonable agreement.

V.2 ${\rm J}/\psi$ meson

The inclusive ${\rm J}/\psi$ production consists of prompt component (includes direct production and feed-down from the radiative decays of higher charmonium states) and non-prompt component (includes feed-down from the weak decays of bottom hadrons). Fraction of the non-prompt component, denoted usually by $f_{B}$, is negligible at $\sqrt{s}<$ 100 GeV but rises monotonically with $\sqrt{s}$ and $p_{\rm T}$. For LHC energies it reaches values about 0.1 at low $p_{\rm T}$ and larger values at high $p_{\rm T}$ (see Fig. 11). The kinematic distributions of prompt and direct ${\rm J}/\psi$ are similar and can be described by the same values of parameters in Eq. (13). Only the normalization constants $\widetilde{V}$ will differ. The non-prompt ${\rm J}/\psi$ has a significantly harder $p_{\rm T}$ spectrum and narrower $y$ spectrum. Its proper description would require the use of Eq. (13) for the production of bottom mesons and baryons which have several decay channels into ${\rm J}/\psi$. To avoid such a complex computation for a rather small fraction of data, we have chosen a simpler approach. Namely, for non-prompt ${\rm J}/\psi$, we use Eq. (13) with the same $T$, $n$, $\mu$ and $v_{s}$ parameters as for prompt ${\rm J}/\psi$. To describe the harder $p_{\rm T}$-spectrum of non-prompt ${\rm J}/\psi$, we assume that the mass in $m_{\rm T}$ in the

corresponding Eq. (13) is larger the ${\rm J}/\psi$ mass by some factor $c_{m}$ and the normalization grows with $p_{\rm T}$ according to the parametrization

[TABLE]

where $\widetilde{V}_{NP}$, $c_{1}$, $c_{2}$ and $c_{m}$ are fit parameters. Also, to ensure the narrowness of the non-prompt $y$-spectrum, we multiply the corresponding $\eta_{max}$ in Eq. (17) by another fit parameter $c_{\eta}<1$. We have performed a combined fit of the available prompt and non-prompt or inclusive ${\rm J}/\psi$ production $p_{\rm T}$-spectra NA3 ; Phen7 ; Star5 ; CDFpsi1 ; D0jpsi ; CDFjpsi1 ; CDFjpsi2 ; CMS3QQ ; Alice3jpsi ; LHCb3jpsi ; Alice513psi ; CMS7jpsi ; Atlas7jpsi ; LHCb7jpsi ; Alice7jpsi ; CMS7psi1 ; Alice7QQ ; CMS7psi2 ; LHCb8QQ ; Alice8QQ ; Atlas78psi ; LHCb13jpsi measured at energies $\sqrt{s}$ from 19.4 GeV NA3 to 13 TeV Alice513psi ; LHCb13jpsi in $pp$ collisions and at 1.8 TeV CDFpsi1 ; D0jpsi ; CDFjpsi1 and 1.96 TeV CDFjpsi2 in $p\overline{p}$ collisions. The $p_{\rm T}$-spectrum from UA1jpsi is not included in the fit since its normalization is a factor of 2.5 lower than expected within our model, which however describes well the shape of this spectrum. The fit gives, in addition to the values for the model parameters, $\chi^{2}$ and $NDF$ listed in Table 1, the following values for the non-prompt component parameters: $\widetilde{V}_{NP}=82.1\pm 0.7~{}{\rm GeV}^{-3},\,c_{1}=2.1,\,c_{2}=26.3~{}{\rm GeV},\,c_{m}=1.4,\,c_{\eta}=0.82$.

We illustrate then some results of the fit. Inclusive ${\rm J}/\psi$ $p_{\rm T}$-spectra at different $\sqrt{s}$ values are shown in Fig. 9 for midrapidity (left) and forward rapidity of $y=3.25$

(right). Data points of ATLAS Atlas78psi in Fig. 9 correspond to $0.25<|y|<0.5$. Figure 10 shows the prompt ${\rm J}/\psi$ $p_{\rm T}$ spectra for different rapidity values at $\sqrt{s}=7$ TeV. In Fig. 11, our predictions for the $p_{\rm T}$ dependence of the ${\rm J}/\psi$ non-prompt fraction at different $\sqrt{s}$ are compared with existing published data. We demonstrate also a good agreement of our predictions with the available data on the inclusive ${\rm J}/\psi$ meson $p_{\rm T}$-integrated cross section (for $p_{\rm T}>0$) as a function of $y$ at three $\sqrt{s}$ values (Fig. 12) and as a function of $\sqrt{s}$ at mid-rapidity and forward rapidity (Fig. 13).

V.3 $\psi({\rm 2S})$ meson

As $\psi({\rm 2S})$ is a charmonium state, similar to ${\rm J}/\psi$, its production prompt and non-prompt components can be described similarly using Eqs. (13) and (23) (where $m$ now is the $\psi({\rm 2S})$ mass) and parameters $c_{m},\,c_{\eta}$. We have performed a combined fit of the available prompt and non-prompt or inclusive $\psi({\rm 2S})$ production data measured in $pp$ collisions at $\sqrt{s}=$ 200 GeV Phen7 , 7–13 TeV Alice513psi ; CMS7psi1 ; Alice7QQ ; CMS7psi2 ; Alice8QQ ; Atlas78psi ; LHCb7psi ; Atlas7psi and in $p\overline{p}$ collisions at $\sqrt{s}=$ 1.8 TeV CDFpsi1 , 1.96 TeV CDFpsi2 . The resulting fit parameter values are listed in Table 1 (some of them are fixed to the corresponding values of ${\rm J}/\psi$). Additional parameters for the non-prompt $\psi({\rm 2S})$ are: $\widetilde{V}_{NP}=15.9\pm 0.1~{}{\rm GeV}^{-3},\,c_{m}=1.3$,  and $c_{1},\,c_{2},\,c_{\eta}$ coincide with the ones of ${\rm J}/\psi$. Examples of the $p_{\rm T}$-spectra fits are shown in Fig. 14 for mid-rapidity and different $\sqrt{s}$ values (left) and for different values of $y$ at $\sqrt{s}=7$ TeV (right). Data points of ATLAS Atlas78psi in the left panel correspond to

$0.25<|y|<0.5$. As an example of our predictions, Fig. 15 presents the inclusive $\psi({\rm 2S})$ meson $p_{\rm T}$-integrated (for $p_{\rm T}>0$) cross section in $pp$ (also high energy $p\overline{p}$) collisions versus $\sqrt{s}$ at mid-rapidity and at forward rapidity of $y=3.25$. Comparison of calculations with the available data shows a reasonable agreement.

V.4 $\rm\Upsilon(1S)$, $\rm\Upsilon(2S)$, $\rm\Upsilon(3S)$ mesons

Here, we discuss the fits of $\rm\Upsilon$-family mesons inclusive production data CMS3QQ ; Alice7QQ ; LHCb8QQ ; Alice8QQ ; E866upsi ; Phen8 ; Star6 ; Phen9 ; CDFupsi ; D0upsi ; LHCb3upsi ; LHCb7upsi ; Atlas7upsi ; CMS7upsi1 ; CMS7upsi2 ; LHCb78upsi measured at energies $\sqrt{s}$ from 38.8 GeV E866upsi to 8 TeV LHCb8QQ ; Alice8QQ ; LHCb78upsi in $pp$ collisions and at 1.8 TeV CDFupsi and 1.96 TeV D0upsi in $p\overline{p}$ collisions. First, a combined fit of more copious $\rm\Upsilon(1S)$ data was done, and the resulting parameter values are given in Table 1. Then, separate combined fits were performed for $\rm\Upsilon(2S)$ and $\rm\Upsilon(3S)$ data with free parameters $\widetilde{V}$ and $p_{3}$–$p_{6}$, fixing all other parameters to the corresponding values of the $\rm\Upsilon(1S)$ fit. The results for $\rm\Upsilon(2S)$ ($\rm\Upsilon(3S)$) are $\widetilde{V}=0.093\pm 0.001\,(0.034\pm 0.001)~{}{\rm GeV}^{-3},\;p_{3}=0.6\,(0.4),\;p_{4}=7.5\,(5.8),\;p_{5}=-8.0\,(-8.0),\;p_{6}=0.0548\,(0.0531),\;\chi^{2}/NDF=3443/491\,(3005/478)$. Note that such large ratios $\chi^{2}/NDF$ for $\rm\Upsilon$ family are mostly due to the somewhat poor match between results of different LHC experiments. Moreover, two different measurements of the LHCb Collaboration LHCb7upsi ; LHCb78upsi at

$\sqrt{s}=7$ TeV do not agree well, and the data at $\sqrt{s}=2.76$ TeV LHCb3upsi seem too high with respect to the model predictions (see Figs. 16 and 17).

Examples of inclusive $\rm\Upsilon(1S)$ meson $p_{\rm T}$-spectra fits at different $\sqrt{s}$ values are shown in Fig. 16 for mid-rapidity (left) and forward rapidity of $y=3.25$ (right). Fig. 17 presents the $\rm\Upsilon(1S)$ $p_{\rm T}$-integrated cross section (for $p_{\rm T}>0$) as a function of $y$ and comparison with the existing measurements at $\sqrt{s}=$ 2.76 and 7 TeV. Our prediction for $\sqrt{s}=13$ TeV is also given. The predictions for this cross section dependence on $\sqrt{s}$ at mid-rapidity and at forward rapidity of $y=3.25$, together with the available data, are shown in Fig. 18.

To illustrate the fits for $\rm\Upsilon(2S)$ and $\rm\Upsilon(3S)$ mesons, we consider the ratios of their inclusive production cross sections times their dimuon branching fractions to the same quantity for $\rm\Upsilon(1S)$, denoted usually as $R^{21}$ and $R^{31}$, respectively. Figure 19 presents the fit results for the $p_{\rm T}$ dependence of these ratios at different values of $\sqrt{s}$ and $y$. Lastly, Fig. 20 demonstrates the model description of the $\rm\Upsilon(1S+2S+3S)$ production $p_{\rm T}$-integrated cross section times the dimuon branching fraction as a function of $y$, measured in $pp$ collisions at $\sqrt{s}=200$ GeV Phen8 ; Star6 ; Phen9 .

VI Conclusion

Thus, we presented a thermal model of a flowing hadronic fireball, based on the TD and BWM, which describes well almost all available data (except diffractive processes at large $y$ values and some other data sets) on the pion and quarkonia production yields $d^{2}N/dp_{\rm T}dy$ in $pp$ collisions at $\sqrt{s}\geq 5$ GeV and in $p\overline{p}$ collisions at $\sqrt{s}>500$ GeV (where the difference between these two collision types can be neglected). Note that the longitudinal boost invariance is broken in the model due to the used fireball geometry.

One of the distinct features of our model is the assumption that the kinetic freeze-out temperature $T$ is the same for all particle types (while their chemical freeze-out temperatures can differ). $T$ is almost constant at $\sqrt{s}>500$ GeV, increases with decreasing energy and reaches its maximum at $\sqrt{s}\sim 10$ GeV (see Fig. 2). In this energy region, the parameter $n$ goes to infinity (see Fig. 3) and TD reduces to BGD. Another feature of the model is that the particle chemical potential $\mu$ is proportional to its mass and vanishes with increase of $\sqrt{s}$. The nonzero $\mu$ can be interpreted as a measure of the chemical non-equilibrium. Also, we provide parametrizations for the $\sqrt{s}$ dependence of the model parameters allowing predictions for the pion and quarkonia yields in $pp$ or $p\overline{p}$ collisions at new energies of the existing and future accelerators. An example script is given in my1 , showing how one can use our model to compute these yields at any $\sqrt{s}$, $p_{\rm T}$ and $y$ in ROOT ROOT .

In our model, the correlation between the parameters $T$ and $q$ (or $n$) and radial flow velocity ($v_{s}$) has similar behavior as in other models (see, e.g., BW1 ; Tang ). Namely, $T$ and $q$ increase with decreasing $v_{s}$. A combined fit of the pion data with $v_{s}=0$ gives about 10% larger $T$ and from 10% to 70% larger $q-1$ when moving from the LHC energies down to $\sqrt{s}\sim 20$ GeV. $\chi^{2}$ of this fit is about 50% larger than the one given in Table 1 for pions. Quarkonia fits also give similar parameter changes. It can be seen in Eq. (13) that the effect of the radial flow diminishes with increasing rapidity. Owing to this feature, our model describes the experimental fact (see, e.g., BW2 ; LHCb7jpsi ; Brahms2 ) that the $p_{\rm T}$ spectrum of a given particle becomes softer ($\langle p_{\rm T}\rangle$ becomes smaller) with the increase of its rapidity.

Finally, since the model includes all the ingredients of the thermal source (fireball), it can be applied for the pion Bose-Einstein correlation studies using, e.g., the methods of Bialas14 .

Acknowledgements

I thank P. Dupieux, C. Hadjidakis, A. Parvan, O. Teryaev, M. Tokarev and G. Wilk for interest and helpful discussions. I would also like to thank the anonymous referee for important comments and the suggestion to discuss the correlation between the model parameters.