Multicomponent meson superfluids in chiral perturbation theory

We show that the multicomponent meson systems can be described by chiral perturbation theory. We chiefly focus on a system of two pion gases at different isospin chemical potential, deriving the general expression of the chiral Lagrangian, the ground state properties and the spectrum of the low-energy excitations. We consider two different kinds of interactions between the two meson gases: one which does not lock the two chiral symmetry groups and one which does lock them. The former is a kind of interaction that has already been discussed in mutlicomponent superfluids. The latter is perhaps more interesting, because seems to be related to an instability. Although the pressure of the system does not show any instability, we find that for sufficiently strong locking, the spectrum of one Bogolyubov mode becomes tachyonic. This unstable branch seems to indicate a transition to an inhomogeneous phase.


I. INTRODUCTION
Cold hadronic matter is an interesting playground for a deep understanding of the properties of the strong interaction. At asymptotic baryonic densities the liberated quarks [1] should pair forming a color superconductor, see [2][3][4] for reviews. At large isospin densities a different kind of collective phenomenon happens, with mesons forming a Bose-Einstein condensate (BEC) [5][6][7][8][9][10]. In general, the matter density of the system is controlled by the baryonic chemical potential, µ B , while the isospin chemical potential, µ I , is associated to its degree of isospin asymmetry, e. g. indicating that the number of neutrons differs from the number of protons. The properties of matter as a function of µ I have been the subject of intensive investigation for a number of reasons. Systems with large isospin asymmetry exist in Nature; in particular neutron stars [11] are believed to be compact stellar objects with a large isospin asymmetry. Recently, the possible existence of pion stars has also been proposed [12][13][14]. Regarding the microscopic properties of matter, the inclusion of µ I can lead to a better understanding of quantum chromodynamics (QCD) in a regime in which lattice QCD simulation are doable [15][16][17][18][19][20][21][22]. Remarkably, the lattice QCD simulations of meson gases with vanishing baryonic density are not affected by the sign problem and can be implemented for not too high values of µ I . These simulations are steadily improving [20][21][22], reaching increasingly precise results on the thermodynamic properties of the system and thus offering powerful tests for the methods and models developed for the effective description of the strong interaction.
Among the various proposed models, it is worth mentioning the Nambu-Jona Lasinio (NJL) model [23][24][25][26][27][28][29][30][31][32] and the quarkmeson model [24,[33][34][35], which can be used in a wide range of values of µ I . Although these models are useful tools for exploring the properties of hadronic matter, they are based on a number of parameters that have to be phenomenologically fixed. Thus, they lead to results which depend on the choice * llepori81@gmail.com † correspondence at: massimo@lngs.infn.it of these parameters and on the number of degrees of freedom used. Moreover, the obtained results cannot be systematically improved because no expansion parameter can be identified. A systematic analysis of hadronic matter can be obtained by effective field theories [36][37][38], which are based on an expansion in a control parameter. Here we focus on chiral perturbation theory (χPT), which is an effective theory designed to describe the low-energy properties of QCD [39][40][41][42][43]. The χPT Lagrangian is derived by the global symmetries of QCD, basically integrating out the high-energy part. The effect of the isospin chemical potential is conveniently included in covariant derivatives, see [40,41,44] for a general discussion. This approach leads to systematic results, which can be improved including higher orders in the χPT expansion [41,45].
The thermodynamic and low-energy properties of mesons at nonvanishing µ I have been studied using the χPT in many different works [12,[46][47][48][49]. In particular, it has been confirmed that the pion condensed phase first discussed in [5][6][7][8][9] sets in at µ I = m π , where m π is the pion mass. Remarkably, χPT can also be used to study different gauge theories with isospin asymmetry, including 2 color QCD with different flavors [50][51][52][53][54][55][56] In the present paper we study the multicomponent meson systems in which each component is characterized by a global symmetry. In general, for each component, the spontaneous breaking of a global symmetry should lead to the formation of a superfluid. Multicomponent superfluids can be realized in He 3 -He 4 mixtures, see [57,58] or in ultracold atoms experiments [59][60][61][62][63]. In the compact star interior neutrons and protons are believed to simultaneously condense [11] and if deconfined quark matter is formed, the color-flavor locked phase [64] supplemented by kaon condensation [65,66] is a phase with two bosonic superfluids. Here we examine the effect of the possible intra-species interactions on multicomponent superfluidity. We focus on the meson condensed phase, employing the χPT framework for deriving the relevant lowenergy Lagrangian. We identify two very different types of interactions: those that lock the two global symmetries and those that do not lock them. Remarkably, at the leading order (LO) in χPT, only the former type of interactions are possible. This kind of interaction is not typically considered in ultracold gases, because in these systems the number of parti-cles of the two species are separately conserved. In our work, we assume that this not the case and we find that the strength of the locking term plays a prominent role. Increasing the locking, we obtain that the transition to the broken phase is favored. Moreover, for sufficiently large couplings the system becomes unstable. Analyzing the dispersion laws of the lowenergy degrees of freedom, we find that the instability can be interpreted as a transition toward an inhomogeneous phase.
Including the next-to-leading order (NLO) χPT corrections, it is possible to include interactions that do not lock the two chiral groups. This type of interaction is akin to the one typically discussed in ultracold atoms systems and indeed in this case we obtain results similar to those of multicomponent Bose gas [67].
The present paper is organized as follows. In Sec. II we report known results for meson systems in χPT. This is useful to fix the notation and for comparison with the multicomponent meson system. In Sec. III we generalize the χPT Lagrangian to two meson gases, introducing the leading interaction terms. In Sec. IV we analyze the effect of one of the possible interaction term leading to chiral locking. In Sec. V we consider the χPT term that does not lock the two chiral groups. We conclude in Sec. VI. A number of results are collected in the Appendices. In the Appendix A, we report the low-energy excitations of a single-component pion gas. In the Appendix B, we discuss the low-energy corrections to the mean-field thermodynamic quantities arising from the vacuum energy of the Bogolyubov modes.

II. SINGLE MESON GAS
The χPT description of the single meson gas is based on the global symmetries of massless QCD, with N f the number of flavors. The meson fields are collected in the Σ field, transforming under G as where L ∈ SU (N f ) L and R ∈ SU (N f ) R . The leading O(p 2 ) χPT Lorenz-invariant Lagrangian [41,42,47] is given by where the mass matrix, M , and the so-called pion decay constant, f π , are the low energy constants (LECs) that cannot be fixed by the symmetry group G and must be determined in some other way. The χPT Lagrangian is constructed assuming that the mass term does not break the global symmetries, thus that M transforms as Σ. Then, the locking of the chiral rotations to the vector SU (N f ) V group is induced by the vev of M , see for example the discussion in [39,42]. The covariant derivative in Eq. (3) allows us to take into account the coupling of the meson fields with the gauge fields and/or with external currents and/or the effect of different chemical potentials [40,41,44]. In the present work, we will only consider the effect of the isospin chemical potential and we will restrict the analysis to pions, corresponding to the N f = 2 case. Thus, we consider the covariant derivative where the isospin chemical potential, µ I , is introduced as the time component of a vector field. Note that the covariant derivative does not include the baryonic chemical potential, µ B , because mesons do not have a baryonic charge. A useful parameterization is where the radial field, ρ, and the unit vector field,φ, encode in a nontrivial way the three pion fields. By this parameterization, the LO χPT low-energy Lagrangian takes the form obtained in [12] where is the potential and the control parameter is γ = µ I /m π . For |γ| > 1, the pion condensed phase is favored [5-10, 46, 47] and in the present parametrization it corresponds to a radial field vev,ρ, satisfying Therefore, in the broken phase the meson field vev is given bȳ where n is a unit vector associated to the residual O(2) symmetry of the vacuum. The pressure and the isospin number density in the broken phase are respectively given by [46,47,49] leading to the O(p 2 ) equation of state [49] ǫ(P ) = −P + 2 P (2f 2 π m 2 π + P ) .
Close to the phase transition point, γ 1, the system is dilute and it is possible to expand the pressure P and the energy density ǫ, as a function of the isospin number density n I . If we define the adimensional isospin density n = n I /(f 2 π m π ), we can expand the control parameter as which is meaningful expansion for n ≪ 1. The pressure can then be expanded as follows where the leading term is the mean-field expression of the pressure of a boson system with coupling g 0 = 1/4f 2 π . This is indeed the correct expression of the coupling close to the phase transition, see Eq. (B2) and the discussion in the Appendix A. The energy density is instead given by ǫ = m π n I + g 0 n 2 which takes into account the energy associated to the mass of the pions. Note that the above expressions are obtained in the mean-field approximation, meaning that the low-energy fluctuations have not been included. Indeed, the order n 3 I corrections are determined by the χPT Lagrangian and not by the contribution of the Bogolyubov modes. The vacuum contribution of the Bogolyubov modes is considered in the Appendix B, and is much smaller than the leading mean-field contribution. However, it can play an important role in a multicomponent gas, as we will see below.

III. SYSTEM OF TWO MESON GASES
We now generalize the discussion of the previous Section to a system with two mesonic gases. In the second quantization formalism we assume that two meson systems with densities n 1 and n 2 are described by the fields Σ 1 and Σ 2 . As for the single meson gas discussed in the previous Section, we use the global symmetries for constructing the χPT Lagrangian.
As a starting point we consider the noninteracting case with symmetry group where is the chiral group of the Σ a field. For simplicity, we will mainly treat the system in which the two meson gases correspond to two fictitious pion systems, paving the way for the discussion of the simultaneous condensation of pions and kaons. In other words, we assume that in the noninteracting case the fields Σ 1 and Σ 2 transform independently under two chiral groups as where L a ∈ SU (2) L,a and R a ∈ SU (2) R,a with a = 1, 2.
The most general O(p 2 ) chiral Lagrangian invariant under these symmetries is where f 1π and f 2π , as well as the matrices M 1 and M 2 , are the low energy constants (LECs) of the system. As for a single meson system described by the Lagrangian in Eq. (3), we have constructed this Lagrangian assuming that the mass terms do not break the global symmetries, which means that M a transforms as Σ a . The covariant derivative D a ν takes into account the interaction of the mesons of the a system with the external fields. If the two meson systems have different isospin chemical potentials, µ 1 and µ 2 , respectively, this can be encoded in the two covariant derivatives for a = 1, 2.
We now introduce the interaction between the two gases. Before doing that, let us first recall that under G a the covariant derivative transforms as the Σ a fields, that is and therefore the two covariant derivatives are independently rotated. Let us now consider the possible interaction terms. If we add to the noninteracting Lagrangian the term it locks the two chiral groups, leaving only the diagonal chiral rotation unbroken. In principle, the k coefficient is a number that depends on the interaction strength between the two chiral fields and, as any LEC, it is independent of the isospin chemical potentials.
Remarkably, the interaction Lagrangian in Eq. (21) is the only O(p 2 ) meaningful coupling leaving the G D group unbroken. One may think to add a Lagrangian term of the type which indeed locks chiral rotations. However, if one of the two fields vanishes, from Eq. (5) we have that say Σ 1 ≡ I. Then the term in Eq. (23) acts as a mass term for the Σ 2 field, breaking G 2 down to the vector subgroup. Therefore, this kind of term or any term of the type with n > 0 is not allowed. For a similar reason the mass-like terms are not allowed, unless n = 0.
If one wants to preserve the G group, then one has to consider the O(p 4 ) terms. At this order, there are only two derivative terms coupling the two meson systems that do not lock the two chiral groups: whereL 1 andL 2 are two LECs analogous to the standard L 1 and L 2 of O(p 4 ) χPT [42]. When including these contributions, one should consistently include the standard O(p 4 ) chiral terms, as well. However, as was shown in [12], the effect of the standard NLO terms on the thermodynamic properties of the system is extremely small and can be accounted for by a renormalization of the LO LECs.
As an aside, we note that in principle one may consider more complicated intra-species interaction terms, like with k µν a Lorentz tensor and a G singlet. This kind of interaction term somehow generalizes Eq. (21) and Eq. (26), however it is not obvious how to fix the values of the k µν components in general.
In the following, we will discuss the effect of the interaction terms in Eq. (21) and in Eq. (26), separately, focusing on the pion system.

IV. CHIRAL LOCKING
To gain insight on the system described by Eqs. (18) and (21), let us first assume that we are making a partition of an ensemble of undistinguishable pions, dividing it in two (interacting) subsets, in such a way that the Σ 1 field describes the pions of the first subset and Σ 2 field the pions of the second subset. Let us first focus on the kinetic terms at vanishing isospin chemical potentials. Since the pions are indistinguishable, one may naively think that the most general O(p 2 ) Lagrangian is where the first term, respectively the second term, describes the propagation and self-interactions of the fields of the subset 1, respectively 2. The third term mixes the two fields and induces the locking between the two subsets. If it were absent, that is for k = 0, there would be no interactions between the two sets. For subsets made of identical particles there must exist a way of reshuffling them. Since Σ 1 Σ † 1 +Σ 2 Σ † 2 = 2, any reshuffling can only correspond to a rotation transforming the Lagrangian in Eq. (28) in  21); k = 0 corresponds to two noninteracting systems, while k = 1 to a system of one single type of particles. For k > 1 the system is unstable.
To maintain the Lagrangian invariant we have to take k = 0 or, more interestingly, k = 1. Indeed, in the latter case wheref 2 1π = f 2 π (1 − sin 2θ),f 2 2π = f 2 π (1 + sin 2θ), and thereforef 1πf2π = f 2 π cos 2θ. Note that one cannot identifŷ f aπ with the pion decay constant of the pions in the subset a, because the fields are mixed by the interaction terms.
If one takes k = 1, the O(2) symmetry in Eq. (29) does not hold and the coefficient of the interaction term cannot be expressed asf 1πf2π , meaning that if one makes the rotation, this term would depend on the rotation angle. In the Lagrangian in Eq. (31) it is possible to eliminate the dependence on the unphysical angle θ in the quadratic terms by writinĝ which is a generalization of the standard nonlinear expression of the pion fields. Therefore, the expression in Eq. (31), where k = 1 is set, is the most general χPT Lagrangian for two gases of undistinguishable pions. We can easily generalize it to N undistinguishable pion gas, writing where f a are a generalization of the pion decay constant. Note that the propagating degrees of freedom are obtained by diagonalizing the quadratic Lagrangian. Including the mass terms, formally considering the vevs of the fields M a in Eq. (18), we can write the total Lagrangian of the system as follows where we have assumed that the two fields have different mass parameters, m 1π and m 2π . These parameters have to be interpreted as LECs for the coupled system and correspond to the pion masses only in the k = 0 case. The actual masses can be obtained by the dispersion laws where the masses are given by for k = 1, and equal to the "reduced mass" for k = 1. From the above expressions it is clear that the interaction term in Eq. (21) induces a mass splitting. For clarity we report the behavior of the meson masses as a function of k in Fig. 1. We remind that k = 0 corresponds to two noninteracting gases, while k = 1 corresponds to two identical pion gases. For k < 1, the mass splitting induced by the locking term is similar to the one induced by µ I between the charged pions, see for example [48]. However, the system is unstable for k > 1. The instability is signaled by the divergent mass of one mode as k → 1 − , which becomes imaginary for k > 1. In the context of ultracold atoms physics, where boson condensates are mostly considered, the latter feature is generally related to the appearance of spatially inhomogeneous phases, see e.g. [68,69] and references therein. We stress, however, that here we are in the presence of a completely different instability. Indeed, in ultracold atoms, the instability is triggered by a sufficiently large coupling between the two systems [68,70,71] (a similar phenomenon is known also for fermions, called Stoner instability, see e.g. [72]). Instead, in the present case, the locking plays the game: indeed, as k varies, the repulsion from the locking term remains fixed and reads In spite of this relevant difference and considering that the locked theory in Eq. (34) is quadratic, it is still quite natural to postulate that the same theory with k > 1 cannot exist with the two species coexisting in the same space domain.
To elucidate the mechanism underlying the locking instability, and its possible resolution, let us consider a simple system consisting of two scalar bosons with a locking term with a manifest discrete Z 2 × Z 2 symmetry for k = 0. This symmetry corresponds to the transformations φ 1 → −φ 1 and φ 2 → −φ 2 , separately. For k = 0 the two discrete symmetries are locked, with the only remaining Z 2 symmetry corresponding to φ 1 → −φ 1 and φ 2 → −φ 2 , simultaneously. This simple system becomes unstable for k > 1, because one of the two eigenmodes has an imaginary mass. One possible solution of the instability corresponds to the realization of an inhomogeneous phase. Let us give an heuristic argument in favor of the inhomogeneous phase. If we assume that one component is realized in the volume V 1 and the other in the volume V 2 , then the action can be written as where S a with a = 1, 2 are the actions of the free scalar fields. The effect of the interaction term is only relevant at the interface, S 12 , of the two volumes. In other words, in the inhomogeneous phase the interaction Lagrangian L int has only support at the interface and therefore the dispersion laws of the field φ 1 , respectively φ 2 , in the volumes V 1 , respectively V 2 , are not tachyonic.

A. Two pion gases at different isospin chemical potentials
We now consider the effect of the isospin chemical potentials for the two pion gases. Including them, the Lagrangian reads where the covariant derivatives are given in Eq. (19). Since the two fields can have different vevs, we generalize Eq. (9) to Σ a = cos ρ a + in a · σ sin ρ i a = 1, 2 , where ρ a are the two radial fields and n a are two unit vectors. Upon substituting Eq. (42) in Eq. (41), we obtain the tree-level potential where γ a = µ a /m aπ and the last term on the right hand side originates from the locking term, which explicitly breaks the G symmetry to the diagonal group, G D The interesting aspect is that the potential depends on the relative angle between n 1 and n 2 . In the ground state the two unit vectors are locked to be aligned, if the isospin chemical potentials have equal signs, or anti-aligned, if the isospin chemical potentials have opposite signs. We can clearly restrict the analysis to the case in which both isospin chemical potentials are positive and aligned. Since the vevs of the two fields are not independent but tend to align, it is clear that the condensation of one field favors the condensation of the other; we will discuss this effect in detail below. From the above expression it is also clear that the system has two NGBs for k = 0, corresponding to the two independent oscillations of the unit vectors, but only one NGB for k = 0, corresponding to the locked oscillations of the two fields. The second mode is massive and corresponds to a pseudo NGB.

B. Phase diagram of the locked pion gases
At the transition to the broken phase, where both gases condense, we can expand with ǫ a ≪ 1. Upon replacing this expression in the stationary condition for the potential, we obtain signaling that the condensation of one gas is deeply related to the condensation of the other: as soon as, say, ǫ 1 > 0, it follows that ǫ 2 > 0. The formation of one superfluid necessarily makes the other gas superfluid by a simultaneous condensation (SCO) mechanism. Upon solving the above system of equations, we easily obtain that the SCO happens for corresponding to the curve, C on the (γ 1 , γ 2 ) plane depicted in Fig. 2 for various values of k. The existence of this curve makes explicit that the onset of one condensate induces the condensation of the other, a manifestation of the interaction between the two. A remarkable aspect is that the SCO happens for any nonvanishing value of k. Clearly, the larger is k, the larger is the effect of one condensate on the other. Moreover, with increasing values of k, the normal phase region shrinks. To better understand this process, let us focus on the For every considered value of k, the broken phase is the region outside the corresponding curve. It corresponds to a system in which there is the simultaneous condensation of both fluids and is indicated with SCO. The only region where the SCO does not happen is along the axes, where γ1 = 0 and γ2 > 1 or γ2 = 0 and γ1 > 1; along these lines only one component is superfluid. The analysis of the low-energy excitations shows that for k > 1 one of the low-energy modes becomes tachyonic, meaning that in this case the mean-field results reported in this figure are not valid. γ 1 = γ 2 = γ case. Since the two isospin chemical potentials are equal, it follows has thatρ 1 =ρ 2 =ρ, cosρ = 1 and the transition happens for γ 2 = 1/(k + 1). Therefore, with increasing values of k, the transition to the SCO phase happens at lower values of γ. One may naively think that increasing k would lead to a system that becomes superfluid for arbitrary values of the isospin chemical potential. As we will see below, this is not the case, because for k > 1 an instability in the low-energy spectrum is triggered. In general, close to the transition curve, C, one can expand the pressure as whereγ a ∈ C and are the susceptibilities. Upon expressing the isospin chemical potential in terms of the number densities, we obtain where the coupling constants are given by where D = L 11 L 22 − L 2 12 , with L ab > 0. It turns out that the equality corresponding to the case γ 1 = γ 2 = 1/ √ 2. For non relativistic distinguishable and dilute superfluid bosons, the equality in Eq. (52) corresponds to the stability threshold against collapse or turn into an inhomogeneous phase (depending on the sign of g 12 ) [68,70,71]. By a similar reasoning, one could expect that, because of the relation in Eq. (52), the two-pion locked system at nonvanishing isospin density is stable. More in detail, the expression in Eq. (52) relies on the mean-field approximation. Instead, in condensed matter system it is known that the inclusion of the vacuum energy contribution of the Bogolyubov modes can only turn a collapsing system into an inhomogeneous one, made of droplets of coexisting gases [71]. Anyway, in the present case the condition in Eq. (52) is not violated, the mean-field pressure is well defined, and the system could be expected to be homogeneous and stable. However, for k > 1, we found that in the normal phase there exists a tachyonic mode. It is therefore important to analyze the low-energy spectrum of the system to figure out what is the fate of the tachyonic mode in the SCO phase.

C. Low-energy excitations
The low-energy excitations of the multicomponent system can be determined studying the fluctuations of the radial com-ponent and of the Bogolyubov modes. We shall employ the same formalism developed in [12] and briefly discussed in the Appendix A, extending it to the two-component pion system.

Radial excitations
In the broken phase, the system has two radial excitations, χ 1 and χ 2 , corresponding to the fluctuations around the corresponding vevs: where it is assumed that χ a ≪ρ a . Upon substituting the above expression in Eq. (41) and restricting to the quadratic order in the fields, we obtain the Lagrangian where c 12 = cos(ρ 1 −ρ 2 ) s 12 = sin(ρ 1 −ρ 2 ) The corresponding dispersion laws are given by thus the two modes have non-negative masses and are stable for any value of k. On the transition region to the BEC phase M 12 = M 1 M 2 , and one of the radial modes becomes massless. The stability in the radial modes for any value of k is clearly a manifestation of the result obtained in the previous Section, that the pressure close to the transition region is positive defined.

Bogolyubov modes
Neglecting the radial excitations, thus taking ρ 1 ≡ρ 1 and ρ 2 ≡ρ 2 , one has the following low-energy Lagrangian where and with L 2 given by a similar expression, while (59) stems from the locking term. The unit vectors fieldsφ 1 and ϕ 2 describe the two angular fluctuations of the condensates and can be parametrized as followŝ ϕ 1 = (cos α, sin α) andφ 2 = (cos θ, sin θ) , (60) which generalize the expression in Eq. (A7). Upon substituting the above expression in the low-energy Lagrangian, we obtain where we have not included the terms and leading to interactions and total derivatives. The Lagrangian in Eq. (61) describes two coupled modes. We restrict to the case µ 1 µ 2 > 0; the other case can be treated in a similar way. The potential term is minimized for α = θ, thus expanding in (α − θ) and keeping only the quadratic terms, we obtain the dispersion laws corresponding to the massless NGB and the massive pseudo NGB, respectively. The propagation velocity of the NGB is equal to 1, however integrating out the radial oscillations would lead to a propagation velocity equal to the speed of sound, see [12] and the discussion in the Appendix A. For k = 0 the mass of the pseudo NGB vanishes and thus the system has two NGBs describing the independent fluctuations of the two decoupled superfluids. We notice that for k → 1 − the mass of the pseudo-NGB diverges and only one low-energy mode exists, which is consistent with the fact that for k = 1 the system is equivalent to a single superfluid. For k > 1 the mass of the pseudo-NGB becomes imaginary, signaling an instability. This is the same instability we previously discussed in Fig 1 in the unbroken phase. Thus, the unstable modes is still present in the SCO phase, now appearing as a pseudo NGB with a tachyonic mass. The presence of this mode indicates that the meanfield approximation breaks down. Therefore, the expression of the pressure in Eq. (50) is incorrect for k > 1. This result is discussed in more detail in the Appendix B, where it is shown that the beyond mean-field contributions are ill-defined for k > 1.

V. INDEPENDENT CHIRAL ROTATIONS
We now consider the interaction terms that do not lock the two chiral groups. Upon expanding the Lagrangian given by Eqs. (18) and (26), we obtain the potential where we have assumed the two gases have unequal masses and decay constant parameters. Unlike the locked case in the previous Sections, now the tree-level potential is independent of the the relative orientation of the two condensates, indeed it does not depend on n 1 · n 2 . In other words, the potential does not break the degeneracy of the two vacua and the two condensates vectors n 1 and n 2 can independently rotate. This is a manifestation of the fact that the interaction term does not lock the two chiral groups and thus the system has two NGBs. ConsideringL 1 +L 2 ∼ 10 −3 , as typical for O(p 4 ) corrections (see for example [42]) the interaction term has a small impact on the favored ground state. In particular, the onset of the simultaneous condensation is for γ 1 1 and γ 2 1. In the following we will consider |L 1 +L 2 | = 10 −2 − 10 −3 , also taking into account possible negative values of (L 1 +L 2 ).
In Fig. 3, we report the phase diagrams obtained with positive (left panel) and negative (right panel) values ofL 1 +L 2 . The behavior with the strength of the intra-species interaction is very similar to the one obtained for a coupled two-fluid system in [67]. TheL 1 +L 2 parameter has the same effect on the phase diagram of the entrainment parameter of [67]: a positive value ofL 1 +L 2 favors the SCO, while a negative value disfa- vors it. In [67] it was also discussed the instability generated by coupled superfluid flows. Although a similar phenomenon might emerge in our model, we postpone its analysis to future work.
In order to infer the effect of one superfluid on the other, we consider the case in which one of the two superfluids is formed, say the superfluid 2, and we seek the critical value γ 1,c for the onset of the condensation of the superfluid 1. At the leading order in the intra-species interaction, we find that the condensation onset for the first species obeys the equation which is depicted in Fig. 4 forL 1 +L 2 = 10 −3 . In principle, for large values of γ 2 it suffices a small µ 1 isospin chemical potential to drive the system 1 in the condensed phase. However, for reasonable values of the NLO LECs, the influence of one condensate on the other is extremely small. The low-energy spectrum in the broken phase consists of two NGBs which have a very small mixing. The system does not show any instability in the spectrum of the Bogolyubov modes.

VI. CONCLUSIONS
We have discussed multicomponent meson superfluids in the χPT framework. We have derived the relevant χPT Lagrangian restricting most of the analysis to the global symmetry group given in Eqs. (15) and (16) with N f = 2, corresponding to two fictitious pion gases with different masses and decay constants. In the noninteracting case, if one of the two isospin chemical potentials exceeds the corresponding pion mass the system becomes superfluid. Turning on the interactions the two condensates influence each other. We have considered two possible interaction terms, one that locks the two chiral groups and one that does not lock them.
The Lagrangian term in in Eq. (21) leads to the tree-level potential in Eq. (43), with the peculiar interaction term between the phases of the two condensates. Minimizing the potential we have obtained the phase diagram reported in Fig. 2. With increasing locking parameter k, the region in which the simultaneous condensation is realized becomes larger. It seems that one can arbitrarily shrink the normal phase region by increasing the value of k. However, the locking turns one low-energy mode becomes in a pseudo NGB with dispersion law given in Eq. (64). For k > 1 the mass of the pseudo NGB becomes imaginary and therefore an instability is triggered. The unusual aspect is that even for k > 1 the potential has a well defined minimum, indeed the low-energy radial excitations studied in Sec. A 1 have a well-defined mass. Since no other homogenous phase is energetically favored, this suggests that there exists an energetically favored inhomogeneous phase, where the two gases do not coexist any longer. Though not rigorously proved, this seems an educated guess, also because of the analysis of the simplified model discussed in Sect IV. It is not obvious to us that this inhomogeneous phase can be treated by a Ginzburg-Landau expansion [73], or any other improved version [74], because in these approaches one expects the appearence of an inhomogeneous phases when the mean-field analysis indicates a first-order phase transition. Instead, in the present case, the tree-level analysis does not show any phase transition or any instability: the only sign of an odd behavior is in the spectrum of the pseudo NGB mode.
The Lagrangian term in Eq. (26), which does not lock the two global symmetries, is also interesting, because it induces a nontrivial interaction between the two condensates. However, in χPT this term can only arise at the NLO in the chiral expansion, thus we expect that it is strongly suppressed. The tree-level interaction potential is reported in Eq. (65): since it is independent of n 1 and n 2 , it is clear that in this case the two condensates are free to oscillate and are not locked. The lowenergy modes consist of two NGBs which do not show any singular behavior. Upon minimizing the potential in Eq. (65) we obtain the phase diagrams reported in Fig. 3.
The present work can be extended in different ways. As already anticipated, it paves the way for the discussion of a two-component system of pions and kaons. We plan to develop this study shortly. It would also be interesting to realize the locking instability in two-component ultracold atoms system.

Radial field
Expanding the radial field around the stationary value as ρ =ρ + χ and neglecting the angular fluctuations we obtain from Eq. (6) where the O(χ 5 ) terms and higher have been suppressed. It is convenient to rescale the field with χ → χ/f π to put the kinetic term in the canonical form, obtaining where the mass and self-couplings are given by We notice that the only nonvanishing term at the phase transition point is the one proportional to χ 4 . Actually, it can be easily proven that any term proportional to χ 2n+1 vanishes at γ = 1, because in the unbroken phase the system is symmetric for ρ → −ρ. Close to the phase transition point, the radial fluctuations can be considered as a self-interacting system of bosons with vanishing mass and cubic interaction but nonvanishing quartic interaction. This Lagrangian for the radial fluctuation is valid in the whole broken phase. For the angular field the situation looks different.

Bogolyubov mode
The Lagrangian of the angular field is given by withf π = f π sinρ playing the role of an effective decay constant. Sinceφ is a unit vector, we can parameterize it by a Bogolyubov mode α as follows: which is the Lagrangian of a free scalar field, α. It can be cast in the canonical form by α → α/f π . The Bogolyubov field can only feel the medium effect by the interactions with the χ field, as will be discussed below. We note that the NLO chiral terms would be proportional to higher powers of momentum, therefore this is the relevant Lagrangian only for p 2 /f 2 π ≪ 1. For this reason, this low-energy expansion is not valid close to the phase transition point, corresponding to γ = 1, wheref π vanishes and thus all the terms of the effective Lagrangian are equally important. Since the momentum scale is dictated by the temperature of the system, one has to consider the T /f π ≪ 1 case.

Mixed terms and dispersion laws
The mixed terms can be obtained from the interaction terms in Eq. (6) and considering that upon substituting Eq. (A7) we have the compact expression in terms of the Bogolyubov field α. Thus, up to the fourth order in the fields, the mixed interaction terms are L χα I = − g 2,1 χ∂ 0 α + g 3,1 χ 2 ∂ 0 α + g 3,2 χ∂ µ α∂ µ α (A10) with the couplings given by: where the first subscript indicates the total number of fields and the second one the number of α fields. The quadratic Lagrangian can be written as where the mixing term allows oscillations between the radial and the angular fields. Integrating out the radial fluctuations one obtains the massless, phonon-like, dispersion law where describes the pressure oscillations propagating at the sound speed. Alternatively, one can diagonalize the quadratic Lagrangian, obtaining the dispersion laws where which agree with the expressions reported in [48]. In conclusion, the low-energy modes of a single-component pion gas correspond to a NGB with dispersion law in Eq. (A16) (in the limit p/m χ → 0) and to a massive mode with mass m eff .

Appendix B: LHY correction
Close to the phase transition to the broken phase, the pressure of the single-component pion gas can be approximated with the expression in Eq. (10). Therefore, χPT analysis gives a correction to the mean-field value proportional to n 3 . However, in the context of condensed matter physics, an additional contribution, due to the vacuum energy of the NGBs, is known to play an important role in certain regimes. This contribution is known as the Lee-Huang-Yang (LHY) term, first evaluated for a hard sphere Bose gas in [75]. The LHY term is proportional to n 5/2 and is the leading correction to the mean-field results, close to the phase transition point. For a single-component pion gas, one can easily obtain the LHY correction using the mapping developed in [49] between the condensed pion gas in χPT and the Gross-Pitaevskii Hamiltonian where M = µ I , and g = 4γ 2 − 1 12f 2 π γ 2 = g 0 1 + where g 0 = 1/(4f 2 π ) is the coupling constant at the phase transition point. The LHY correction to the pressure close to the phase transition point is given by ǫ GP,LHY = M 3/2 15π 2 (g n I ) 5/2 ∝ m 4 π n 5/2 , with the particular dependence on n indicating that this is a nontrivial effect beyond mean field. The LHY contribution is the first one in the series expansion na 3 , where a = gM/(4π) is the s-wave scattering length. Close to the transition point and using the values of the coupling constant and of the mass of the GP expansion, we find that na 3 ≪ 1 that means the diluiteness condition for any γ ∈ [1,2]. However, the evaluation of the LHY term by Eq. (B3) assumes that the GP expansion is reliable, implying that 1 ≤ γ ≪ 2.
For a general evaluation of he LHY correction in the χPT context, we consider the vacuum contribution of the NGBs where E ph = c s p is the dispersion law of the NGBs obtained integrating out the radial fluctuations, see Eq. (A16). The hard cutoff, Λ, takes into account that the NGBs describe the lowenergy fluctuations below the mass scale, m χ , of the radial field, see Eq. (A3). Taking for simplicity Λ = m χ , considering the expression of the speed of sound in Eq. (A17), and that, close to the phase transition, γ ≈ 1 + n/4, see Eq. (12), we find ǫ LHY ∝ m 4 π n 5/2 , in agreement with Eq. (B3). In Fig. 5 we compare the isospin number density evaluated in χPT (solid red line), with that obtained including the LHY correction (dashed blue line), as well as with the results of the lattice simulations of Refs. [13,21,22] using the same value of their pion mass, m π = 135 MeV, and of the pion decay constant, f π = 133/ √ 2. The χPT results systematically underestimate the number density. With the inclusion of the LHY term the agreement slightly improves. It follows that the χPT + LHY pressure is always larger than the χPT one. However, the difference between the two is extremely small.
Generalizing the previous discussion to the two-component pion gases with the interaction term in Eq. (21), it is clear that there are two relevant low-energy contributions. One from the NGB, and one from the pseudo NGB. Since the latter becomes tachyonic for k > 1, the LHY contribution is illdefined. Again, this is a signal that the mean-field approximation breaks down for k > 1, and thus the evaluation of the pressure of the system given by the expression in Eq. (50) is incorrect.