*Article* **Nonequilibrium Pion Distribution within the Zubarev Approach**

#### **David Blaschke 1,2,3,***<sup>∗</sup>* **, Gerd Röpke 3,4, Dmitry N. Voskresensky 2,3 and Vladimir G. Morozov <sup>5</sup>**


Received: 11 April 2020; Accepted: 30 April 2020; Published: 3 May 2020

**Abstract:** We discuss how the non-equilibrium process of pion production within the Zubarev approach of the non-equilibrium statistical operator leads to a theoretical foundation for the appearance of a non-equilibrium pion chemical potential for the pion distribution function for which there is experimental evidence in experiments at the CERN LHC.

**Keywords:** Zubarev formalism; non-equilibrium statistical operator; pion chemical potential; LHC

#### **1. Introduction**

The thermal statistical model [1–6] for chemical freeze-out of hadron species gives a successful description of a set of particle ratios produced in heavy-ion collisions (HIC) at different center of mass energies ranging from the energies provided by the Schwerionensynchrotron (SIS-18) at GSI Darmstadt over those of the Alternating Gradient Synchrotron (AGS) at BNL Brookhaven and the Super Proton Synchrotron (SPS) at CERN Geneva up to the highest energies at the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron Collider (LHC) at CERN. It came therefore as a surprise that for LHC at <sup>√</sup>*<sup>s</sup>* <sup>=</sup> 2.76 TeV the measured proton abundances [7,8] do not agree with the most common version of the thermal mode (the inclusion of resonance formation due to (multi-)pion-nucleon interaction and further correlations in the continuum within the Beth–Uhlenbeck approach [9,10] improves the agreement with the experiment) based on the grand canonical ensemble [3,4]. As a possible explanation of this effect it has been suggested that the freeze-out may take place off chemical equilibrium [11–13]. Hereby, a key feature is the enhancement of low-transverse momentum pion spectra above the expectation from equilibrium statistical models which was seen already at lower energies in pion spectra at SPS and clearly seen in the RHIC and the LHC data. The effect can be seen as a precursor of pion Bose–Einstein condensation due to high phase space occupation at low momenta and has consequently been parametrized by adopting a pion chemical potential very close to the pion mass [14,15]. This concept is based on the assumption that the total pion number is dynamically fixed on a time scale between the pion chemical freeze out *tπ*,cfo and the thermal freeze-out (or simply freeze-out) *<sup>t</sup>*fo, *<sup>t</sup>π*,cfo <sup>&</sup>lt; *<sup>t</sup>* <sup>&</sup>lt; *<sup>t</sup>*fo, where at *<sup>t</sup>π*,cfo the pion number becomes frozen and at *t*fo the momentum distributions stop to change [16]. Thereby, for pions, we assume dominance of elastic collisions over inelastic ones. We will assume that for pions the time typical for absorptive processes *<sup>t</sup>π*,abs, which change the pion particle number, is *<sup>t</sup>π*,abs <sup>&</sup>gt; *<sup>t</sup>*fo.

In the present paper, we want to elucidate the non-equilibrium evolution of the initial fireball and the emergence of a non-equilibrium chemical potential for hadrons, in particular for the pions. The theoretical background shall be provided by the Zubarev formalism of the non-equilibrium statistical operator (NSO) [17] which introduces a generalization of the thermodynamical Gibbs ensemble by including non-equilibrium observables into the derivation of the statistical operator of the non-equilibrium, see also [18]. This is facilitated by extending the set of Lagrangian multipliers by the additional non-equilibrium chemical potentials for the hadrons will appear in the NSO. If the non-equilibrium chemical potential for the pions coincides with the pion effective mass, Bose–Einstein condensation will occur, and strong effects are expected on the measured pion spectra [19,20].

Many formulations have been given for the approach to equilibrium [21,22] and the kinetics of Bose–Einstein condensation, see for instance Refs. [23–27]. New ansätze have been developed, e.g., in References [28,29] with a state that is a fixed point and the evolution towards it is universal. Via this fixed point, the system develops then dynamically.

Such a behavior is long known in the context of the Zubarev formalism, which is able to describe, for example, the transition from the kinetic stage to the hydrodynamic stage. The short-time evolution goes to a relevant statistical operator with local time-dependent thermodynamic parameters. The long-time scale evolution is given by the time dependence of these thermodynamic parameters approaching thermodynamic equilibrium.

#### **2. The Nonequilibrium Statistical Operator Method**

HIC at ultrarelativistic energies are violent non-equilibrium processes which need, in principle, a genuine non-equilibrium approach. At present, simple approximations are used such as transport models based on the kinetic equations for single-particle distribution functions. Transport codes based on the relativistic Boltzmann–Uehling–Uhlenbeck (BUU) or Vlasov–Uehling–Uhlenbeck (VUU) equations have been worked out [23–25,30]. However, a non-equilibrium single-particle distribution is not sufficient to describe correlations in the evolving system. As example, cluster formation in an expanding fireball requires the inclusion of higher order correlation functions to describe bound states like hadrons or nuclei. Alternatively, the freeze-out concept assumes nuclear statistical equilibrium (NSE) during the expansion of the fireball, which is justified if the time *τ*therm for the relaxation to thermodynamic equilibrium is small compared to the variation time *τ*exp = *q*/*q*˙ of a parameter *q* describing the thermodynamic state of the expanding system. The treatment of thermodynamic equilibrium is able to include all equilibrium correlations, in particular cluster formation. At freeze-out, *<sup>t</sup>* <sup>=</sup> *<sup>t</sup>*fo, collision processes that change the composition and the distribution die out. For *<sup>t</sup>* > *<sup>t</sup>*fo, baryon distributions evolve according to the mean-field description of the expansion. Note that, at least under conditions of the LHC and highest RHIC energies for soft pions (with energies and momenta smaller than *mπ*), the time scales characterizing elastic and inelastic (absorptive) processes are such that *τπ*,el *<sup>t</sup>*fo <sup>&</sup>lt; *τπ*,abs. Thereby, we may speak about the evolution of *μπ*(*t*) till thermal freeze-out. Although the expanding fireball approach with time dependent temperature *T*(*t*) and chemical potentials *μc*(*t*) for all hadrons including pions is rather reasonable, it is just an approximation to a more sophisticated many-body non-equilibrium approach given below.

A systematic general approach to non-equilibrium is given by the NSO *ρ*(*t*) which is the solution of the von Neumann equation with a given initial state [17],

$$\rho(t) = \lim\_{\varepsilon \to 0} \varepsilon \int\_{-\infty}^{t} dt' e^{\varepsilon(t - t')} \mathcal{U}(t, t') \rho\_{\text{rel}}(t') \mathcal{U}^{\dagger}(t, t'), \tag{1}$$

where for a Hamiltonian *H* which is not time-dependent holds *U*(*t*, *t* ) = exp[(*i*/¯*h*)*H*(*t* − *t* )], and the limit *ε* → 0 has to be taken after the thermodynamic limit. Instead of a distribution *ρ*initial(*t*0) at an initial time *t*0, according to Zubarev, a relevant statistical operator *ρ*rel(*t* ) has been introduced that contains all relevant information of the system in the past *t* < *<sup>t</sup>*. This relevant information characterizes the state of the system in non-equilibrium and will be discussed below. The missing irrelevant correlations *ρ*irrel(*t*) = *ρ*(*t*) − *ρ*rel(*t*) are assumed to be formed dynamically by the action of the time evolution operator *U*(*t*, *t* ).

The non-equilibrium state of the system is characterized by observables *Bn* in addition to the conserved observables such as total energy and particle numbers. The relevant information is given by the averages *Bn<sup>t</sup>* <sup>=</sup> Tr{*ρ*(*t*)*Bn*} of the relevant observables *Bn*. The maximum of the information entropy *S*inf = −Tr{*ρ*rel(*t*)ln *ρ*rel(*t*)} at given averages is the generalized Gibbs distribution

$$\rho\_{\rm rel}(t) = \frac{1}{Z\_{\rm rel}(t)} e^{-\sum\_{n} F\_{n}(t) B\_{n}}; \qquad Z\_{\rm rel}(t) = \text{Tr} e^{-\sum\_{n} F\_{n}(t) B\_{n}} \tag{2}$$

where the Lagrange multipliers *Fn*(*t*) are determined by the self-consistency conditions

$$
\langle \langle B\_{\rm n} \rangle\_{\rm rel}^{t} = \langle B\_{\rm n} \rangle^{t} \tag{3}
$$

with *Bn<sup>t</sup>* rel ≡ Tr{*ρ*rel(*t*)*Bn*}. With respect to the selection of the set of relevant observables, *Bn*, it should contain all conserved quantities which cannot be changed owing to the dynamics of the system. In an ergodic system, all other correlations are produced dynamically. We can include any set of observables to the relevant ones. In particular, if we include all slowly varying observables in the set of relevant observables {*Bn*}, we can expect that a shorter time is necessary to produce the missing non-equilibrium correlations. This means that Eq. (2) already gives a good approximation for *ρ*(*t*) at finite *ε* so that memory effects are less important.

According to the NSO method, the equations of evolution (generalized kinetic equations) are obtained from

$$\frac{d}{dt}\langle B\_n\rangle^t = \lim\_{\varepsilon \to 0} \frac{i\varepsilon}{\hbar} \int\_{-\infty}^t dt' e^{\varepsilon(t'-t)} \text{Tr}\left\{\rho\_{\text{rel}}(t') e^{iH(t'-t)/\hbar} [H, B\_n] e^{iH(t-t')/\hbar}\right\},\tag{4}$$

where we inserted the time derivative of the NSO (1) and used the self-consistency conditions (3).

The correct reproduction of the relevant information in the past gives the possibility to form the irrelevant correlations very fast so that a perturbation expansion is possible. Although the expression (1) is correct for any choice of relevant observables after performing the limit *ε* → 0, an appropriate choice of the set of relevant observables {*Bn*} allows us to make expansions quickly convergent so that, for instance, the Markov approximation can be performed (We speak here about memory effects associated with dying initial correlations. There are memory effects related to processes described by diagrams with more than two vertices in the non-equilibrium Greens function technique. These effects can be neglected only for dilute gases, but, in general, they are important. They result in the famous *T*<sup>3</sup> ln *T* correction to the specific heat of 3He, cf. [31,32]. Since this correction is quite comparable (numerically) to the leading term in the specific heat (∝ *T*), one may claim that liquid 3He is a liquid with quite strong memory effects from the point of view of kinetics). In Section 4, we discuss this issue in detail. It is our main goal to show that the optimal path for the non-equilibrium evolution (i.e., a sufficient broad choice of relevant observables *Bn*) must be found to provide us with a precise description already in low orders of perturbation theory.

A correct description of the evolution is also necessary for the expanding fireball if the freeze-out approximation is used. We cannot assume that the system at freeze-out is strictly in thermodynamic equilibrium. A simple relaxation-time ansatz is not sufficient, but a more detailed description of the time evolution is necessary. An important feature of the evolution is that a perturbation expansion in powers of the interaction which is often used cannot give the formation of bound states or quantum condensates in any finite order of perturbation theory. The introduction of the relevant NSO provides us with the description of those phenomena. In particular, we show that the formation of an intermediate pion Bose–Einstein condensate can be described in our approach. It is an advantage of the Zubarev approach that initial correlations are included via the relevant statistical operator *ρ*rel(*t*), in contrast to the kinetic theory which is based on the single-particle distribution function.

#### **3. Model for Pions in Heavy-Ion Collisions at Ultra High Energies**

To explain our approach denoted as the NSO method, we discuss the pion production in heavy-ion collisions. After an initial stage where hadrons are formed, we consider a fireball consisting of pions, nucleons *N*, Δ resonances, and further hadronic states such as *N*∗ and *K* mesons containing strangeness degrees of freedom. In occupation number representation, *a*<sup>+</sup> **<sup>p</sup>**,*c*, *a***p**,*<sup>c</sup>* are the creation/annihilation operators of a particle in the quantum state **p**, *c* given by the species *c* (including spin) and momentum **p**. The Hamiltonian *H* = *H*<sup>0</sup> + *H* contains the kinetic part

$$H^0 = \sum\_{\mathbf{p}, \mathcal{c}} \sqrt{p^2 + m\_{\mathcal{c}}^2} a\_{\mathbf{p}, \mathcal{c}}^+ a\_{\mathbf{p}, \mathcal{c}} = \sum\_{\mathbf{p}, \mathcal{c}} E\_{\mathbf{p}, \mathcal{c}} a\_{\mathbf{p}, \mathcal{c}}^+ a\_{\mathbf{p}, \mathcal{c}} \tag{5}$$

and the interaction part *H* . At highest RHIC and LHC energies, the fireball is dominated by pions. Approximately, one may speak of a pure pion gas. Thus, further within our model, we consider only special processes, which concern the pion distribution, *c*, that runs over pion species. We split the interaction into an elastic part describing collisions which conserve the pion number, *H*col, and an inelastic part describing reactions, *H*reac, where the pion number is changed, *H* = *H*col + *H*reac.

Expressions for the meson–meson interaction are found in the literature [16,30]. *π* − *π* interactions at high energies are predominantly elastic implying that at low density the number of pions is effectively conserved. We assume elastic *π* − *π* scattering of the form

$$H\_{\rm col} = \frac{1}{2} \sum\_{\mathbf{p}\_1, \mathbf{p}\_2, \mathbf{p}\_1', \mathbf{p}\_2'; x, d} \lambda\_{c,d}(\mathbf{p}\_1, \mathbf{p}\_2; \mathbf{p}\_{\mathbf{p}\_1}', \mathbf{p}\_2') a\_{\mathbf{p}\_1,c}^\dagger a\_{\mathbf{p}\_2,d}^\dagger a\_{\mathbf{p}\_2',d} a\_{\mathbf{p}\_1',c} \tag{6}$$

We assume that at *<sup>t</sup>* <sup>&</sup>lt; *<sup>t</sup>π*,cfo a state overpopulated by soft pions is formed, for *<sup>t</sup>* <sup>&</sup>gt; *<sup>t</sup>π*,cfo the collisions conserve the particle number, but evolve the distribution function to a thermal equilibrium distribution with a corresponding short relaxation time *τ*col. These collision processes may happen via a virtual states such as the *ρ* and *σ* mesons or other resonances. We note that the matrix element *λπ*,*π*(**p1**, **p2**; **p <sup>1</sup>**, **p <sup>2</sup>**) of the *π* − *π* interaction can be taken in a separable form [33,34] so that the Bethe–Goldstone equation for the T-matrix of the *π* − *π* scattering in the pion gas can be solved straightforwardly, resulting in in-medium scattering phase shifts and cross sections with resonances [34–36] as well as the corresponding equation of state [37,38].

Interactions of pions with baryons can also be described as particle number conserving 2 → 2 processes with a Hamiltonian of the form (6). See, for example, the recent work [39] on the ANL-Osaka model which provides an excellent description of existing *π* − *N* scattering data. These processes contain, in particular, the Δ and *N* resonances that become very important for heavy-ion collision experiments with lower c.m. energies at SPS, SIS-18, and the future FAIR and NICA experiments. For our discussion of results from the LHC experiment in the present work, the processes involving baryons are not important and will not be treated explicitly here.

If the particle number is fixed, neglecting processes described by *H*reac, a pion gas in thermodynamic equilibrium may form a Bose–Einstein condensate at high phase space occupation densities and sufficiently low temperatures. It is possible that the expanding fireball will meet such parameter values of the pion phase space during the non-equilibrium evolution.

There are also processes which change the particle numbers of the different species which contribute to the interaction part of the Hamiltonian *<sup>H</sup>*reac. As an example, we have *<sup>π</sup>* <sup>+</sup> *<sup>π</sup>* <sup>4</sup>*<sup>π</sup>* [40] or the formation of other mesons such as *<sup>π</sup>* <sup>+</sup> *<sup>π</sup>* <sup>→</sup> *<sup>K</sup>*¯ <sup>+</sup> *<sup>K</sup>* which decay in other channels, see Lin and Ko [30]. As shown there, because of the threshold for these reactions and the small cross sections compared to the elastic collisions, the corresponding relaxation time *τ*reac to establish chemical equilibrium is large. This is a slow process not of relevance for the time scales considered here. Other reactions involving resonant correlations such as <sup>Δ</sup> *<sup>N</sup>* <sup>+</sup> *<sup>π</sup>* contribute to collision processes via virtual states and conserve the particle number, but also have a small branching ratio for a number of non-conserving processes. These reactive collisions which change the pion number are assumed to be weak in comparison to the quasielastic collisions and can be discarded for the short-time evolution, but are relevant for the evolution on long time scales to produce chemical equilibrium.

The *π* − *π* collision term was treated in Boltzmann equation calculations by Welke and Bertsch [30]. Strong pion interaction processes which change the pion number have been considered, e.g., in Reference [40].

In addition, other work has been performed using transport codes to describe the time evolution of the pion distribution function and to solve the low-*pT* enhancement puzzle. We discuss here the more general NSO approach to give an approach which goes beyond the single-particle distribution function considered in the transport codes which are based on kinetic equations.

#### **4. The Relevant Statistical Operator**

Our main point is the selection of the set of relevant observables *Bn* which determines the convergence and the accuracy of the non-equilibrium description. We discuss three examples, the Kubo case considering only conserved quantities, the kinetic theory considering the single-particle occupation numbers, and the formation of a condensate where amplitudes are added. In principle, all three choices for the set of relevant observables should give the same results if the limit *ε* → 0 is correctly performed. However, because we use perturbation expansions and Fermi's Golden rule, these approximations lead to different results.

#### *4.1. Kubo Case*

Within the NSO method, there is no prescription for the choice of the set {*Bn*} of relevant observables. Only conserved observables have to be included because their averages cannot be changed dynamically.

A minimum set of relevant observables of the pion–nucleon system (Kubo case) is the energy *H* that is conserved. The number of pions is not strictly conserved. Because the pion number is not prescribed, in equilibrium, no corresponding chemical potential *μπ* can be introduced. Formally, one takes *μπ* = 0 similar to the photon system. Thus, in the Kubo case, one supposes that *τπ*,abs *t*fo (contrary to the case studied by us in this paper). The number of baryons *Nb* is conserved. In addition, the charge has to be considered as conserved quantity. With this selection of relevant observables, we obtain from the maximum principle for the relevant entropy the grand canonical distribution

$$\rho\_{\rm rel}^{(0)}(t) = \frac{1}{Z\_{\rm rel}^{(0)}(t)} e^{-\beta(t)\left(H - \sum\_{a} \mu\_{a}(t)N\_{a}\right)}.\tag{7}$$

Because the system is expanding, the thermodynamic averages are also changing with time and also the corresponding Lagrange parameters. We can adopt the blast wave model [41–45] to describe the expansion of the fireball. If we assume that the velocity is proportional to the distance from the center, the density is decreasing with time. Assuming adiabatic expansion, the entropy is constant, but the temperature is also decreasing. This hydrodynamical model may serve as an approximation to describe the time dependence of the average density and energy and, according to the self-consistency conditions (3), *μα*(*t*) and *β*(*t*).

Starting from a non-equilibrium state, we consider the relaxation to the local thermodynamic equilibrium. The time behavior of the observables *Nc* and *H*MF *<sup>c</sup>* of the particle number of the species *c* and its mean-field energy (see below) are given by

$$\frac{d}{dt}\langle\mathcal{N}\_{\mathbf{c}}\rangle^{t} = \frac{i}{\hbar}\langle[H,\mathcal{N}\_{\mathbf{c}}]\rangle\_{\text{rel}}^{t} - \frac{1}{\hbar^{2}}\int\_{-\infty}^{0} dt' \epsilon^{\mathbf{c}\mathbf{t'}} \text{Tr}\left\{\rho\_{\text{rel}}(t)[H(\mathbf{t'}),[H,\mathcal{N}\_{\mathbf{c}}]] \right\},\tag{8}$$

$$\frac{d}{dt}\langle H\_{\rm c}^{\rm MF}\rangle^{t} = \frac{i}{\hbar}\langle [H, H\_{\rm c}^{\rm MF}]\rangle\_{\rm rel}^{t} - \frac{1}{\hbar^{2}}\int\_{-\infty}^{0} dt' e^{\varepsilon t'} \text{Tr}\{\rho\_{\rm rel}(t)[H(t'), [H, H\_{\rm c}^{\rm MF}]]\}\,. \tag{9}$$

*Particles* **2020**, *3*

Evaluating the correlation functions in Born approximation for the pion system, we observe a behavior different from the one that would be consistent with the assumption of *μπ* = 0. Because the particle numbers *N<sup>π</sup>* are conserved with respect to the elastic collisions *H*col, only the inelastic collisions *H*reac contribute. This makes the time derivative (8) small. In contrast, the thermalization process (9) is dominated by *H*col so that the exchange of energy between the different components *c* and the momentum states is a fast process.

We can calculate the corresponding relaxation times *τ*−<sup>1</sup> *<sup>i</sup>* <sup>=</sup> <sup>−</sup>[*dBi<sup>t</sup>* /*dt*]/*Bi<sup>t</sup>* for the observables of local thermodynamic equilibrium and compare it with the expansion time scale *<sup>τ</sup>*−<sup>1</sup> exp <sup>=</sup> *∂μu<sup>μ</sup>* ≈ −*n*−<sup>1</sup> *<sup>b</sup>* [*dnb*/*dt*]. The freeze-out time *t*fo is given by the condition that the increasing *τi*(*t*) becomes equal to *τ*exp.

It is evident that this is a very global approach. We cannot assume that at any freeze-out time *t*fo the system is well approximated by the equilibrium distribution (7). A more detailed description of the non-equilibrium state is necessary, in particular if there exist long-living correlations. Indeed, the relaxation to thermodynamic equilibrium (7) also implies the achievement of the total pion number in equilibrium which is determined only by *T* and *μπ* = 0 in thermodynamic equilibrium. Because the processes which change the pion number are weak under conditions at RHIC and LHC which we focus on, the corresponding relaxation times are long and the Kubo case is not valid for our considerations. To have an appropriate description of the non-equilibrium process, these slow modes should be included in *ρ*rel(*t*). At freeze-out, we therefore expect that not the thermodynamic equilibrium but a more general non-equilibrium distribution is seen.

#### *4.2. Pion Number as a Relevant Observable*

Long-living correlations have to be implemented in the set of relevant observables to improve the convergence of the perturbation expansion, and to apply the Markov approximation. For the pion system considered here, we have elastic collisions which conserve the particle numbers and in general inelastic reactions where the particle numbers of the constituents *c* are changed. Because the conserving interaction *H*col leads to cross sections which are large compared with cross sections for the non-conserving interaction *H*reac, the pion number *N<sup>π</sup>* is an observable which changes slowly with time and should be included in the set of relevant observables so that the index *α* in Equation (7) goes over all pion species *c*. The new condition

$$
\langle N\_{\pi} \rangle\_{\text{rel}}^{t} = \langle N\_{\pi} \rangle^{t} \tag{10}
$$

is not given by the external condition of the expanding fireball but must be calculated self-consistently solving the corresponding equation of evolution (8).

A new feature of the relevant distribution including the pion number is the possibility of a singularity when the self-consistency conditions (3) are solved. As well known from the ideal Bose gas, we have to treat the occupation of the ground state separately so that below a critical temperature we have *Nπ* <sup>=</sup> *N*norm *<sup>π</sup>* <sup>+</sup> *N*cond *<sup>π</sup>* with the normal component *<sup>N</sup>*norm *<sup>π</sup>* <sup>=</sup> <sup>∑</sup>*p*><sup>0</sup> *<sup>a</sup>*<sup>+</sup> **<sup>p</sup>**,*πa***p**,*<sup>π</sup>* and the condensate *N*cond *<sup>π</sup>* = *<sup>a</sup>*<sup>+</sup> 0,*πa*0,*π*. The corresponding relevant operator reads

$$\rho\_{\rm rel}^{\rm Rease}(t) = \frac{1}{Z\_{\rm rel}^{\rm Rease}(t)} e^{-\beta(t)[H - \sum\_{\ell} \mu\_{\ell}(t)N\_{\ell}^{\rm norm}] - F\_{0,\rm r}(t)a\_{0,\rm r}^{+}a\_{0,\rm r}}.\tag{11}$$

The new Lagrange parameter *F*0,*π*(*t*) follows from the self-consistency relation (10) in perturbation expansion with respect to *H* as

$$\langle N\_{\pi}^{\text{cond}} \rangle\_{\text{rel}}^{t} = \frac{1}{e^{\beta(t)[E\_{0,\pi} - \mu\_{\pi}(t)] + F\_{0,\pi}(t)} \dots} \, \, \tag{12}$$

which is a macroscopic number if the temperature is below the critical one. The normal component *N*norm *<sup>π</sup> <sup>t</sup>* rel is only a function of *β*(*t*) as is well known; we have *E*0,*<sup>π</sup>* − *μπ*(*t*) = 0 in the condensate state. Then, the Bose condensate *N*cond *<sup>π</sup> <sup>t</sup>* rel is a macroscopic number so that *F*0,*π*(*t*) is infinitesimally small. Below, we improve this shortcoming introducing coherent states.

We can also solve the evolution equations (8) with the relevant statistical operator (11). The corresponding relaxation times are given by the interaction part *H* . However, they are very different. As before, the elastic collisions thermalize the kinetic energies of the components (9) leading to the relaxation time *τ*therm.

To describe chemical relaxation (8), we can simplify the von Neumann equation as

$$\frac{\partial}{\partial t}\rho(t) = \frac{1}{i\hbar}\left[ (H^0 + H\_{\text{reca}})\_\prime \rho(t) \right] - \frac{1}{\tau\_{\text{therm}}} \left( \rho(t) - \rho\_{\text{rel}}^{\text{Rose}}(t) \right). \tag{13}$$

This is possible if the thermalization is very fast compared to the formation of the chemical equilibrium. The solution is given by Equation (1) after replacing  by 1/*τ*therm and *H* by *H*reac.

With respect to the evolution of the fireball, we have the result that full thermodynamic equilibrium must not necessarily occur at the freeze-out time. The reaction rates become small during the expansion so that the relevant distribution (11) at *t*fo is seen. In contrast to the thermodynamic equilibrium (7), a Bose–Einstein condensate of pions is possible.

#### *4.3. Kinetic Equations*

We can further improve the relevant statistical operator considering the occupation numbers *n***p**,*<sup>c</sup>* = *a*<sup>+</sup> **<sup>p</sup>**,*ca***p**,*<sup>c</sup>* of the single-particle states as relevant observables. Formally, instead of *Nc*, each single-particle state is taken similar to a new species. The time dependence of the mean occupation of this state leads to the kinetic equations. For details of the derivation, see, e.g., Reference [46], Equation (4.100).

We consider the time evolution of the pion distribution function *n***p**,*c<sup>t</sup>* as the diagonal part of the Wigner distribution function [46]. The relevant statistical operator has the form

$$\rho\_{\rm rel}^{\rm kin}(t) = \frac{1}{Z\_{\rm rel}^{\rm kin}(t)} e^{-\sum\_{\mathbf{p},\epsilon} s\_{\mathbf{p},\epsilon}(t) a\_{\mathbf{p},\epsilon}^{+} a\_{\mathbf{p},\epsilon}} \tag{14}$$

with the corresponding self-consistency conditions to eliminate the Lagrange parameters *s***p**,*c*(*t*),

$$
\langle n\_{\mathbf{p},\varepsilon} \rangle^t = \frac{1}{e^{s\_{\mathbf{p},\varepsilon}(t)} - 1}. \tag{15}
$$

For the time evolution, the following kinetic equation is obtained from (4) after integration by parts and using (3), see also [46]

$$\frac{d}{dt}\langle n\_{\mathbf{p},c}\rangle^t = \frac{1}{\hbar^2} \int\_{-\infty}^0 dt' e^{\mathbf{t'}t'} \text{Tr}\left\{ [H\_\prime n\_{\mathbf{p},c}] e^{(i/\hbar)Ht'} [H\_\prime \rho\_{\text{rel}}^{\text{kin}}(t)] e^{-(i/\hbar)Ht'} \right\} \tag{16}$$

if we neglect the explicit time dependence of *ρ*kin rel (*t*). In the approximation of binary collisions, we get the quantum statistical Boltzmann equation

$$\begin{array}{rcl} \frac{d}{dt} \langle u\_{\mathbb{P}^1} \rangle\_{\mathrm{coll}}^{\mathbb{I}} &=& \frac{2\pi}{\hbar} \sum\_{\mathbf{p}\_1 \mathbf{p}\_1', \mathbf{p}\_2'} \delta \left( E\_{\mathbf{p}\_1} + E\_{\mathbf{p}\_2} - E\_{\mathbf{p}\_1'} - E\_{\mathbf{p}\_2'} \right) \delta\_{\mathbf{p}\_1 + \mathbf{p}\_2 - \mathbf{p}\_1' - \mathbf{p}\_2'} |t(\mathbf{p}\_1 \mathbf{p}\_2, \mathbf{p}\_1' \mathbf{p}\_2') + t(\mathbf{p}\_1 \mathbf{p}\_2, \mathbf{p}\_2' \mathbf{p}\_1')|^2 \\ & \times \left\{ \langle n\_{\mathbf{p}\_1'} \rangle^{\mathbb{I}} \langle n\_{\mathbf{p}\_2'} \rangle^{\mathbb{I}} [1 + \langle n\_{\mathbf{p}\_1} \rangle^{\mathbb{I}}] [1 + \langle n\_{\mathbf{p}\_2} \rangle^{\mathbb{I}}] - \langle n\_{\mathbf{p}\_1} \rangle^{\mathbb{I}} \langle n\_{\mathbf{p}\_2} \rangle^{\mathbb{I}} [1 + \langle n\_{\mathbf{p}\_1'} \rangle^{\mathbb{I}}] [1 + \langle n\_{\mathbf{p}\_1'} \rangle^{\mathbb{I}}] \right\} \end{array} \tag{17}$$

where the two-particle T-matrix is given by the interaction potential in Born approximation [46]. Near thermodynamic equilibrium,

$$
\langle n\_{\mathbf{p},\varepsilon} \rangle\_{\text{eq}} = \frac{1}{e^{E\_{\mathbf{p},\varepsilon}/T - \mu\_{\varepsilon}/T} - 1} \,\,\,\,\tag{18}
$$

we can approximate the Boltzmann equation in relaxation time approximation as:

$$\frac{d}{dt}\langle n\_{\mathbf{p},c}\rangle^t = -\frac{1}{\tau\_{\mathbf{p},c}}(\langle n\_{\mathbf{p},c}\rangle^t - \langle n\_{\mathbf{p},c}\rangle\_{\text{eq}}) \tag{19}$$

where the relaxation time *τ***p**,*<sup>c</sup>* is calculated from a microscopic collision process. This approach of a relaxation time *τ*col of collisions as the thermal average over *τ***p**,*<sup>c</sup>* [46] or the corresponding collision frequency is used in the relevant literature (see, e.g., [30]). Note that the relaxation time ansatz (19) with **p**-dependent relaxation time does not obey, in general, the conservation of particle number. According to Mermin [47], this defect is removed if the relaxation occurs not to the equilibrium state but to a relevant operator which accounts for the conservation of particle number, see also [48].

Thermal freeze-out is obtained at the time *t*fo when *τ*col(*t*fo) = *τ*exp(*t*fo), where the scattering time scale for a given particle species *c*, working in favor of equilibration, can be computed locally from the local densities *nd*, thermal (relative) velocities *vcd*, and total scattering cross sections *σcd* between the particles *c* and *d* [42–44] after momentum average

$$\frac{1}{\pi\_{\rm col}} = \sum\_{d} \langle v\_{\rm cd} v\_{\rm cd} \rangle n\_d \,. \tag{20}$$

#### *4.4. Nonequilibrium State with Condensate Formation*

Alternatively, we can construct another relevant operator with the single-particle occupation numbers *n***p**,*c*, but containing also non-diagonal parts and, in addition, also single construction operators *a***p**,*<sup>c</sup>* and *a*<sup>+</sup> **<sup>p</sup>**,*c*. Corresponding expressions are known from the theory of coherent states which are of interest to describe Bose–Einstein condensates. We can construct the relevant entropy with arbitrary powers of the creation and annihilation operators, corresponding to a very general expansion of the entropy operator in occupation number representation.

As a simple case, we construct the relevant statistical operator

$$\rho\_{\rm rel}^{\rm coh}(t) = \frac{1}{Z\_{\rm rel}^{\rm coh}(t)} e^{\sum\_{\mathbf{p},\varepsilon} \left[ F\_{\mathbf{p},\varepsilon}^{\circ}(t) a\_{\mathbf{p},\varepsilon} + F\_{\mathbf{p},\varepsilon}(t) a\_{\mathbf{p},\varepsilon}^{+} - s\_{\mathbf{p},\varepsilon}(t) a\_{\mathbf{p},\varepsilon}^{+} a\_{\mathbf{p},\varepsilon} \right]} = e^{-S^{(2)}(t)} \tag{21}$$

with the corresponding expression for the partition function *Z*coh rel (*t*). Higher order contributions such as *a*<sup>+</sup> **p**,*ca*<sup>+</sup> **<sup>p</sup>**,*<sup>c</sup>* are also possible, as well as non-diagonal terms (describing systems which are not homogeneous in space) but will not be considered here. A similar approach has been used for superfluidity in strongly coupled fermion systems [49]. Note that this bilinear form of the entropy *S*(2)(*t*) may be extended including higher than second order terms in *a*<sup>+</sup> **<sup>p</sup>**,*<sup>c</sup>* and *a***p**,*c*. This is necessary to describe, e.g., total energy conservation or the formation of bound states.

We have to eliminate the Lagrange multipliers *F*∗ **<sup>p</sup>**,*c*(*t*), *F***p**,*c*(*t*),*s***p**,*c*(*t*) using the self-consistency conditions (3). The evaluation of correlation functions becomes quite simple if the relevant statistical operator (7) is diagonal in the occupation number representation. We transform *a***p**,*<sup>c</sup>* = *b***p**,*<sup>c</sup>* + *B***p**,*c*(*t*) where *b***p**,*<sup>c</sup>* obey the usual commutation relations for bosons, and *B***p**,*c*(*t*) = *F***p**,*c*(*t*)/*s***p**,*c*(*t*) is a c-number. We obtain the diagonal form

$$S^{(2)}(t) = \sum\_{\mathbf{p}, \mathcal{L}} [s\_{\mathbf{p}, \mathcal{L}}(t) b\_{\mathbf{p}, \mathcal{L}}^{+} b\_{\mathbf{p}, \mathcal{L}} - |F\_{\mathbf{p}, \mathcal{L}}(t)|^2 / s\_{\mathbf{p}, \mathcal{L}}(t)] \tag{22}$$

for the bilinear relevant entropy *S*(2)(*t*), the c-number term can be canceled with *Z*coh rel (*t*). Then, the evaluation of *n***p**,*c*rel is quite simple and yields the well-known result

$$\langle \langle b\_{\mathbf{p},\varepsilon}^{+} b\_{\mathbf{p},\varepsilon} \rangle\_{\mathrm{rel}}^{t} \rangle\_{\mathrm{rel}}^{t} = \frac{1}{e^{s\_{\mathbf{p},\varepsilon}(t)} - 1} = f\_{\mathbf{p},\varepsilon}(t), \qquad \langle b\_{\mathbf{p},\varepsilon}^{+} \rangle\_{\mathrm{rel}}^{t} = \langle b\_{\mathbf{p},\varepsilon} \rangle\_{\mathrm{rel}}^{t} = 0,\tag{23}$$

*Particles* **2020**, *3*

which is the Bose distribution for the ideal quantum gas, but with non-equilibrium parameter *s***p**,*c*(*t*) ≥ 0 which are determined by the given averages. The mean occupation numbers follow as *n***p**,*c<sup>t</sup>* rel = [*es***p**,*<sup>c</sup>* (*t*) <sup>−</sup> <sup>1</sup>] <sup>−</sup><sup>1</sup> <sup>+</sup> <sup>|</sup>*B***p**,*c*(*t*)<sup>|</sup> 2. In addition, we find *a***p**,*c<sup>t</sup>* rel = *B***p**,*c*(*t*). With these relations, the Lagrange parameters in Equation (21) can be eliminated.

As before, the dynamics of the many-particle system is described by the Hamiltonian *H* = *H*<sup>0</sup> + *H*col + *H*reac, defined in Equations (5) and (6). We extract the mean-field terms (MF) from the interaction (1 = {**p**, *c*})

$$H = \sum\_{1} E^{\rm MF}(1, t) a\_1^+ a\_1 + \frac{1}{2} \sum\_{12} \Delta\_{\rm pair}^{\rm MF}(12, t) a\_1^+ a\_2^+ + \text{c.c.} + \frac{1}{2} \sum\_{12 \text{V} 2'} V(12, 1' 2') a\_1^+ a\_2^+ a\_{2'} a\_{1'} - (\text{MF}) \tag{24}$$

with *<sup>E</sup>*MF(1, *<sup>t</sup>*) = *<sup>E</sup>*(1) + <sup>∑</sup><sup>2</sup> *<sup>V</sup>*(12, 12)ex*n*2*<sup>t</sup>* and <sup>Δ</sup>MF pair(12, *t*) = ∑12 *V*(12, 1 2 )*a*2 *<sup>a</sup>*1*<sup>t</sup>* . We will not consider pairing so that ΔMF pair(12, *t*) = 0. In the case of fermions, pairing was considered in Reference [49], which can also transformed to the bilinear form (22) applying the Bogoliubov transformation. In the case of a Bose gas considered here, the mean-field terms contain also averages *a*+ <sup>2</sup> *<sup>a</sup>*2 *<sup>a</sup>*1*<sup>t</sup>* of the condensate mode so that

$$H^{\rm MF}(t) = \sum\_{1} E^{\rm MF}(1, t) a\_1^+ a\_1 + \frac{1}{2} \Delta\_{\rm cond}^{\rm MF}(1, t) a\_1^+ + \text{c.c.} \tag{25}$$

with ΔMF cond(1, *t*) = ∑212 *V*(12, 1 2 )*a*<sup>+</sup> <sup>2</sup> *<sup>a</sup>*2 *<sup>a</sup>*1*<sup>t</sup>* .

According to the NSO method, the kinetic equations are obtained from the equations of evolution (4) for the relevant observables

$$\frac{d}{dt}\langle B\_n\rangle^t = \lim\_{\varepsilon \to 0} \frac{i\varepsilon}{\hbar} \int\_{-\infty}^t dt' e^{\varepsilon(t'-t)} \text{Tr}\left\{ \varepsilon^{-S(2)} (t') \, e^{iH(t'-t)/\hbar} [H\_\prime B\_n] \varepsilon^{-iH(t-t')/\hbar} \right\}.\tag{26}$$

We apply perturbation theory with respect to the deviation Δ*H* from the mean-field expression which can be incorporated into *s***p**,*c*(*t*) = *βc*(*t*)[*E*MF **<sup>p</sup>**,*<sup>c</sup>* (*t*) <sup>−</sup> *<sup>μ</sup>c*(*t*)] + *<sup>δ</sup> <sup>f</sup>***p**,*c*(*t*) = *<sup>f</sup>* MF **<sup>p</sup>**,*<sup>c</sup>* (*t*) + *δ f***p**,*c*(*t*). The new Lagrange parameters *βc*(*t*), *μc*(*t*) are introduced to describe the total particle number *Nc* and mean-field energy *H*MF *c* of the species *c*. The perturbation expansion is performed with respect to Δ*S*(*t*) = *S*(2) (*t*) <sup>−</sup> *<sup>S</sup>*0(*t*) with *<sup>S</sup>*0(*t*) = *<sup>β</sup>*(*t*)[*H*MF(*t*) <sup>−</sup> <sup>∑</sup>*<sup>c</sup> <sup>μ</sup>c*(*t*)*Nc*], where *<sup>β</sup>*(*t*) is determined by the average of the total energy *<sup>H</sup>*. We have (*S*(2) (*t*) commutes with *S*0(*t*) in the lowest order of perturbation theory)

$$\frac{d}{dt}\langle a\_1^+\rangle^t = \lim\_{\varepsilon \to 0} \frac{i\varepsilon}{\hbar} \int\_{-\infty}^0 dt' e^{\varepsilon t'} e^{-i\mu\_1 t'/\hbar} \text{Tr} \left\{ e^{-S^{(2)}(t'+t)} \left[ \mathcal{E}(1)a\_1^+ + \frac{1}{2} \sum\_{1'22'} V(1'2', 12)a\_1^+ a\_2^+ a\_2^+ \right] \right\} \tag{27}$$

and the corresponding equations of evolution for the other relevant observables *Nc*, *H*MF *<sup>c</sup>* , *a***p**,*c*. To evaluate the trace, we perform the transformation of the relevant statistical operator (7) to the diagonal form (22) and have

$$\frac{d}{dt}\langle a\_1^+\rangle^t = \frac{i}{\hbar} E^{\text{MF}}(1, t) \lim\_{\varepsilon \to 0} \varepsilon \int\_{-\infty}^0 dt' e^{\varepsilon t'} e^{-i\mu\_1 t'/\hbar} F\_1^\*(t + t'),\tag{28}$$

where it was assumed that the mean-field energy *E*MF(1, *t*) = *E*(1) + ∑<sup>2</sup> *V*(12, 12)ex *f*2(*t*) depends only weakly on time so that it can be extracted from the integral which is determined by the collisions. In addition, we suppose that a condensate mode is only in the state **p**1, *c*<sup>1</sup> and *V*(11, 11) = 0 with 1 denoting the state of lowest mean field energy. In the stationary state, we assume a periodic dependence on time, *F*1(*t*) = *F*<sup>0</sup> <sup>1</sup> *<sup>e</sup>iω<sup>t</sup>* . Then, the integral can be performed with the result *ω* = *μ*1/¯*h*. The amplitude *a*1*<sup>t</sup>* <sup>=</sup> *<sup>F</sup>*(*t*) depends periodically on time. We obtain the condition *<sup>h</sup>*¯ *<sup>ω</sup>* <sup>=</sup> *<sup>μ</sup>*<sup>1</sup> <sup>=</sup> *<sup>E</sup>*MF(1) for a stationary solution, considering the lowest order (mean-field approximation). Similar results hold for *a*<sup>+</sup> 1 *t* . For the other relevant observables, the time derivative vanishes in the lowest order of perturbation expansion with respect to the interaction.

To obtain the evolution of the averages, one has to consider higher orders of the interaction. It is convenient to use the following expression for the NSO obtained from (1) after integration by parts

$$\rho(t) = \rho\_{\rm rel}(t) - \lim\_{\varepsilon \to 0} \int\_{-\infty}^{t} \varepsilon^{\varepsilon(t'-t)} \mathcal{U}(t, t') \left\{ \frac{i}{\hbar} [H\_{\prime} \rho\_{\rm rel}(t')] + \frac{\partial}{\partial t'} \rho\_{\rm rel}(t'), \right\} \mathcal{U}(t', t) dt'. \tag{29}$$

so that, for the averages,

$$\begin{split} \frac{d}{dt} \langle B\_n \rangle^t &= \frac{i}{\hbar} \text{Tr} \{ \rho\_{\text{rel}}(t) [H, B\_n] \} \\ + \frac{1}{\hbar^2} \int\_{-\infty}^0 dt' e^{\varepsilon t'} \text{Tr} \left\{ [H, B\_n] e^{iHt'/\hbar} \left( [H, \rho\_{\text{rel}}(t')] + \frac{\hbar \partial}{i \partial t'} \rho\_{\text{rel}}(t') \right) e^{-iHt'/\hbar} \right\} \,. \end{split} \tag{30}$$

In Markov approximation, the time dependence of *ρ*rel(*t* ) is neglected, and we have the Boltzmann-like form of the equations of evolution; see also (8), (9), which is obtained after cyclic permutation within the trace. In our case, we cannot assume that the time dependence of *F*1(*t*) ∝ exp(*iωt*) is slow; only after the transformation to the diagonal form may the remaining *s*1(*t*) be slow.

Using Wick's theorem, the evaluation of the first term of the r.h.s. of (30) is immediately done if we transform to the diagonal form of the relevant statistical operator. The second term describes the collision between the pions and gives the relaxation to the intermediate relevant state showing the condensate distribution. The evaluation of the time dependence of occupation numbers *n***p**,*c<sup>t</sup>* coincides with the expression (17) but replacing *n***p**,*c<sup>t</sup>* by *<sup>f</sup>***p**,*c*(*t*) + <sup>|</sup>*B***p**,*c*(*t*)<sup>|</sup> 2.

The time evolution of the amplitude follows in Born approximation as

$$\begin{array}{rcl} \frac{d}{dt} \langle a\_1^+ \rangle^t &=& \frac{i}{\hbar} E\_1^{\text{MF}} B\_1^\*(t) + \frac{\pi}{2\hbar} B\_1^\*(t) \sum\_{l' 2 \mathbf{2}'} |V\_{l\mathbf{k}}(12, l' 2')|^2 \delta(E\_{\mathbf{p}\_l} + E\_{\mathbf{p}\_2} - E\_{\mathbf{p}\_l'} - E\_{\mathbf{p}\_2'}) \delta\_{\mathbf{p}\_1 + \mathbf{p}\_2 - \mathbf{p}\_1' - \mathbf{p}\_2'} \\ & \times \left\{ f\_{\mathbf{l}} f\_{\mathbf{2}'} (1 + f\_2) - f\_2 (1 + f\_{\mathbf{2}'}) (1 + f\_{\mathbf{l}'}) \right\}. \end{array} \tag{31}$$

Stationary solution is the grand canonical distribution with the chemical potential given by the pion number density. If the pion chemical potential approaches the lowest pion energy state, a quantum condensate will be formed which is described by a coherent state. The time evolution of the condensate is characterized by the collision time *τ*cond.

#### *4.5. Quantum Master Equation*

The non-equilibrium evolution of the pion system can also be treated considering it as an open system coupled to a bath. We can consider the gluon system as the bath in the stage of evolution where the pions are formed from the hot quark–gluon plasma, or we can consider the interaction between pions (e.g., via *ρ* mesons) as a bath. We also can consider the pion Bose–Einstein condensate as the relevant subsystem interacting with the normal pion gas. A quantum master equation is derived that contains a Lindblad term [17], and a solution can be performed using coherent states. We will not discuss this interesting approach here but mention that the treatment of open many-particle systems is also possible within the Zubarev NSO method, leading to quantum master equations [17]. See, for instance, Reference [50] for the treatment of heavy quarkonia kinetics in a quark–gluon plasma.

#### **5. Discussion**

We refer to the central 200 AGeV 16O+Au data of the NA35 Collaboration, see [30]. For more recent data, see also the discussion of Reference [51]. Assuming hadronization (formation time) at *t*<sup>0</sup> ∼ 1 fm/c and temperature of about 160 MeV, the equilibrium pion density at *μπ* <sup>=</sup> 0 is *<sup>n</sup>π*,normal <sup>≈</sup> 0.15 fm−3, in contrast to the observed density *<sup>n</sup>π*(*t*0) <sup>∼</sup> 1 fm−3. This motivates the consideration of a strong macroscopic occupation of the lowest momentum state, forming a pion Bose–Einstein condensate.

We use the hydrodynamical expansion under conservation of the initial entropy *S*<sup>0</sup> which determines the temperature evolution according to *s*(*T*)*V*(*τ*) = *S*<sup>0</sup> = const., where for the entropy density we use a fit to lattice QCD data from Reference [52].

For the expansion of the fireball, we can adopt a Bjorken like picture where in the first stage (*<sup>t</sup>* < 10 fm/c) we have one-dimensional expansion with *<sup>n</sup>π*(*t*) <sup>∼</sup> *<sup>n</sup>π*(*t*0)*t*0/*t*, and afterwards three-dimensional expansion where *nπ*(*t*) ∝ *t* <sup>−</sup>3. The expansion rate *<sup>τ</sup>*−<sup>1</sup> exp <sup>=</sup> <sup>|</sup>*n*˙ <sup>|</sup>/*<sup>n</sup>* <sup>∼</sup> 1/*<sup>t</sup>* drops down as 1/*t*.

The relaxation *τ*col = 1/[*σvn*(*t*)] is estimated by a thermal average of the elastic *π* − *π* cross section *σ* ≈ 23 mb [30], determined by scattering phase shift data. The thermal average *σv* ∼ 7 − −10 mb is nearly not depending on time, so that the collision time drops down proportional to 1/*t* in the first stage, but proportional to 1/*t* <sup>3</sup> in the later stage which induces the freeze-out which occurs at *<sup>t</sup>* <sup>≈</sup> 10 fm/c, where this transition occurs. The relaxation of the condensate mode *a*<sup>+</sup> <sup>1</sup> differs from that of the normal modes. A slower relaxation entails that this mode freezes out while the normal part is further thermalizing.

For a recent discussion of the chemical freeze-out in the phase diagram on the basis of a kinetic criterion, see Reference [53] and references therein.

The experimental data are well reproduced by the fit of a pion (and kaon) distribution with a non-equilibrium chemical potential as a Lagrangian multiplier in CERN SPS [14]. The thermal freeze-out process of pions and kaons at LHC conditions is characterized by just two parameters, the freeze-out time *τ*fo and the transverse size *r*max, whereby the shape of the transverse momentum spectra is described with only one parameter, *r*max/*τ*fo because the volume at freeze-out *V* = *πτ*fo*r*<sup>2</sup> max fixes the overall normalization [13]. For an excellent simultaneous fit of pion, kaon, and proton spectra in most central Pb+Pb collisions at <sup>√</sup>*<sup>s</sup>* <sup>=</sup> 2.76 TeV, including the low-momentum enhancement of pions, a non-equilibrium chemical potential of pions *μπ* = 134.9 MeV is required which is very close to the neutral pion mass *m<sup>π</sup>* = 134.98 MeV. The other parameters are *T*kin = 138 MeV, *τ*fo = 7.68 fm/c and *r*max = 11.7 fm, according to [13]. A scenario as described by the Zubarev approach to the non-equilibrium statistical operator, with a fast relaxation to an intermediate relevant operator describing a Bose–Einstein condensate of pions and the slow relaxation to full thermodynamic equilibrium seems to be realistic with these estimates of time scales.

A detailed numerical calculation within the presented approach, e.g., on the basis of a separable model Hamiltonian for *π* − *π* scattering as discussed in Section 3 is the subject of ongoing work that shall be reported in a subsequent publication.

#### **6. Conclusions**

There are different models to describe the low momentum enhancement of pions observed in HIC at SPS, RHIC, and LHC energies. We discuss this effect as a signature of a quantum condensate of the high-density pion gas. As an origin for the high phase space density of pions, one may think of an initial state in the form of a color glass condensate state (gluon saturation) which gets converted to a pion gas by particle number conserving process as described, e.g., [54,55].

After the hadronization time, a hadron gas is formed, which, under LHC conditions, mainly consists of pions. The time evolution of the fireball is governed by particle-conserving binary collisions; processes that change the pion number are weak and influence only the long-time evolution. In the oversaturated pion gas, the cross sections of pion rescattering processes are relatively large. As a consequence, the pion distribution function quickly relaxes to local thermodynamic equilibrium (here denoted as relevant distribution) which slowly evolves to full equilibrium.

The expansion of the fireball produced in HIC changes not only the parameter of the local thermodynamic equilibrium but influences also the relaxation time, and freeze-out happens if the relaxation rate becomes smaller than the expansion rate. A general description of this non-equilibrium process is given here within the Zubarev method of the non-equilibrium statistical operator. As discussed in the literature, it appears that one can capture the essence of the effect with fixed

point dynamics. The relevant statistical operator may be considered as a transient distribution proposed also recently from other approaches [28,29]. Here, we formulate this behavior using the Zubarev concept of a relevant statistical operator. The system quickly relaxes to a relevant distribution (pre-equilibrium state) which evolves slowly to equilibrium, but is frozen out at the freeze-out time. This relevant distribution at *t*fo describes the composition to be observed in the experiments.

We show different possibilities to introduce a relevant statistical operator to derive the corresponding equations of evolution of the state. Treating the binary collisions in relaxation time approximation, a quantum condensate may appear in the relevant statistical operator. After freeze-out, where this relevant distribution stops evolving further, the non-equilibrium evolution of the pion system is described by kinetic equations with initial condition at freeze-out time for the distribution function, which is approximated by the relevant statistical operator at *t*fo. To obtain a optimum description already in lowest order perturbation theory (Markov approximation), the relevant statistical operator should contain already all relevant correlations, in particular the formation of quantum condensates.

The method of the Zubarev NSO as presented here for the application to the pion production in heavy-ion collision experiments considers not only a systematic description of the collision processes but also a simultaneous treatment of the hydrodynamical evolution as well as the evolution of the condensate.

**Author Contributions:** Conceptualization, G.R. and D.B.; methodology, G.R.; validation, D.B., D.N.V. and V.G.M.; formal analysis, G.R.; investigation, D.B.; writing–original draft preparation, G.R.; writing–review and editing, D.B. and D.N.V.; funding acquisition, D.B. and G.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the MEPhI Academic Excellence Program under contract number 02.a03.21.0005.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
