1 Introduction

Studying the evolution of quark–gluon plasma (QGP), a hot and dense medium created in the high-energy collision experiments at the LHC, CERN and RHIC, BNL or the evolution of the high-energy particles passing through it, has been a subject matter of intense research. To study the evolution of QGP, which is a short lived medium that expands very fast, hydrodynamic equations have been employed for a long time [1,2,3,4,5,6,7]. In addition to the vigorous change of the background medium, a bunch of high-energy particles (jets) passing through the bath also display sufficient modification in their distribution. These particles are formed much earlier than the formation of QGP, and if they happen to pass through the medium, deposit energy inside it. The evolution of their phase-space distribution has been studied microscopically using the Fokker–Planck equation [8,9,10,11] or the Boltzmann transport equation [12,13,14] in many articles.

Apart from the ones mentioned above, there are other types of studies which have utilized the Boltzmann transport equation [15] or the hydrodynamic equation to investigate the evolution of the energy density perturbations in nucleus-nucleus collision [16,17,18,19,20] using, in many cases, different equations of state inspired from the MIT bag model, the mean field treatment etc.

Energy density perturbations can originate due to high-energy particles depositing energy inside the QGP medium. The initial perturbations generated in QGP propagate through the medium showing nonlinear features and undergo modification due to medium effects. Now, if the perturbation does not break, it undergoes particle production. Otherwise, it is completely absorbed inside the bath. For example, when a high energy particle enters quark–gluon plasma, it produces energy density perturbation which may maintain its shape and is detected as the same particle. However, a lighter, low energy particle may not be able to maintain the shape of its perturbation and may be absorbed inside the medium simulating opaqueness.

Studying nonlinear waves in nuclear matter started with the derivation of Korteweg de-Vries (KdV) equation describing the propagation of baryon density pulses in proton-nucleus collisions [16]. Zero and non-zero temperature QGP is considered in [17] where breaking wave solution is found using the MIT bag equation of state. This kind of nonlinear wave structures are also found in [20] for hot QGP in two spatial dimensions with cylindrical symmetry. However, the above calculations were done using the Boltzmann–Gibbs statistics which does not consider a fluctuating ambiance that appears in the high energy nuclear collisions.

Fluctuations of the positions of the nucleons inside the colliding nuclei may lead to fluctuation in energy density. Fluctuation inside the medium through which the perturbation travels, impacts the distribution of the particles which are created after the medium freezes-out. This can be understood from the studies of the particle spectra originating from the p–p or Pb–Pb collisions taking place in high-energy collision experiments. In the analyses of those data, experimentally obtained particle distributions have been seen to be describable by the Tsallis-like power law distributions [21,22,23,24,25,26,27,28,29,30] (also called the ‘Tsallis non-extensive distributions’) that arise in a fluctuating environment [31,32,33,34] which is very well the case for the QGP medium [35].

In the present study, we examine the evolution of the energy density perturbation inside a hot QGP medium with the help of the Euler’s equation in which the Tsallis-MIT bag equation of state [36, 37] has been utilized to take into account the fluctuating ambiance. We observe that the non-extensive bag equation of state gives rise to a breaking wave solution for the first order energy density perturbation. However, the appearance of the breaking waves is influenced by the factors like the Tsallis q parameter, the Tsallis temperature and the parameters describing the initial condition. While carrying out our analysis, we obtain analytical closed form expressions for the Tsallis thermodynamic variables for gases comprising of very light (massless and almost massless) quantum particles (Eqs. 9–12). As far as our knowledge goes, these expressions were not used in any other studies, and our finding is expected to have applications in studies pertaining to thermodynamics of quantum Tsallis gas.

This article is organized as follows. In the next section, we describe the mathematical model containing some discussions about the relativistic Euler’s equation, the conservation equation, the non-extensive equation of state, and the derivation of the thermodynamic variables to be used in the equation of state. Section 3 contains a detailed account of the evolution equation and its solution. Section 4 highlights the results, and we summarize as well as conclude in Sect. 5.

2 Mathematical model

Since here we derive the nonlinear evolution equation for energy density perturbation in hot QGP using Tsallis statistics and relativistic hydrodynamics, we need to evaluate the required thermodynamic variables and the equation(s) of state. As the basic dynamical equations of the system we use the relativistic Euler’s equation and continuity equation for entropy density. We consider hot QGP produced at the LHC, where the baryon number density is zero at the central rapidity region and evaluate the energy density and pressure appearing in Euler’s equation using the Tsallis-MIT bag equations of state.

2.1 Relativistic Euler equation

Throughout this article we will follow the natural units i.e., \(c = 1, \hbar = 1\), \(k_B = 1\). As discussed in [17], the correct description of nonlinear wave structure in QGP including quantum effects should be given by relativistic hydrodynamics of colored superfluids which is too complicated and still in an initial stage. Admitting this fact, we stick to the classical relativistic hydrodynamics to study the nonlinear wave structure in hot QGP following Tsallis distribution.

The velocity field is given by, \(\mathbf {v}=v(x,t){\hat{x}}\) with a magnitude v in the x-direction. We use the one-dimensional relativistic Euler’s equation [38, 39] which is written as,

$$\begin{aligned} \frac{\partial v}{\partial t} + v \frac{\partial v}{\partial x} = \frac{v^2-1}{\epsilon +P} \left( \frac{\partial P}{\partial x}+ v \frac{\partial P}{\partial t}\right) , \end{aligned}$$
(1)

where \(\epsilon \equiv \epsilon (x,t)\) is the energy density and \(P \equiv P(x,t)\) is the pressure.

Our calculation will be aimed at quark–gluon plasma produced at the LHC energies where at the central rapidity region the net baryonic number density is zero, and hence we encounter vanishing chemical potential. This leads to the fact that the particle and the anti-particle distributions are identical.

Now, defining the entropy four-current as \(s^{\mu } = s u^{\mu }\), where \(u^{\mu } \equiv (\gamma , \gamma \mathbf {v})\) is the four-velocity, the continuity equation for the entropy density s can be derived as (for details see [17]),

$$\begin{aligned} v s \left( \frac{\partial v}{\partial t} + v \frac{\partial v}{\partial x}\right) +(1-v^2)\left( \frac{\partial s}{\partial t} + s \frac{\partial v}{\partial x}+v\frac{\partial s}{\partial x}\right) = 0, \nonumber \\ \end{aligned}$$
(2)

where we have used the Lorentz factor as \(\gamma =1/\sqrt{1-v^2}\).

2.2 Non-extensive MIT bag equation of state

From Eqs. (1) and (2), we observe that we have four unknown quantities appearing in two equations. Hence, we require two additional equations which may be provided by the equations of state. In the present work we consider the non-extensive MIT bag model [36, 37] equations of state which enable us to express the pressure and the entropy density in terms of the energy density. In this model, quarks and gluons follow the quantum Tsallis fermionic (‘f’) and bosonic (‘b’) single particle distributions [40] given by,

$$\begin{aligned} n_f=\frac{1}{\left[ 1+(q-1)\frac{E_p-\mu }{T}\right] ^{\frac{q}{q-1}}+1}\nonumber \\ n_b=\frac{1}{\left[ 1+(q-1)\frac{E_p-\mu }{T}\right] ^{\frac{q}{q-1}}-1} , \end{aligned}$$
(3)

where \(E_p=\sqrt{p^2+m^2}\) is the single particle energy of a particle of mass m, \(\mu \) is the chemical potential, q is the Tsallis parameter, and T is the Tsallis temperature.

Tsallis single particle distributions can be obtained from the Tsallis entropy defined by [41],

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

where \(p_i\) are the probabilities of the micro-states. This entropic form is non-additive, i.e. the total entropy of the system is not a summation of those of its sub-parts. This situation may appear due non-ideal plasma effects, fluctuation etc [42, 43]. Extremizing the Tsallis thermodynamic potential with respect to \(p_i\) one arrives at the expression for the Tsallis probabilities and the single particle distribution [44]. It is shown in Ref. [44] that the transverse momentum spectra obtained from the single particle distribution is expressed in terms of an infinite summation and the zeroth term coincides with the phenomenological classical Tsallis distribution widely used in the literature. The forms of the Tsallis quantum distributions used in literature are similar to (but not always exactly the same as) the zeroth order approximated term when one uses the factorization approximation used in Ref. [40]. In this work we restrict ourselves within this approximated version of the Tsallis quantum distribution. We reserve the analysis with the exact distributions for future.

Assuming an ideal gas of quarks and gluons, we write down the expressions for the thermodynamic variables like the net baryonic number density \(\rho _{\mathrm {B}}\), energy density \(\epsilon \) and pressure P.

$$\begin{aligned} \rho _{\mathrm {B,bag}}= & {} \rho _f-{\bar{\rho }}_f, \end{aligned}$$
(5)
$$\begin{aligned} \epsilon _{\mathrm {bag}}= & {} {\mathcal {B}}+\epsilon _{b}+2\epsilon _f, \end{aligned}$$
(6)
$$\begin{aligned} P_{\mathrm {bag}}= & {} -{\mathcal {B}}+P_{b}+2P_f, \end{aligned}$$
(7)

where \({\mathcal {B}}\) is the bag constant which represents the pressure of the vacuum [45], and the subscript ‘b’(‘f’) signifies the bosonic (fermionic) contribution of the thermodynamic variables. The factor 2 in front of the fermionic parts lets us consider both particles and anti-particles. According to the bag model, quarks and gluons are assumed to stay in a spherical cavity (the ‘bag’) in QCD vacuum and the constant \({\mathcal {B}}\) takes care of the confinement property. \(\rho _f\) (\({\bar{\rho }}_f\)) signifies the particle (anti-particle) number density and as already mentioned in the previous section, the net baryonic number at the central rapidity region at the LHC energies is zero. Hence, we have zero baryonic chemical potential which leads to \(\rho _f={\bar{\rho }}_f\).

For the energy density and pressure, we identify that the bag variables have contributions from the massless bosons (gluons) and from the massive (where mass is much less than temperature) fermions (up and down quarks) and they can be expressed in terms of \(n_b\) and \(n_f\) in the following way:

$$\begin{aligned} \epsilon _{i}= g\int \frac{d^3p}{(2\pi )^3}~E_p~n_i, P_i=g\int \frac{d^3p}{(2\pi )^3}\frac{p^{2}}{3E_p}~n_i, \end{aligned}$$
(8)

where g is the degeneracy factor and \(i=f,b\).

2.2.1 Energy density and pressure for massless bosons

The closed analytic expressions of the energy density and pressure of a non-extensive massless bosonic gas are given by,

$$\begin{aligned} \epsilon _b= & {} \frac{g T^4}{2 \pi ^2 (q-1)^3 q} \left[ 3 \psi ^{(0)}\left( \frac{3}{q}-2\right) + \psi ^{(0)}\left( \frac{1}{q}\right) \right. \nonumber \\&\left. - 3 \psi ^{(0)}\left( \frac{2}{q}-1\right) - \psi ^{(0)}\left( \frac{4}{q}-3\right) \right] , \end{aligned}$$
(9)
$$\begin{aligned} P_b= & {} \frac{g T^4}{6 \pi ^2 (q-1)^3 q} \left[ 3 \psi ^{(0)}\left( \frac{3}{q}-2\right) + \psi ^{(0)}\left( \frac{1}{q}\right) \right. \nonumber \\&\left. - 3 \psi ^{(0)}\left( \frac{2}{q}-1\right) - \psi ^{(0)}\left( \frac{4}{q}-3\right) \right] , \end{aligned}$$
(10)

and they are related by \(\epsilon _b=3P_b\). In Eqs. (9) and (10) \(\psi ^{(0)}\) is the digamma function [47]. We defer the detailed calculations till the appendix of the paper.

2.2.2 Energy density and pressure for massive fermions

We evaluate the fermionic thermodynamic variables up to \({\mathcal {O}}(m^2T^2)\). We consider the plasma to be consisted of the up and down quarks having masses of 5–10 MeV and we have checked that for a wide range of the q and T parameter values appearing in the phenomenological studies of high-energy collisions, the \({\mathcal {O}}(m^2T^2)\) approximation works very well when the mass of the particle is around 10 MeV. In some papers it has been suggested that the \({\mathcal {O}}(m^2T^2)\) contribution helps one account for the non-perturbative effects [46] in QGP. The closed analytic expressions for the thermodynamic variables for a non-extensive gas of massive fermions up to \({\mathcal {O}}(m^2T^2)\) are given by the following equations:

$$\begin{aligned} \epsilon _f= & {} \frac{g T^4}{2 \pi ^2 (q-1)^3 q} \left[ 3 \varPhi \left( -1,1,\frac{2}{q}-1\right) \right. \nonumber \\&\left. - 3 \varPhi \left( -1,1,\frac{3}{q}-2\right) \right. \nonumber \\&\left. +\varPhi \left( -1,1,\frac{4}{q}-3\right) -\varPhi \left( -1,1,\frac{1}{q}\right) \right] \nonumber \\&-\frac{g m^2 T^2}{4 \pi ^2 (q-1) q} \left[ \varPhi \left( -1,1,\frac{2}{q}-1\right) -\varPhi \left( -1,1,\frac{1}{q}\right) \right] , \nonumber \\ \end{aligned}$$
(11)
$$\begin{aligned} P_f= & {} \frac{g T^4}{6 \pi ^2 (q-1)^3 q} \left[ 3 \varPhi \left( -1,1,\frac{2}{q}-1\right) \right. \nonumber \\&\left. -3 \varPhi \left( -1,1,\frac{3}{q}-2\right) \right. \nonumber \\&\left. +\varPhi \left( -1,1,\frac{4}{q}-3\right) -\varPhi \left( -1,1,\frac{1}{q}\right) \right] \nonumber \\&-\frac{g m^2 T^2}{4 \pi ^2 (q-1) q} \left[ \varPhi \left( -1,1,\frac{2}{q}-1\right) -\varPhi \left( -1,1,\frac{1}{q}\right) \right] , \nonumber \\ \end{aligned}$$
(12)

where \(\varPhi \) is the Lerch’s transcendent [47]. For this section also, we defer the detailed calculations till the appendix.

2.2.3 Calculation of \(\epsilon _{\mathrm {bag}}\) and \(P_{\mathrm {bag}}\)

Using Eqs. (9)–(12), we write down the expressions for the energy density and pressure in the non-extensive bag model,

$$\begin{aligned} \epsilon _{\mathrm {bag}} (x,t)= & {} {\mathcal {B}} + \epsilon _{b,1}T^4 + 2 (\epsilon _{f,1}T^4+\epsilon _{f,2}m^2T^2), \nonumber \\= & {} {\mathcal {B}} + \epsilon _{\mathrm {bag},1}T^4 + \epsilon _{\mathrm {bag},2}m^2T^2 \end{aligned}$$
(13)
$$\begin{aligned} P_{\mathrm {bag}} (x,t)= & {} -{\mathcal {B}} + P_{b,1}T^4 + 2 (P_{f,1}T^4+P_{f,2}m^2T^2) \nonumber \\= & {} -{\mathcal {B}} + P_{\mathrm {bag},1}T^4 + P_{\mathrm {bag},2}m^2T^2. \end{aligned}$$
(14)

When written out long-hand, \(\epsilon _{i,\ell }\) and \(P_{i,\ell }\) (\(i=f,b\); \(\ell =1,2\)) are given by,

$$\begin{aligned} \epsilon _{b,1}= & {} \frac{g}{2 \pi ^2 (q-1)^3 q} \left[ 3 \psi ^{(0)}\left( \frac{3}{q}-2\right) + \psi ^{(0)}\left( \frac{1}{q}\right) \right. \nonumber \\&\left. - 3 \psi ^{(0)}\left( \frac{2}{q}-1\right) - \psi ^{(0)}\left( \frac{4}{q}-3\right) \right] , \end{aligned}$$
(15)
$$\begin{aligned} \epsilon _{f,1}= & {} \frac{g}{2 \pi ^2 (q-1)^3 q} \left[ 3 \varPhi \left( -1,1,\frac{2}{q}-1\right) \right. \nonumber \\&\left. -3 \varPhi \left( -1,1,\frac{3}{q}-2\right) +\varPhi \left( -1,1,\frac{4}{q}-3\right) \right. \nonumber \\&\left. -\varPhi \left( -1,1,\frac{1}{q}\right) \right] , \end{aligned}$$
(16)
$$\begin{aligned} \epsilon _{f,2}= & {} -\frac{g}{4 \pi ^2 (q-1) q} \left[ \varPhi \left( -1,1,\frac{2}{q}-1\right) \right. \nonumber \\&\left. -\varPhi \left( -1,1,\frac{1}{q}\right) \right] , \end{aligned}$$
(17)
$$\begin{aligned} P_{b,1}= & {} \frac{\epsilon _{b,1}}{3},~P_{f,1} = \frac{\epsilon _{f,1}}{3},~ P_{f,2} = \epsilon _{f,2} . \end{aligned}$$
(18)

In terms of the above variables, \(\epsilon _{\mathrm{{bag}},\ell }\) and \(P_{\mathrm{{bag}},\ell }\) (\(\ell =1,2\)) are given by,

$$\begin{aligned} \epsilon _{\mathrm {bag},1}= & {} \epsilon _{b,1} + 2\epsilon _{f,1}, ~P_{\mathrm {bag},1} = P_{b,1} + 2P_{f,1}, \nonumber \\ ~\epsilon _{\mathrm {bag},2}= & {} 2\epsilon _{f,2},~P_{\mathrm {bag},2} = 2P_{f,2}. \end{aligned}$$
(19)

Once we get the pressure, the entropy density is given by the partial derivative of the pressure with respect to the temperature T, i.e. \(s_{\mathrm{{bag}}}=\partial P_{\mathrm{{bag}}}/\partial T\). Hence we obtain,

$$\begin{aligned} s_{\mathrm {bag}} = 4 P_{\mathrm {bag},1}T^3 + 2 P_{\mathrm {bag},2}m^2T. \end{aligned}$$
(20)

It has been verified that in absence of the baryonic chemical potential, the pressure, energy density and the entropy density obey the following relationship,

$$\begin{aligned} \epsilon _{\mathrm {bag}} + P_{\mathrm {bag}} = T s_{\mathrm {bag}}. \end{aligned}$$
(21)

2.2.4 Equations of state

In order to find out the equations of state, we express the pressure P and the entropy density s as functions of the energy density \(\epsilon \). We solve Eq. (13) for the temperature T in terms of the bag model energy density \(\epsilon _{\mathrm{{bag}}}\), and put the solution in Eqs. (14), and (20) to express \(P_{\mathrm{{bag}}}\), and \(s_{\mathrm{{bag}}}\) as functions of \(\epsilon _{\mathrm{{bag}}}\). Solving Eq. (13) for a real and positive value of T and denoting the solution with \(T_{\mathrm{{sol}}}\) we obtain,

$$\begin{aligned} T_{\mathrm {sol}} = \left( \frac{-m^2 \epsilon _{\mathrm{{bag}},2} + {\mathcal {R}}(\epsilon _{\mathrm{{bag}}})}{2 \epsilon _{\mathrm{{bag}},2}} \right) ^{\frac{1}{2}}, \end{aligned}$$
(22)

where \({\mathcal {R}}(\epsilon _{\mathrm{{bag}}}) = \sqrt{m^4 \epsilon _{\mathrm{{bag}},2}^2 + 4 \epsilon _{\mathrm{{bag}},1}(\epsilon _{\mathrm{{bag}}} - {\mathcal {B}})}\). The expressions of \(P_{\mathrm{{bag}}}(\epsilon _{\mathrm{{bag}}})\) and \(s_{\mathrm{{bag}}}(\epsilon _{\mathrm{{bag}}})\) are given by,

$$\begin{aligned} P_{\mathrm{{bag}}}= & {} -{\mathcal {B}} + P_{\mathrm {bag,1}}\left( \frac{-m^2 \epsilon _{\mathrm{{bag}},2} + {\mathcal {R}}(\epsilon _{\mathrm{{bag}}})}{2 \epsilon _{\mathrm{{bag}},2}} \right) ^{2} \nonumber \\&+ P_{\mathrm{{bag,2}}} m^2 \left( \frac{-m^2 \epsilon _{\mathrm{{bag,2}}}+ R(\epsilon _{\mathrm{{bag}}})}{2\epsilon _{\mathrm{{bag,1}}}} \right) , \nonumber \\ s_{\mathrm{{bag}}}= & {} 4 P_{\mathrm {bag,1}} \left( \frac{-m^2 \epsilon _{\mathrm{{bag}},2} + {\mathcal {R}}(\epsilon _{\mathrm{{bag}}})}{2 \epsilon _{\mathrm{{bag}},2}} \right) ^{\frac{3}{2}} \nonumber \\&+2 P_{\mathrm{{bag,2}}} m^2 \left( \frac{-m^2 \epsilon _{\mathrm{{bag}},2} + {\mathcal {R}}(\epsilon _{\mathrm{{bag}}})}{2 \epsilon _{\mathrm{{bag}},2}} \right) ^{\frac{1}{2}}. \end{aligned}$$
(23)

3 Nonlinear evolution equation of energy density perturbation

In this section, we derive the evolution equation of energy density perturbation of the QGP system using the well known Reductive Perturbation Theory (RPT) [48] which helps one to deal with the perturbation which can not be neglected with respect to the mean value. For this procedure, we consider the two dynamical equations i.e, the relativistic Euler’s equation (Eq. 1) and the entropy conservation equation (Eq. 2). We expand the dependent variables in terms of a perturbation parameter \(\sigma \) following the RPT. Finally combining Eqs. (1) and (2) and solving the system of equations order by order we arrive at the space-time evolution of a perturbation of the energy density in quark–gluon plasma. Now, we put \(P_{\mathrm{{bag}}}(\epsilon _{\mathrm{{bag}}})\), and \(s_{\mathrm{{bag}}}(\epsilon _{\mathrm{{bag}}})\) in Eqs. (1) and (2) and in the resulting equations express (\(\epsilon _{\mathrm{{bag}}}-{\mathcal {B}}\)) in terms of temperature using Eq. (13). Hence we obtain,

$$\begin{aligned}&\frac{\partial v}{\partial t} + v \frac{\partial v}{\partial x} + {\mathcal {E}}_1(1-v^2) \left( \frac{\partial \epsilon _{\mathrm {bag}}}{\partial x} + v \frac{\partial \epsilon _{\mathrm {bag}}}{\partial t} \right) = 0 \end{aligned}$$
(24)
$$\begin{aligned}&{\mathcal {C}}_1 (1-v^2) \left( v \frac{\partial \epsilon _{\mathrm {bag}}}{\partial x} + \frac{\partial \epsilon _{\mathrm {bag}}}{\partial t} \right) + {\mathcal {C}}_2 \left( \frac{\partial v}{\partial x}+v\frac{\partial v}{\partial t} \right) =0 , \nonumber \\ \end{aligned}$$
(25)

where \({\mathcal {E}}_1,~{\mathcal {C}}_1\) and \({\mathcal {C}}_2\) are given by,

$$\begin{aligned} {\mathcal {C}}_1= & {} m^2 \epsilon _{\mathrm {bag},2}+2 T^2 \epsilon _{\mathrm {bag},1} \end{aligned}$$
(26)
$$\begin{aligned} {\mathcal {C}}_2= & {} \left( m^2 \epsilon _{\mathrm {bag},2}+\frac{2}{3} \epsilon _{\mathrm {bag},1}T^2\right) \nonumber \\&\times \left( 4T^4\epsilon _{\mathrm {bag},1}+2m^2T^2 \epsilon _{\mathrm {bag},2}\right) \end{aligned}$$
(27)
$$\begin{aligned} {\mathcal {E}}_1= & {} \frac{\left( m^2 \epsilon _{\mathrm {bag},2}+\frac{2}{3} \epsilon _{\mathrm {bag},1}T^2\right) }{\left( m^2\epsilon _{\mathrm {bag},2}+2T^2\epsilon _{\mathrm {bag},1}\right) \left( \frac{4}{3}T^4 \epsilon _{\mathrm {bag},1} + 2m^2T^2 \epsilon _{\mathrm {bag},2} \right) }. \nonumber \\ \end{aligned}$$
(28)

Now we express Eqs. (24), and (25) in terms of the following dimensionless variables:

$$\begin{aligned} {\hat{\epsilon }} = \frac{\epsilon _{\mathrm{{bag}}}}{\epsilon _0}, ~ {\hat{v}} = \frac{v}{c_s}, \end{aligned}$$
(29)

where \(\epsilon _0\) is a reference energy density, and \(c_s\) is the velocity of sound. We express \({\hat{\epsilon }}\) and \({\hat{v}}\) in terms of a small parameter \(\sigma \) following RPT as,

$$\begin{aligned} {\hat{\epsilon }}= & {} 1+ \sigma \epsilon _1 + \sigma ^2\epsilon _2 + {\mathcal {O}}(\sigma ^3) \nonumber \\ {\hat{v}}= & {} \sigma v_1 + \sigma ^2 v_2 + {\mathcal {O}}(\sigma ^3). \end{aligned}$$
(30)

We also change the independent variables from (x, t) to (\(\xi \), \(\tau \)) which are related by,

$$\begin{aligned} \xi = \sigma ^{\frac{1}{2}} \frac{(x-c_s t)}{R}, ~ \tau = \sigma ^{\frac{3}{2}} \frac{c_s t}{R}, \end{aligned}$$
(31)

where R is the characteristic length scale of the problem.

The stretched coordinates used in Eq. (31) is a part of the ‘reductive perturbation technique’ (RPT) where the small parameter \(\sigma \) signifies the smallness of the perturbed quantities relative to the corresponding equilibrium quantities. Eq. (31) involves two time scales in order to explain fast dynamics at the linear limit and slower dynamics occurring at the nonlinear level. This means, at short time scale (at the lowest limit of \(\sigma \)) the perturbation obeys linear equation and travels with velocity \(c_s\). Over a longer time scale, the wave form is influenced by nonlinearity giving rise to the breaking wave structure. The particular scaling with \(\sigma \) used in Eq. (31) comes from the idea of two time scales for long waves. The scaling must satisfy the required invariance and compatibility condition as discussed in [49, 50] in order to keep the perturbation scheme valid. A similar scaling is also used in [17] for studying nonlinear waves in cold and hot QGP.

Using the expansion (30) in Eqs. (24) and (25) and equating the coefficients of different powers of \(\sigma \) to zero, we get the system of differential equations for the perturbations.

3.1 Calculation at \({\mathcal {O}}(\sigma )\)

At the lowest order, i.e at \({\mathcal {O}}(\sigma )\) of Eqs. (24), (25), we get the following relations,

$$\begin{aligned} c_s^2 \frac{\partial v_1}{\partial \xi } = \epsilon _0 {\mathcal {E}}_1 \frac{\partial \epsilon _1}{\partial \xi }, \ \ {\mathcal {C}}_2 \frac{\partial v_1}{\partial \xi } = \epsilon _0 {\mathcal {C}}_1 \frac{\partial \epsilon _1}{\partial \xi }. \end{aligned}$$
(32)

Solving (32) we can determine the speed of sound \(c_s\) as,

$$\begin{aligned} c_s^2 = \frac{{\mathcal {C}}_2 {\mathcal {E}}_1}{{\mathcal {C}}_1}. \end{aligned}$$
(33)

3.2 Calculation at \(O(\sigma ^2)\)

At the next order i.e, at \({\mathcal {O}}(\sigma ^2)\) we get the nonlinear evolution equation for the energy density perturbation \(\epsilon _1\). Using the relations between \(v_1\) and \(\epsilon _1\) from Eqs. (32) and  (33) we finally get,

$$\begin{aligned}&\frac{\partial \epsilon _1}{\partial \tau } + \frac{2 \epsilon _0 \epsilon _{\mathrm {bag},1} \epsilon _1}{3 m^4 \epsilon _{\mathrm {bag},2}^2+8 m^2 T^2 \epsilon _{\mathrm {bag},1} \epsilon _{\mathrm {bag},2}+4 T^4 \epsilon _{\mathrm {bag},1}^2} \frac{\partial \epsilon _1}{\partial \xi } \nonumber \\&\quad =0, \nonumber \\ \end{aligned}$$
(34)

which is the main result of our work which estimates the evolution of the first order energy density perturbation.

Coming back to the x-t space, we get,

$$\begin{aligned}&\frac{\partial {\hat{\epsilon }}_1}{\partial t} + c_s \frac{\partial {\hat{\epsilon }}_1}{\partial x} \nonumber \\&\quad + \frac{2 c_s \epsilon _0 \epsilon _{\mathrm {bag},1} {\hat{\epsilon }}_1}{3 m^4 \epsilon _{\mathrm {bag},2}^2+ 8 m^2 T^2 \epsilon _{\mathrm {bag},1} \epsilon _{\mathrm {bag},2}+4 T^4 \epsilon _{\mathrm {bag},1}^2} \frac{\partial {\hat{\epsilon }}_1}{\partial x} = 0, \nonumber \\ \end{aligned}$$
(35)

where \({\hat{\epsilon }}_1=\sigma \epsilon _1\) is the first order perturbation term in scaled energy density. This differential equation is similar to the form obtained in Ref. [17], differing only in the coefficient of the last non-linear term. For a constant \(T = T_{\mathrm{{con}}}\), the equation (35) turns to be of the form of inviscid Burger’s equation [51] the details of which is discussed in the next section. If we put, \(m=0, \epsilon _{\mathrm {bag},2}=0\), and \(\epsilon _{\mathrm {bag},1}=37\pi ^2/30\), as used in Ref. [17] for the Boltzmann–Gibbs statistics of massless fermions and bosons, we exactly get back the equation derived in there.

3.3 Breaking wave solutions

For a constant \(T = T_{\mathrm{{con}}}\), the coefficient of the nonlinear term of Eq. (34) can be written as,

$$\begin{aligned} B_c = \frac{2 \epsilon _0 \epsilon _{\mathrm {bag},1}}{3 m^4 \epsilon _{\mathrm {bag},2}^2+8 m^2 T_{\mathrm{{con}}}^2 \epsilon _{\mathrm {bag},1} \epsilon _{\mathrm {bag},2}+4 T_{\mathrm{{con}}}^4 \epsilon _{\mathrm {bag},1}^2} . \end{aligned}$$
(36)

Hence Eq. (34) becomes,

$$\begin{aligned} \frac{\partial \epsilon _1}{\partial \tau } + B_c \epsilon _1 \frac{\partial \epsilon _1}{\partial \xi } =0, \end{aligned}$$
(37)

which is nothing but the inviscid Burgers equation [51]. For a linear initial condition, S. Chandrasekhar found the explicit exact solution [52] of (37) as,

$$\begin{aligned} \epsilon _1 (\xi , \tau ) = \frac{1}{B_c} \left( \frac{a \xi + b}{a \tau + 1}\right) , \end{aligned}$$
(38)

where a, b are arbitrary free parameters. Explicit solutions of Eq. (37) for other relevant initial conditions are not known as far as our knowledge goes. In the x-t frame the form of the exact solution (38) becomes,

$$\begin{aligned} \epsilon _1 (x, t) = \frac{1}{B_c} \left( \frac{ a \sigma ^{\frac{1}{2}} (x-c_s t) + b R}{a \sigma ^{\frac{3}{2}} c_s t + R} \right) , \end{aligned}$$
(39)

where \(\sigma , R\) are the small parameter and characteristic length of the system respectively ( defined in the previous section) and \(c_s\) is given by Eq. (33). The solution (39) behaves as a breaking wave the energy of which dissipates quickly with distance unlike solitons.

For a more general initial condition, the solution of

$$\begin{aligned} \frac{\partial u(x,t)}{\partial t} + f(u) \frac{\partial u(x,t)}{\partial x} = 0, \end{aligned}$$
(40)

is given by \(u\left( \chi \equiv x-f(u)t\right) \) [53]. Assuming an initial (\(t=0\)) energy perturbation distribution of

$$\begin{aligned} {\hat{\epsilon }}_1 (x) = A ~\mathrm {Sech}^2\left( \frac{x}{B}\right) , \end{aligned}$$
(41)

the solution of Eq. (35) can be written as,

$$\begin{aligned} {\hat{\epsilon }}_1 (\chi ) = A ~\mathrm {Sech}^2\left( \frac{\chi }{B}\right) , \end{aligned}$$
(42)

where

$$\begin{aligned} \chi = x-c_s(1+B_c {\hat{\epsilon }}_1 (\chi ))t \end{aligned}$$
(43)

4 Results and discussion

For a fixed value of the reference energy density (\(\epsilon _0=0.16\) GeV \(fm^{-3}\)) and mass \(m=10\) MeV (which is of the order of the down quark mass), the final solution depends on the amplitude A of the initial wave, its width B, the Tsallis entropic parameter q, and the Tsallis temperature T.

In Fig. 1, we plot the solutions of Eq. (35) on the x-t plane. For the sake of a better understanding, we present Fig. 2 which displays a two-dimensional version of Fig. 1. For both Figs. 1 and 2, \(A=2.5\), \(B=0.5\) fm, \(q=1.08\), and \(T=140\) MeV. It is observed from Fig. 2 that at around \(t=10\) fm, the solutions start becoming multiple-valued functions thereby implying the breaking of the waves.

Fig. 1
figure 1

Evolution of the energy density perturbation with time t and space coordinate x. \(A=2.5\), \(B=0.5\) fm, \(q=1.08\), and \(T=140\) MeV

Full size image
Fig. 2
figure 2

Solutions of Fig. 1 at different times on a two-dimensional plot. \(A=2.5\), \(B=0.5\) fm, \(q=1.08\), and \(T=140\) MeV

Full size image
Fig. 3
figure 3

Solutions at \(t=15\) fm for different amplitudes of the initial wave. \(B=0.5\) fm, \(q=1.08\), and \(T=140\) MeV

Full size image
Fig. 4
figure 4

Solutions at \(t=15\) fm for different widths of the initial wave. \(A=0.5\), \(q=1.08\), and \(T=140\) MeV

Full size image
Fig. 5
figure 5

Solutions at \(t=15\) fm for different Tsallis temperature values. \(A=2.5\), \(B=0.5\) fm, and \(q=1.08\)

Full size image
Fig. 6
figure 6

Solutions at \(t=15\) fm for different Tsallis parameter values. \(A=2.5\), \(B=0.5\) fm, and \(T=140\) MeV

Full size image

The solutions, however, depend on the width and the amplitude of the waves as well. It is apparent from Fig. 3, that beyond \(A=1\) (when \(t=15\) fm, \(B=0.5\) fm, \(q=1.08\), and \(T=140\) MeV), the solutions start becoming multiple-valued. This means, the more the amplitude of the initial energy density perturbation is, the more it becomes prone to be a breaking wave. On the other hand, a narrower initial energy density perturbation is more likely to give rise to a breaking wave solution (as seen in Fig. 4).

The temperature T, and the Tsallis parameter q also influence the solutions so that moving from a lower to a higher temperature value makes functions single-valued. Hence, the lower the Tsallis temperature is, the more likely it is to get a breaking wave solution (as seen in Fig. 5, where \(t=15\) fm, \(A=2.5\), \(B=0.5\) fm, and \(q=1.08\)). Similarly, a lower value of the q parameter indicates a higher likelihood of obtaining breaking waves (as seen in Fig. 6, where \(t=15\) fm, \(A=2.5\), \(B=0.5\) fm, and \(T=140\) MeV).

5 Summary and conclusions

In summary, we have studied the evolution of the first order energy density perturbation in hot, ideal, and non-extensive quark–gluon plasma. To take into account the fluctuating ambiance, we have used a non-extensive version of the MIT bag model and find a breaking wave solution. However, we observe that the breaking of the energy density perturbation is influenced by the temperature and the Tsallis q parameter which according to Ref. [31] is related to the relative variance in a scale variable i.e. temperature. We observe that the breaking is favoured by decreasing values of both temperature and the q parameter. This may have implications in the LHC phenomenology as the quark–gluon plasma medium formed in this region may have higher initial temperature and q value in comparison with that formed at the experiments having lower beam energies. So, the resulting wave may behave as a stable one even after a long time.

We have also verified that the at a particular instant a wave with higher amplitude and a smaller width is more prone to breaking which is intuitively understandable. During the analysis we have also provided the closed analytical formulae for the Tsallis thermodynamic variables of an ideal gas of massless bosons and very light fermions. This can be seen as an extension of an earlier work for the classical Tsallis distribution by one of us [54]. These results can be used in the studies pertaining to the thermodynamics of quantum Tsallis gases. However, an unapproximated expression for the thermodynamic variables may be used in the future studies to examine whether that affects the present conclusions. In addition to that, the present work relies on some simplifications like one-dimensional expansion. An extension of the present work considering expansion in higher dimensions should be done. Within the domain of the Boltzmann–Gibbs statistics, whether one obtains a breaking wave or a localized wave depends on the equation of state. For example, in [18] an equation of state inspired by the mean field QCD model resulted in the Korteveg-de Vries (KdV) soliton. But, no similar study has been reported in the field of the Tsallis statistics. It will also be interesting to explore this problem. We reserve these studies for future.