On the effective lagrangian at nonzero isospin chemical potential

We revisit the most general effective lagrangian within Chiral Perturbation Theory at nonzero isospin chemical potential. In addition to the contributions already considered in the literature, we discuss the effects of new terms allowed by the symmetries, derived within the external source method including spurion fields, as well as of linear-field corrections. We study the influence of those new contributions on the energy density at zero temperature and observables derived from it, such as the pion and quark condensates and the isospin density. Corrections are shown to be compatible with lattice results, which favor a nonzero value for the only undetermined low-energy constant (LEC) to leading order ${\cal O}(p^2)$, rendering in particular a shift of the critical value for Bose-Einstein condensation. To ${\cal O}(p^4)$ we study the physical constraints on the new LEC, which renormalize the energy density and whose numerical effect is estimated within natural values. The new ${\cal O}(p^4)$ corrections give rise to more significant deviations than those previously considered and remain compatible with lattice results.


I. INTRODUCTION
The study of the QCD phase diagram has experienced a notable boost over the last decade, due to the advance of both lattice field theory at nonzero temperature and chemical potentials [1][2][3][4][5][6] and experimental data of Relativistic Heavy-Ion Collisions (RHIC) within the so called Beam Energy Scan (BES) program probing the transition at chemical freeze-out [7,8].
Regarding chemical potentials, the main interest has been focused on the baryon chemical potential µ B , motivated mostly by the possible existence of a critical point separating the µ B = 0 crossover transition from a first-order one [5,6]. However, the difficulties for the lattice analyses at µ B = 0 associated to the sign problem [4,[9][10][11][12] have motivated the study of QCD when chemical potentials associated to additional, physically relevant, charges are present, which do not present such sign problem. This is the case of isospin µ I and strangeness µ S chemical potentials, which are actually relevant phenomenological quantities at RHIC at chemical freeze-out, related to electric charge and strangeness conservation [8,13]. Another example is the chemical potential µ 5 associated to chiral imbalance, which has been explored recently both in the lattice [14][15][16] and within effective theories [17,18]. The latter has been mainly motivated by the possible existence of local P -breaking regions within RHIC and related phenomena such as the Chiral Magnetic Effect [19].
The case of isospin chemical potential µ I has become particularly interesting. As it was first showed in [20] within a leading-order (LO) Chiral Perturbation Theory (ChPT) effective theory approach, QCD at increasing µ I exhibits a second-order phase transition from the normal vacuum phase to a Bose-Einstein Condensation (BEC) phase for the charged pion modes. The transition point µ c predicted by ChPT analyses is at the physical pion mass, including next-to-leading order (NLO) corrections [21,22]. For high enough µ I , the system would enter a BCS-like phase [20] where effective theory approaches based on the low-density regime, such as ChPT, are expected to break down [23,24].
Other recent analyses at µ I = 0 have included finite temperature corrections, allowing to study the effects on the QCD transition [22,[25][26][27][28]. Temperature effects are expected to smooth out the BEC transition, increasing the value of the critical µ I value for BEC and setting an upper temperature limit for which the BEC phase does no longer take place. In addition, the QCD transition temperature T c is expected to drop with increasing µ I . Other effects considered in this context have been µ B = 0 [29] and µ S = 0 [30,32] corrections. A recent analysis within the NJL model confirms the ChPT results, showing some deviations for µ I > ∼ 2M π [31].
Lattice works have confirmed most of the previous findings. In particular, in [24] the energy density as a function of the isospin density n I , as well as n I versus µ I were investigated at a small but nonzero temperature of around T ∼ 20 MeV. Significant deviations from the ChPT prediction take place for µ I > ∼ 2M π and the ratio of the energy density with the Stefan-Boltzmann energy density corresponding to free quarks exhibits a maximum around µ I 1.3 M π , which was interpreted in [24] as the onset of the BEC phase. Such peak behaviour is well reproduced both within the ChPT and NJL approaches [31]. The apparent plateau of that ratio above µ I = 3M π might indicate a BEC-BCS crossover.
On the other hand, lattice simulations at finite temperature and nonzero isospin chemical potential confirm the decreasing behaviour of T c with µ I as well as the increase of the BEC onset with temperature [33,34]. Recent analyses provide results for the pion and quark condensates from T = 113 MeV up to the QCD transition temperature, confirming those trends for the T − µ I phase diagram [35,36]. Lattice results for n I as a function of µ I for low temperatures are provided in [37]. Recent comparisons between lattice data and the ChPT and NJL model approaches can be found respectively in [38] and [39].
In this work we will discuss some new relevant aspects related to the formulation of the effective ChPT lagrangian at nonzero isospin density. The usual approach to construct the effective lagrangian follows from considering the isospin contribution in the QCD lagrangian as an external constant vector source in the τ 3 direction and using the well-know external source method [40,41] so that the chemical potential enters through covariant derivatives rendering the theory locally invariant under the chiral symmetry. However, since the µ I term in QCD is not invariant under isospin rotations in the τ 1,2 directions, the most general effective lagrangian may include additional terms with such symmetry-breaking patttern, which would be multiplied by additional low-energy constants (LEC). A systematic procedure to account for all the possible terms of such type is provided by the so-called spurion method, which is an extension of the external source method developed originally to include correctly the effects of the electromagnetic field in the chiral lagrangian [42][43][44][45]. Following similar ideas, the most general chiral lagrangian for µ 5 = 0 has been derived up to O(p 4 ) in [18].
Thus, the plan of the paper is as follows. In section II we will revisit the construction of the most general µ I = 0 effective lagrangian up to O(p 4 ) including possible new terms of the type commented above, which include both LO and NLO contributions, as well as a new correction coming from the linear terms in the effective lagrangian, which arise to NLO. The effect on different observables steming from the energy density and available in the lattice will be carefully examined to LO in section III, where we will perform a fit analysis to lattice results using the only undetermined LEC to that order. NLO corrections will be discussed in section IV, including the effect of new LEC associated to additional terms allowed by the symmetries, as well as the linear correction. In this paper we will work for SU (2) at T = 0 to obtain a first glance of those new effects, leaving the finite temperature and strangeness corrections for future work.

II. CHIRAL LAGRANGIAN INCLUDING EXPLICIT ISOSPIN-BREAKING OPERATORS
We start from the QCD lagrangian including a nonzero isospin chemical potential, corresponding to the grandcanonical partition function. The fermionic part of the lagrangian reads with q T = (u, d), D is the covariant derivative corresponding to the gluon field, M is the quark mass matrix, which here we take as M = m1 with m = m u = m d , and τ k are the Pauli matrices. The µ I = 0 lagrangian is chiral invariant SU (2) L × SU (2) R for M = 0 (chiral limit) which reduces to the isospin symmetry SU (2) V with V = L = R for nonzero quark masses. However, for nonzero µ I , the latter symmetry breaks down to U (1) I3 since the µ I term is only invariant under vector transformations in the τ 3 direction. Note also that Lorentz invariance is also broken by the inclusion of the µ I term, as a consequence of the preferred reference frame of the thermal bath at rest. In addition the isospin chemical potential term breaks C invariance, since it is essentially a charge operator. Therefore, the low-energy effective lagrangian must share the above symmetry requirements in the most general way. That is, at a given order in the chiral (low-energy) power counting, one must include all possible operators compatible with the symmetries and their breaking. This can be ensured by following the external source method, where the scalar, pseudoscalar, vector and axial vector sources are space-time dependent to ensure local chiral invariance, as well as Lorentz, P , C invariance, of the QCD lagrangian. The effective lagrangian is then constructed out of the most general set of operators satisfying such invariance at a given order in the generic momentum scale p, which accounts for meson masses and derivatives, as well as µ I = O(p) which limits the ChPT analysis to low and moderate values of µ I , as discussed below.
An important comment is in order here. The building blocks for constructing such operators are in principle the Goldstone Boson (GB) field operator codified in a SU (2) matrix field U , its covariant derivative d µ U including the external sources, as well as the flavor matrices entering those external sources, which in the case given by (1) are the mass matrix M and µ I 2 τ 3 . The latter point, i.e, the fact that one can have additional operators including µ I 2 τ 3 in a compatible way with the symmetry pattern, is actually one of the main novelties of our present work. One can easily understand this by considering all possible operators to orden O(p 2 ). In addition to the standard terms the following independent term is also allowed since it breaks chiral symmetry but preserves U (1) I3 (see details below). The philosophy to include such additional terms is the same as when introducing electromagnetic corrections to the chiral lagrangian [42] and the systematic procedure to account for all those terms consists in introducing the so-called "spurion" fields Q L,R (x) as additional space-dependent external sources transforming suitably under chiral transformations [18,[42][43][44][45].
Let us explain this procedure in more detail for our present case. For that purpose, we cast the lagrangian (1) in terms of external sources as where v µ ∈ SU (2). The above lagrangian corresponds to the choice relevant for this work, Q L = Q R = Q, and we have set the axial source a µ = 0 with respect to the general source lagrangian considered in [18,[43][44][45], where we have kept a nonzero pseudoscalar source p(x) in order to derive expectation values of pionic fields, such as the pion condensate.
The lagrangian in (1) corresponds to the particular choice v µ + A µ Q = µ I 2 τ 3 , s = M and p = 0. Thus, after the general effective lagrangian is constructed, we will choose, without loss of generality, consistently with our power counting as long as µ I Λ. The choice of the parameter Λ is irrelevant, as it should, since it will be absorbed in the new independent LEC involved, which will have to be determined. The important point is to include µ I in the Q contribution, consistently with the counting µ I = O(p). In addition, as customary we will write p(x) = jτ 1 so that thaking derivatives with respect to j we reproduce expectation values of the π 1 field, which is the direction we have chosen for the condensed field.
The lagrangian (4) can be made locally invariant under chiral SU (2) L × SU (2) R rotations of the quark and source fields. Here, it is enough to restrict to the vector subgroup g L = g R = g ∈ SU (2), under which the lagrangian in (4) can be made invariant by considering the following transformations with with M the tree level pion mass. In addition, one has to consider the C and P transformations of those fields, given in [18,[43][44][45], leading to a C, P invariant lagrangian.
The key point here is that the v µ and A µ Q sources can be transformed independently (we choose here to leave A µ invariant for simplicity, which does not affect our arguments). As we are about to see, those "spurion" transformations are essential to ensure that all possible operators are accounted for. Now, the chiral lagrangian can be constructed order by order in the chiral expansion. The building blocks and their chiral counting are the pseudo-GB field U ∈ SU (2), which is O(1), its covariant derivative which is O(p), and the external fields, which in this case are χ = O(p 2 ) and v µ = O(p), A µ Q = O(p), where we follow the same convention as in [18,[43][44][45] assigning the counting A µ = O(1), Q = O(p) as commented above. In addition, the following covariant derivative of the Q field which is O(p 2 ), has to be formally considered, although in practice it will not enter our effective lagrangian at the order considered here. The chiral transformations of the U field, i.e., U → g R U g † L , as well as those of the external sources and the covariant derivatives, and their C, P transformations, allow to construct the most general effective lagrangian which is locally chiral and C, P invariant for our present case. In particular, regarding the external sources, in addition to the usual corrections coming from the covariant derivative (7), which depends only on v µ + A µ Q, such as the first term in (2), additional terms are allowed, such as (3). That term comes from the invariant operator Tr QU QU † which will be therefore multiplied by µ 2 I times an independent LEC. Let us now proceed order by order within this scheme, following the discussion in [18] which considers all possible terms coming from arbitrary constant Q L,R fields. At the lowest O(p 0 ) order, only trivial µ I -independent terms can be constructed out of the U fields. At O(p) the only allowed structure is Tr(Q), which according to (5) would give rise to a term linear in µ I in the energy density and therefore would imply a nonzero isospin density for µ I = 0, which is ruled out by the lattice [24,37]. Therefore, demanding n I (0) = 0 we will fix to zero the LEC multiplying that term and ignore it in what follows.
To O(p 2 ), in addition to the structures (2) and (3), the terms Tr(Q 2 ) and (TrQ) 2 are also allowed, the latter vanishing with our choice (5). Therefore, the most general O(p 2 ) lagrangian at nonzero µ I is given by where F is the tree-level pion decay constant. a 1 and a 2 are new low-energy constants to be determined below and we have included a F 2 factor in front of the new terms to keep easy track of the contributions of different chiral orders to the observables, so that the chiral expansion can be parametrically tracked in powers of F 2 . Note also that the equations of motion (eom) derived from the above lagrangian are modified with respect to the µ I = 0 case as The above O(p 2 ) eom can be used in the construction of the O(p 4 ) lagrangian to eliminate some operators in favor of a minimal set [43] together with standard SU (2) operator identities.
In principle, the following O(p 3 ) operators are also allowed [18]: The five operators in (11) vanish trivially with Q in (5). In addition, one can check that the operators in (12) also vanish for d µ U and U in SU (2) and Q in (5) as long as we remain in the isospin limit m u = m d , as we will do here consistently with isospin conservation. The last term in (12) is actually proportional to Tr( As for the O(p 4 ) lagrangian, following the derivation in [18] we have in our present case with Q in (5). The usual terms considered in the literature are included in the L 0 4 lagrangian, for which we have used the basis in [46] 1 , whereas L Q 4 contains all possible operators including explicitly the Q field and whose influence on different observables will be analyzed in detail in section IV D. With respect to the basis considered in [21,22,28], the h 1 constant considered in those papers corresponds to h 1 − l 4 here.
In addition, following also the derivation in [20][21][22]28], we allow for a nontrivial vacuum configuration of the GB field, which will parametrize the pion condensed phase. Namely, we write with π a the pion fields, so that the condensed phase is characterized by a nonzero value of the α angle, which is determined by minimizing the vacuum energy density at a given chiral order, from which we will actually obtain the main properties of interest here. In the present work we are interested in the vacuum energy density at zero temperature, i.e., where Z(T, µ 5 ) is the euclidean partition function constructed out of the effective lagrangian for nonzero µ I and n = O(p n ). The above definition of the vacuum energy density is equivalent to the effective potential considered in [20][21][22]28].
The observables of interest we can calculate in both phases are the quark and pion condensates and the isospin density, which are given respectively by

III. LEADING ORDER RESULTS
The lowest order of the energy density, 2 , parametrically of O(F 2 ), is given just by minus the constant part of the L 2 lagrangian in (9): The value α LO 0 minimizing the above expression with respect to α is given for j = 0 by 1 In the basis used in [18] there is a typo in the LEC multiplying the Tr χ † U − χU † 2 operator in eq.(4.1) in that paper, which should read −(l 4 + l 7 )/16. (21) Therefore, at this order, the constant a 1 displaces the critical BEC value µ c from the squared pion mass. We will see below that this is perfectly compatible with lattice results, which will constrain the a 1 value and its uncertainty. Note that the upper bound a 1 < 1 arises here from the very existence of a BEC phase.
The above result is also consistent with the modifications of the leading order pion dispersion relations steming from the µ I -dependent lagrangian above. In fact, let us consider the linear and quadratic terms in the pion field, from the L 2 lagrangian in (9): where, following the notation in [21], we have Now, we follow the same steps as in [21] to obtain the pion dispersion relation in terms of the parameters m 12 , m 2 1,2,3 . Note that to leading order, we can replace α = α LO 0 , which in particular cancels the contribution proportional to π 1 (x) in (22) so that such linear lagrangian becomes a total derivative in the condensed phase and vanishes in the normal phase. Therefore, at this order, the linear lagrangian does not play any role for the dispersion relation nor in the vacuum energy density. As we will discuss in section IV C, this does no longer hold to NLO, so that term will have to be included.
We then get for the dispersion relation of charged and neutral pions to leading order, respectively, where p ≡ | p|. Now, setting α = α LO 0 given in (21) and p = 0 we get then the dependence of the static pion masses on the isospin chemical potential, now including the correction from the a 1 term, which for j = 0 read Thus, as in the a 1 = 0 case, the vanishing mass of one of the charged pions signals the onset of BEC condensation, as a Goldstone mode corresponding to the U (1) I3 spontaneous symmetry breaking of the vacuum in that phase, i.e, with a nonzero pion expectation value. The above mass dependence is plotted in Fig. 1 for a sample value of a 1 = −0.1, compared to the a 1 = 0 case, the qualitative dependence with µ I being quite similar, although note that a 1 introduces a non-polynomial dependence below µ c . From the energy density in (20) we get the LO to the quark and pion condensates and the isospin density. Namely, from (17), (18), (19), we get, replacing α = α 0 and for j = 0, Note that the only dependence with a 2 shows up in the pion density, which should remain zero below the BEC point. Therefore, with that physical requirement we fix up to higher order corrections, which leaves us only with one free parameter at this order.

A. LO fits to lattice
Here, we will discuss more quantitatively the effects of the new terms to leading order, by contrasting the results for different observables with those obtained in the lattice. For that purpose, we have considered the most recent lattice results. In particular, we will use the lattice points provided in [22] at T = 0, coming from the collaboration [37], for the quark and pion condensates. The latter are given for a finite pionic source. We take j = 0.00517054M π , one of the two values considered in [22]. As for the isospin density, we compare with the results quoted in [37] for j = 0 and T small enough to provide an accurate enough description of the T = 0 case. We also fix the numerical values of M π = 131 MeV, F π = 90.51 MeV and m = 3.47 MeV for an easier comparison with the results in [22,37]. In addition, we assign a characteristic uncertainty of 5% to the lattice data, following again [22,37].
As discussed above, to LO the new contributions considered here are parametrized only in the constant a 1 . The relevant question whether a 1 = 0 is preferred for lattice results over the standard choice a 1 = 0. Recall that a 1 = 0 implies in particular a displacement in the critical value µ c with respect to the pion mass. Actually, it is worth mentioning here that in some recent lattice analysis [36] the condition µ c = M π is actually imposed, not from a physical requirement but from the ChPT result, i.e., a sufficiently small variation of µ c = M π could be equally assumed. In fact, as we are about to see, fits to lattice points favor a negative nonzero value for a 1 , which according to (21) yield µ c < M π to LO, consistently for instance with the position of the lowest lattice point in the isospin density given in [37].
Thus, we show here the results of three different fits, which are summarized in Table I and Figs. 2, 3 and 4. The uncertainty bands in the figures and the uncertainty range for the a 1 parameter correspond to the 95% confidence level.
Those results lead to the following partial conclusions. First, the quark and pion condensates can be reasonably fitted with a 1 = 0 within the uncertaintites considered (fit 1). Fitting only the quark (pion) condensate favors a positive (negative) a 1 value, still compatible with zero. However, the conclusion is very different when fitting the isospin density in fit 2, for which we have included only the lattice points with µ I > M π . As explained in [37], the lattice uncertainty for µ I ∼ M π is actually much higher, although we have still plotted the first lattice point below µ I and given in [37]. As mentioned above, that point lies below M π with a nonzero isospin density central value, which is actually favored by our present analysis with negative a 1 . In fact, as can be seen from the results in Table I, the lattice results for n I are much better fitted with a negative value for a 1 than for a 1 = 0, as the values of the corresponding χ 2 /dof clearly show. The same conclusion is reached when performing a combined fit of the three observables (fit 3), with significant corrections from the a 1 = 0 case. Note also that the fits for n I reproduce quite well the two ends of the lattice points, i.e., the closest points to µ I = M π and those for larger µ I > M π , improving over the a 1 = 0 case.   To end this section, we show the results for α LO 0 , i.e., the solution of the equation ∂ 2 ∂α = 0, for the value of a 1 obtained in the joint fit 3, both for j = 0 and for the value of j used here. The explicit expression of α LO 0 is given in (21). While for j = 0, the a 1 contribution changes the transition point to the BEC phase, as already comented, for j = 0 the transition is a crossover and a 1 merely modifies the inflection point of the α LO 0 function. To summarize this section, lattice results are perfectly compatible with a coefficient a 1 = 0, which actually improves the description of the isospin density and the overall description of n I and the two condensates.

IV. NEXT TO LEADING ORDER RESULTS
The next to leading corrections 4 to the energy density come from three different sources: • Loop corrections coming from the quadratic field terms (23) in the L 2 lagrangian, which can be obtained in a similar fashion as to the standard free-field contributions to the partition function in vacuum, including now the µ I -corrections to the pion dispersion relation [21].
• The constant terms coming from the L 4 lagrangian in (13)- (14). The LEC coming from this contribution, including the new ones coming from L Q 4 , will be renormalized to absorb the loop divergences. • The linear field terms coming from L 2 and given in (22) will also contribute as long as the minimizing angle α N LO 0 = α LO 0 . This contribution has not been considered in previous works on this subject and is discussed below.

A. Loop contributions
Following the same steps as in previous works [21] we can write the one-loop contribution to the energy density as with E ± (p) and E 0 (p) in (28)-(29) and, following the notation in [21], with D = 4 − 2 and → 0 + and µ the Dimensional Regularization (DR) scale.  The treatment of the above integrals, separating their UV divergent contribution in DR follows the same steps as in [21] except for the modifications proportional to a 1 in equations (24)- (27). The contribution from the charged pions to the one-loop energy density can be written as and where we have introduced the quantities E 1,2 (p) given by Separating in this way the divergent part of the charged contribution, div 4,+− in (39) has the form of the neutral part in (36) and is easier to handle in DR. Actually, note that the large-p behavior of E 1 + E 2 is the same as the sum E + + E − . The finite contribution (subtraction integral) f in 4,+− can be calculated numerically. Thus, in DR we get for the full divergent part with N = 1 − γ + log(4π) and γ the Euler constant.

B. Fourth order lagrangian and renormalization
The contributions to the energy density coming from the constant part in the fourth order lagrangian (13) and (14) are respectively given by 40 4 Comparing (44) with (43), we see that new terms introduce µ I -dependent corrections as follows:q 1 shifts the (l 1 − l 2 )µ 4 I contribution,q 2,3,4,5 modify respectively the four l 4 µ 2 I terms, whereasq 6,7 introduce new µ 6 I terms. Of those new seven independent LECs appearing in the energy density,q 4 andq 5 will not contribute for j = 0 and thê q 7 term is independent of α so it does not contribute to the energy density minimum.
The new LECs will precisely absorb the new loop divergences dependent on a 1 included in (42). Thus, the renormalized energy density at NLO resulting of adding all the contributions mentioned before can be written as 45) where the renormalized and scale-independentl i ,h i are the standard ones given in [22,40] while the new LECs are renormalized asq with where, as usual, the µ-dependence of the renormalized LECs cancels with that of the loops, rendering the energy density finite and scale independent. In the following sections we will analyze the effect of these new LECs on the µ I -dependence of the different observables obtained from the energy density.
The numerical values we will use for thel i as well as for the physical pion mass and decay constant will be the same as the previous works on this subject [22,38] for an easier comparison. Thus, we will take the central values of which have been used to study the quark and pion condensates in [22] and the isospin density in [38] at zero temperature.
As for the physical M π , F π , as commented above, we are considering those used in the lattice simulations, which using the one-loop standard ChPT expressions [40] give rise for the following values of the tree-level M, F appearing in our previous expressions for the energy density: Recall that thel i ,h i , although scale-independent are mass-dependent and therefore there might be slight numerical variations from the values (48) when taking the masses in (49). Those variations are logarithmic and therefore numerically negligible, so for our purposes of comparing with previous works we will still use the values in (48).

C. Linear terms
The solution for the angle minimizing the energy density to leading order, given in (21) for j = 0, is such that the linear term in (22) proportional to π 1 (x) vanishes so the linear terms can be ignored since the derivative term in (22) does not contribute to the vacuum energy density. However, this is not neccessarily true to higher orders. That is, if we consider a perturbative deviation of the LO minimizing angle, where the dots denote derivative terms and, according to (22), The above linear contribution to NLO can be reabsorbed into a redefinition of the π 1 field, over which we are integrating to get the energy density. Namely, The above shift eliminates the linear term at this order and completing squares generates the following additional O(1) perturbative contribution to the NLO energy density: with m 2 1 in (25), which added to the contributions given in (36), (43) and (44) has to be minimized with respect to α, around α LO 0 , to find α N LO 0 . Moreover, it must be taken into account to calculate different observables because of its m, j and µ I dependences: To have a more quantitative idea of the effect of this correction, we have plotted in Fig. 6 the result for the minimizing angle α 0 , comparing the LO contribution with the NLO with and without including the linear term contribution (54). For easier comparison with previous works, we have not included in that plot the new contributions coming from the a 1 and q i terms and we have used the same LEC as in [38]. We consider both the j = 0 and j = 0 situations. As we can see in that figure, the inclusion of the linear term may generate sizable differences between the NLO and LO results. Actually, for some values of the constants involved, those linear corrections can be such that the effective potential stops having a minimum above a certain µ I value. We can actually see this behaviour in the plot showed in Fig. 6 for which that limiting value is µ I 300 MeV for j = 0 and µ I 340 MeV for the j = 0 value considered. The deviations with respect to the LO indicate that we are reaching the borderline of the ChPT validity limit where in particular the very same approximation followed in (50) and (51) would fail. This is actually consistent with lattice analyses showing that deviations from ChPT around those µ I values signal the onset of the BCS phase [24]. Nevertheless, it should be taken into account that the actual value of such validity limit depends also on the rest of the LEC involved, as we will discuss in detail below. . Therefore, taking the derivative with respect to µ I and setting j = α = 0 in our previous expressions for the energy density, we get, where f 0 = f (j = 0) with f in (52) and we have followed [21] for the calculation of the loop integrals. Thus, since the µ I dependence of the n 1 function above (coming from the linear term) is non-polynomical, the only way to ensure that n I vanishes for all µ I below µ N LO c is that theq r i satisfy the constraint On the other hand, the condition that the n 0 function in (58) vanishes, together with (35), implies at this order: The constraints (61) and (62) allow us to eliminate the dependence on two of the LEC involved in terms of the remaining ones. In fact, note that equation (61) is actually the explicit expression for the O(1/F 2 ) in (35). We also remark that the condition (60) is fully consistent with having obtained a numerical value a 1 1 in our fit study in section III A. The situation is similar to the EM corrections in ChPT, where there are operators which are formally O(p 2 ), such as (3), but multiplied by e 2 which is a numerically small parameter.
On the other hand, we must keep in mind the constraint (59) when considering the possible numerical range for theq r i in the following sections.

The critical BEC value at NLO
Now, let us evaluate the effect of the NLO corrections to the critical value µ c . Since ∂ N LO ∂α | α=0 = 0, µ c can be determined as the value for which α = 0 flips from a local minimum to a local maximum. Therefore, expanding the energy density around α = 0, we have From the energy density discussed above, we find where M π and F π are the NLO pion decay constant and mass respectively [40] and where we have made use of the condition (60) so that the a 1 dependence in the NLO has been ignored. From the previous expression, we see that onlyq r 2,3,6 and a 1 modify the critical BEC value. We show in Table II the value of µ c expected within the range of natural values for those constants, defined within a typical ChPT uncertainty range [44,45] as where we have set the scale as µ = M ρ 770 MeV and we have taken a 1 = −0.019, the central value of the LO fit 3 to lattice results in Section III A. In addition, we have highlighted in the table the values for which the condition (59) holds. As explained,q r i values not satisfying that condition are not acceptable, since they give rise to a nonzero isospin density below the critical value.
Note that, as happened at LO, our results support µ c < M π . In turn, we remark that the above analysis would allow in principle to fix µ c = M π , as in previous ChPT studies [21,22], by imposing an additional constraint relatinĝ q r 2,3,6 and a 1 from (65), namely wherel 3 arises from the M 2 π − M 2 difference [40] and we have replaced M by M π in the NLO when the difference is of higher order.
Nevertheless, as explained in previous sections, lattice results are actually compatible with µ c < M π and in principle there is not a physical reason to impose the above constraint, unlike those coming from the vanishing of the isospin density and discussed in section IV D 1, so we will allow theq r i to fluctuate within natural values without imposing (67) when evaluating the different observables in the following section.

NLO Results for different observables
We will consider now the NLO evolution of chiral observables for nonzero µ I , in particular regarding the role of theq r i LEC and the comparison with lattice analyses. From the results in the previous sections, we see that the NLO energy density depends on seven independent new LECs, namely a 1 ,q r 1−6 whose numerical values will then influence the µ I dependence of the observables. Note that all the observables considered depend onq r 1−6 since, in addition to the explicitq r i dependence one must consider that in α N LO 0 . We will represent our results for the range of natural values (66) at the scale µ = M ρ and setting a 1 = −0.019 (the mean value obtained in Fit 3 in section III A). Thus, we have 3 6 points for each µ I , corresponding to the three values 0, ±1/(16π 2 ). In doing so, we discard thoseq r i violating the condition (59) and we calculate the mean square error of the results which provides a dispersion estimation. Note that, as commented in section IV C, due to the linear term, the minimum could disappear above a given µ I . Thus, the uncertainty bands in the following figures for a given µ I > µ c correspond to thoseq r i combinations for which the minimum exists. The effect of the linear term will be actually showed separately in the figures in order to calibrate better its effect.
First, in Fig.. 7, we show α N LO 0 , i.e., the angle minimizing the energy density including NLO corrections. We see that including all the corrections discussed here implies sizable deviations above µ c with respect to the LO, larger than in previous analyses [21,22]. Note in particular the tail below µ I = M π coming from the reduction in the numerical value of µ N LO c as Table II shows. We plot in Fig. 8 the quark and pion condensate deviations (as defined in [22]) at NLO for natural values ofq r i with and without linear term, comparing with lattice results. The NLO corrections are again significant and remain within the uncertainties of lattice points, taking into account that we are not performing a NLO fit so there would be still room for improvement. The NLO corrections actually improve over the LO fit for the pion condensate for high µ I values, whereas in the case of the quark condensate, with the inclusion of the linear term, the theoretical curve seems to depart from the lattice points with respect to the LO. Finally, we study the isospin density at NLO for j = 0. The result, including the uncertainty bands ofq r i natural values, are showed in Fig. 9. As for previous observables, the NLO with the new terms considered in this work provides significant deviations from the LO as µ I increases, still accommodating the lattice results.
As explained above, for the NLO analysis we have fixed a 1 to the LO fit and consider the newq r i within natural values. In view of the results in Figs. 8 and 9, we would surely have obtained a much better description of lattice results by keeping a 1 and theq r i as fit parameters with the NLO curves, but such precision analysis is outside the scope of this work.

V. CONCLUSIONS
In the present work, we have analyzed the most general effective chiral lagrangian for nonzero isospin chemical potential for two light flavours up to fourth order. We have followed the technique of external sources including spurion fields, which allows to account for all possible operators respecting the symmetry breaking problem at hand. We have analyzed in detail the effect of new terms with respect to previous analyses.
In particular, at leading order, two new-independent terms in the L 2 lagrangian have to be considered, whose corresponding low-energy constants a 1 , a 2 are related from the constraint of vanishing isospin density below the critical BEC value µ c . The constant a 1 generates a shift in µ c with respect to the pion mass at that order. To estimate the preferred value of a 1 , we have performed several fits to lattice results for the quark and pion condensates (for nonzero pionic source j) as well as for the isospin density (for j = 0). When comparing to the a 1 = 0 results, a small nonzero negative value for a 1 is favored, improving the description of lattice results, specially for the isospin density.
To NLO, the fourth order lagrangian L 4 contains seven new terms multiplied by new low-energy constantq i , which we have consistently renormalized to absorb the divergences coming from loops with vertices of the new L 2 terms. Apart from those terms, we have shown that to NLO one has to consider an additional contribution coming from a term in L 2 linear in the pion fields. That term does not contribute to LO but it does to NLO due to the corrections in the angle minimizing the energy density. The effect of the linear term is qualitatively important, since it eventually makes the minimum of the energy density disappear, which sets a natural limit of validity for the ChPT framework, consistently with lattice analyses.
Imposing that the isospin density vanishes below µ c to NLO gives rise to additional constraints, involving now the a 1 , a 2 ,q r i LEC. Those constraints allow on the one hand to understand the small value for a 1 obtained in the lattice fits and on the other hand, to eliminate one of the seven new LEC. In addition, the critical BEC value to NLO must remain below the LO one, which restricts further the admisible values for theq r i . A further constraint could be obtained by demanding µ c = M π . We have estimated the effect of the new LEC and the linear term to NLO by keeping theq r i within natural values. The results for the different observables show again consistency with the lattice points, leaving room for improvement with respect to the LO.
In summary, our present analysis establishes systematically the most general way to describe low-energy QCD at nonzero isospin density, consistently with lattice results and complementing previous theoretical analyses. We believe that this work will be useful towards a better understanding of the QCD phase diagram and we leave for future works its extension to three flavours and finite temperature.