Critical point and Bose-Einstein condensation in pion matter

The Bose-Einstein condensation and the liquid-gas first order phase transition are studied in the interacting pion matter. Two phenomenological models are used: the mean-field model and the hybrid model. Free model parameters are fixed by fitting the lattice QCD data on the pion Bose condensate density at zero temperature. In spite of some minor differences the two model demonstrate an identical qualitative and very close quantitative behavior for the thermodynamic functions and electric charge fluctuations. A peculiar property of the considered models is an intersection of the Bose-Einstein condensation line and the line of the first order phase transition at the critical end point.


I. INTRODUCTION
The Bose statistics [1] and Bose-Einstein condensation (BEC) phenomenon [2] were predicted almost hundred years ago. The BEC can take place in equilibrium systems of noninteracting bosons when a macroscopic part of all particles begins to occupy a single zero-momentum state. This was approximately realized experimentally in atomic systems at very small temperatures and particle number densities [3][4][5][6].
The BEC can also happen in condensed matter, nuclear physics, astrophysics, and cosmology (see, e.g., Refs. [7][8][9][10][11][12][13][14][15]).In most of these situations, particle interactions should be taken into account. If both the attractive and repulsive interactions between particles are taken into account, the system reveals the first-order liquidgas phase transition (FOPT) and the critical endpoint (CP). Therefore, the two phenomena BEC and CPcan be expected for interacting bosons.
Pions are three pseudo-scalar mesons, π + , π 0 , π − , that obey the Bose statistics, thus, an emergence of the BEC of pions is possible. The BEC can take place during the cooling of the early Universe [16], in the gravitationally bound pion stars [17][18][19], or as a non-equilibrium phenomenon in heavy-ion collisions [8,9,20]. It has been also predicted to occur at large isospin (electric) chemical potentials [21,22]. This suggestion has been studied in Quantum Chromodynamics (QCD), a theory of strong interactions. A presence of the BEC was supported by recent first principle lattice QCD simulations [23,24]. The Bose condensate (BC) of pions was observed at low temperature, T < m, and large electric (isospin) chemical potential m < µ < 2m, where m is the pion mass (in what follows we use an approximate value m = 140 MeV neglecting a small difference between the masses of neutral and charged pions). In this specific region of the (µ, T ) phase diagram the QCD matter is expected in a form of the interacting pions. The heavier hadrons and/or quarkgluon degrees of freedom are expected to be suppressed at T m and not too large µ. Various approaches to a description of the pion matter were developed: chiral perturbation theory [25,26], linear sigma model [27,28], Nambu-Jona-Lasinio model [29] or Polyakov-loop extended quark meson model [30,31], functional renormalisation group [32,33], hard thermal loops [34], self-interacting mean field theory [35], etc.
Influence of particle interactions on the thermodynamic properties was also considered within the Smatrix formulation of statistical mechanics [36]. In particular, the attractive and repulsive interaction from hadron-hadron scatterings were discussed in a number of works [37][38][39][40][41]. A possibility of the BEC in the pion system with zero chemical potential was also discussed using a Skyrme-like model including both attractive and repulsive interaction [14,35,42].
A description of the repulsive and attractive interaction in statistical systems of hadrons is often performed in terms of the following phenomenological approaches: Mean-field approximation, effective-mass model, excluded volume approximation, etc. The mixtures of these approaches are also used. For example, a famous Walecka model [43] for nuclear matter describes repulsive interactions in terms of the mean field U (n) being a linear function of nucleon number density, and attractive interactions in terms of the effective nucleon mass. These different phenomenological models belong to the same universality class, the so-called classical models [44] (the names mean-field models and van der Waals-type models are also used). These classical models lead to very similar description of the FOPT and CP (see, e.g., Ref. [45]). The models from the same universality class for the CP can however lead to rather different consequences for the BEC. This was discussed in Ref. [39] where only the repulsive parts of interactions were taken into account.
In Ref. [46] the effective mass model with φ 4 attractive and φ 6 repulsive interaction was studied. It was found that both phenomena, FOPT and BEC, take place. An additional peculiar feature of the model was an observation that the CP belongs to the line of the BEC. How do these features of the pion matter depend on the specific model used in Ref. [46]? This question motivates the present studies. We discuss two other phenomenological descriptions of the pion matter. The first model considers both attractive and repulsive interactions in terms of the mean field U (n) depending on the pions number density. The second "hybrid" model treats the repulsive interactions as the mean field, and the attractive interactions in terms of the effective pion mass m * (T, µ) < m. We first fit the lattice QCD (lQCD) data at zero temperature and finite isospin chemical potential µ to fix the model parameters. Then, thermodynamic functions and electric (isospin) charge fluctuations are calculated in the (µ, T ) plane.
Intensive measures of the electric charge fluctuations, the scaled variance, skewness, and kurtosis, appear to be very sensitive to a presence of the CP and the BEC. The paper is organized as follows. The ideal Bose gas of pions is considered in Sec. II. The two phenomenological models of the interacting pion matter are discussed and compared in Sec. III. Summary in Sec. IV closes the paper.

II. IDEAL PION GAS
The ideal pion gas (IdPG) is described in the grand canonical ensemble by the pressure function [47] p id (T, µ) = where integration by parts was used. The chemical potential µ corresponds to the electric charge conservation, and is the Bose momentum distribution. The particle number densities of the i-th sort of pions equals to The electric charge density is calculated as and the number density of all pions as where K 2 is the modified Bessel function.
The inequality |µ| ≤ m should be satisfied in the IdPG. An onset of the BEC occurs at |µ| → m − 0. The condition n Q = n Q (T, |µ| = m) defines then a line T BC = T BC (n Q ) denoted further as the BEC line. Under this line at T < T BC (n Q ) there is a region with the nonzero BC. In what follows we consider µ ≥ 0. It leads to n + ≥ n 0 ≥ n − and n Q = n + − n − ≥ 0, where n + , n 0 , and n − denote the number densities of π + , π 0 , and π − , respectively. The results for µ ≤ 0 can be obtained by interchanging π + and π − . A line of the BEC T BC (n Q ) for the IdPG corresponds to µ → m − 0. It is shown in the (n Q , T ) plane in Fig. 1 (a) 1 . At this line there is an onset of the π + BEC. Under this line nonzero values of the π + BC are formed. The total pion number density (6) should be then modified at T < T BC (n Q ), where n + BC ≥ 0 corresponds to the BC of π + . At T → 0 all thermal densities (4) vanish. Thus, n = n Q = n + BC at T = 0, i.e., at zero temperature the pion system consists from the pure BC of π + . A dashed line in Fig. 1 (a) shows the BEC line for only one sort of pions, i.e., for π + It gives a good approximation of the general case with all three sorts of pions. This is because of n − n + on the BEC line T BC (n Q ). Therefore, an analytic behavior of the T BC (n Q ) at T BC /m 1corresponds approximately to the well-known textbook result for one sort of nonrelativistic bosons [44]: where ζ(3/2) ∼ = 2.612 is the Riemann zeta function. Corrections to Eq. (8) from π − are small as seen from Fig. 1 (a). Relativistic corrections to Eq. (8) at T /m ≥ 1 are considered in Ref. [9]. The ratios n 0 /n + and n − /n + along the BEC line are shown in Fig. 1 (b). The inequalities n + n 0 n − remain valid on the BEC also for the interacting pion matter. An intensive measure for fluctuations of the electric charge Q = N + − N − is the scaled variance where N + and N − are the total numbers of π + and π − , respectively, and . . . denotes the grand canonical averaging. In the IdPG it takes the following form becomes thus divergent at the lower limit as k −2 . Therefore, the scaled variance ω Q becomes infinite on the BEC line.
At µ → m − 0, the electric charge fluctuations show an anomalous behavior ω Q ∝ V 1/2 in the finite volume V . This leads to infinite values of ω Q on the BEC line [8,9]: It happens due to the BEC of π + , and is rather similar to a behavior of the ideal Bose gas for one sort of particles. Under the BEC line an additional contribution ω C from the BC n + BC > 0 should be added to ω Q (10). It behaves as [9] in the large finite system and goes to infinity in the thermodynamic limit V → ∞.
It should be noted that the scaled variance (9) can not be applied at µ → 0. It becomes meaningless in this limit as Q → 0. At N + − N − = 0 the other fluctuation measure (N + − N − ) 2 / N + + N − is usually used to describe the charge fluctuations.

III. INTERACTING PION GAS
In this section, two phenomenological models for the repulsive and attractive interactions in the pion matter are considered.

A. Mean-field model
The mean field model for the system of interacting pions is given by the following set of self-consistent equations: where U (n) is the mean field that describes pion interactions, µ i in Eq. (14) is given by Eq. (2), and n + BC in Eq. (13) is the BC density of π + that can attains nonzero values when the BEC condition, µ * = m, is fulfilled. A second term in the right hand side of Eq. (12) corresponds to the so-called excess pressure that makes Eqs. (12)- (14) to be thermodynamically self-consistent (see, e.g., Ref. [48]). The mean field potential U (n) will be taken in the following form: where constants A and B correspond to the attractive and repulsive interactions, respectively. In the model with mean field the BEC can take place at both µ < m, for U (n) < 0, and µ > m, for U (n) > 0. For the poten-tial U (n) given by Eq. (15) these conditions correspond to small and large values of n, respectively.
At T = 0 the pion system can only exist in the form of the π + BC. A condition of the BEC, µ − U (n + BC ) = m, leads to the following solution for n + BC : The system pressure at T = 0 is given by a second term in the right hand side of Eq. (12): At small µ the system exists in the "gaseous" phase with n = 0 and p = 0. At some µ = µ 0 there is the FOPT. According to the Gibbs criteria the pion BC density jumps to the "liquid" phase with n + BC = n 0 > 0 and the pressure p = 0. This FOPT takes place at where the liquid pressure (17) equals to zero, and the chemical potential µ 0 found from Eq. (16) is which defines also the ground state binding energy per At µ > µ 0 the BC n + BC as a function of µ is defined by Eq. (16). In Fig. 2 (a) the fit of the Monte Carlo lattice data for n + BC at T = 0 with Eq. (16) is presented. The fitting parameters A and B are fixed as With these parameters one finds: This ground state of the pion BC looks rather rarefied and weakly bounded when it is compared with the ground state of the nuclear matter, n nuc 0 ∼ = 0.16 fm −3 and W nuc ∼ = −16 MeV.
In Fig. 2 (b) the ratio of the pressure of interacting pion gas to that of the IdPG is shown by a solid line as a function of T at µ = 0. The effects of the repulsive interaction suppress the pion pressure, p < p id , and they are seen at large T . The repulsive effects are still rather moderate and correspond approximately to the excluded volume corrections with rather small hard-core pion radius r ∼ = 0.13 fm. Tiny effects of the attractive interactions are only seen at small T where p > p id . Note that our modelling concerns the non-resonance part of the pion-pion interactions. These interactions contribute to the n + BC (µ) at T = 0 presented in Fig. 2 (a). At µ = 0 a contribution of these non-resonance residual pion interactions to the thermodynamic functions are small. There are almost no chances to observe them in hadron statistical equilibrium models using the data on heavy-ion collisions. This is because of |µ| m in the equilibrium systems created in heavy ion collisions. Previous suggestions to observe the BEC in heavy ion collisions assumed a large values of µ ≈ m for all three types of pions due to chemical non-equilibrium effects [8,9] .
The mean-field model with free parameters A and B was discussed in Ref. [35] for µ = 0. Several interesting phenomena including the BEC can take place in this case due to large attractive interactions, A ≥ A cr = 2(Bm) 1/2 . The required large values of the parameter A are however fully unrealistic ones, A cr is about 45 times larger than the A value given by Eq. (21). These large values A ≥ A cr are in a strong contradiction with lattice data at T = 0.
In contrast to µ = 0 case, at large µ ∼ = m and µ > m the small pion interactions are however crucially important. A vertical dotted line µ = m in Fig. 2 (a) presents a behavior µ = m in the IdPG. As seen from Fig. 2 (a) the IdPG behavior is far away from the lattice data. A dashed line in Fig. 2 (a) presents the model results at A = 0, i.e., when only repulsive interactions are included. The both fits, with A > 0 and A = 0, are of the similar quality with χ 2 per degree of freedom (dof) ∼ = 2. This fact means a dominant role of the repulsive interactions. The lattice data at T = 0 presented in Fig. 2 (a) can not give indisputable estimates of the parameter A. Similar conclusions were made in Ref. [46].
The lattice data shown in Fig. 2 (a) defines the value of parameter B rather accurately, but give only some restrictions from above for A. Thus, a presence or absence of the FOPT remains as an open question. In Fig. 2 (a) a dashed-dotted line presents the π + BC calculated in the leading-order chiral perturbation theory χ PT in Ref. [21]. This model being in a good agreement with the lattice data for π + BC at T = 0 does not include the FOPT and predicts a second order phase transition.
In the mean-field model the BEC line is defined by the condition µ = m + U (n). At U (n) > 0, one thus observes µ > m. The BEC line T BC (n) is approximately a universal function of particle number density 2 and it remains the same as for the IdPG where U (n) = 0. This is valid for U (n) > 0 in the region with no FOPT. At A = 0 there is no FOPT and the BEC line T BC (n Q ) look almost the same as for the IdPG shown in Fig. 1  (a). A behavior of the electric charge fluctuations ω Q (10) is however drastically changed on the BEC line and under this line due to a presence of the repulsive meanfield interactions. As will be seen below, in contrast to the IdPG, the scaled variance ω Q becomes finite at all n Q > 0 due to the repulsive interactions.
The phase diagram of the interacting pion gas is shown in Fig. 3. For A > 0 and B > 0 the considered meanfield model reveals the FOPT. The mixed phase region is constructed by the standard Gibbs procedure: when several solutions exist for the pressure function p at the same T and µ values, the physical solution corresponds to that with larger value of p. The line of the FOPT in the (µ, T ) plane corresponds to the equal pressures of gaseous small density and liquid large density solutions.  the left and right binodal curves that are shown in Fig. 3 (a). At T < T c , the pion system inside the mixed phase is an inhomogeneous mixture of a rarefied gas and a dense liquid. The electric charge density n Q of the pion matter is then given by a linear combination of the gaseous phase with n g Q < n Q lying on the left binodal, and the liquid phase with n l Q > n Q lying on the right binodal. The right liquid binodal of the mixed state includes n + BC > 0 while the conditions for the BEC are not fulfilled in the left gaseous binodal. At T → 0, one obtains n g Q → 0 and n l Q = n + BC , i.e., the gaseous phase is vanished, and the liquid phase consists of the π + BC.
A peculiar property of the model is a position of the CP lying on the BEC line. This is similar to the results of Ref. [46] and is a consequence of smallness of the attractive forces in the pion matter. We checked that for A larger than some critical value the BEC line intersects a line of the FOPT at the triple point T = T tr < T c . Such a behavior was found for the interacting α-matter in Ref. [7].
We compare the mean-field model with the lattice QCD results [49,50] at T > 0. Figures 4 (a) and (b) present, respectively, a position of the BEC line in the (µ, T )-plane and n Q as a function of µ at T = 124 MeV. One observes a good agreement with the lattice QCD data in Fig. 4 (a) at not too large T . At T > 100 MeV meson resonances give a substantial contribution to the thermodynamic observables. These degrees of freedom are absent in the present version of the mean field-model. This point will be a subject for the future studies. To estimate the resonance contribution to n Q we calculate the n ρ Q value that comes from a presence of non-interacting ρ ± mesons in the pion system: where g ρ = 3 and m ρ ∼ = 775 MeV. A dashed line in Fig. 4 (b) presents a sum of n Q value in the interacting pion system and n ρ Q contribution (24) from the ideal ρ-meson gas.
We introduce the quantities (i = +, 0, −) that are equal to the scaled variances of the IdPG ω id i , but with shifted chemical potentials µ i → µ * i . One finds an explicit expression for ω Q (see the Appendix): where n i ≡ n id i (T, µ * i ). Approaching the BEC line, i.e., at µ * + → m − 0, the ω i values demonstrate the following behavior: ω + → ∞ while both ω 0 and ω − remain finite. This leads to the result on the BEC line: If A = 0 the value of dU/dn is always positive at n > 0. Thus, in contrast to the IdPG, the scaled variance ω Q on the BEC line becomes finite due to the repulsive interactions. This conclusion remains also valid for A > 0 at T BC (n Q ) > T c where dU/dn > 0. It can be also shown that the scaled variance ω Q is continuous across the BEC line.
For A > 0 the CP becomes the end point of both the FOPT line and the BEC line. When T → T c along the BEC line the first term in the right hand side of Eq. (27) goes to infinity, and, thus, the second finite term gives a negligible contribution to ω Q . A relative contribution of this second term increases monotonously with increasing T BC (n Q ) along the BEC line. It still remains small for the discussed region of the system temperature. Note that ω 0 ∼ = 1, ω − ∼ = 1, and n 0 < n + , n − n + . At T c < T BC (n) < 0.8 m the relative contributions to ω Q (27) from each terms, n 0 ω 0 /n Q and 4n − ω − /n Q , are smaller than 3%. Neglecting these small contributions to ω Q one obtains an expression It corresponds to the mean field model result for one particle species n = n + (see, e.g., Ref. [39]). At the CP, dU/dn = 0 and ω Q according to (27) goes to infinity. For U (n) given by Eq. (15) the condition dU/dn = 0 gives the pion number density n = n c and the chemical potential µ = µ c at the CP: The critical temperature is calculated numerically as T c ∼ = 0.21 m ∼ = 28 MeV. At the CP one also finds n + . n 0 and n + n − . Therefore, the electric charge density at the CP is approximately A behavior of ω Q on the (µ, T ) plane is shown in Fig. 5. Approaching the CP by any path, one observes ω Q → ∞. At n Q > 0 this is the only point on the phase diagram with infinite value of the scaled variance ω Q . Most interesting regions of the phase diagram are µ ∼ = m and µ > m. These are regions of the FOPT and the BEC. For these two phenomena both the repulsive and attractive interactions between pions play a crucial role.

B. Hybrid model
Another phenomenological model considered in our paper is constructed by combining the two frameworks: the mean field U (n) = Bn 2 to describe the repulsive interactions (the same as in the mean-field model considered in the previous subsection) and effective pion mass m * (T, µ) < m for attractive ones. It resembles the Walecka model [43] for the symmetric nuclear matter. A principal difference however is the Bose statistics for interacting pions instead of Fermi statistics for interacting nucleons. There are also some technical differences. For example, we consider quadratic repulsive mean field Bn 2 instead of the linear function of the nucleon number density in Walecka model. The considered model for interacting pion matter will be called the hybrid model, and it is defined by the following equations: where the scalar density At T = 0 one obtains: Equations (35)- (37) give the functions µ = µ(n + BC ) and p = p(n + BC ) that are identical to the corresponding expressions for the mean-field model with U (n) given by (15) and considered in the previous subsections. Therefore, these two models are identical at T = 0, and free model parameters of the hybrid model will be taken as in Eq. (21). At T > 0 the models are not identical. However, all conclusions concerning the FOPT and BEC are the same for the both models. Even more, the numerical values of T c and µ c in the hybrid model are also approximately equal to those in the mean-field model. This is shown in Fig. 6 (a). Because of these reasons we do not present n Q (T, µ) and ω Q (T, µ) for the hybrid model. These results are very similar to those in Figs. 2 and 3 for the mean-field model.
The ratio m * /m is shown on the (µ, T ) plane for the hybrid model with parameters (21) in Fig. 6 (b). This ratio is slightly smaller than unity, but rather close to it. This corresponds to the small attractive effects in the hy- brid model. In a nonrelativistic approximation one finds n s ∼ = n and observes a straightforward correspondence of the effective mass m * to the attractive part of the mean field potential −An: Nevertheless some differences between these two models can be found. For example, the ratio of the pressure function in the hybrid model to pressure of the IdPG at µ = 0 is shown as a function of T by a dashed line in Fig. 2 (b). The two models with the same A and B parameters lead to different p = p(T ) functions.

C. Fluctuations of higher orders
Electric (isospin) charge susceptibilities χ j in the grand canonical ensemble defined as (j = 1, 2, . . .): Susceptibilities are derivatives of the thermodynamic potential and give, therefore, additional important information of the equation of state. Particularly their values are very sensitive to both the CP and BEC phenomena. Some ratios of the susceptibilities (40) are well known and used to quantify the fluctuations of conserved charges.
Most familiar of these measures are the scaled variance ω, skewness Sσ, and kurtosis κσ 2 (see, e.g., Ref. [51]): In this subsection the results for Sσ and κσ 2 are presented within the hybrid model with parameters (21). These fluctuation measures are presented in Figs. 7 (a) and (b), respectively. The results of the mean-field model are essentially the same. At small chemical potential µ m and not too large T the pion densities are small. In this case both the pion interaction and Bose statistics effects can be neglected. It gives Sσ ∼ = 1 and κσ 2 ∼ = 1 that corresponds to the ideal classical gas limit. Both these measure strongly deviate from these baseline values of the ideal Boltzmann gas at µ ∼ = m and µ > m. This is due to a presence of the FOPT and the BEC effects, respectively.
Skewness. The skewness Sσ is presented in Fig. 7  (a). This measure attains both positive and negative values on the (µ, T ) plane. The positive values correspond to those regions of the phase diagram where n + BC = 0. The skewness Sσ has a discontinuity along the BEC line and jumps to negative values in the phase with the BC n + BC > 0. At the CP, Sσ shows the singular behavior: it can go to both −∞ and +∞ depending on the path of approaching to the CP. When crossing the FOPT there is a discontinuity of Sσ from positive values in the gaseous phase to the negative ones in the liquid phase.
Kurtosis. The kurtosis κσ 2 in the (µ, T ) plane is shown in Fig. 7 (b). The considered models demonstrate only positive values κσ 2 > 0 for the whole (µ, T ) plane. This is in contrast to the universal behavior of fluctuations in the Ising model, as well as in various phenomenological model calculations (see, e.g., Ref. [35] and references therein), where negative values of κσ 2 are observed in the so-called analytic crossover region above the critical temperature. The large but finite values of the kurtosis are generally obtained in a vicinity of the BEC line with a discontinuity to smaller values on the FOPT and BEC line at increasing µ. The singular behavior with κσ 2 → ∞ takes place at the CP. The values of the kurtosis remain large with κσ 2 1 even far away from the CP. This is due to a large sensitivity of the higher order fluctuations to the CP.

IV. SUMMARY
Thermodynamic properties of the interacting pion matter are studied in the two phenomenological models: the mean-field model with the potential U (n) and the hybrid model. The potential is chosen as a function of the pion number density U (n) = −An + Bn 2 , and it includes both the repulsive Bn 2 and attractive −An interactions. The hybrid model assumes the same repulsive mean field potential Bn 2 , while the attractive interactions are described by the effective mass with the excess pressure −(m − m * ) 2 /(2A) similar to that in the Walecka model for interacting nucleons [43]. Model parameters A > 0 and B > 0 are fixed by fitting the lattice QCD data on the BC pion density at T = 0 as a function of the chemical potential µ. At zero temperature, the both considered model become identical to each other. Thus, the fitting procedure at T = 0 leads to the same set (21) of the A and B parameters. The two phenomena -the BEC and the FOPT with the CP -take place in the pion matter. In spite of some minor differences the two model demonstrate an identical qualitative and very close quantitative behavior for the thermodynamic functions and electric charge fluctuations in the whole (µ, T ) plane. Note that the qualitative features found in these two models are also in agreement with the results obtained in Ref. [46].
The interaction parameters A and B (21) found from fitting the lattice data correspond to rather moderate interactions in the pion matter. At µ = 0 these interactions are completely unimportant in the pion thermodynamics. It should be emphasized that the mesonic resonances as a part of the pion-pion interactions are not included in our consideration. The residual non-resonance pion-pion interactions are however crucially important at µ ∼ = m and µ > m.
If A > 0 and B > 0 the both models demonstrate the FOPT with a position of the CP at µ c ∼ = m and T c ∼ = 28 MeV. The BEC line merges to the CP. At T < T c only the liquid (dense) pion phase includes the BC n + BC > 0, while n + BC = 0 in the gaseous (rarefied) phase. In the ideal pion gas, the scaled variance of electric charge fluctuations becomes infinite on the BEC line and under this line. In contrast to this ideal gas behavior, a presence of the repulsive interactions makes ω Q to be finite and continuous function. The only point of anoma-lous electric charge fluctuations is the CP. At the CP both the scaled variance ω Q → ∞ and kurtosis κσ 2 → ∞. The skewness Sσ has a more complicated behavior. It can go to both +∞ and −∞ depending on the way of approaching the CP. A special feature of the considered models is an absence of negative values of the kurtosis. The negative values κσ 2 < 0 are usually happen in a crossover region near the CP. These negative values are absent in the considered models of the pion matter. This is because of the fact that two phenomena -an onset of the BEC and the CP -takes place at the same point.
Both the FOPT and BEC are mainly defined by π + mesons. A presence of π − and π 0 mesons give only moderate numerical corrections and does not change the qualitative properties of the pion matter at µ ∼ = m and µ > m.
The critical point in the system of interacting pions can be searched in the lQCD. The lattice simulations should focus on direct computations of the charge fluctuations measures and include some "small" temperatures between 0 and 50 MeV, i.e., in a vicinity of the hypothetical critical point.