*Article* **Generalizing Bogoliubov–Zubarev Theorem to Account for Pressure Fluctuations: Application to Relativistic Gas**

#### **Yuri G. Rudoy \* and Yuri P. Rybakov**

Institute for Physical Research and Technologies, RUDN University, Miklukho-Maklay Str. 6, Moscow 117198, Russia; soliton4@mail.ru

**\*** Correspondence: rudikar@mail.ru

Received: 18 January 2019; Accepted: 17 March 2019; Published: 21 March 2019

**Abstract:** The problem of pressure fluctuations in the thermal equilibrium state of some objects is discussed, its solution being suggested via generalizing the Bogoliubov–Zubarev theorem. This theorem relates the thermodynamic pressure with the Hamilton function and its derivatives describing the object in question. It is shown that unlike to other thermodynamic quantities (e.g., the energy or the volume) the pressure fluctuations are described not only by a purely thermodynamic quantity (namely, the corresponding thermodynamic susceptibility) but also by some non-thermodynamic quantities. The attempt is made to apply these results to the relativistic ideal gases, with some numerical results being valid for the limiting ultra-relativistic or high-temperature case.

**Keywords:** Gibbs equilibrium statistical mechanics; Bogoliubov's quasi-averages; pressure fluctuations; relativistic ideal gas

**PACS:** 05.70.-a; 05.30.-d; 05.40.-a

#### **1. Introduction**

The long-standing and rather non-trivial problem of calculating pressure fluctuations in the Gibbs equilibrium statistical mechanics is revised. The previous attempts are critically analyzed and it is shown that the application of the Bogoliubov's ideas gives the full and unambiguous solution to this problem. The crucial role plays the Bogoliubov's idea of quasi-average (in our case—quasi-dynamic) quantities—specifically, the pressure *P* and the dynamic compressibility Ψ. Following the Bogoliubov's idea of spontaneous symmetry breaking, we introduce the virtual conjugate field, which appears to be the singular potential ε*U* of the container impenetrable walls. The translational invariance of the Hamilton function *H* being broken, finally we consider the limiting case ε→0. General relations expressing *P* and Ψ in terms of the derivatives of *H* are presented and some examples are studied. In particular, we consider the cases when the Hamilton function can be expressed as the sum of uniform functions (in the Euler sense).

In our case the virtual conjugate field, which in the limit ε→0 breaks the translational invariance of the Hamilton function *H*, appears to be the singular potential ε*U* of the container impenetrable walls. The general relations expressing *P* and Ψ in terms of the derivatives of *H* are presented and some examples are studied—i.e., the cases of the ideal vs. non-ideal as well as those of uniform vs. non- and quasi-uniform (in the Euler sense) Hamilton function *H* describing the system (here—the gas, presumably in the classical regime).

The problem of the equilibrium pressure fluctuations is one of the oldest and most difficult problems in classical statistical mechanics. In 1902, Gibbs [1] in Ch.VII wrote down the appropriate expression, which included the quantity named by him the "dynamical compressibility" (see also Fowler [2], Hill [3], Kittel [4], Terletzky [5]), the latter one being important for the problem of thermodynamic stability (more precise definitions are given below).

Many unsuccessful attempts were undertaken for calculating this quantity even for the simplest case of 'ordinary' non-relativistic ideal gas. One can find the details in the works by Fowler [2], Wergeland [6], Münster [7], and M. Klein [8].

A 'pessimistic' point of view was expressed in [3], Ch. 4, §19 (see also [5], §69), where the calculation of the pressure fluctuations was connected with the detailed knowledge of the kind of forces acting between the gas particles and the container walls. Some attempts [2,3] to follow this route brought physically unsatisfactory—i.e., divergent—results. This appears to be unphysical because, due to Maxwell, the gas in the container should relax to the thermal equilibrium state quite independently on physical properties of the walls.

Finally, it was even claimed (see, e.g., [4], Ch. 11) that the solution to the problem of Gibbs' pressure fluctuations is generally outside the scope of the equilibrium theory, so all these failures were sometimes considered as the inconsistency of the Gibbs's approach as a whole. The situation becomes even more involved by noting that some of the physically acceptable results for the pressure fluctuations obtained earlier [6–8] refer in fact not to the Gibbs's approach itself, but to the Einstein's one, which is called "quasi-thermodynamic" by Landau and Lifshitz [9]. These approaches differ significantly by the choice of the thermodynamic variables fixed by calculation. For example, in the case of the pressure fluctuations it is the entropy in the Einstein's approach, though it is the volume in the Gibbs's ensemble approach (more details relating to these approaches are given in [10]). The goal of our paper is to show the efficiency of the Gibbs approach.

Concerning the problem of application, it is not obvious that just Gibbs approach (and thus BZ and our results) should be more useful than Einstein one. To decide this one should analyze the specific experimental situation but it is outside scope of our paper.

The way out concerning only the Gibbs approach was actually outlined in 1946 by Bogoliubov [11], who used the coordinate scale transformation in order to connect the thermodynamic pressure with the dynamical quantities—namely, with the first derivatives of the Hamilton function and the particles' pair distribution function. Later on in 1971, Zubarev [12] obtained the analogous expression for the dynamic (yet not somehow thermally-averaged) pressure as the function defined only in the phase space.

In fact, in [12] there was implicitly used the idea of quasi-averages—or, in our case, quasi-dynamic quantities—which was also formulated by Bogoliubov [13] in 1961. From the computational point of view, the Zubarev's result became possible by virtue of the explicit usage of generalized functions (in this case—the singular potential of the container walls) following the lines of Vladimirov [14].

However, it was only in 2000 one of the present authors with Sukhanov [10] succeeded to extend the Bogoliubov–Zubarev approach and obtained for the first time the general expression for the Gibbs 'dynamical compressibility' in terms of the second derivatives of the Hamilton function. This generalized form of the Bogoliubov–Zubarev theorem is valid for any reasonable kind of the kinetic energy and the interaction potential, but only the non-relativistic Maxwell gas was considered in [10] as an example.

Later on in [15–17], these results were extended to the ideal gas with any uniform (in the Euler sense) dependence of the Hamilton function upon the (quasi)particle momentum. The most general case of the non-uniform Hamilton function—i.e., that of Lorentz as well as Lorentz-violated form [16]—was considered for the classical ideal gas. In the present paper the pressure fluctuation problem is considered for the more complicated case—the non-uniform gas in the ultra-relativistic limiting case—both from dynamic as well as thermodynamic points of view. We stress that the logical and computational completeness of the Gibbs statistical mechanics, which was sometimes brought to doubt—especially in connection with the problem of the pressure fluctuations—is fully restored by means of Bogoliubov's seminal ideas.

This paper is organized as follows. In Section 2, the problem of pressure fluctuations is formulated and in Section 3 the solution to this problem is given in general form. Sections 4 and 5 are devoted to some illustrative examples, whereas Section 6 contains the main result for the thermodynamic properties of the ideal gas in the ultra-relativistic approximation. Section 7 provides the short summary. Appendix A clears some computational problems encountered in Section 6.

#### **2. Rigorous Formulation of the Pressure Fluctuations Problem**

The equilibrium pressure fluctuations <(Δ*P*(Г))2> are defined in a standard way as <(Δ*P*(Г))2> = <(*P*(Г))2> <sup>−</sup> <sup>&</sup>lt;*P*(Г)>2, where <sup>Г</sup> = {*q*,*p*} is the phase space of coordinates *<sup>q</sup>* and momenta *<sup>p</sup>* and < ... > denotes the canonical averaging for the system in the isothermal-isochoric ensemble with fixed values of the inverse temperature *β* and the volume *V*. The value of *β* is introduced by the canonical distribution function whereas *V*—by the restriction of the region of Г.

Following Gibbs [1], if the Hamilton function *H*(Г) for the dynamic system is given, then

$$\lambda < \dots > = Z^{-1}(\pounds V) \int d\Gamma \exp[-\beta H(\Gamma)] \langle \dots \dots \rangle,\\ Z(\beta\_\prime V) = \int d\Gamma \exp[-\beta H(\Gamma)] \lrcorner \Phi(\beta\_\prime V) = \ln Z(\beta\_\prime V), \tag{1}$$

where the partition function *Z*(*β*,*V*)is supposed to be finite and strictly positive, so that the Massieu—Planck thermodynamic potential *Φ*(*β*,*V*) does exist. The latter one is usually a smooth function of *β*and *V*, so there exist also the relevant thermodynamic derivatives, in particular the equilibrium (i.e., isothermal) pressure *P*(*β*,*V*) and the compressibility χ(*β*,*V*)<0

$$P(\beta, V) = (1/\beta)[\partial \Phi(\beta, V) / \partial V], \chi(\beta, V) \equiv \exists P(\beta, V) / \partial V = (1/\beta)[\partial^2 \Phi(\beta, V) / \partial V^2];\tag{2}$$

the expressions (2) being known in thermodynamics as thermic equations of state.

According to the Gibbs lemma [1] (Ch. VII, Equations (252) and (255), see also [3,5]), the equilibrium pressure fluctuations <(Δ*P*(Г))2> are given by the expression

$$\beta \cdot \left(\triangle P(\Gamma)\right)^2 \rhd = \chi(\beta, V) + \Psi(\beta, V),\\\chi(\beta, V) = \partial \lhd P(\Gamma) \rhd / \partial V,\\\Psi(\beta, V) = -\varsigma \partial P(\Gamma) / \partial V >,\tag{3}$$

or, following Gibbs and introducing the additional dynamic quantity Ψ(Г),

$$\Psi(\beta, V) = <\Psi(\Gamma)>,\ \Psi(\Gamma) = -\partial P(\Gamma)/\partial V = \partial^2 H(\Gamma)/\partial V^2, P(\Gamma) = -\partial H(\Gamma)/\partial V. \tag{4}$$

Gibbs called the quantity Ψ(Г) dynamic compressibility, but gave no example of its calculation; in general, calculation of quantities in (3) and (4) consists of two stages: firstly, the adequate definition of *P*(Г) and Ψ(Г), and secondly—their correct averaging according to (1).

Note that for the pressure P the first stage may be in fact bypassed due to the first of Equation (2) along with definitions in (1), and thus the pressure *P*(*β*,*V*) is called the 'thermodynamic' average. On the contrary, though according to (3) Ψ(*β*,*V*) also belongs to the set of Gibbs's averages, it is a 'non-thermodynamic' one because it needs a direct calculation according to (1).

Furthermore, in order to satisfy the conditions of thermodynamic stability relative to the external mechanical disturbance, it is necessary for<(Δ*P*(Г))2> to be *positive*. It requires Ψ(*β*,*V*) not only to be positive but also to exceed −χ(*β*,*V*). Hence Ψ(*β*,*V*) cannot be equal to −χ(*β*,*V*), this fact implying (with the account for (3) and (4)) that the operation < ... > is in general not permutable with the operation *∂*/*∂V*—just this circumstance is of decisive significance for further presentation.

It is worthwhile to note, that due to the Gibbs lemma the expressions analogous to (3) hold also for the equilibrium thermal fluctuations of other (thermo)dynamic quantities—e.g., the energy *H* or the generalized force *A* = −*∂H*/*∂a*. In all cases the relevant derivatives in the Gibbs lemma refer to the (thermo)dynamically conjugate variables (for *H* and *A* those are the inverse temperature *β* = 1/kB*T* and the relevant generalized parameter *a* respectively), but in the cases with *H* and *A* the terms *∂H*(Г)/*∂β* and <sup>−</sup>*∂A*/*∂<sup>a</sup>* <sup>=</sup> *<sup>∂</sup>*2*H*/*∂a*<sup>2</sup> fully disappear and thus no difficulties arise at all. Indeed, the energy *<sup>H</sup>*(Г) is a pure dynamic variable and so—by definition—does not depend upon the thermal parameter *β*, whereas variables *A* and *a* are mutually independent and enter *H*(Г) in the bilinear form (−*Aa*).

Quite a different situation takes place for the pair of relevant (thermo)dynamically conjugate variables—the pressure *P* and the volume *V*. Strictly speaking, all the derivatives of the energy *H*(Г) with respect to the volume *V*, entering the definitions (4), are identically zero by definition. Indeed, all the quantities *H*(Г), *P*(Г), and Ψ(Г) in (4) have pure dynamical origin and do not contain the kinematical parameter *V*. All these quantities are defined in the whole phase space Г for the 'free' system without any 'walls' (therefore, for *V*→∞), while the finite value *V* enters only at the final stage, that is, after the averaging procedure.

Let us recall here, that the method of quasi-averages was created by Bogoliubov [13] just in order to cope those frequently encountered problems, when the symmetry of the Hamiltonian *H* of the physical system is higher than that of the ground state or of the state of thermal equilibrium. In these cases, the formal calculations in accordance with (1) prescribed by the Gibbs approach [1] (as well as by its quantum generalization) give unphysical zero results for ordinary average values. It was shown by Bogoliubov [13], that the reason lies in the existence of some kind of degeneration in the system's energy, so the notion of quasi-averages was suggested in order to obtain physically meaningful results.

Note that the term 'degeneration' is fully deprived here of any 'quantum' sense and is used only to designate the presence of some additional symmetry in *H* (e.g., with respect to translations or rotations in the configuration part of the phase space Γ). The ingenious—though almost 'obvious'—Bogoliubov's idea was to remove this 'degeneration' by means of relevant conjugate (in wide sense) infinitesimal 'external field' before the averaging procedure is carried out and then, after all calculations are made, fully eliminate the field. Spoken figuratively, the quasi-averages are alike to such fictitious personages as the Moor of Venice or the Cheshire Cat.

To be specific, in our case the Hamilton function *H*(*p*,*q*) ≡ *H*(Γ) describing the energy of the system of particles in classical regime, is translation invariant and does not distinguish between the 'interior' and the 'exterior' of some 'container'. Therefore, *H*(Γ) does not depend upon the volume *V* of this container and thus both quantities *P* and Ψ are formally identically equal to zero. But in fact, any system in thermal equilibrium should be confined in space, so the system's energy should depend upon the value of the volume *V*—in the opposite case no pressure and no dynamic compressibility may be formally defined at all.

#### **3. Solution to the Problem of Pressure Fluctuations**

In order to overcome this contradiction and to obtain the adequate definitions of *P* and Ψ, we act in the spirit of the Bogoliubov's method [13] and, following partly to Zubarev [12], violate (may be virtually) the translational symmetry of 'free' Hamilton function *H*(Г). To this end, one can simply add to *H*(Г) the singular repulsive potential *U*V(*q*)

$$H\_V^{(\varepsilon)}(\Gamma) = H(\Gamma) + \varepsilon \mathcal{U}\_V(q); \mathcal{U}\_V(q) = \begin{cases} \ 0, q \notin \mathcal{S}\_{V\_{\varepsilon}} \\ \ \infty, q \in \mathcal{S}\_V. \end{cases} \tag{5}$$

The potential *U*V(*q*) is called also the 'contact delta-like', or the 'hard core' potential, which describes dynamically the container of the volume *V* and the surrounding surface *S*<sup>V</sup> with the idealized 'impenetrable' walls. Evidently, *U*V(*q*) should not depend on the form of any actually present 'wall–particle' interaction. Its role reduces to introducing the dependence of the *ε*-deformed Hamilton function (5) on the volume *V*.

By virtue of the suggested properties (5) of *U*V(*q*), the configuration part of Γ is divided into the 'interior' and the 'exterior' parts (relative to the container). Therefore, the potential *U*V(*q*) acquires the properties of the generalized function (in particular, see [10]). Possibly, just these circumstances have led to the failure of perturbation approaches in papers [2,6].

Taking into account these definitions, it is natural to define the quantities *P* and Ψ in the proper and unambiguous way as the 'quasi-dynamical' variables in the following 'limiting' sense

$$P\_{\mathcal{V}}(\Gamma) \equiv \lim [\neg \partial H\_{\mathcal{V}}^{\{\varepsilon\}}(\Gamma) / \partial V] , \: \mathbb{V}\_{\mathcal{V}}(\Gamma) = \lim [\partial^2 H\_{\mathcal{V}}^{\{\varepsilon\}}(\Gamma) / \partial V^2] \ (\varepsilon \to 0); \: \mathbb{V}\_{\mathcal{V}}(\Gamma) \neq -\partial P\_{\mathcal{V}}(\Gamma) / \partial V. \tag{6}$$

Note that the mathematical hallmark of the Bogoliubov's method of quasi-averages [13] consists in their non-analytic dependence upon the infinitesimal parameter ε, and this is just the main reason why the results of the limiting procedure (6) differ drastically from the identically zero results for *P*(Г) and Ψ(Г), when ε is taken equal to zero from the very beginning.

It can be shown (details of calculation see in [10], App. 6) that in accordance with the definition (6), *P*V(Г) coincides exactly with the previously known result of Zubarev [12], whereas ΨV(Г) in the form (6) was presented in [10] explicitly for the first time. It is worthwhile to note that, in [10], the quantum generalization of these results was also obtained based on the well-known Hellman—Feynman theorem for the operator's parameter differentiation.

Finally, the averaging of ΨV(Г) according to (1) gives Ψ(*β*,*V*) and thus allows one to obtain in quite general way the solution of the long standing and rather controversial problem (see [2–8]) of thermal equilibrium pressure fluctuations (3) in the isothermal-isochoric Gibbs ensemble for the non-ideal systems of particles in classical regime.

The key mathematical device for obtaining the quasi-dynamical equations of state (6) is the equality of volume derivatives of the n-th order for the two types of functionals, namely ZV(β) = <sup>d</sup>Гexp[−βH(Г)] and ZV(ε)(β) = dГexp[−βHV(ε)(Г). In the first case the integral is taken over the kinematic-confined coordinate subspace of Г with the volume V. On the contrary, in the second case, the integral is taken over the whole coordinate subspace of Г and only after this the limiting procedure ε→0 is performed.

For *n* = 0 the equality is quite obvious because the dynamical factor ΔV(Г) = exp[−β*U*V(*q*)] acts as the projection operator onto the relevant coordinate subspace of Г. Indeed, according to the definition (5) of the external wall potential *U*V(*q*), ΔV(Г) = 1, if *U*V(*q*) = 0—i.e., when *q* belongs to the interior of the container, and ΔV(Г) = 0, if *U*V(*q*)→∞, when *q* belongs to the exterior of the container or even to its walls. Details of calculations for *n* = 1 and *n* = 2, which give the constructive realization of the definitions (6), can be found in [10], App. 7, and were repeated in the recent paper [18]. The main result is the following.

Suppose that a macroscopic dynamic system is confined within the finite volume *V* and is described by the Hamilton function of the form (5). Then the explicit expressions for *PV*(Г) and Ψ*V*(Г) are determined only by the 'free' part *H*(Г) of the Hamilton function (5) and are quite independent upon the specific form of the "wall potential" *UV*(Г). This result is quite typical when one works with the generalized functions.

Furthermore, the canonical scale transformation in the phase space Г = (*q*,*p*)→Г<sup>λ</sup> = (λ*q*,*p*/λ) is performed, and thus we obtain the following expressions for *P*V(Г)and ΨV(Г) in terms of the partial derivatives of the Hamilton function *H*(Гλ) for the "free", or unconfined, system but with λ-deformed phase space (*f* is the degree of freedom)

$$\begin{aligned} P\_V(\Gamma) &= -(1/fV)[\mathcal{D}\_\lambda H(\Gamma\_\lambda)] \mid\_{\lambda=1\prime} \Psi\_V(\Gamma) \equiv (1/V)P\_V(\Gamma) + \Delta \Psi\_V(\Gamma), \\ &\Delta \Psi\_V(\Gamma) = (1/fV)^2[\mathcal{D}\_\lambda (1 + \mathcal{D}\_\lambda)H(\Gamma\_\lambda)] \mid\_{\lambda=1}. \end{aligned} \tag{7}$$

Here D<sup>λ</sup> <sup>≡</sup> d/dλ, and 1 <sup>≡</sup> <sup>D</sup>λ<sup>0</sup> is the symbolic designation of the unity operator in the operator family, {Dλn} (*<sup>n</sup>* <sup>≥</sup> 0—integer) is the *<sup>n</sup>*-fold differentiation with respect to <sup>λ</sup>. Finally, one should put everywhere λ = 1. Expressions (7) are well defined for sufficiently smooth Hamilton function *H*(*q*,*p*)—i.e., twice differentiable with respect to the arguments *p*и*q*, while this operation does not yield the dependence of *PV*(Г) and Ψ*V*(Г) upon *V*.

The auxiliary variable λ establishes the connection between the change of the volume *V* and the equivalent change of the coordinates *q*, where the condition of canonicity requires also the corresponding change of the momenta *p*. In other words, λ is a parameter of canonical scaling transformation preserving the Liouville dynamic measure—i.e., the volume element of phase space *d*Г: clearly, *d*Г<sup>λ</sup> = (λd*q*)(d*p*/λ) = (d*q*d*p*) = *d*Г.

Note that terms of different order in D<sup>λ</sup> entering ΔΨ*V*(Г) may give contributions of the same order; e.g., in the case (7) and (8) (see below) the terms in (6) take the form—e.g., for *H*k(*p*) (just the same expressions will be valid for the contribution of *H*p(*q*)into (6), with replacing *m* with *l*.):

$$[\mathrm{D}\_{\lambda}H\_{\mathrm{k}}(p/\lambda)]\mid\_{\lambda=1} = -mH\_{\mathrm{k}}(p), \;[\mathrm{D}\_{\lambda}^{\;2}H\_{\mathrm{k}}(p/\lambda)]\mid\_{\lambda=1} = m(m+1)H\_{\mathrm{k}}(p),$$

$$P\_V(p) = (1/fV)mH\_{\mathrm{k}}(p), \;\Delta\Psi\_V(p) = (1/fV)^2m^2H\_{\mathrm{k}}(p).$$

The expression for *PV*(Г) in (7) is usually cited as the Bogoliubov–Zubarev theorem [11,12], whereas the expression for Ψ*V*(Г) for the first time was obtained in the paper byRudoy and Sukhanov [10]. It is natural to call the expressions (7) (quasi)-dynamical equations of state, because they connect the (quasi)-dynamic quantities—the pressure *P*(Г) and the compressibility Ψ(Г) with the main characteristic of the dynamic system—the Hamilton function *H*(Г).

It is essential that thermodynamic equations of state (7) do not include external thermal parameter—the temperature *T*, but they explicitly depend on the external mechanical parameter—the volume *V*. It is evident that the dynamic functions *H*, *P* and Ψ are defined in the system's phase space Г. Note also that all functions entering (7) are usually (but not always!) additive, so their average values are proportional to particle's number *N*. Moreover, functions entering (7) possess various – but universal for all dynamical systems – kinds of behavior relative to deformations of the volume *V*, namely *H*(Г) = O(*V*0), *P*V(Г) = O(*V*−1), ΨV(Г) = O(*V*−2). Indeed, the external parameter *V* enters the right-hand parts of (6) only as entire negative (or zero) powers, so in the limit *V*→∞ (i.e., for the case of 'free' system) quantities *P*V(Г) and ΨV(Г) really tend to zero in full accord with (2), whereas *H*V(Г) rests in this limit invariable.

For most non-ideal macroscopic non-relativistic systems the Hamilton functions *H*(*q*,*p*) appear to have additive and separable nature in *q* and *p*, so they can be represented as sums of three terms: the constant rest energy *E*0, the kinetic energy *H*k(*p*) and the potential energy *H*p(*q*). These energies usually are also additive relative to all particles (*E*<sup>0</sup> and *H*k(*p*)) and to their pairs (*H*p(*q*). Evidently, the energy *E*<sup>0</sup> gives no contribution to the Equation (7) for the pressure *P* and the compressibility Ψ.

#### **4. Uniform Ideal and Non-Ideal Systems**

#### *4.1. Uniform Non-Ideal Case*

In [10] the particular case was considered, where both energies *Hκ*(*p*) and *H*p(*q*) are *uniform* (in the Euler's sense) functions of their arguments with exponents *m* and *l* respectively. This means that

$$H\_{\mathbb{P}}(\lambda q) = \lambda^l H\_{\mathbb{P}}(q),\ H\_{\mathbb{k}}(\lambda^{-1}p) = \lambda^{-m} H\_{\mathbb{k}}(p),\tag{8}$$

so the expressions (7) can be represented as

$$P\_{\nabla}(q,p) = (1/fV)[mH\_{\mathbf{k}}(p) - lH\_{\mathbf{p}}(q)], \,\Delta\Psi\_{\nabla}(q,p) = (1/fV)^2[m^2H\_{\mathbf{k}}(p) + l^2H\_{\mathbf{p}}(q)].\tag{9}$$

It should be remarked that the expression for *P*V(*q*,*p*) in (9) includes the quantity (−*lH*p(*q*)) = *q*F(*q*), where *F*(*q*) = −*∂H*p(*q*)/*∂q*, which is in fact the *Clausius force virial*. Therefore, after the Gibbs averaging the resulting expression is nothing else as the virial theorem. The "uniform" expressions (8) and (9) possess the following useful properties. Note, that in this approach it is not necessary to invoke the dynamical equations of motion with the additional assumptions of their stationary behavior relative to the time averaging: here we operate only with the phase space variables without using the time variable.

1. For any non-zero exponents *m* and *l* in (8) both energies *H*k(*p*) and *H*p(*q*) enter the right-hand parts of (8) *linearly*, every differentiation with respect to λ increasing by unity both indices *m* and *l*.

2. Physical dimension for the pressure in (9) corresponds to the energy volume density, whereas that for the dynamic compressibility being the pressure volume density, and therefore every differentiation with respect to λ increases by unity the power of the factor 1/*V*.

3. There exist conditions when the pressure *P*V(*q*,*p*) as well as the compressibility ΨV(*q*,*p*) are proportional to the total energy *H*(*q*,*p*) = *H*k(*p*) + *H*p(*q*). Under these conditions, according to (9), the average value <Ψ> is proportional to <*P*> and/or <*H*>; thus, <Ψ> is the usual thermodynamic average and its calculation does not amount to any additional problem. Clearly, these conditions can be realized only in two cases: at *m* = −*l* or at *l* = 0, while *m* may be arbitrary.

In [15–17,19], we concentrated on the case of an *ideal* dynamic system, where the coordinatedependent potential energy *H*p(*q*) of the inter-particle interaction vanishes, so *l* = 0. The total energy *H*(*q*,*p*) for this system is given by the sum of the constant term *E*<sup>0</sup> and of the kinetic energy *H*k(*p*) which depends only on the particle's momentum:

$$H\_{\mathcal{P}}(q) = 0,\\ H(\Gamma) \equiv H(q, p) = E\_0 + H\_{\mathcal{k}}(p). \tag{10}$$

#### *4.2. Uniform Ideal Gas*

In the case when *H*k(*p*) is a uniform function (in the Euler's sense) with the exponent *m*, the expressions (8) obtain the following simple form

$$P\_{\nabla}(p) = \mu[H\_{\mathbf{k}}(p)/V]\_{\prime} \,\Delta^{\prime} \mathbb{V}\_{\nabla}(p) = (1/V)\mu P\_{\nabla}(p) = (1/V)\mu^2[H\_{\mathbf{k}}(p)/V]\_{\prime} \,\mu \equiv m/f.$$

Note that both expressions (11) contain the constant *μ* = *m*/*f*, which is the ratio of the uniformity exponent *m* to the number of degrees of freedom *f*. The ratio *μ* characterizes the given dynamic system in the course of its dynamic (and also thermodynamic) description in both classic and quantum regimes; thus, *μ* represents some kind of 'similarity index' and specifies the whole class of dynamic systems.

For the given values of *f* = 1,2,3 typical values of index *μ* may vary from *μ*nr = *m*nr/*f* = 2/*f* up to *μ*ur = *m*ur/*f* = 1/*f*, where the subscripts "nr" and "ur" correspond to the non- and ultra-relativistic limiting expressions for the kinetic energy *H*k(*p*)

$$H\_{\mathbf{k}}\,^{\text{nr}}(\mathbf{p}) \approx \text{(cp)}^2 / 2E\_0\,\text{(cp/}E\_0\,\text{\*1)},\\H\_{\mathbf{k}}\,^{\text{nr}}(\mathbf{p}) \approx \text{cp}\,\text{ (cp/}E\_0\,\text{\*1)}.\tag{11}$$

Note that in the particular case of massless particles (e.g., photons) with *E*<sup>0</sup> = 0 the expression for *H*<sup>k</sup> ur(*p*) becomes exact. Obviously, for both limiting cases in (12) the kinetic energy has the form

$$H\_{\rm k}(p) = \mathfrak{a}\_{\rm m} p^{\rm m}, \; m\_{\rm ur} = 1, \; \mathfrak{a}\_{1} \equiv \mathfrak{a}\_{\rm ur} = \mathfrak{c}; \; m\_{\rm ur} = 2, \; \mathfrak{a}\_{2} \equiv \mathfrak{a}\_{\rm nr} = (\mathfrak{a}\_{1})^{2} / 2E\_{0}. \tag{12}$$

which is the exponential—and thus uniform (in the Euler's sense)—function of the momentum *p* with the uniformity exponent *m* equal to 2 and 1, respectively.

In more general situations, for any possible values 1 <sup>≤</sup> *<sup>f</sup>* <sup>≤</sup> 3 and 1 <sup>≤</sup> *<sup>m</sup>* <sup>≤</sup> 2, one obtains 1/3 <sup>≤</sup> *<sup>μ</sup>* <sup>≤</sup> 2, but in some models of the 'ideal gas' (e.g., used in modern cosmology) the ranges of the parameters *m*, *f* and *μ = m*/*f* may differ in magnitude (and sometimes also in sign); nevertheless, the expressions (11) preserve the applicability for these cases too.

Note that if the energy density is positive, the pressure fluctuations are also positive for any sign of *μ*. This fact means that the system may be mechanically stable (ΔΨ *>* 0) even if the pressure is negative (*P* < 0), and this is just the case (if *μ* < 0) due to the unusual value *m* < 0 (e.g., for the Chaplygin gas). As is easily seen, the condition *f* > 0 is always fulfilled by definition.

#### **5. Non-Uniform and Quasi-Uniform Ideal Gas**

#### *5.1. Non-Uniform Ideal Gas*

Rather general case of the ideal gas is that of the free isotropic relativistic particles with *the non-uniform* Hamilton function *H*(*p*) consisting of the rest energy *E*<sup>0</sup> ≡ *H*(0) and the kinetic energy *H*k(*p*), where *H*k(0) = 0. Note that study of this model from the point of view of statistical mechanics was started by Jüttner [20,21] and Glaser [22] many years ago; however, the problem of pressure fluctuations was not even mentioned in these papers.

The expression for *H*(*p*) is given by the Lorentz–Einstein equation

$$H(p) \equiv E\_0 + H\_\mathbf{k}(p) = [E\_0^\ 2 + (cp)^2]^{\frac{1}{2}},\\ H(p) = E\_0 h(p), \\ h(p) = 1 + h\_\mathbf{k}(p);\tag{13}$$

which can be rewritten in the dimensionless form

$$h(\boldsymbol{\xi}) = 1 + h\_{\mathbf{k}}(\boldsymbol{\xi}) = (1 + \boldsymbol{\xi}^2)^{\frac{1}{2}}; \\ h\_{\mathbf{k}}(\boldsymbol{\xi}) = h(\boldsymbol{\xi}) f \mu^{\{-\}}(\boldsymbol{\xi}), \ \boldsymbol{\xi} = \boldsymbol{cp} / E\_0 \ (E\_0 \neq 0). \tag{14}$$

Here *c* is the velocity of light in vacuum, *h* and *h*<sup>k</sup> being the dimensionless energies (total and kinetic). In the case *E*<sup>0</sup> = 0, the ultra-relativistic limit becomes an exact one: it is the uniform case with *m* = 1 (see Equation (13)).

The dynamic equations of state follow immediately from (7) but differ noticeably from (11). Using the dimensionless variable ξ = *cp*/*E*0, we obtain instead of (11) the following exact dynamical equations of state

$$P\_V(\xi) = (E\_0/fV)[[h^2(\xi) - 1]/h(\xi)] = [H(\xi)/V] \mathbf{v}^{(-)}(\xi) = [H\_k(\xi)/V] \mathbf{u}^{(+)}(\xi),\tag{15}$$

$$\Delta\Psi\_V(p) = E\_0(1/fV)^2 \langle [h^4(\xi) - 1]/h^3(\xi) \rangle = (1/V)P\_V(\xi)\mathbf{v}^{(+)}(\xi) = (1/V)[H\_k(\xi)/V]\boldsymbol{\upmu}^{(+)}(\xi)\mathbf{v}^{(+)}(\xi). \tag{16}$$

It is evident that the non-uniform expressions (14)–(16) are much more complicated as compared to their uniform counterparts (11). In particular, instead of the unique and constant 'similarity index' *μ* in (11) one obtains in (14)–(16) a whole family of the variable dimensionless factors *ν*(±) (ξ) and μ(±) (ξ). These factors have the meaning of the generalized 'similarity indices' and depend (though weakly enough) on ξ through the function *h*(ξ)

$$f\mu^{(\pm)}(\xi) = 1 \pm \left[h(\xi)\right]^{-1}, f\nu^{(\pm)}(\xi) = 1 \pm \left[h(\xi)\right]^{-2}, f\kappa^{(\pm)}(\xi) = 1 \pm \left[h(\xi)\right]^{-4};$$

$$f\mu^{(+)}(\xi)\mu^{(-)}(\xi) = \nu^{(-)}(\xi), f\nu^{(-)}(\xi)\nu^{(+)}(\xi) = f\kappa^{(-)}(\xi), h\_{\mathbf{k}}(\xi) = h(\xi)f\mu^{(-)}(\xi). \tag{17}$$

The system of exact Equations (14)–(17) is rather complicated, but in practice only their approximate forms are of real interest, namely, the two limiting cases: the non-relativistic (nr) (ξ→0) and the ultra-relativistic (ur) (ξ→∞) one. The lowest order corrections to the functions *<sup>h</sup>κ*(ξ) and 1/*hκ*(ξ), as compared to their 'uniform' analogs (12) and (13), have the form

$$h\_{\mathbf{k}}(\boldsymbol{\xi}) \approx h\_{\mathbf{k}} \,^{\text{nr}}(\boldsymbol{\xi}) [1 - \frac{1}{4} \boldsymbol{\xi}^2] = h\_{\mathbf{k}} \,^{\text{nr}}(\boldsymbol{\xi}) [1 - \frac{1}{2} h^{\text{nr}}(\boldsymbol{\xi})] \, h\_{\mathbf{k}} \,^{\text{nr}}(\boldsymbol{\xi}) = \frac{1}{2} \boldsymbol{\xi}^2 \, (\boldsymbol{\xi} \to 0), \tag{18}$$

$$\ln 1/h\_{\mathbf{k}}(\xi) = [h\_{\mathbf{k}}\,\mathrm{ur}(\xi)]^{-1} \langle 1 - \frac{1}{2}\xi^{2} \rangle = [h\_{\mathbf{k}}\,\mathrm{ur}(\xi)]^{-1} \langle 1 - \frac{1}{2}[h\_{\mathbf{k}}\,\mathrm{ur}(\xi)]^{-2} \rangle,\\ [h\_{\mathbf{k}}\,\mathrm{ur}(\xi)]^{-1} = \xi^{-1} \text{ (}\xi \to \infty\text{)}. \tag{19}$$

Note that *h*<sup>k</sup> nr(0) = 1/*h*<sup>k</sup> ur(∞) = 0, this fact enabling one to consider the quantities *hκ*(ξ) and 1/*hκ*(ξ) as small in corresponding ranges of the variable ξ.

In some physical problems there may be of interest to obtain the corrections to the limiting 'uniform' Equations (11) and (12), which are stipulated by the variable nature of the functions *κ*(±) (ξ) and *μ*(±) (ξ) entering the dynamic equations of state (16) and (17) for *P*(ξ;*V*) and Ψ(ξ;*V*). In order to

carry out some perturbation procedure at small values ξ«1 in the non-relativistic (nr) limit and at large values ξ»1 in the ultra-relativistic (ur) limit, it is convenient to use in Equations (17)–(19) as the corresponding small parameters not ξ and 1/ξ, but the quantities *hκ*(ξ) and 1/*hκ*(ξ).

Omitting simple but lengthy calculations, one obtains the following approximate results

$$P\_{\rm V}(\xi) \approx \mu\_{\rm tr} [H\_{\rm k}(\xi)/V] [1 - \frac{1}{2} h\_{\rm k}(\xi)], \Delta \Psi\_{\rm V}(\xi) \approx \mu^2 \,\_{\rm tr} [H\_{\rm k}(\xi)/V] [1 - (3/2) h\_{\rm k}(\xi)] \tag{20} (\xi \to 0);\tag{21}$$

$$P\_{\rm V}(\xi) \approx \mu\_{\rm ur} [H\_{\rm k}(\xi)/V] [1 + [h\_{\rm k}(\xi)]^{-1}], \,\Delta\Psi\_{\rm V}(\xi) \approx \mu^2 \,\_{\rm ur} [H\_{\rm k}(\xi)/V] \,\left(\xi \to \infty\right). \tag{21}$$

These expressions reveal the tendency of 'sloping' the dependence upon ξ both for kinematical (*κ*(+), μ(+)) as well as dynamical (*h*k, *P*, ΔΨ) quantities: at small (but finite) ξ all these quantities become smaller than their 'uniform' limits at ξ = 0, whereas at large (but finite) ξ, on the contrary, they become larger than their 'uniform' limits at 1/ξ = 0.

#### *5.2. Quasi-Uniform Ideal Gas*

Evidently, the most general case of ideal gas includes the Hamilton function *H*(*p*) with the non-uniform dependence upon *p*. However, in practice only certain limiting cases (e.g., non- or ultra-relativistic ones) are of interest, where *H*(*p*) (and hence also its derivatives) may be presented as an expansion in integer powers *m* of *p* with *m* > 0 or *m* < 0 (i.e., in 1/*p*), where *H*0(*p*) = *h*<sup>0</sup> ≡ *E*<sup>0</sup> = const, *m*<sup>0</sup> ≡ 0, but *m*<sup>i</sup> and *h*<sup>i</sup> at *i* = 1, 2, . . . may have both signs

$$\begin{array}{c} H(p) = \sum\_{i=0}^{n} H\_{i}(p) = \sum\_{i=0}^{n} h\_{i} p^{m\_{i}}, P\_{V}(p) = (1/fV)\sum\_{i=1}^{n} m\_{i} H\_{i}(p) = (1/fV)\sum\_{i=1}^{n} m\_{i} h\_{i} p^{m\_{i}}\\ \Psi\_{V}(p) = (1/fV)^{2} \sum\_{i=1}^{n} m\_{i}^{2} H\_{i}(p) = (1/fV)^{2} \sum\_{i=1}^{n} m\_{i}^{2} h\_{i} p^{m\_{i}}. \end{array} \tag{22}$$

Evidently, every term in (22) is a uniform one, whereas the whole expression (22) is not; so it can be considered as a quasi-uniform one and characterized not by the single uniformity exponent but by the whole discrete set of them. The examples can be found in [15–17]

$$\text{INR-limit: } m\_1 = 2, h\_1 > 0; m\_2 = 4, h\_2 = -\frac{1}{4}h\_1 < 0; \\ \text{UR-limit: } m\_1 = 1, h\_1 > 0; m\_2 = -1, h\_2 = \frac{1}{2}h\_1 > 0. \quad \text{(23)}$$

Obviously, the final sign of the quantities presented in (22) is determined by the non-trivial interplay of the coefficients *h*<sup>i</sup> and *m*i. As a rule, *h*<sup>i</sup> contains some small parameter and decreases in magnitude with increasing *i*, whereas *m*i, on the contrary, increases with *i* in magnitude.

Note that nowadays the Lorentz–Einstein expression (14) seems to be not the uniquely possible one and therefore in [16] the scheme outlined in that paper was carried out for this more general case. In particular, it appears, that in the Lorentz-violated case the UR-limit in (23) is supplemented by the third term with *h*<sup>3</sup> > 0, *m*<sup>3</sup> = 2, which has a typical NR-form. This term enters (22) due to the appearance in the Lorentz-violated case of the parameter *H*(*p*)/*E*Pl (here *E*Pl is the Planck energy) which is always small—even in the extreme UR-situation when *H*(*p*)/*E*<sup>0</sup> is large. In other words, the ratio *E*0/*E*Pl is always very small for any reasonable choice of particles constituting the system. The analysis of relevant expressions shows the existence of some critical value *p*\* defined as *cp*\*~(*E*<sup>0</sup> <sup>2</sup>*E*Pl) 1/3 . When the particle's momentum takes the value *p*\*, then the usual Lorentz behavior breaks and the velocity *v*(*p*)=d*H*(*p*)/d*p* exceeds the critical value *c* (details are given in [16]).

As it was mentioned in Section 2, the calculation of the equilibrium pressure fluctuations (3) in terms of *β* and *V* (and, may be, *N*) will be completed after averaging the quasi-dynamic quantities obtained in Sections 4–6. This procedure is much more traditional but far from being simple, so we give here only its general outline for the ideal system in the case (14) of the non-uniform kinetic Hamilton function *H*k(*p*). The partition function is of the multiplicative form

$$Z\_N(\beta, V) = [Vz(\beta)]^N, z(\beta) = \exp(-\beta E\_0) z\_k(\beta), z\_k(\beta) = \int d\Gamma(p) \exp[-\beta H\_k(p)];\tag{24}$$

$$\rho\_{\mathbf{k}}(p) = [z\_{\mathbf{k}}(\beta)]^{-1} \exp[-\beta H\_{\mathbf{k}}(p)], \\ d\Gamma(p) = A\_{\mathbf{f}} p^{f-1} dp \ (A\_1 = 1, A\_2 = 2\pi, A\_3 = 4\pi).$$

#### **6. Thermodynamic Equations of State—Relativistic Ideal Classical Gas**

Let us turn now to the derivation of the thermodynamic equations of state for the general, i.e., non-uniform, case of the relativistic ideal classical gas with the Hamilton function *h*(ξ), defined in (14). Here and below the dimensionless energetic (*h* = *H*/*E*0) and momentum (*ξ* = *cp*/*E*0) units are used, and all the extensive (i.e., proportional to the particle number *N*) quantities are given per particle.

#### *6.1. Representations for the Partition Function and Some Moments*

Using the definition (24) for the 'small' partition function *z*(*β*), we express it in the dimensionless temperature units *a* = *βE*<sup>0</sup> = *T*0/*T*, *a* ≥ 0, putting so far *E*<sup>0</sup> = 0; where*T*<sup>0</sup> ≡ *E*0/kB is the characteristic temperature and *p*<sup>0</sup> = *E*0/*c* = *T*0(kB/*c*) is the characteristic momentum for the given sort of particles

$$z(a) = \int d\Gamma\_{\xi} \exp[-ah(\xi)],\\d\Gamma\_{\xi} = d\Gamma\_{\mathbb{P}(\xi)} = A\_{\mathbb{P}(\xi)} ^ {\xi} \xi^{\frac{\xi}{\xi} - 1} d\xi,\\z(a) = \zeta(\infty \mu) - \zeta(0; \mu). \tag{25}$$

Here *A*<sup>1</sup> = 1, *A*<sup>2</sup> = 2π, *A*<sup>3</sup> = 4π, the integration limits on ξ in (25) (as well as on *p* in (24)) being equal to 0 and ∞ respectively; ζ(ξ;*a*) is the indefinite Riemann integral in the left part of (25).

Clearly, at any *a* > 0 the convergence of the integral (25) is ensured and improved with the growth of *a*, however the limiting value *a* = 0 should be excluded. Physically, the limit *a* = *T*0/*T*→0 corresponds to the high-temperature approximation *T*→∞ or to the ultra-relativistic case *E*<sup>0</sup> = *T*<sup>0</sup> = 0. Therefore, the representation (25) is convenient at large values *a*»1 (when *T*→0 and/or *T*0→∞) in order to obtain the low-temperature (LT) and/or the non-relativistic (NR) expansions, so it is natural to call it the *LT*/*NR-representation* of the partition function *z*(*a*) for the classical relativistic gas. The inclusion of the point *a* = 0 implying the high-temperature approximation can be realized through the change of variables (see Equation (26)).

The exclusion of the point *a* = 0 for the LT/NR-representation is stipulated by the fact that the quantity ζ(ξ;0) = *<sup>d</sup>*Гξ~ξ<sup>f</sup> at any *<sup>f</sup>* > 0 diverges on the upper limit at <sup>ξ</sup>→∞; the same conclusion follows from the asymptotic behavior *z*(*a*)~ *<sup>d</sup>*Гξe−*a*ξ~*a*−<sup>f</sup> at *<sup>a</sup>*→0, where the property *<sup>h</sup>*(ξ) <sup>≈</sup> <sup>ξ</sup> at large values of ξ is used. For the possibility of considering small values *a*«1, including *a* = 0 (when *T*→∞ and/or *T*0→0), i.e., to obtain the high-temperature (HT) and/or the ultra-relativistic (UR) expansions, it is necessary to go over from the *LT/NR*-*representation* to the *HT/UR-representation* for *z*(*a*).The latter one can be of interest in the case of the hot dense quark–gluon–plasma (QGP).

To this end, one should carry out in (25) the change of variable *ah*(ξ) = *η*, where *h*(ξ) = (1 + ξ2) 1 <sup>2</sup> ≥ 1, so that *η* ≥ *a*. Moreover, it is convenient to introduce the denotation *p*<sup>T</sup> = *T*(kB/*c*) for the characteristic thermal momentum of gas particles, thus obtaining

$$z(a) = \int\_{a}^{\infty} d\Gamma\_{\eta} [1 - (a/\eta)^2]^{(f-2)},\\ d\mathcal{G}\_{\eta} = A\_{\mathbf{f}}(p\_{\mathbf{T}})^{\mathbf{f}} \mathbf{e}^{-\eta} \eta^{f-1} d\eta\_{\mathbf{}}(a) = \zeta(\infty; a) - \zeta(0; a). \tag{26}$$

However, the structure in (26) contrasts with that in the integral (25), since the variable *a* enters not only into the integrand, but also into the lower limit of the integral (26). Moreover, according to (26) the quantity *z*(*a*)→0 as *a*→∞, whereas at *a* = 0 the quantity *z*(0)/*A*f(*p*T) <sup>f</sup> takes its finite limiting value equal to Г(*f*).

According to the results of Section 4, all the thermodynamic quantities of the relativistic ideal classical gas can be expressed through the ordinary (not central) moments of the partition function; these moments being defined in the following way (the quantities *zκ*(*a*) and *hκ*(ξ) all will be analogous, replaced with the *h*(*n*) (*a*) with *hκ* (*n*) (*a*))

$$h^{(n)}(a) \equiv \int\_0^\infty d\Gamma\_{\vec{\xi}}[h(\xi)]^n \exp[-ah(\xi)], \\ h^{(0)}(a) = z(a), \\ h^{(n)}(a) = \zeta^{(n)}(\infty; a) - \zeta^{(n)}(0; a). \tag{27}$$

After changing the variable <sup>ξ</sup>→<sup>η</sup> the quantities *<sup>h</sup>*(*n*) (*a*) will have the form

$$h^{\langle \eta \rangle}(a) = (1/a)^{\eta} \int d\Gamma\_{\eta} \eta^{\eta} \left[1 - (a/\eta)^{2}\right] \overline{\hat{\}}^{(f)} \,. \tag{28}$$

It is natural to call the quantities *h*(*n*) (*a*) (or *hκ* (*n*) (*a*)) the Jüttner integrals for the total and kinetic energies (in analogy with Maxwell, Bose, Fermi, and other similar integrals in statistical mechanics). Indeed, the definition of the canonical averages (9) reads <[*h*(ξ)]*n*> = *h*(*n*) (*a*)/*h*(0)(*a*), so for the caloric quantities—i.e., the internal energy and its fluctuations—we obtain immediately from their definitions

$$H(a) = E\_0[h^{(1)}(a)/h^{(0)}(a)],<\text{(\Delta H)}^2> = E\_0^{-2}[h^{(2)}(a)/h^{(0)}(a)]-[H(a)]^2.\tag{29}$$

For the thermal quantities—i.e., the pressure and its fluctuations—one obtains

$$P(a,V) = (E\_0/fV)[h^{(1)}(a) - h^{(-1)}(a)]/h^{(0)}(a) = E\_0(1/V)(1/a),\tag{30}$$

$$<\langle \Delta P \rangle^2 > = (1/\beta)\Delta \Psi(a, V) = E\_0^2 (1/a)(1/fV)^2 [h^{(1)}(a) - h^{(-3)}(a)] / h^{(0)}(a). \tag{31}$$

#### *6.2. Perturbation Expansion for the HT/UR-Representation*

Obviously, the exact expressions for the partition function in the representation (26) and the corresponding Jüttner integrals (28) are not available, and so we are not able to construct the corresponding asymptotic behavior as *a*→0. Therefore, we obtain only approximate expressions for *h*(*n*) (*a*) (*n*is any integer including zero), namely, the expansion in degrees of *hκ* (*m*) (*a*) with positive *m*. Expansions of this kind for all thermodynamic quantities arise in the limit of large values of the parameter 1/*a=T*/*T*0»1 (including the value *a* = 0 at *E*<sup>0</sup> = 0). Physically, this corresponds to the smallness of the ratio *E*0/*Hκ*(*T*), i.e., to the high temperature case (if the rest energy *E*<sup>0</sup> is fixed) or, on the contrary, to the small values of *E*<sup>0</sup> (at fixed temperature *T*). As can be seen from (19), *E*0/*H<sup>κ</sup>* ur(*T*) = *κ*ur(*T*0/*T*) = *aκ*ur, where *κ*ur = 1/*f* is the factor of the order unity.

In order to obtain the desired expansions, we use the binomial power series at ν ≥ 0

$$\mathbb{E}\left[1-(a/\eta)^2\right]^{\nu-} = \sum\_{m=0}^{\infty} \left[\mathbf{v},m\right] a^{2m} \eta^{-2m},\\ \left[\mathbf{v},m\right] = (-1)^m (2^m m!)^{-1} \prod\_{l=0}^{m-1} \left\{2\mathbf{v} - (2l-1)\right\}.\tag{32}$$

The coefficients in (32) satisfy the recurrence relation [ν,*<sup>m</sup>* + 1] = <sup>−</sup><sup>1</sup> <sup>2</sup> [ν,*m*]{2<sup>ν</sup> <sup>−</sup> (2*<sup>m</sup>* + 1)}(*<sup>m</sup>* + 1)<sup>−</sup>1, so that [*ν*,0] <sup>≡</sup> 1, [*ν*,1] = <sup>−</sup><sup>1</sup> <sup>2</sup> (2*<sup>ν</sup>* <sup>−</sup> 1), [*ν*,2] = [*ν*,1](−<sup>1</sup> <sup>4</sup> )(2*ν* − 3) = (1/8)(2*ν* − 1)(2*ν* − 3). Furthermore, it is more convenient to designate these coefficients [*f*,*m*], going over from the variable *<sup>ν</sup>* <sup>≡</sup> <sup>1</sup> <sup>2</sup> (*f* − 1) (*ν* ≥ 0) to the variable *f* = 2*ν* +1(*f* ≥ 1).

Let us substitute the expansion (32) into the integrand of the Jüttner integral (26) and introduce the special denotation for the combined exponent *k*(*m*)

$$k(m) \equiv 2\nu + n - 2m + 1 = k(0) - 2m,\\ k(0) = f + n \text{ (all } f, n, m, k \text{ being integers)} \tag{33}$$

This exponent at given values of the number of particle's degrees of freedom *f* = 1, 2, 3 as well as the order of the Jüttner integral *n* = 0, ±1, ±2, ... depends only upon the value *m* ≥ 0. Then for the *h*(*n*) (*f*;*a*) one obtains the infinite sum of the following integrals

$$h^{(\mathbf{n})}(f;a) = A\_{\mathbf{f}}(p\_{\mathbf{f}})^{\mathbf{f}}a^{-\mathbf{n}}\sum\_{m=0}^{\infty}[f,m]a^{2\mathbf{m}}\Gamma[k(m);a],\\\Gamma[k(m);a] \equiv \int\_{a}^{\infty}d\eta \mathbf{e}^{-\eta}\eta^{\mathbf{k}(m)-1}.\tag{34}$$

The quantity Г[*k*(*m*);*a*] is the incomplete gamma-function (related to the integral exponential function, see, e.g., [23]), and its expansion into the power series in *a* (at fixed value of *k*(*m*)) depends significantly upon the sign of *k*(*m*). According to the definition (33), the quantity *k*(*m*) decreases linearly with the increase of *m* and changes at the critical value *m* = *m*0, where

$$m\_0 = \frac{1}{2}k(0) \text{ (k(0) > 0 \text{ even)}, } m\_0 = \frac{1}{2}(k(0) + 1) \text{ (k(0) > 0 \text{ odd)}, m\_0 = 0 \text{ (k(0) \le 0)}.\tag{35}$$

Therefore, the infinite sum (34) is appropriate to be presented in the following form

$$h^{(n)}(f;a) = A\_{\mathbf{f}}(p\_{\Gamma})^{\mathbf{f}}a^{-n}\left\{\sum\_{0}^{m0-1}[f,m]\,a^{2m}\Gamma[k(m)>0;a]+\sum\_{m=m0}^{\infty}[f,m]\,a^{2m}\Gamma[k(m)\le0;a]\right\},\tag{36}$$

where the desired power expansions in *a* for Г[*k*(*m*);*a*] at *k*(*m*) > 0 and *k*(*m*) ≤ 0 are qualitatively different and should be considered separately (all the definitions are given in Appendix A, Equations (A5)–(A7)).

Finally, the expression (36) for *h*(*n*) (*f*;*a*) with the account for only lowest corrections in degrees of *a* may be written in the following form (recall that *p*<sup>T</sup> = *T*(*k*B/*c*), *a* = *T*0/*T* = *E*0/*k*B*T*)

$$h^{\{\eta\}}(f;\mu) = A\_{\mathbf{f}}(k\_{\mathbf{B}}/c)^{\mathbf{f}}T^{\mathbf{f}}a^{-n}\langle\Sigma(f;\nu;\mu) + S(f;\nu;0)a^{\mathbf{f}}\rangle \,. \tag{37}$$

Clearly, however, that if the values of the parameters *f* and *n* (just their sum defines *k*(0) in (33)) are such that *k*(0) < 0 and *m*<sup>0</sup> = 0, then the first summand on the right-hand side of (36) vanishes. In the second summand, pole divergences arise in of the form *a*−|*k*(0)| = (*T*/*T*0) *<sup>|</sup>*f+n| which are now not compensated, so the corresponding Jüttner integral *h*(n)(*f*;*a*) exists only at finite values of *a*, the same being valid for *T*0.

Using this fact, let us consider qualitatively the problem of thermodynamic stability of the so-called *Wien gas*, or the ideal gas of massless particles (*E*<sup>0</sup> = *k*B*T*<sup>0</sup> = 0, *a* = 0), in the context of its dependence upon the dimension *f*. In order to ensure such a stability, it is necessary that in the limit *a* = 0 the corresponding Jüttner integrals *h*(*n*) (*f*;*a*) should exist, since according to Section 4 they determine the main thermodynamic quantities and their fluctuations.

Recall that for the partition function the similarity index *n* = 0, for the average energy *n* = 1 and for its fluctuation (specific heat) *n* = 2. However, for the average pressure it is necessary to choose the values *n* = 1 and *n* = −1, whereas for the pressure fluctuations (compressibility)—the values *n* = 1 и *n* = −3. Note that in the structure of the perturbation theory expansions there appear two specific dependences: upon the dimensionality *f* (i.e., upon the particle's spatial degrees of freedom) as well as upon the order *n* of the moment (i.e., the average value of the *n*-th power of the particle's energy).

That is why in the HT/UR-representation it is impossible to write down general expressions for the coefficients of expansions in (34). Thus, one should enumerate all the terms, considering consequently different combinations of integer values of *f* and *n*.

Note that all the thermodynamic parameters and their fluctuations are determined by the dimensionless quantity χ(*n*) (*f*;*a*) = *h*(*n*) (*f*;*a*)/*h*(0)(*f*;*a*). Indeed, we obtain for the average energy, the specific heat, the pressure and the compressibility the following expressions

$$H(f; \boldsymbol{\mu}) = \mathbb{E}\_0 \chi^{(1)}(f; \boldsymbol{\mu}), \mathbb{C}\_V(f; \boldsymbol{\mu}) = \mathbb{E}\_0^{-2} \{ \chi^{(2)}(f; \boldsymbol{\mu}) - \|\chi^{(1)}(f; \boldsymbol{\mu}) \|^2 \},$$

$$P(f; \boldsymbol{\mu}) = (\mathbb{E}\_0 / fV)[\chi^{(1)}(f; \boldsymbol{\mu}) - \chi^{(-1)}(f; \boldsymbol{\mu})] = k\_\mathbb{B} T / V,\tag{38}$$

$$\Delta \Psi(f; \boldsymbol{\mu}) = (\mathbb{E}\_0 / fV)^2 (1/a)[\chi^{(1)}(f; \boldsymbol{\mu}) - \chi^{(-3)}(f; \boldsymbol{\mu})].$$

$$\chi^{(n)}(f; \boldsymbol{\mu}) = a^{-n}[\Sigma(f; n; a) + S(f; n)a^{f + n}] / [\Sigma(f; 0; a) + S(f; 0)a^f].$$

The quantity χ(*n*) (*f*;*a*) below should be approximated in the spirit of the perturbation theory with the accepted accuracy in *a* in the following way

$$\chi^{(n)}(f;\mathfrak{a}) - \chi^{(n)}(f;\mathfrak{a}) = a^{-n} \left( \left[ \Sigma(f;\mathfrak{a}\mathfrak{a}) + S(f;\mathfrak{a})a^{f+n} \right] - a^{n-n} \left[ \Sigma(f;\mathfrak{a}';\mathfrak{a}) + S(f;\mathfrak{a}')a^{f+n} \right] \right) / \left[ \Sigma(f;\mathfrak{a}\mathfrak{a}) + S(f;\mathfrak{a}')a^{f+n} \right] = 0$$

while for the cases *n* = 1 at *n*' = −1 and *n*' = −3 the contribution of the second summand in braces contains additional small factors *a*<sup>2</sup> and *a*4. Note that for the quantities *H* and *C*<sup>V</sup> this accuracy (within the scope of the applied here direct moments method) can be superfluous in comparison to the usual method, when the quantities *H* and *C*<sup>V</sup> are expressed through the derivatives with respect to *a* of the function ln*h*(0)(*f*;*a*). In the latter case the final accuracy is confined by that of calculating the partition function *h*(0)(*f*;*a*).

Consider now the HT/UR-expansions for the thermodynamic quantities of the ideal gas of particles with various numbers *f* of degrees of freedom. We start with the dimensionless moments *h* \*(*n*)(*f*;*a*) <sup>≡</sup> *<sup>h</sup>*(*n*) (*f*;*a*)/*A*f(*k*B/*c*) *f Tf* , which determine these thermodynamic quantities according to (38)

*n* **= 0**. *h̃* (0)(1;*a*) = 1 + *S*(1,0)*a*; *h̃* (0)(2;*a*) = 1 + *S*(2,0)*a*2; *h̃* (0)(32;*a*) = 2 —½*a*2 + *S*(3,0)*a*3.

*n* **= 1.** *h̃* (1)(1;*a*)*a =* 1 + *S*(1,1)*a*2; *h̃* (1)(2;*a*)*a* = 2 + *S*(2,1)*a*3; *h̃* (1)(3;*a*)*a* = 6 —½*a*2 + *S*(3,1)*a*4.


$$H = -\mathbf{3} . \qquad h^{\vec{\tau} - 3}(f\_{\prime}^{\prime}a) = S(f\_{\prime}^{\prime} - \mathbf{3})a\prime.$$

Then we write down the quantities (38), with accounting for the lowest (in *a*) correction terms. Average energy *H*(*f*;*a*); *Hκ* ur(*T*) = *H*(*f*;0) = *fk*B*T*.

$$H(1;a) = k a T \{1 - S(1,0) + [S(1,0) + S(1,-1)] \} a^2; \ H(2;a) = 2k a T \{1 + S(2,0) a^2\};\tag{39}$$

$$H(3;a) = 3k a T \{1 + \mathbb{W}a^2\}.$$

Specific heat *C*V(*f*;*a*); *C***V**ur = *C*V(*f*;0) = *fk*B.

$$\text{Cov}(1; \mu) = k\mu \{ 1 - [\text{CS}(1, 0) + 2S(1, 1)]a^2 \} \text{. Cv}(2; \mu) = 2k\mu \{ 1 + S(2, 0)a^2 \} \text{. Cv}(2; \mu) = 3k\mu \{ 1 - \text{Vs} \cdot a^2 \} \tag{40}$$

Note that due to the multiplication of *a*−<sup>1</sup> by the 'small' factor *E*<sup>0</sup> the quantity *Hur*(*T*) = *H*(*f*;0) proves to be not more 'large' and coincides with the first of the expressions (33) for the average (kinetic) energy of the UR *Wien gas* (*k* = 1, κur = 1/*f*). The correction within the second order of smallness in *a* = *T*0/*T*«1 for the expression (39) is stipulated by the account for the corresponding correction in ξ to 1/*h*(ξ) in the second of the expressions (26). This correction in (39) is positive, and physically it corresponds to the increase of the average energy with that of the rest energy.

Analogously, the correction to the specific heat in (40) at *f* = 2 and 3 is also within the second order of smallness in *a* and differs from the corresponding correction to the average energy only in sign. At *f* = 1 this tendency also takes place (because *S*(1,0) < 0), but the connection between the coefficients looks more intricate due to the fact that the linear in *a* correction to *C*<sup>V</sup> disappears. One can easily see that this property always takes place and does not depend upon the specific value of the linear in *a* term for the average energy. Naturally, *C***V**ur coincides with the second expression in (23) for the specific heat of the UR Wien gas with *k* = 1 and does not depend upon the temperature *T*. It is worth-while to note that in this case, just as before, both *C*V(*T*) and d*C*V(*T*)/d*T* are positive, so the thermodynamic stability is guaranteed.

Pressure vs. temperature *P*(*f*;*a*) = *k*B*T*/*V* at all *f* and *a*.

Pressure vs. kinetic energy *P*(*f*;*Hκ*); *P*ur(*f*;*H<sup>κ</sup>* ur)=(*κ*ur/*V*)*H<sup>κ</sup>* ur.

$$P(1;H\_k) = 1 - S(1,-1) \text{[E}\_0/H\_\text{x}\text{ur]};\\P(2;H\_k) = 1 - 2 \text{[E}\_0/H\_\text{x}\text{ur]};\\P(3;H\_k) = 1 - (3/2) \text{[E}\_0/H\_\text{x}\text{ur]}\tag{41}$$

Compressibility ΔΨ(*f*;*a*) is equal

$$
\Delta\Psi(f;\mu) = \Delta\Psi^{\rm sur}\{1 - S(f,-\Im)\mu^f\}.\tag{42}
$$

Formally, the correction for the *κ*ur = 1/*f* can be obtained by using the second of the Equation (39), with the result of the lowest order *h* (−1)(*f*;*a*) <sup>≈</sup> *<sup>a</sup>*. Taking into account that the small parameter reads *a* = *T*0/*T* = (1/*κ*ur)(*E*0/*H<sup>κ</sup>* ur) = *f*(*E*0/*H<sup>κ</sup>* ue), one obtains (41). Finally, the corrections to the limiting

UR-value of the compressibility ΔΨur = ΔΨ(*f*;0) = (*κ*ur/*V*)*P*ur at *a* = 0 start with *af* and in full analogy with the pressure they are negative.

Therefore, it is obvious that the HT/UR-corrections do not violate the thermodynamic stability of the system because these corrections cannot change the sign of fluctuations for the energy (39) and the pressure (41).

#### **7. Conclusions**

In this paper, we have revised the long-standing problem of equilibrium pressure fluctuations and showed that its solution can be obtained on the grounds of generalizing the Bogoliubov–Zubarev theorem by using the method of quasi-averages (applied to the introduction of the volume) as well as that of scale transformation in the phase space of a physical object in question. Besides general formulation for the proof of the theorem (which can be found in Refs. [10,18]), we have presented some numerical results for the thermodynamic quantities of the relativistic gases. We hope that these results could be partly applied to the description of the hot quark-gluon plasma within the scope of thermodynamics as well as of statistical mechanics (in this connection see, e.g., papers [24,25]). However, for the moment the thermal equations of state for the pressure are formulated mostly within the phenomenological approach on the grounds of QCD thermodynamics, whereas the application of the generalized Bogoliubov–Zubarev theorem needs some dynamical description in the object's phase space (e.g., if possible, for the Mott–Hagedorn resonance gas described in [24,25]).

**Author Contributions:** Conceptualization, Y.G.R. and Y.P.R.; Methodology, Y.G.R. and Y.P.R.; Software, Y.G.R.; Validation, Y.P.R.; Formal Analysis Y.P.R.; Investigation, Y.G.R. and Y.P.R.; Resources, Y.G.R.; Data Curation, Y.G.R.; Writing-Original Draft Preparation, Y.G.R.; Writing-Review & Editing, Y.P.R.; Visualization, Y.P.R.; Supervision, Y.P.R.; Project Administration Y.P.R.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are grateful to D. Blaschke for invitation to the Symposium "Nonequilibrium Phenomena in Strongly Correlated Systems" held at Dubna, April 2018. We are also indebted to V.G. Morozov for the useful discussion of the manuscript. The paper is devoted to the memory of Dmitrii Nikolaevich Zubarev on the occasion of his 100th birthday.

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

#### **Appendix A Details of Calculating Sum (36) in Section 6**

#### **1. Case** *k***(***m***)>0**

At *k*(*m*) > 0 (0 ≤ *m* < *m*0) in the expansion for Г[*k*(*m*) > 0;*a*] only positive degrees of *a* appear

$$\Gamma[k(m)>0;a] = \Gamma[k(m)>0;0] - \int\_0^t d\eta \, e^{-\eta} \eta^{l(m)-1} = \Gamma[k(m)>0] - \sum\_{l=0}^{\infty} (-1)^l (l!)^{-1} d^{l(m)+l} (k(m)+l)^{-1}, \tag{A1}$$

where Г[*k*(*m*)>0;0] is the ordinary (i.e., complete) gamma-function Г[*k*(*m*)] (see, e.g., [22]); for integer values *k*(*m*) = 1, 2, . . . it possesses the most simple form [*k*(*m*) − 1]!

In order to obtain the expansion (A1), it is sufficient to expand the exponent e−η, which enters the integrand in the definition (34), in the Taylor series and then to integrate over *η* the relevant convergent series. In this course no singularities in (A1) in the limit *a*→0 arise, because they would appear in any of the terms on the right-hand side of (A1) only in the case of violating the condition *k*(*m*) > 0. Indeed, the function Г[*k*(*m*)] in this case would be not well defined and some of the denominators *k*(*m*) + *l* might take zero values.

Taking into account the form of the product *a*<sup>2</sup>*mak*(*m*) = *ak*(0), it can be easily seen that the first term in the braces on the right-hand side of Equation (36) contains in general case two groups of expansion terms: one running even degrees of *a* (starting with *a*0, *a*2, ... ), and another running all degrees of *a*, starting with *ak*(0) and taking subsequent values *ak*(0) + *<sup>l</sup>* , *l* = 1, 2, ... Clearly, the lowest order contribution "surviving" in the limit *a* = 0 is of the form

$$[f0] \Gamma[k(0) \ge 1] \mu^0 = [k(0) - 1]! = (f + n - 1)!,$$

with the values of *f* and *n* satisfying the aforementioned condition *k*(0) ≥ 1.

If the quantity *k*(0) = *f* + *n* has the *minimal* possible (for the case in question) value *k*(0) = 1, the lowest order contribution (linear in *a*) will be given by the first term of the second group *a<sup>k</sup>*(0). The next order contribution will be given by the second term of the first group, which is quadratic in *a*. However, if *k*(0) = 2, the terms mentioned will give the contribution of one and the same order in *a*, and only at *k*(0) = 3 the contributions of the second group (starting with *a*3) will follow, the two first terms of the first group joining the battle.

#### **2. Case** *k***(***m***)** *≤* **0**

At *k*(*m*) ≤ 0 (i.e., at *m* ≥ *m*<sup>0</sup> ≥ 0), in contrast with the case *k*(*m*) > 0, only *negative* degrees of *a* enter into the expansion for Г[*k*(*m*) ≤ 0;*a*]. This fact implies the arising of the pole singularities of all orders from 1 till |*k*(*m*)|, as well as also the logarithmic singularity in *a*, the latter singularity being the only "surviving" one even in the limiting case *k*(*m*0) = 0. However, as can be seen, these singularities do not become apparent in the final result for *h*(n)(*f*;*a*), since they are fully suppressed by the factor *a*2m appearing in every order at *m* ≥ *m*0.

In order to obtain the expansion for Г[*k*(*m*) ≤ 0;*a*] in degrees of *a*, it is worth-while to note that this quantity is defined by the integral in (34) and the condition of its convergence at the upper limit for any value of *k*(*m*) (independently on the sign) is guaranteed by the factor *e*−η. However, at the lower limit the convergence condition is violated already for the maximally possible in our case value *k*(*m*)=0 implying the logarithmic singularity. Moreover, with the decrease of *k*(*m*) (i.e., the increase of |*k*(*m*)|) there arise pole singularities of maximal order |*k*(*m*)|.

So it is appropriate to use for Г[*k*(*m*) ≤ 0;*a*] the recurrence relation enabling one to increase by the unity the value *k*(*m*) (and respectively to decrease the value |*k*(*m*)|), thus selecting the pole singularities. The relation of this kind can be easily found through the integration by parts of the original integral in (34), with the result reading

$$\Gamma[k(m)\le 0\&] = \mathbf{e}^{-a} a^{-\lfloor k(m)\rfloor} - (1/\lfloor k(m)\rfloor)\Gamma[k(m)+1\le 0\&].\tag{A2}$$

Finally, the relation (A2) permits one to express Г[*k*(*m*) ≤ 0;*a*] with an arbitrary value *k*(*m*) ≤ 0 as a function of Г[*k*(*m*) = 0;*a*]

$$\Gamma[k(m)\le 0;a] = \mathbf{e}^{-a} \sum\_{l=0}^{k(m)-1} (-1)^l ^{l+1} [|k(m)|\dots(|k(m)|-l)]^{-1} a^{-|k(m)|+l} + (-1)^{|k(m)|} (|k(m)|!)^{-1} \tag{A3}$$

The finite sum entering the right-hand side of (A3) is different from zero only under the condition *k*(*m*) < 0. Otherwise (at *k*(*m*) = 0) the relation (A3) reduces to the identity. In particular, just this sum contains all the pole singularities mentioned above.

The quantity Г[*k*(*m*) = 0;*a*] is the limiting one for all possible values *k*(*m*) ≤ 0 and coincides (up to the sign) with the integral exponent function Ei(−*a*) (see, e.g., [23])

$$\Gamma[k(m) = 0; a] = \int\_{a}^{\infty} d\eta e^{-\eta} \eta^{-1} \equiv -\text{Ei}(-a), \text{Ei}(-a) = \text{C} + \ln a + \sum\_{l=0}^{\infty} \left( l! l \right)^{-1} \text{a}^{l}, \tag{A4}$$

where *C* ≈ 0,577 is the Euler constant. The power series on the right-hand side of (A4) converges for all finite real values of *a*, but the term ln*a* possesses an obvious singularity at the limiting value *a* = 0, corresponding to the case of massless particles with *E*<sup>0</sup> = 0.

It can be shown that all the singularities mentioned above of the quantity Г[*k*(*m*) ≤ 0;*a*] disappear, as was expected, after its substitution into the second term in the braces on the right-hand side of Equation (36) due to its multiplication by the factor *a*2*<sup>m</sup>* in every order of the infinite sum over the index *m* ≥ *m*<sup>0</sup> > 0. It is quite clear for the logarithmic singularity (and also for the constant term) entering (A4). As to the pole singularities entering the finite sum in (A3) at *k*(*m*) ≡ *k*(0) − 2*m* ≤ 0, one obtains <sup>|</sup>*k*(*m*)| = <sup>−</sup>*k*(*m*)=2*<sup>m</sup>* <sup>−</sup> *<sup>k</sup>*(0) and <sup>−</sup>|*k*(*m*)| = <sup>−</sup>2*<sup>m</sup>* <sup>+</sup> *<sup>k</sup>*(0), so *<sup>a</sup>*2m*a*−|*k*(*m*)| <sup>=</sup> *ak*(0).

Thus, the two cases 2 and 1, which look on the first glance as quite different, appear to be in sufficiently complete accordance one with another. Indeed, in the case 2 the two groups of the expansion terms in degrees of *a* prove to appear: those with even degrees and also with all degrees, starting with *ak*(0) and taking subsequent values *ak*(0) + *<sup>l</sup>* (*l* = 1, 2, ... ). Note that the first group of terms in the case 2 starts not with *a*<sup>0</sup> (with the coefficient Г[*k*(0);0]), as in the case 1, but with the term *am*<sup>0</sup> (with the coefficient *C*), where according to Equation (35) the value *m*0, in general differs from zero. Otherwise, just this term proves to be the starting one for the whole expansion (36), so that the case 1 cannot be realized.

If the case 1 is nevertheless realized, the first group of terms may be represented as

$$\Sigma(f, n; a) \equiv \sum\_{m=0}^{m0-1} [f, m][k(m) - 1]! a^{2m} = [(f+n) - 1]! + [f, 1][(f+n-2) - 1]! a^2 + \mathcal{O}(a^4), \tag{A5}$$

where it was taken into account that [*f*,0] <sup>≡</sup> 1, *<sup>k</sup>*(0) = *f+n*( <sup>≥</sup> 1) and [*f*,1] = <sup>−</sup><sup>1</sup> <sup>2</sup> (*f* − 2), *k*(1) = *k*(0) − 2. The number of terms in (A5) depends upon the value of the index *m*0, which according to (33) и (35) depends in turn upon the values *f* and *n*.

As for the second groups of terms in both cases 1 and 2, it follows that they should be unified, so that the resulting contribution into the right-hand side of Equation (36) takes the form *ak*(0)*S*(*f*;*a*). Here *<sup>S</sup>*(*f*,*n*;*a*) = *<sup>S</sup>*<(*f*,*n*;*a*) + *<sup>S</sup>*≥(*f*,*n*;*a*) is the expansion in *<sup>a</sup>*, including all the degrees (starting with *<sup>a</sup>*0), and the quantities *S*<(*a*) (with *m* < *m*0) and *S*≥(*a*) (with *m* ≥ *m*0) are the following double sums

$$\begin{array}{c} \mathbf{S}\_{<}(f,n;a) = -\sum\_{m=0}^{m0-1} [f,m][f,m] \sum\_{l=0}^{\infty} (-1)^{l} (l!)^{-1} a^{l} (k(m)+l)^{-1},\\ \mathbf{S}\_{\geq}(f,n;a) = \mathbf{e}^{-a} \sum\_{m=m0}^{\infty} [f,m] \sum\_{m=m0}^{k(m)-1} [f,m] \sum\_{l} (-1)^{l+1} [k(m)|\dots(|k(m)|-l)]^{-1} a^{l}. \end{array} \tag{A6}$$

It is necessary to underline that at the point *a* = 0 in the "inner" sums over the index *l* only the first term with *l* = 0 remains. We do not study here the infinite sum (42), but in virtue of definitions (32) and (33) for [*f*,*m*] and *k*(*m*) it is seen that the general term of this sum with alternating signs is of the form (−1)*m*(2*mm*!)−<sup>1</sup> and even in the worst (in the sense of convergence) case *<sup>m</sup>*<sup>0</sup> = 0 the series (42) converges, with *S*(0) taking the finite value

$$S(f, n; 0) = -\sum\_{m=0}^{m0-1} s(f, m) + \sum\_{m=0}^{\infty} s(f, m), \\ s(f, m) \equiv [f, m](k(m))^{-1}. \tag{A7}$$

Here the quantity *S*(*f*,*n*;0) ≡ *S*(*f*,*n*), like Σ(*f*,*n*;*a*), depends on *n* through *m*0, which is determined by the relations (33) and (35).

#### **References**


© 2019 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/).

*Article*
