Now, the phenomenological single-particle Tsallis distribution [1, 2] is used for the description of the experimental data on the transverse momentum distributions of hadrons created in the proton-proton and heavy-ion collisions at RHIC and LHC energies [3,4,5,6,7,8,9,10,11,12,13,14]. The phenomenological Tsallis transverse momentum distribution for the Maxwell–Boltzmann statistics of particles introduced in Refs. [1, 2] has the form

$$\begin{aligned} \frac{d^{2}N}{dp_{T}dy}= & {} \frac{gV}{(2\pi )^{2}} p_{T} m_{T} \cosh y \nonumber \\&\times \left[ 1- (1-q) \frac{m_{T} \cosh y-\mu }{T} \right] ^{\frac{q}{1-q}}, \end{aligned}$$
(1)

where \(m_{T}\) is the transverse mass of particles. Initially, the phenomenological Tsallis transverse momentum distribution was proposed in another format (see Refs. [15, 16])

$$\begin{aligned} \frac{1}{\sigma }\frac{d\sigma }{dp_{T}}= & {} c p_{T} \int \limits _{0}^{\infty } dp_{L} \left[ 1- (1-q) \frac{\sqrt{p_{L}^{2}+m_{T}^{2}}}{T} \right] ^{\frac{q}{1-q}}, \end{aligned}$$
(2)

where \(p_{L}\) is the longitudinal momentum. However, further we will discuss only the function (1) as it is widely used in high-energy physics.

In high-energy physics, the phenomenological transverse momentum distribution (1) is associated with the Tsallis statistics based on the Tsallis entropy [17]

$$\begin{aligned} S = \sum \limits _{i} \frac{p_{i}^{q}-p_{i}}{1-q}, \quad \sum \limits _{i} p_{i}=1, \end{aligned}$$
(3)

where \(p_{i}\) is the probability of the i-th microscopic state of the system and \(q\in {\mathbf {R}}\) is a real parameter taking values \(0<q<\infty \). Note that there is also another supplementary definition of the statistical entropy named as the escort entropy, which is written in the form as [18, 19]

$$\begin{aligned} S_{E} = \frac{1-\left( \sum _{i} p_{i}^{1/q}\right) ^{-q}}{q-1}. \end{aligned}$$
(4)

However, in the present study, we do not use it, because the escort entropy (4) is not suitable for the search of the origin of the transverse momentum distribution (1).

The Tsallis generalized statistical mechanics is constructed on the base of the statistical entropy (3), different definitions of the expectation values of thermodynamic quantities and the second law of thermodynamics [17]. There is at least three classification schemes of the Tsallis statistical mechanics depending on the choice of the definition of the expectation values of the thermodynamic quantities [18]. The Tsallis-1 statistics deals with the standard linear definition of the expectation values, \(\langle A \rangle = \sum _{i} p_{i} A_{i}\). The Tsallis-2 statistics is defined on the base of the generalized expectation values of the thermodynamic quantities, \(\langle A \rangle = \sum _{i} p_{i}^{q} A_{i}\), and the Tsallis-3 statistics is associated with the expectation values with the escort probabilities as \(\langle A \rangle = \sum _{i} p_{i}^{q} A_{i}/\sum _{i} p_{i}^{q}\).

In spite of the success of the function (1) in the applications its connection with the basis of the statistical mechanics is still under the question. The phenomenological single-particle Tsallis distribution (1) was obtained using the method of the maximization of generalized entropy of the ideal gas instead of the general Tsallis entropy (3). This method is correct only for the usual Boltzmann-Gibbs statistics, but for the Tsallis statistics it is under the question [20]. Nevertheless, the correct connection of the distribution function (1) with the Tsallis entropy (3) was established in Ref. [20]. In that paper, it was rigorously demonstrated that the phenomenological transverse momentum distribution (1) corresponds to the zeroth term approximation of the Tsallis-2 statistics. In the case of massive particles, this result was confirmed in Ref. [21]. However, the Tsallis-2 statistics is improperly defined as the generalized expectation values of the thermodynamic quantities, \(\langle A \rangle = \sum _{i} p_{i}^{q} A_{i}\), in this formalism are not consistent with probability normalization condition, \(\sum _{i} p_{i}=1\) (see Ref. [18]). This means that the Tsallis-2 statistics disagrees with the probability theory and the second law of thermodynamics. Thus considering the phenomenological transverse momentum distribution (1) as a single-particle distribution function of the Tsallis statistics based on the Tsallis entropy (3) is a serious lack, because in this case the distribution (1) is obtained on the base of the erroneous general premises.

In the present study we demonstrate that this problem of connection of the phenomenological transverse momentum distribution (1) with the basis of the statistical mechanics can be uniquely solved by introducing the new nonextensive statistics on the base of the q-dual entropy instead of the Tsallis entropy (3).

q-dual statistics The q-dual entropy is obtained from the Tsallis entropy (3) by the multiplicative transformation \(q\rightarrow 1/q\). Then, we have

$$\begin{aligned} S = q\sum \limits _{i} \frac{p_{i}^{1/q}-p_{i}}{q-1}, \end{aligned}$$
(5)

where \(q\in {\mathbf {R}}\) is a real parameter taking values \(0<q<\infty \). In the Gibbs limit \(q\rightarrow 1\), the entropy (5) recovers the Boltzmann-Gibbs-Shannon entropy, \(S=-\sum _{i} p_{i} \ln p_{i}\). The properties of the q-dual entropy are the same as those of the Tsallis entropy (3) but at different values of the parameter q. The values of the Tsallis entropy for \(q<1\) correspond to the values of entropy (5) for \(q>1\) and, inversely, the values of entropy (3) for \(q>1\) correspond to the values of the q-dual entropy for \(q<1\). Therefore, for the entropy (5) the properties of the Tsallis entropy (3) are completely preserved with respect to the replacement of the ranges of \(0<q<1\) to \(1<q<\infty \) under the rule \(q\rightarrow 1/q\). The comprehensive description of the properties of the Tsallis entropy (3) can be found, for example, in Ref. [19].

The statistical mechanics based on this entropy is defined by Eq. (5) with the probabilities \(p_{i}\) of the microstates of the system normalized to unity

$$\begin{aligned} \varphi =\sum \limits _{i} p_{i} - 1 = 0 \end{aligned}$$
(6)

and by the standard expectation values

$$\begin{aligned} \langle A \rangle = \sum \limits _{i} p_{i} A_{i}. \end{aligned}$$
(7)

Here, the definition of the expectation values (7) is consistent with probability normalization constraint (6). Let us consider the grand canonical ensemble. The thermodynamic potential \(\varOmega \) of the grand canonical ensemble is

$$\begin{aligned} \varOmega= & {} \langle H \rangle -TS-\mu \langle N \rangle \nonumber \\= & {} \sum \limits _{i} p_{i} \left[ E_{i}-\mu N_{i} - T q \frac{p_{i}^{\frac{1}{q}-1}-1}{q-1}\right] , \end{aligned}$$
(8)

where \(\langle H \rangle =\sum _{i} p_{i} E_{i}\) is the mean energy of the system, \(\langle N \rangle =\sum _{i} p_{i} N_{i}\) is the mean number of particles of the system, and \(E_{i}\) and \(N_{i}\) are the energy and the number of particles, respectively, in the ith microscopic state of the system.

The unknown probabilities \(\{p_{i}\}\) are obtained from the second law of thermodynamics (the principle of maximum entropy). In the grand canonical ensemble they are found from the constrained local extrema of the thermodynamic potential (8) by the method of the Lagrange multipliers (see, for example, Refs. [22,23,24]):

$$\begin{aligned} \varPhi= & {} \varOmega - \lambda \varphi , \end{aligned}$$
(9)
$$\begin{aligned} \frac{\partial \varPhi }{\partial p_{i}}= & {} 0, \end{aligned}$$
(10)

where \(\varPhi \) is the Lagrange function and \(\lambda \) is an arbitrary real constant. Substituting Eqs. (6) and (8) into Eqs. (9), (10) and using Eq. (6) again, we obtain the normalized equilibrium probabilities of the grand canonical ensemble of the q-dual statistics as

$$\begin{aligned} p_{i} = \left[ 1+(1-q)\frac{\Lambda -E_{i}+\mu N_{i}}{T}\right] ^{\frac{q}{1-q}} \end{aligned}$$
(11)

and

$$\begin{aligned} \sum \limits _{i} \left[ 1+(1-q)\frac{\Lambda -E_{i}+\mu N_{i}}{T}\right] ^{\frac{q}{1-q}}=1, \end{aligned}$$
(12)

where \(\Lambda \equiv \lambda -T\) and \(\partial E_{i}/\partial p_{i}=\partial N_{i}/\partial p_{i}=0\). In the Gibbs limit \(q\rightarrow 1\), the probability \(p_{i}=\exp [(\Lambda -E_{i}+\mu N_{i})/T]\), where \(\Lambda =-T\ln Z\) is the thermodynamic potential of grand canonical ensemble and \(Z=\sum _{i} \exp [-(E_{i} - \mu N_{i})/T]\) is the partition function.

Substituting Eq. (11) into (7), we obtain the statistical averages of the q-dual statistics in the grand canonical ensemble as

$$\begin{aligned} \langle A \rangle = \sum \limits _{i} A_{i} \left[ 1+(1-q)\frac{\Lambda -E_{i}+\mu N_{i}}{T}\right] ^{\frac{q}{1-q}}, \end{aligned}$$
(13)

where the norm function \(\Lambda \) is the solution of Eq. (12). Note that the quantities (11)–(13) of the q-dual statistics exactly recover the corresponding quantities of the Tsallis-1 statistics under the multiplicative transformation of the entropic parameter \(q\rightarrow 1/q\) (see Ref. [23]).

Let us rewrite the norm equation (12) and the statistical averages (13) in the integral representation in the case when the parameter \(q>1\). Using the integral representation of the Gamma-function [25], we can rewrite the norm equation (12) for \(q>1\) in the following form

$$\begin{aligned} \frac{1}{\varGamma \left( \frac{q}{q-1}\right) } \int \limits _{0}^{\infty } t^{\frac{1}{q-1}} e^{-t\left[ 1+(1-q)\frac{\Lambda -\varOmega _{G}\left( \beta '\right) }{T}\right] } dt =1 \end{aligned}$$
(14)

or

$$\begin{aligned}&\sum \limits _{n=0}^{\infty } \frac{1}{n!\varGamma \left( \frac{q}{q-1}\right) } \int \limits _{0}^{\infty } t^{\frac{1}{q-1}} e^{-t\left[ 1+(1-q)\frac{\Lambda }{T}\right] } \nonumber \\&\quad \left( -\beta '\varOmega _{G}\left( \beta '\right) \right) ^{n} dt = 1, \end{aligned}$$
(15)

where \(\beta '=t(q-1)/T\) and

$$\begin{aligned} \varOmega _{G}\left( \beta '\right)= & {} -\frac{1}{\beta '} \ln Z_{G}\left( \beta '\right) , \end{aligned}$$
(16)
$$\begin{aligned} Z_{G}\left( \beta '\right)= & {} \sum \limits _{i} e^{-\beta '(E_{i}-\mu N_{i})}. \end{aligned}$$
(17)

The statistical averages (13) in the integral representation for \(q>1\) can be rewritten as

$$\begin{aligned} \langle A \rangle= & {} \frac{1}{\varGamma \left( \frac{q}{q-1}\right) } \int \limits _{0}^{\infty } t^{\frac{1}{q-1}} e^{-t\left[ 1+(1-q)\frac{\Lambda -\varOmega _{G}\left( \beta '\right) }{T}\right] } \langle A \rangle _{G}\left( \beta '\right) dt \nonumber \\= & {} \sum \limits _{n=0}^{\infty } \frac{1}{n!\varGamma \left( \frac{q}{q-1}\right) } \int \limits _{0}^{\infty } t^{\frac{1}{q-1}} e^{-t\left[ 1+(1-q)\frac{\Lambda }{T}\right] } \nonumber \\&\times \left( -\beta '\varOmega _{G}\left( \beta '\right) \right) ^{n} \langle A \rangle _{G}\left( \beta '\right) dt, \end{aligned}$$
(18)

where

$$\begin{aligned} \langle A \rangle _{G}\left( \beta '\right) =\frac{1}{Z_{G}\left( \beta '\right) } \sum \limits _{i} A_{i} e^{-\beta '(E_{i}-\mu N_{i})}. \end{aligned}$$
(19)

Equation (18) links the statistical averages of the q-dual statistics with the corresponding statistical averages (19) of the Boltzmann-Gibbs statistics.

Transverse momentum distribution Let us consider the relativistic ideal gas with the Fermi–Dirac \((\eta =1)\), Bose–Einstein \((\eta =-1)\) and Maxwell–Boltzmann \((\eta =0)\) statistics of particles in the grand canonical ensemble. The transverse momentum distribution of particles is related to the mean occupation numbers of the ideal gas as

$$\begin{aligned} \frac{d^{2}N}{dp_{T}dy} = \frac{V}{(2\pi )^{3}} \int \limits _{0}^{2\pi } d\varphi p_{T} \varepsilon _{{\mathbf {p}}} \ \sum \limits _{\sigma } \langle n_{{\mathbf {p}}\sigma }\rangle , \end{aligned}$$
(20)

where \(\varepsilon _{{\mathbf {p}}}=m_{T} \cosh y\) is the single-particle energy, \(p_{T}\) and y are the transverse momentum and rapidity variables, respectively, and \(m_{T}=\sqrt{p_{T}^{2}+m^{2}}\). The mean occupation numbers \(\langle n_{{\mathbf {p}}\sigma }\rangle \) are calculated from Eq. (18) using the mean occupation numbers of the ideal gas of the Boltzmann-Gibbs statistics:

$$\begin{aligned} \langle n_{{\mathbf {p}}\sigma } \rangle _{G}\left( \beta '\right) = \frac{1}{e^{\beta ' (\varepsilon _{{\mathbf {p}}}-\mu )}+\eta }. \end{aligned}$$
(21)

Considering Eqs. (18), (20) and (21), we obtain

$$\begin{aligned} \frac{d^{2}N}{dp_{T}dy}= & {} \frac{gV}{(2\pi )^{2}} p_{T} m_{T} \cosh y \sum \limits _{n=0}^{\infty } \frac{1}{n!\varGamma \left( \frac{q}{q-1}\right) } \nonumber \\&\times \int \limits _{0}^{\infty } t^{\frac{1}{q-1}} e^{-t\left[ 1+(1-q)\frac{\Lambda }{T}\right] } \frac{\left( -\beta '\varOmega _{G}\left( \beta '\right) \right) ^{n}}{e^{\beta ' (m_{T} \cosh y-\mu )}+\eta } dt, \end{aligned}$$
(22)

where

$$\begin{aligned} -\beta '\varOmega _{G}\left( \beta '\right) = \sum \limits _{{\mathbf {p}},\sigma } \ln \left[ 1+\eta e^{-\beta ' (\varepsilon _{{\mathbf {p}}}-\mu )} \right] ^{\frac{1}{\eta }} \end{aligned}$$
(23)

is the thermodynamic potential for the relativistic ideal gas of the Boltzmann–Gibbs statistics of microstates and \(\varepsilon _{{\mathbf {p}}}=\sqrt{{\mathbf {p}}^2+m^{2}}\). The norm function \(\Lambda \) for \(\eta =-1,0,1\) is calculated from the norm equation (15) using Eq. (23).

Maxwell–Boltzmann statistics of particles. In the limit \(\eta \rightarrow 0\), the thermodynamic potential of the ideal gas of the Boltzmann-Gibbs statistics (23) can be written as

$$\begin{aligned} \varOmega _{G}\left( \beta '\right) = - \frac{gV}{2\pi ^{2}} \frac{m^{2}}{\beta '^{2}} e^{\beta '\mu } K_{2}\left( \beta 'm\right) , \end{aligned}$$
(24)

where \(K_{\nu }(z)\) is the modified Bessel function of the second kind. Substituting Eq. (24) into (15), we obtain the norm equation for the Maxwell-Boltzmann statistics of particles as

$$\begin{aligned}&\sum \limits _{n=0}^{\infty } \frac{\omega ^{n}}{n!} \frac{1}{\varGamma \left( \frac{q}{q-1}\right) } \int \limits _{0}^{\infty } t^{\frac{1}{q-1}-n} e^{-t\left[ 1+(1-q)\frac{\Lambda +\mu n}{T}\right] } \nonumber \\&\quad \left( K_{2}\left( \frac{t(q-1)m}{T} \right) \right) ^{n} dt = 1, \end{aligned}$$
(25)

where

$$\begin{aligned} \omega =\frac{gVTm^{2}}{2\pi ^{2}} \frac{1}{q-1}. \end{aligned}$$
(26)

Substituting Eq. (24) into (22) and taking the limit \(\eta =0\), we obtain the transverse momentum distribution for the Maxwell–Boltzmann statistics of particles as

$$\begin{aligned} \frac{d^{2}N}{dp_{T}dy}= & {} \frac{gV}{(2\pi )^{2}} p_{T} m_{T} \cosh y \sum \limits _{n=0}^{\infty } \frac{\omega ^{n}}{n!} \frac{1}{\varGamma \left( \frac{q}{q-1}\right) } \nonumber \\&\times \int \limits _{0}^{\infty } t^{\frac{1}{q-1}-n} e^{-t\left[ 1+(1-q)\frac{\Lambda -m_{T} \cosh y+\mu (n+1)}{T}\right] } \nonumber \\&\times \left( K_{2}\left( \frac{t(q-1)m}{T} \right) \right) ^{n} dt. \end{aligned}$$
(27)

Zeroth term approximation Let us find the transverse momentum distribution of particles in the zeroth term approximation, which was introduced in Ref. [20]. In this approximation we retain only the zeroth term \((n=0)\) in the series expansion. Taking \(n=0\) in Eq. (15), we obtain that the norm function \(\Lambda =0\). Substituting \(\Lambda =0\) into Eq. (22) and considering only the zeroth term and the equation

$$\begin{aligned} \frac{1}{e^{x}+\eta } = \sum \limits _{k=0}^{\infty } (-\eta )^{k} e^{-x(k+1)}, \end{aligned}$$
(28)

where \(|e^{-x}|<1\), we obtain the transverse momentum distribution in the zeroth term approximation as

$$\begin{aligned} \frac{d^{2}N}{dp_{T}dy}= & {} \frac{gV}{(2\pi )^{2}} p_{T} m_{T} \cosh y \nonumber \\&\times \sum \limits _{k=0}^{\infty } (-\eta )^{k} \left[ 1- (k+1) (1-q) \frac{m_{T} \cosh y-\mu }{T} \right] ^{\frac{q}{1-q}} \nonumber \\&\mathrm {for} \quad \eta =-1,0,1 \end{aligned}$$
(29)

or

$$\begin{aligned} \frac{d^{2}N}{dp_{T}dy}= & {} \frac{gV}{(2\pi )^{2}} p_{T} m_{T} \cosh y \nonumber \\&\times \left[ 1 - (1-q) \frac{m_{T} \cosh y-\mu }{T} \right] ^{\frac{q}{1-q}} \nonumber \\&\mathrm {for} \quad \eta =0. \end{aligned}$$
(30)

The transverse momentum distribution (30) exactly coincides with the phenomenological Tsallis distribution (1) introduced in [1, 2] (see Eq. (56) in Ref. [1]). Thus, the phenomenological Tsallis distribution (1) for the Maxwell–Boltzmann statistics of particles exactly corresponds to the zeroth term approximation of the q-dual nonextensive statistics based on the entropy (5). Moreover, the transverse momentum distribution (30) of the q-dual nonextensive statistics recovers the transverse momentum distribution of the Tsallis-1 statistics in the zeroth term approximation under the transformation of the entropic parameter \(q\rightarrow 1/q\) [20, 21]. However, the phenomenological Tsallis distributions for the quantum (Fermi–Dirac and Bose–Einstein) statistics of particles introduced in Ref. [1] do not recover the transverse momentum distributions in the zeroth term approximation of both the q-dual nonextensive statistics (cf. Eqs. (31) and (33) of Ref. [1] along with Eq. (29) of the present paper) and the Tsallis statistics [21].

To summarize, we have introduced the q-dual entropy which is defined from the Tsallis entropy by the multiplicative transformation of the entropic parameter \(q\rightarrow 1/q\). The q-dual nonextensive statistics based on this statistical entropy and consistently defined expectation values of thermodynamic quantities was formulated in the framework of the grand canonical ensemble. The thermodynamic potential was obtained from the fundamental thermodynamic potential by the Legendre transform. The probabilities of microstates of the system were derived from the second law of thermodynamics (the Jaynes maximum entropy principle) by the constrained local extrema of the statistical thermodynamic potential. The norm equation and statistical averages were expressed analytically in terms of the corresponding Boltzmann–Gibbs quantities using the integral representation and the exponential function containing the thermodynamic potential was expanded into the series. We have obtained the exact analytical formulae for the transverse momentum distributions of the relativistic massive particles following the Bose–Einstein, Fermi–Dirac and Maxwell–Boltzmann statistics in the framework of q-dual statistics in the grand canonical ensemble.

In the present study, we have found that the phenomenological Tsallis distribution for the Maxwell–Boltzmann statistics of particles introduced in Refs. [1, 2] is equivalent to the transverse momentum distribution of the q-dual nonextensive statistics in the zeroth term approximation. This demonstrates that the correct link between the phenomenological Tsallis distribution and the fundamental theory of the statistical mechanics is provided by the q-dual entropy instead of the Tsallis one. The q-dual statistics is properly defined in comparison with the Tsallis-2 statistics as the expectation values of this formalism are consistent with probability normalization condition. Thus we have found that the well-known phenomenological single-particle Tsallis distribution is thermodynamically founded only if it belongs to the q-dual nonextensive statistics instead of the Tsallis one.