Deforming Gibbs Factor Using Tsallis q -Exponential with a Complex Parameter: An Ideal Bose Gas Case

: The paper presents a study of a non-standard model of fractional statistics. The exponential of the Gibbs factor in the expression for the occupation numbers of ideal bosons is substituted with the Tsallis q -exponential and the parameter q = 1 − α is considered complex. Such an approach predicts quantum critical phenomena, which might be associated with PT -symmetry breaking. Thermodynamic functions are calculated for this system. Analysis is made both numerically and analytically. Singularities in the temperature dependence of fugacity and specific heat are revealed. The critical temperature is defined by non-analyticities in the expressions for the occupation numbers. Due to essentially transcendental nature of the respective equations, only numerical estimations are reported for several values of parameters. In the limit of α → 0 some simplifications are obtained in equations defining the temperature dependence of fugacity and relations defining the critical temperature. Applications of the proposed model are expected in physical problems with energy dissipation and inderdisciplinarily in effective description of complex systems to describe phenomena with non-monotonic dependencies.


Introduction
The concept of space-time reflection symmetry commonly referred to as PT -symmetry (parity-time symmetry) [1] has recently entered a vast class of problems in classical and quantum physics, with both theoretical and experimental domains ranging from acoustics and mostly optics to topological insulators and metamaterials, see [2] and references therein. Being primarily associated with non-Hermitian Hamiltonians containing complex-valued potentials [3][4][5][6], PT -symmetry breaking in non-conservative quantum systems is linked to quantum critical phenomena of a special sort [7]. This brings a broader context to such structures as appearance of complex-valued quantities in physical problems means dissipative or decay processes, from a classical example of complex refractive index [8] (Chap. 9.4) to Bose-Einstein condensation in leaking optical lattices [9]. Applications of complex thermodynamic quantities are exemplified by temperature [10], chemical potential [11][12][13], energy [14] or magnetic field [15].
Another concept considered in the present work is nonextensive and nonadditive distributions, which originated from information theory almost sixty years ago [16,17]. They were brought into physics by Tsallis [18] and gained multidisciplinary applications in studies of complex systems in the domain of biology, climatology, economics, linguistics, and many other fields [19]. Tsallis-like generalizations are usually applied to systems with significantly pronounced non-Markovian (memory) effects [20].
Various approaches are known to nonextenside and nonadditive generalizations of quantum gases, in particular both ideal [21][22][23][24][25][26] and interacting [27] Bose gas. In the present work, the approach from paper [25] is extended to complex values of the nonextensivity parameter q in the Tsallis exponential [28]. Curiously, generalizations of the Tsallis entropy to the complex values domain were not popular until last decade [29][30][31][32][33]. Applications of the complex nonextensivity parameter include interpretation of data on high-energy particle collisions [32,34], studies of directed networks [35], analysis of complex neural networks [36], and even are related to the income distribution [37].
The paper is organized as follows. Section 2 contains basic expressions and the description of the calculation procedure. These are further detailed in Section 3. Numerical results and analytical high-temperature behavior of thermodynamic functions are presented in Section 4. Critical point is analyzed in Section 4 and some limiting analytical results are obtained in Section 6. Brief discussion in Section 7 followed by a short Materials and Methods section conclude the paper.

Starting Points
For the sake of completeness, let us briefly summarize the procedure to calculate thermodynamic functions, cf. [25]. Mean occupation numbers n(ε j , T, z) of the jth level with energy ε j are functions of temperature T and fugacity z. Their sum over all the levels with degeneracies g j yields the total number of particles in the system: Solving this equation implicit for z, one obtains z = z(T, N). This solution has to be inserted in the definition of the total energy leading to E = E(T, N). Subsequent calculations are quite straightforward. For instance, specific heat, which is used in the paper to demonstrate peculiarities in the behavior of thermodynamic functions, is just The expression for the mean occupation numbers is chosen in the following form: where the exponential in the Gibbs factor e ε j /T is substituted with the so called Tsallis q-exponential [28] This is a phenomenologically introduced generalization of the well-known expression for bosons [38] (p. 82), and thus the model considered in the present work corresponds to a fractional quantum statistics. This notion is understood here in a wide sense, as a statistics differing from standard Bose-Einstein or Fermi-Dirac one, cf. [39]. Such types of statistics might be also referred to as "intermediate" or "exotic".
To make a step further comparing to previous studies [25], the parameter q is considered a complex number being represented in the form Note that complex parameters in various types of deformed statistics can be used, in particular, to effectively account for maximal level occupancy [40], interparticle interactions [41] or small dissipative branch in elementary excitation spectra [42,43].
To simplify the mathematical treatment of the problem, the summation over the levels in (1) and (2) is changed to the integration over energies with the density of states g(ε): and Such a change is justified primarily for interlevel separations much smaller comparing to the temperature T. Infinitely small energy level separations also appear in the thermodynamic limit: for instance, for free particles in a box with side L the quantized values of energies scale as 1/L at L → ∞ and for D-dimensional harmonic oscillators the frequencies scale as The density of states is considered in the form [44] where the constant A depends on parameters of the system under study, e.g., particle mass or oscillator frequency. Some particular values of s include s = D/2 for free particles in a D-dimensional box or s = D for harmonic oscillators in D dimensions [45]. A more general dispersion relation ε = ap b , yields the dependence g(ε) ∝ e D/b−1 [39] (p. 150). The explicitly written factor of N simplifies analysis in the thermodynamic limit [44].

Calculations of Fugacity and Energy
The equation defining fugacity z is For complex α, the function (1 + αε/T) 1/α is multivalued; the principal branch is considered in the calculations throughout this work. Making a simple change of variables x = ε/T, we can write two equations-separately for real and imaginary parts-as follows: They define the temperature dependence of the complex-valued fugacity z = z + iz . The total energy per particle is then given by Note that complex energies are usually associated with dissipation processes [12,42].
The convergence of the integrals in (8) or (11) is ensured if We thus immediately see that purely imaginary parameter α = iα is not applicable. From intermediate transformations, it becomes clear that negative values of α should not be considered as well. So, after simple manipulations the relation between s and α is as follows: Moreover, one should also take into consideration that the calculation of energy requires integration [see (9) or (13)] with s + 1. Finally, we arrive at Consider, for instance, s = 2 and α = 0.1 + 0.15i being one of the values attempted by the author. This yields 0.975 < 1, although satisfying the condition but making the integral for the calculation of energy very slowly convergent. To avoid such numerical issues, further in the paper the values of α are chosen to ensure the right part of (17) well apart from unity.
The integrations in (12) and (13) can be carried out similarly to ordinary Bose gas problem [46] [Chap. V], by expanding the integrands into series (at least for |z| < 1) [25]. The resulting expressions, e.g., are written using a generalization of the polylogarithm function, which is represented by the following series for |z| ≤ 1 where B(x, y) is the beta function, and defined by the integral in (18) otherwise. In the limit of α → 0 it coincides with the ordinary polylogarithm (also called Jonquière's function or Bose-Einstein integral) The first (integral) definition in (20) applies for all complex z except for purely real z with Re z ≥ 1. The second one is valid for complex z with |z| < 1. Note that if z = 1 the polylogarithm reduces to Riemann's zeta function, Li s (1) = ζ(s).
Technical difficulties in solving the equation for fugacity are caused by problems in inverting the Li α,s (z) function, as the respective analytical expressions are not even known for ordinary polylogarithms, cf. [43]. We thus will rely on numerical analysis and only obtain some limiting dependencies analytically.

Numerical Results and High-Temperature Behavior
In Figures 1 and 2, the results of numerical calculations for fugacity and specific heat are shown for several values of α and s. The units of temperature and energy are A −1/s . Except for s = 1, singularities in the temperature dependencies of the fugacity are observed at certain temperatures T c . This critical point is discussed in more detail in the next section. For T < T c , the values of z are fixed equal to their values at the critical temperature z c = z c + iz c ≡ z (T c ) + iz c (T c ). This is done by analogy with both ordinary ideal Bose gas and its nonadditive generalization [25], where the fugacity equals unity below the critical temperature. The low-temperature behavior is just entirely guided by the exponent s, yielding in particular the specific heat C/N ∝ T s . From Equation (11) defining fugacity, it is immediately clear that the following symmetry holds: the sign change in α yields the sign change in z , as one can observe in the right panel of Figure 1. Indeed, these sign changes just correspond to complex conjugation, while the right-hand side of the respective equation contains a real number N. Obviously, real part z remains unchanged at α → −α . We can see in particular that above the critical point a maximum is observed on the specific heat curve for complex-valued α, but not for the real one (α = 0). A similar picture was discovered for another fractional statistics with a complex parameter [43].
At high temperatures, specific heat tends to a constant limit, which can be obtained exactly as in [25], namely: The values of high-temperature specific heat together with critical temperatures and the respective critical fugacities are summarized in Table 1. Note that for α = 0, the expected classical value of the specific heat C/N = s is recovered.  For s = 1, the zero value of the critical temperature is written formally: this means that no critical point exists for such a system. The same situation is known for the ordinary ideal Bose gas of particles in a box in two dimensions or two-dimensional harmonic oscillators, both having the density of states g(ε) = const.
As one can see, the values of the critical temperature decrease with the imaginary part α increasing. This observation correlates, in particular, with the results of [47], where the critical behavior of an interacting Bose gas in an optical lattice was studied.
Small values of z c at s = 2 (see Table 1) could have misled one to the condition defining the critical temperature in the form z (T c ) = 0 similar to that obtained for the Polychronakos statistics with a complex parameter [42,43]. More details on the critical temperature are given in the next section.

Towards the Critical Point
The critical point is achieved when singularities occur under integrals in (12) or (13). This in turn corresponds to Therefore, the following set of four equations defines four variables at the critical point: T c being the critical temperature, the fugacity value at the critical temperature z c + iz c , and the coordinate x 0 where the singularity occurs.
Consider small deviations of q from unity, i.e., α , α → 0. In this limit, from (23) we obtain: Keeping only terms up to linear in α and α , we have thus yielding two equations (for the real and imaginary parts): From the second equation we obtain (linear in α and α ) The leading order of (27) defines the singular coordinate x 0 as follows: The link between the real and imaginary parts of the fugacity is thus Thus, the critical temperature T c corresponds to the point when the solutions z c and z c of first two equations in (24) satisfy relation (31). Note that the solution with z c > 1 should be taken, see (30).
We can now discuss what happens at temperatures below T c . It is worth recalling the ordinary ideal Bose gas with z c = z(T c ) = 1. From Equation (30) we obtain the singular coordinate x 0 = 0 corresponding to the ground state ε = 0 for all temperatures T ≤ T c , and, consequently, a macroscopic (mathematically infinite) number of particles in the ground state, see (6). This phenomenon is known as the Bose-Einstein condensation and n 0 /N is the condensate fraction. As the contribution of the single level j = 0 is a macroscopic number, it is not taken into account properly when changing summation over energy levels with integration and should be written explicitly. In the considered model with a complex q-exponential the situation is slightly different. Equation (30) means in this case that singularities occur at different energy levels ε c at different temperatures, so that x 0 = ε c /T. The number of particles on these levels form an analog of the Bose condensate, with the fraction easily obtained in the following form: The respective illustration is shown in Figure 3.

The Limit of α → 0
One can also perform analysis of other expressions in the limit of small α , α . Using the series definition at |z| < 1 given by (19), from (18) one gets: The limiting behavior of the beta function with one parameter fixed and the other tending to infinity is [48] This yields immediately: Keeping only terms linear in α we finally obtain: The same equation can be derived with the integral definition (18) by expanding the Tsallis exponential into series over α: The second integral evaluates to the derivative of the polylogarithm, which can be subsequently transformed using the identity This yields identical to the previously obtained Equation (37) due to the well-known recursion property of the gamma function Γ(s + 1) = sΓ(s).
Let us return to Equation (37) and write polylogarithms as the series over z = z + iz explicitly: Keeping only linear terms in α , α and z (which should also be small at small deviations from the real-valued statistics), we obtain for the real and imaginary parts: so the following relation holds in the limit of α → 0: Further analytical treatment, however, is not possible even in this limit. Series representations of polylogarithms around z = 1 or z = ∞ involve powers of ln(−z) [48] and thus equations remain essentially transcendental.

Discussion
A new type of fractional statistics has been suggested in the paper. It pleaches the nonadditivity and complex nature of the statistics parameter, yielding thus a broader spectrum of properties. Thermodynamics of the ideal Bose gas with the Gibbs factor deformed using the Tsallis q-exponential in case of complex q = 1 − α has been studied. The analysis has been conducted mostly numerically, with the limiting behavior for α → 0 obtained analytically.
The studied system exhibits critical behavior caused by non-analyticities in the expression for occupation numbers. At the critical temperature T c , specific heat has a discontinuity meaning thus a second-order phase transition. The values of T c decrease as the power in the density of states (or space dimensionality) increase. The decrease of T c has been also observed for the imaginary part of the statistics parameter α = Im α increasing while keeping other parameters fixed.
Analytical expressions for thermodynamic functions involve a deformed polylogarithm function Li α,s (z). Detailed studies of its properties in the complex domain would significantly facilitate future analysis of the considered model and allow for deeper understanding of the relation between the α parameter and critical behavior. This include better analytical estimations for the critical temperature T c and clarification of the temperature dependencies for T < T c .
Discontinuities in the specific heat at the critical temperature resemble those known for ordinary bosons, which can be caused by interparticle interactions or effective space dimensionality greater than three. The complex nature of the statistics parameter suggests that such bosonic systems with energy dissipation could be potential fields for the application of the model proposed in the present work. The respective quantum critical phenomena could be connected with the P T -symmetry breaking and the order parameter appearing in such problems is linked to the imaginary part of the fugacity. Finally, interdisciplinary applications of the analyzed statistics is envisaged for phenomena with non-monotonic dependencies, which is related to the complex parameter in the expression for occupation numbers.

Materials and Methods
Numerical computations in the work are made in wxMaxima 13.04.2, a graphical interface for the computer algebra system MAXIMA (http://maxima.sourceforge.net). In particular, integrations are carried out by the Quadpack function quad_qagi, see https://web.csulb.edu/~woollett/mbe8nint.pdf. Numerical equation solutions are obtained by the mnewton procedure from the mnewton package (https://web.csulb.edu/~woollett/mbe4solve.pdf).
Funding: This research was partly funded by the Ministry of Education and Science of Ukraine, grant number 0119U002203.