Phase diagram of interacting pion matter and isospin charge fluctuations

Equation of state and electric (isospin) charge fluctuations are studied for matter composed of interacting pions. The pion matter is described by self interacting scalar fields via a $\phi^4-\phi^6$ type Lagrangian. The mean-field approximation is used, and interaction parameters are fixed by fitting lattice QCD results on the isospin density as a function of the isospin chemical potential at zero temperature. Two scenarios for fixing the model parameters -- with and without the first order phase transition -- are considered, both yielding a satisfactory description of the lattice data. Thermodynamic functions and isospin charge fluctuations are studied and systematically compared for these two scenarios, yielding qualitative differences in the behavior of isospin charge susceptibilities. These differences can be probed by lattice simulations at temperatures $T \lesssim 100$ MeV.

In the present work we study the BEC phenomenon in strongly interacting QCD matter. The effective lowenergy degrees of freedom in QCD are pions -the three pseudo-Goldstone bosons in the confined phase. The pions obey the Bose-Einstein statistics, thus an emergence of the BEC of pions is possible and has been predicted to occur at large isospin chemical potentials, both in effective QCD theories [26,27] and in first-principle lattice QCD simulations [28,29]. In nature, the pion BEC may occur during the cooling of the early Universe [30], in the gravitationally bound pion stars [29,31,32], or as a nonequilibrium phenomenon in heavy-ion collisions [8,9,33]. The hypothetical boson stars [34][35][36] may exist and can be a candidate for the dark matter in the Universe [37][38][39][40][41][42].
Different effective QCD descriptions of the phase diagram of interacting pion matter with a BEC include chiral perturbation theory [43,44], Nambu-Jona-Lasinio model [45], Polyakov-loop extended quark meson model [46,47] etc. Recently, a possibility of the BEC in the pion system at zero chemical potential was considered within a Skyrme-like model including both attractive and repulsive interaction terms [14,48,49]. Effects of repulsive interactions on the BEC of pions were studied in Ref. [50] at non-zero chemical potential. The system of pions at zero chemical potential was described [14] by an effective Lagrangian with the attractive (φ 4 ) and repulsive (φ 6 ) terms of a scalar field φ. In the present paper we extend this model to the finite isospin 1 chemical potential µ I .
The phase diagram on the whole plane of isospin chemical potential µ I and temperature T is investigated. Most macroscopic systems with both repulsive and attractive interactions between constituents display the first order liquid-gas phase transition (FOPT) which is ended by the critical point (CP). Therefore, these phenomena can also be expected for the interacting pions in addition to the BEC.
Lattice QCD results support an existence of the pion BEC at finite isospin chemical potential [28]. We use the recent lattice data at zero temperature to fix the repulsive and attractive interaction parameters of the model. Then, thermodynamic functions and electric (isospin) charge fluctuations up to the fourth order are calculated in the (µ I , T )-plane. Two different scenarios are employed and systematically compared. The first one includes only the repulsive interactions via the φ 6 term, but not the attractive φ 4 term. In this case no FOPT transition is observed, only the BEC transition. The second possibility takes into account both the repulsive and attractive pion-pion interactions. In this case the FOPT is observed at small T and generates a non-trivial interplay between the FOPT and BEC transitions on the phase diagram. The measures of the isospin charge fluctuations -scaled variance, skewness, and kurtosis -appear to be very sensitive to a presence of the CP and BEC phenomena. They are used to differentiate these two scenarios.
The paper is organized as follows: the theoretical description of interacting pion system is presented in Sec. II. The two choices of the interaction potential, the fixing of the model parameters, and the resulting phase diagrams are discussed. Section III is dedicated to fluctuations of the isospin charge, in particular the scaled variance, skewness, and kurtosis are discussed in some detail. Summary in Sec. IV closes the paper.

A. Model formulation
The three pions species, (π + , π − , π 0 ), are represented as a triplet of interacting pseudo-scalar fields φ = (φ 1 , φ 2 , φ 3 ) that are described by an effective relativistic Lagrangian density: where m π is the vacuum pion mass 2 and L int is the interaction part of the Lagrangian. We omit here the electromagnetic interactions. Consider now this system in statistical equilibrium within the grand canonical ensemble (GCE). The independent variables are the temperature T and the isospin chemical potential µ I . The isospin chemical potential couples to the conserved isospin charge, the pion species π + , π − , and π 0 carry the isospin charges of +1, -1, and 0, respectively. To proceed, we apply a relativistic mean-field approximation, i.e. series L int in terms of δσ = φ 2 − σ, where σ = φ 2 is the expectation value of the scalar field and ... denotes the GCE averaging. The effective mean-field Lagrangian can then be represented as [14]: where M (σ) is the effective pion mass and p ex (σ) is the so-called excess pressure, (3) 2 We use the natural units, = c = k = 1, and assume equal masses of all three pion species, mπ = 140 MeV.
The effective Lagrangian form of Eq. (2) implies that the main effect of interactions in our description leads to an appearance of a medium-dependent effective mass M (σ). The excess pressure p ex (σ) -the second term in the right hand side of Eq. (2) -ensures the proper counting of the interaction energy.
The details of the model formulation can be found in Refs. [14,51]. This model was previously used to describe the pion system at zero chemical potential [14] and the system of interacting alpha particles [51]. In the present study we apply this model to the new physical situation and consider the pion system at non-negative values of the isospin chemical potential µ I ≥ 0. The results at µ I ≤ 0 can then be obtained by interchanging π + and π − . Values of µ I > 0 correspond to positive values of the isospin charge density n I ≡ n + − n − > 0, where n + and n − correspond to π + and π − particle number densities, respectively. We consider the possible BEC of the positively charged pions in this regime. The expectation value of the scalar field σ is presented as (see Refs. [14,51] for the derivation details) where σ th i correspond to the contributions of the thermal pions, i = (+, −, 0 ), while σ bc + corresponds to a possible contribution of the Bose condensate (BC) of π + . Here µ + = µ I , µ − = −µ I , µ 0 = 0, and We will use a Skyrme-like parameterization of the interaction term: In Eq. (7), a ≥ 0 and b > 0 are model parameters which define the strength of, respectively, attractive and repulsive interactions between particles. The effective mass and the excess pressure for this choice of the interaction terms are equal to: Inverting Eq. (8) with respect to σ we obtain 3 At given T and µ I we use the system of self-consistent Eqs. (4) and (8) to determine σ and M . The pressure p, the number densities of thermal pions n th i , and the isospin charge density n I can be calculated as Here n i are the total number densities of pions. The n + density may include a contribution from a Bose condensate (BC). The condensation does not occur if µ I < M .
In the case µ I < M only the thermal pions contribute to the total number densities, The BEC of π + occurs when their chemical potential µ + ≡ µ I reaches the value of the effective mass, 4 i.e. µ I = M . In this case the number density n + may receives a contribution n bc + from the BC: The number of densities of the two other pions species are unchanged: n − = n th − , n 0 = n th 0 . The number density of π + in a condensate reads In Eq. (16) the quantities σ and σ th i are calculated according to Eq. (10) and Eq. (5), respectively.
An onset of the BEC takes place when µ I reaches the effective mass M . This condition defines a line in the phase diagram -the BEC line. This line is calculated by substituting M → µ I and σ bc + → 0 into the system of equations (4), (8), and (12) and solving it with respect to µ I at given value of T .
B. Fixing the parameters using lattice data at zero temperature Lattice QCD simulations at finite isospin density provide constraints on the equation of state from first principles [28]. In particular, the isospin density n I (T = 0, µ I ) at zero temperature has been presented in [29]. Here we will use these lattice data to constrain the parameters of our model.
In the limit of zero temperature, T = 0, the thermal pion excitations are absent, i.e. all thermal densities (12) vanish. In this case the system consists solely of the BC of π + -mesons, thus, the isospin density coincides with the number density of the condensed pions, n I (T = 0, µ I ) = n bc + , and the total pressure equals to the excess pressure, p(T = 0, µ I ) = p ex . The explicit expression for n I (T = 0, µ I ) in the the considered model follows from Eqs. (5), (10), and (16): . (17) 1

. Scenario I: Repulsive interactions only
In the first scenario we consider purely repulsive interactions between pions. To achieve this we set a = 0 and b > 0. 5 Equation (17) in this case is reduced to An onset of the BEC occurs at µ I = m π . The isospin density is a continuous function of µ I since n I (T = 0, µ I = m π ) = 0. On the other hand, the µ I -derivative of n I exhibits a discontinuity at µ I = m π . Therefore, the transition between vaccuum and a pion-condensed phase at µ I = m π is a second-order phase transition at T = 0. Qualitatively, this is consistent with predictions of many different theories, including for instance chiral perturbation theory [26,43] or Polyakov-loop extended quark meson model [46,47].
To fix the value of the parameter b we fit the lattice QCD data on n I (T = 0, µ I ) of Ref. [29] in the range of chemical potentials µ I /m π < 2. We obtain b 9.09/m 2 π 5 Another option would be to set b = 0 and take a < 0. The results in such a case are qualitatively similar to a = 0 and b > 0. with χ 2 /dof 1.62. A comparison with the lattice data is shown in Fig. 1 (a) by blue dashed line.

Scenario II: Repulsion + attraction
Let us turn now to the more general case when both the attractive and repulsive interactions are present, a > 0 and b > 0. In this case the system undergoes a first-order phase transition between vacuum (the gaseous phase) and a pion-condensed phase (the liquid phase), with the coexistence point being characterized by a vanishing pressure. To see this consider Eq. (9): this equation have two solutions for p ex = p(T = 0, µ I ) = 0, defining the expectation values of the scalar field at the FOPT boundaries, σ g = 0 and σ l = 3a/4b. This corresponds, via Eq. (8), to the following values of the effective mass in the gaseous and liquid components: The FOPT takes place at µ 0 = M l . The gaseous phase at T = 0 corresponds to the vacuum, thus, n g = 0. The isospin density jumps at µ I = µ 0 from n g = 0 to To fix the numerical values of a and b we again fit the lattice data on n I at T = 0. We obtain a 0.93 and b 11.39/m 2 π with a fit quality of χ 2 /dof 1.35 -a slightly better fit compared to Scenario I.
As discussed above, Scenario II predicts the FOPT. Using the numerical values of the a and b parameters fitted to the lattice data, one obtains the FOPT at µ 0 ≈ 0.993 m π , where the isospin density jumps from n g = 0 to n l > 0. The values of the π + density n I and the binding energy per particle W at T = 0 and µ I = µ 0 , obtained in Scenario II correspond to the ground state of the pion matter. This density is about 7 times smaller than the normal nuclear matter density of n 0 = 0.16 fm −3 , and the binding energy is about 16 times smaller than that in the nuclear ground state.
The behavior of n I in Scenario II at zero temperature is shown in Fig. 1 (a) by the solid red line. A comparison with the lattice data and the predictions of Scenario I are also shown. Overall, the behavior of n I in the both scenarios is similar. Even though a nature of the phase transition differs between the two scenarios, due to the small latent heat of the FOPT in Scenario II it is difficult to distinguish it from the second-order phase transition in Scenario I using the presently available lattice data. In Sec. III we discuss fluctuations as a possibility to make such a distinction.

C. Phase diagram at finite temperatures
Model calculations at finite temperatures are straightforward. Important constraints on the equation of state of pion matter can be obtained at zero chemical potentials and large temperatures, 120 ≤ T ≤ 160 MeV, where the QCD equation of state is known from lattice QCD [52,53]. In this range, the pressure and energy density are reasonably well described by the ideal hadronresonance gas (see e.g., Ref. [54]). This indicates that effects of pion interactions in this regime are small. To verify this we plot in Fig. (b) a ratio of the pressure of interacting pions to the ideal pion gas baseline (i.e., at a = 0 and b = 0) for the two scenarios. For purely repulsive interactions (scenario I) the pressure demonstrates small suppression relative to the ideal gas. If attractive effects are included (scenario II) the pressure at small temperatures is higher than that of the ideal gas of pions. However, at large T the repulsive effects become dominant and the pion pressure is again suppressed. In both scenarios the corrections to the ideal gas pressure at µ I = 0 are small, not exceeding 1%. This is not the case at non-zero values of µ I : the thermodynamics of the interacting pion gas differs drastically from that in the ideal pion gas in the (µ I , T )-region of the phase diagram where a BEC pions is formed, as discussed in the following.
In Scenario I, where the attractive pion interactions are absent (a = 0), there is no FOPT in the pion system. An onset of the π + BEC takes place when µ I reaches the value of the effective mass M . The BEC line can thus be obtained by substituting M → µ I in the system of equations (4) and (8) and solving it with respect to T . The resulting BEC line T bc (µ I ) is shown in Fig. 2 (a) by the dashed line. In the ideal gas limit one would obtain a vertical BEC line, µ I = m π . The deviation from the ideal gas behavior thus becomes evident as T is increased. This is due to large particle number densities, and thus stronger effects of interactions, as the temperature is increased. Note that in the ideal Bose gas, a region of the (µ I , T )-plane with µ I > m π is forbidden, whereas in the interacting system considered here this region is legitimate. It follows from Eq. (10) that the effective mass is always larger than the vacuum mass, M (T, µ I ) > m π , the pure repulsion scenario (a = 0). The (n I , T ) phase diagram in the a = 0 scenario is shown in Fig. 2 (c). The thermodynamic states below the dashed lines in Figs. 2 (a) and (c) correspond to a non-zero density of the BC, i.e. to a macroscopic number of π + -mesons occupying the zero momentum level k = 0.
In Scenario II, with both the repulsive (b > 0) and attractive (a > 0) pion interactions present, the FOPT phase transition takes place in addition to the BEC formation. The (µ I , T ) and (n I , T ) planes are presented for this scenario in Figs. 2 (b) and (d), respectively. The line of the FOPT is shown by a thick solid line in Fig. 2 (b). This line ends in a critical point (CP) at T = T c ≈ 0.369 m π , µ I = µ c ≈ 0.991 m π , and n I = n c ≈ 0.051 m 3 π , which is shown by the green star. The µ c value can be expressed explicitly in terms of model parameters An approximate analytical dependence of T c on a, b, and m π can also be obtained Here ζ(x) is the Riemann zeta function. The relation (24) has been obtained assuming n − n 0 n + as well as the non-relativistic approximation in vicinity of the CP. Using the previously obtained parameters a and b from fitting the lattice data one obtains T c ≈ 0.393 m π . This is within 6% of the numerical result obtained without approximations. Note that the limit a → 0 corresponds to T c → 0 and µ c → m π .
At the FOPT line in the (µ I , T )-plane the pressures of the gaseous and liquid phases are equal to each other. On the other hand, the isospin charge density n I has a discontinuity. The mixed phase shown in Fig. 2 (d) is bounded by the gas-like (left) and liquid-like (right) binodals presented by solid lines that intersect each other at the CP. The pion states inside the mixed phase correspond to linear combinations of the diluted (gaseous) and dense (liquid) states lying on the left and right binodals, respectively. The liquid component of the mixed phase always lies below the BEC line, thus it always contains a non-zero fraction of condensed π + -mesons. The gaseous component, on the other hand, does not contain the BEC.
A remarkable feature of the considered model is that the BEC line enters the mixed phase at the CP. This property of the model is robust with regard to variations in the values of the a and b parameters. Another peculiar property is the non-smooth intersection of the left and right binodals at the CP.
Scenarios I and II provide a similar picture of the phase diagram at T T c . At T T c , on the other hand, the differences are significant. We argue that these differences can be most clearly seen by studying the behav- ior of isospin charge fluctuations. This is discussed in Sec. III.

III. FLUCTUATIONS
The two presented descriptions of the isospin charge density from the lattice results at T = 0 both contain the BEC. Within the second description, the FOPT at T < T c leads to the n I discontinuity. We argue that the difference between the two scenarios can be probed by considering isospin charge fluctuations.
In the GCE, the j-th order susceptibility of the isospin charge is determined by a j-th order partial derivative of the pressure p with respect to the chemical potential µ I : Ratios of susceptibilities given by (25) can be particularly useful as such quantities are intensive in the thermodynamic limit. Some of the most well known such quantities include the scaled variance ω, skewness Sσ, and kurtosis κσ 2 (see, e.g., Ref. [55]): Using Eq. (25) together with Eq. (11) the scaled variance ω can be written as In the ideal pion gas the scaled variance diverges at the BEC line [8], i.e., ω id → ∞ at µ I → m π − 0: Due to the repulsive interactions in the considered model the scaled variance remains finite. On the BEC line, M = µ I , one finds: where µ c is given by Eq. (23). The value of ω remains also finite in a presence of the BC, n bc + > 0. In Scenario II (a > 0), ω exhibits singular behavior at the CP, T = T c , µ I = µ c , where it diverges. A systematic expansion of the thermodynamic functions in a vicinity of the CP allows to obtain the critical exponents. We expect that the critical exponents of the considered system are different from those in the mean-field class universality. This is due to a presence of the two order parameters, n l I −n g I > 0 and n bc + > 0, which disappear simultaneously at the CP (see, e.g., Ref. [56]). A detailed discussion of this subject is however outside of the scope of the present study.
Scaled variance. The behavior of the scaled variance ω in the plane of temperature and isospin chemical po- tential is shown in Fig. 3. In Scenario I (pure repulsion), ω is a continuous function, in particular across the BEC boundary [see Fig. 3 (a)]. In Scenario II (full potential), on the other hand, ω exhibits a jump discontinuity over the FOPT line an becomes divergent at the CP [ Fig. 3 (b)]. It is still a continuous function across the BEC line, however. Note that in the limit a → 0 (Scenario I), the CP approaches T c → 0 and µ c → m π . Therefore, the point µ I = m π at zero temperature retains some of the proper tires of the CP and exhibits large fluctuations in its vicin-ity. One can observe ω of any magnitude in the vicinity of this point, the exact magnitude depending on the path of approach. In particular, approaching this point along the BEC line one finds ω → ∞ at T → 0.
Skewness. The skewness, Sσ, for Scenario I (a = 0) and II (a > 0) is shown in Figs. 4 (a) and (b), respectively. At small µ I m π values, where the pion densities are small, both the pion interactions and Bose statistics effects can be neglected, thus, Sσ ≈ 1. The skewness attains positive values in those regions of the phase diagram where there is no BC. Sσ is discontinuous along the BEC line, jumping from positive values outside the BEC phase to negative values in the phase with a BC. The above observations are valid for both scenarios. In Scenario II (a > 0) Sσ shows singular behavior at the CP. The skewness can reach both −∞ and +∞ at the CP depending on the path of approach. When crossing the FOPT in Scenario II Sσ undergoes a jump discontinuity.
Kurtosis. κσ 2 presented in Fig. 4. In both the scenarios it also always attains positive values everywhere on the phase diagram. Kurtosis can strongly deviate from the baseline κσ 2 = 1 of an ideal Boltzmann gas. This is due to the presence of interactions and Bose statistics.
The largest values of the kurtosis are generally obtained in the vicinity of the BEC line. The kurtosis exhibits a non-monotonic behavior as a function of µ I at both the BEC-line and the FOPT-line, where it jumps down as µ I is increased. κσ 2 is an increasing function of µ I elsewhere on the phase diagram. The values of κσ 2 remain large even far away from the CP and the Bose condensation boundary. This is due to its large sensitivity to interactions in the system. κσ 2 diverges at the CP. The model does not predict negative values of κσ 2 anywhere on the phase diagram. This is in contrast to the universal behavior of fluctuations in the Ising model [57,58], as well as various model calculations [59][60][61][62][63], where negative values of κσ 2 are observed in the so-called analytic crossover region above the critical temperature. In the present work the negative values of κσ 2 are not observed because of the Bose-Einstein condensation. The BECline, which itself corresponds to a phase transition of a higher order, crosses the CP, thus no region in the vicinity of the CP can be identified as an analytic crossover.
We would like note that κσ 2 exhibits a singular behavior also in Scenario I (a = 0), at a point (T = 0, µ I = m π ) [see Fig. 4 (c)]. Approaching this point along the BEC line one finds κσ 2 → ∞ at T → 0.
In the present work we do not discuss the behavior of fluctuations inside the mixed phase of the FOPT. These fluctuations can be addressed using the method developed in Ref. [64] and will be the subject of a future study.

IV. SUMMARY
We studied thermodynamic properties of interacting pion matter in the framework of a mean-field model with a φ 4 -φ 6 type Lagrangian. The phase structure has been studied at non-zero isospin chemical potential µ I that corresponds to the conserved 3rd component of isospin. Parameters of the repulsive and attractive interactions were fixed using lattice QCD data on the isospin density as the function of the chemical potential µ I at zero temperature. The lattice data can be reasonably fitted within the two qualitatively different scenarios: Scenario I with only repulsive interactions, and Scenario II with both repulsive and attractive interactions. In both scenarios a phase with a Bose condensate of π + pions was found to occur at sufficiently large µ I , the transition between ordinary pion matter and matter with a BC taking place along the so-called BEC lines. The presence of the attractive interactions in Scenario II leads, in addition to the BEC, also to a first order liquid-gas phase transition of pions with a CP at T c ≈ 0.369 m π and µ c ≈ 0.991 m π . A notable qualitative feature of the model, present for a broad range of values of parameters a and b. is the fact that the BEC line merges with the FOPT line at the CP. The system is characterized by two order parameters: (i) the difference n l − n g > 0 between the liquid and gas phase densities the corresponds to the liquid-gas transition; (ii) the density n bc + of the Bose condensed π + pions that characterizes the BEC transition. This makes the model qualitatively different from the usual systems with a CP and FOPT where only a single order parameter is present.
The susceptibilities of isospin charge fluctuations up to the 4-th order studied in the paper can serve as a robust observable to distinguish between the two different scenarios. In the both scenarios, the scaled variance ω = χ 2 /χ 1 , skewness Sσ = χ 3 /χ 2 , and kurtosis κσ 2 = χ 4 /χ 2 remain finite on the BEC line. This happens due to the repulsive interactions in the pion system, in contrast to the ideal pion gas where these measures become infinite on the BEC line. All three fluctuation measures demonstrate anomalous properties approaching the CP: ω → ∞, κσ 2 → ∞, and Sσ can reach both +∞ and −∞ depending on a path to the CP. Note a significance of the higher order susceptibilities (e.g., skewness and kurtosis fluctuation measures) which are highly sensitive to a presence of the CP. In the scenario I the CP is absent. In this case the anomalous fluctuations take place in the point T = 0 and µ I = m π . Approaching this point: ω and κσ 2 can reach any value from 0 to ∞, and Sσ can reach any value between −∞ and +∞ depending on a path to the (µ I = m π , T = 0)-point. The susceptibilities χ i can be computed in lattice QCD which are free of sign problem at finite µ I . Analysis of their behavior can be used to establish a point (region) of anomalously large fluctuations. Determining whether this point corresponds to zero or finite temperatures will allow to distinguish Scenarios I (T = 0) and II (T = T c > 0). Besides, there is a qualitative difference in a behavior of the scaled variance ω near T = 0 in Scenario I and near T = T c in Scenario II: ω → ∞ at T → T c in Scenario II, and ω can reach any value from 0 to ∞ depending on the path of approaching to T = 0 in Scenario I.
The results obtained in this paper can be used in systems where the pion densities are large and a BC of pions may occur. This can happen, for example, in heavy-ion or proton-proton collisions where the pion condensation may occur as a chemical non-equilibrium effect. Other possibilities include pion stars as well as the Early Universe which may have passed through a pion-condensed phase if the lepton flavor asymmetries during its evolution were large.