Inhomogeneous chiral condensates in three-flavor quark matter

We investigate the effect of strange-quark degrees of freedom on the formation of inhomogeneous chiral condensates in a three-flavor Nambu–Jona-Lasinio model in mean-field approximation. A GinzburgLandau study complemented by a stability analysis allows us to determine in a general way the location of the critical and Lifshitz points, together with the phase boundary where the (partially) chirally restored phase becomes unstable against developing inhomogeneities, without resorting to specific assumptions on the shape of the chiral condensate. We discuss the resulting phase structure and study the influence of the bare strange-quark mass ms and the axial anomaly on the size and location of the inhomogeneous phase compared to the first-order transition associated with homogeneous matter. We find that, as a consequence of the axial anomaly, the critical and the Lifshitz point split. For realistic strange-quark masses, the effect is, however, very small and becomes sizeable only for small values of ms.


I. INTRODUCTION
Mapping the phase diagram of QCD at nonvanishing temperature T and quark chemical potential μ is one of the major goals in strong interaction physics [1,2]. Lattice-QCD calculations revealed that chiral symmetry, which is spontaneously broken in vacuum, gets approximately restored at high temperature via a smooth crossover transition [3]. At low temperature but nonvanishing chemical potential, standard lattice Monte Carlo techniques are not applicable. In this regime, one therefore has to rely on other approaches to QCD, like Dyson-Schwinger equations [4] or the Functional Renormalization Group (FRG) [5], or on effective models, like the Nambu-Jona-Lasinio (NJL) or the quark-meson model. Assuming spatially uniform chiral condensates, these approaches typically predict a first-order phase transition from the chirally broken to the (approximately) restored phase, ending at a critical end point [4][5][6][7]. In the two-flavor chiral limit, one finds instead a tricritical point where the first-order phase boundary at low T is joined to a second-order one which replaces the crossover in the high-T, low-μ region. For simplicity, we will refer to both the tricritical point and the critical end point as the "critical point" (CP) in the following. Some time ago, it was found that spatially inhomogeneous condensates could be favored over homogeneous ones in certain regions of the phase diagram [8][9][10][11]; see Ref. [12] for a review. In the chiral limit, one then finds a so-called Lifshitz point (LP) where three phases, the homogeneous chirally broken phase, the restored phase, and the inhomogeneous phase, meet. In particular, it was found for the two-flavor NJL model in the mean-field approximation that the inhomogeneous phase completely covers the first-order phase boundary between the homogeneous phases [11] and that the LP exactly coincides with the CP [10]. 1 Recently, we generalized this result to the case away from the chiral limit [14]. A similar picture also emerged from a truncated Dyson-Schwinger approach to QCD, although the coincidence of CP and LP has not yet been shown rigorously [15]. Indications for the presence of an inhomogeneous phase have also been found in a very recent FRG study of the QCD phase diagram [5].
In the present article, we want to investigate the effect of strangeness degrees of freedom (d.o.f.) on the locations of the CP and the LP in the NJL model. To this end, we will employ the three-flavor NJL model and perform a Ginzburg-Landau (GL) analysis. In order to get a better feeling on the form of the phase diagram, this GL study will be supplemented by a stability analysis, which, under certain conditions, can provide us the location of the phase boundary between the inhomogeneous phase and the partially chirally restored phase.
Our motivation for this study is twofold. First, in the twoflavor model, the CP and the LP are typically found in a temperature and density regime where strange-quark effects may not be negligible. Including strangeness d.o.f. will thus be a step toward a more realistic description. A first study of inhomogeneous phases in three-flavor matter is performed in Ref. [16], albeit only for vanishing temperature and employing a simple explicit Ansatz for the spatial dependence of the chiral condensate. The authors found that the coupling of strange and nonstrange quarks, which is related to the axial anomaly, enlarges the inhomogeneous domain. Our investigation will be complementary, as its range of validity will be close to the locations of the CP and the LP at nonvanishing temperature, as well as along the phase boundary between the inhomogeneous and the partially restored phase. Moreover, it does not rely on a certain Ansatz for the inhomogeneity and allows for a more transparent analysis of the role of the axial anomaly. Our second motivation is the fact that QCD with three very light quark flavors is expected to have a first-order chiral phase transition even at μ ¼ 0 [17]. Hence, if in three-flavor QCD the LP coincided with the CP, as it does in the two-flavor NJL model, it would mean that the inhomogeneous phase would also reach down to μ ¼ 0 for small enough quark masses. At least in principle, it could then be studied on the lattice in this case. For a first step, we therefore want to investigate the behavior of CP and LP in the three-flavor NJL model. While it has been known for a long time that the CP indeed moves to lower chemical potentials when lowering the quark masses and eventually reaches the T axis [18], until recently, the behavior of the LP had not yet been studied for the three-flavor case. In particular, it has not been checked whether CP and LP still coincide in that case. 2 Our paper is organized as follows. In Sec. II, we introduce the three-flavor NJL model and derive a general expression for the mean-field thermodynamic potential as a functional of several strange and nonstrange chiral condensates with arbitrary spatial dependence. Based on this, we perform in Sec. III a GL and a stability analysis, focusing on the locations of the CP and the LP as well as the position of the phase transition to the partially restored phase. In Sec. IV, we give a quantitative illustration of the results by numerical evaluation of the GL coefficients before we draw our conclusions in Sec. V.

II. NJL MODEL SETUP
Our starting point is the three-flavor NJL Lagrangian where ψ ¼ ðu; d; sÞ T denotes a quark field with three flavor d.o.f. andm is the matrix containing their current masses. The last two terms describe a Uð3Þ L × Uð3Þ R invariant four-point interaction, and a six-point ["Kobayashi-Maskawa-'t Hooft" (KMT)] interaction, which is SUð3Þ L × SUð3Þ R symmetric but breaks the Uð1Þ A symmetry, thus mimicking the axial anomaly. In the former τ a , a ¼ 1; …; 8, denote Gell-Mann matrices in flavor space, while τ 0 ¼ ffiffiffiffiffiffiffi ffi 2=3 p 1 is proportional to the unit matrix.
Following standard procedures, we perform a mean-field approximation, considering scalar and pseudoscalar flavordiagonal condensates, σ f ðxÞ≡hffiðxÞ; π f ðxÞ≡hfiγ 5 fiðxÞ; f¼fu;d;sg; ð4Þ allowing them to be spatially dependent. The mean-field quark self-energies can then be expressed in terms of the mass operatorŝ where f ≠ g ≠ h ≠ f. We furthermore assume isospin symmetry in the light sector, m u ¼ m d ≡ m l , and thus due to the flavor symmetry of the interaction. For the pseudoscalar condensates, given the isospin structure in the light sector, we identify while we neglect the pseudoscalar condensate in the strange sector, π s ¼ 0. The latter is motivated by the fact that pseudoscalar condensates are usually suppressed by the presence of nonvanishing bare quark masses, which we will always consider in the strange sector.
The mass operators then reduce tô which is a considerable simplification compared to Eq. (5) but still catches the most important structures. For instance, it is compatible with the Ansatz of Ref. [16], which we will briefly discuss at the end of this section. Neglecting fluctuations, the mean-field thermodynamic potential at temperature T and quark chemical potential μ is then given by where V is a quantization volume and Tr denotes a functional trace to be taken in the Euclidean space-time (ðit; xÞ ∈ V 4 ¼ ½0; 1 T × V) and over internal (Dirac, color, and flavor) d.o.f. The inverse dressed quark propagator is given by From the above expressions, we can see that the different flavors only couple through the six-point interaction, and thus for K ¼ 0, we obtain the usual expression for the twoflavor model plus a decoupled s-quark.
Having set up the model, we can now analyze its phase structure by minimizing ΩðT; μÞ with respect to the condensates. When restricting the analysis to homogeneous condensates, this is straightforward and has been done many times in the literature; see, e.g., Ref. [20]. With realistic parameters, one basically obtains the result already mentioned in the Introduction, i.e., a first-order phase transition at low temperatures turning into a crossover at higher T, both of which can mainly be associated with the chiral restoration in the light-flavor sector. Further outside, i.e., at larger values of μ or T, there is another transition (usually a crossover) related to the partial symmetry restoration in the strangequark sector. Although of minor importance for our GL analysis, this second transition will play a role for the behavior of the inhomogeneous phase at higher chemical potentials, as we will briefly discuss at the end of Sec. IV.
The numerical analysis of the model phase diagram allowing for inhomogeneous condensates is obviously much more challenging, as one should in principle allow for arbitrary spatial modulations of σ l , π l , and σ s .
The authors of Ref. [16] have therefore restricted themselves to a relatively simple Ansatz, namely a chiral density wave in the light sector, σ l ðxÞ ¼ σ 0 cosðq · xÞ, π l ðxÞ ¼ σ 0 sinðq · xÞ, m l ¼ 0, embedded into a homogeneous strange-quark background σ s ¼ const. Inserting this into Eq. (8), one finds thatM s ¼ m s − 4Gσ s þ 2Kσ 2 0 is homogeneous as well, whileM u andM d form an isospin doubletM l ¼ M 0 expðiγ 5 τ 3 q · xÞ with amplitude M 0 ¼ −ð4G − 2Kσ s Þσ 0 . This feature allowed the authors to find the eigenvalue spectrum of the corresponding Hamiltonian by adopting the known two-flavor result [21]. It should be noted, however, that this straightforward generalization of a known two-flavor solution was only possible because the space dependence of the σ l and π l contributions toM s cancel each other for this Ansatz. In contrast, the purely scalar "real-kink crystal," which in the two-flavor model is energetically favored over the chiral density wave [11], cannot self-consistently be generalized in this simple way since coupling it to a homogeneous σ s would still lead to an oscillatingM s as long as K is nonzero.
In the remainder of this work, we will not resort to a particular Ansatz for the condensates but consider general functions σ l ðxÞ, π l ðxÞ, and σ s ðxÞ both within a GL as well as a stability analysis of the homogeneous phase. This will allow us, under certain conditions, to locate the position of the Lifshitz point together with the critical point and the phase boundary where inhomogeneous condensates become energetically favored.

III. GINZBURG-LANDAU AND STABILITY ANALYSIS
The GL expansion of the thermodynamic potential allows for a systematic study of the phase structure while requiring only a limited number of assumptions. Generally, it corresponds to an expansion of the free energy in terms of one or more order parameters and their gradients, assuming that these quantities are small. In the context of chiral symmetry breaking, the order parameter is usually identified with the chiral condensate, and the expansion is performed about the chirally restored solution where the condensate vanishes. Of course, this only works in the chiral limit, while in the presence of nonzero bare quark masses, there is no chirally restored solution, and therefore the expansion has to be performed around nonvanishing condensate values [14].
On the other hand, in a realistic calculation, the bare strange-quark mass is much larger than the bare masses of the up and down quarks. We may therefore study a partially simplified problem where we take into account m s ≠ 0 but consider m l ¼ 0. In this case, a two-flavor restored solution σ l ¼ π l ¼ 0 exists, which we take as the expansion point of the GL analysis. In the strange-quark sector, we proceed similarly as in Ref. [14] for the case of light quarks away from the chiral limit and expand around a T-and μ-dependent but spatially constant condensate σ ð0Þ s , corresponding to a stationary point of Ω at σ l ¼ π l ¼ 0. To this end, we write and introduce the order parameters which are proportional to the fluctuations around our expansion point and have the dimension of a mass. We then write and expand ω GL in powers of Δ l and Δ s and gradients thereof.
Note that in Eq. (12) we combined the light-quark scalar and pseudoscalar condensates to a complex order parameter. In fact, as a consequence of the remaining two-flavor chiral symmetry, σ l and π l can be rotated into each other by a symmetry transformation, and therefore ω GL can only depend on the modulus of Δ l but not on its phase. Moreover, since we expand about the two-flavor restored phase, only even powers of jΔ l j are allowed. For Δ s , on the other hand, we can also have odd terms. Finally, the possible gradient terms are restricted by the requirement that the potential must be a scalar under spatial rotations.
Treating initially both light and strange order parameters as well as gradients to be of the same order, the GL functional up to order 4 thus takes the form where the α and β terms correspond to the contributions from the light and strange condensates, respectively, and the γ terms describe the mixing.
Since we are expanding about a homogeneous stationary point, the linear coefficient β 1 associated with Δ s has to vanish, This condition is of course equivalent to a gap equation for σ ð0Þ s at given T and μ.
For vanishing γ i , the GL analysis of the nonstrange sector would be identical to the two-flavor case in the chiral limit, which has been analyzed in Ref. [10]. More specifically, the location of the CP, where the first-order phase transition turns into a second-order one, is given by the point where both the quadratic and quartic coefficients α 2 and α 4;a vanish, while the LP appears where both α 2 and the lowest-order gradient coefficient α 4;b become zero. For the standard two-flavor NJL model, one finds that α 4;a ¼ α 4;b , and as a consequence, the CP and the LP coincide [10].
In order to study how these relations get modified in the presence of a third flavor coupled to the light sector, as reflected by γ i ≠ 0, the first step is to eliminate Δ s by extremizing the thermodynamic potential with respect to this function. To this end, we employ the Euler-Lagrange equations, which, after using the gap equation Here, the orders in Δ l and gradients are treated equally, Oð∇ n Þ ¼ OðΔ n l Þ, as we already did when we wrote down Eq. (14). From Eq. (17), we then see, however, that Δ s is of the order OðΔ 2 l Þ. 3 Inserting Eq. (17) back into Eq. (14) and keeping only terms up to the order OðΔ 4 l Þ, one then obtains We thus find that the quartic term in Δ l gets an additional contribution through the coupling to the strange quarks, while the gradient term does not. So, the relevant equations for locating the CP and the LP become and therefore, even if α 4;a and α 4;b were still equal (as they are in the two-flavor model), the CP and the LP would no longer coincide for γ 3 ≠ 0.

A. GL coefficients
In order to make more quantitative statements on the relative position of CP and LP, we need to determine the GL coefficients explicitly. To this end, we basically follow the procedure of Ref. [10], although the present case is technically somewhat more involved due to the nonvanishing strange-quark mass and the presence of the KMT term.
In a first step, we decompose the mass operator for flavor f into a constant and a fluctuating piece, by inserting Eqs. (11) and (12) into Eq. (8). Explicitly, one finds The mean-field thermodynamic potential (9) can then be written as where S −1 f;0 ðxÞ ¼ iγ μ ∂ μ þ μγ 0 − M f;0 is the inverse bare propagator of a free fermion with mass M f;0 . We have already turned out the trivial color trace (leading to a factor of N c ) and write the flavor trace explicitly as a sum, so that the remaining functional trace is only to be taken in Euclidean space-time and Dirac space.
The condensate part, corresponding to the second term on the right-hand side of Eq. (9), becomes with from which we can straightforwardly read off the condensate contributions to the GL coefficients β 1 , α 2 , β 2 , and γ 3 . In order to obtain the contributions from the Tr log term in Eq. (23) as well, we expand the logarithm into a Taylor series about S −1 0 : The functional traces are given by where the integrals are again over V 4 and tr D indicates a trace in Dirac space. Noting that S f;0 is the standard propagator of a free fermion with mass M f;0 at temperature T and chemical potential μ, the evaluation of the Dirac trace is straightforward using the momentum-space representation of S f;0 in the Matsubara formalism. After performing a gradient expansion of one can also perform the integrations over all space-time variables x i ≠ x 1 and then compare the results with Eq. (14) to read off the GL coefficients. Here, some attention has to be payed to the fact that, as a consequence of the KMT coupling, δM f contains linear and quadratic terms in the fluctuations Δ i and therefore ðS f;0 δM f Þ n contains different orders as well. We find where we have introduced the functions with fermionic Matsubara frequencies ω j ¼ ð2j þ 1ÞπT.
From the stationary condition β 1 ¼ 0, we thus obtain which corresponds to the gap equation for a constant dressed strange-quark mass at temperature T and chemical potential μ in the absence of light condensates, as expected.
With the help of this equation, some of the GL coefficients can be simplified. In particular, we find ; ð36Þ As we have seen earlier, the different flavors are only coupled through the six-point interaction. Indeed, for K ¼ 0 and thus κ ¼ δ ¼ 0, the flavor mixing coefficient γ 3 vanishes, and the α coefficients reduce to the corresponding two-flavor expressions in the chiral limit. In particular, we reproduce again the result of Ref. [10] that CP and LP coincide in this case.
For K ≠ 0, on the other hand, γ 3 ≠ 0, and therefore the CP and the LP would not even coincide if α 4;a and α 4;b were equal, as seen in Eq. (19). Moreover, α 4;a and α 4;b are not equal, and the two effects do not cancel each other. We thus find that CP and LP split for K ≠ 0, i.e., as a consequence of the axial anomaly.

B. Stability analysis
The GL analysis is valid in regimes where both the amplitudes and the gradients of the order parameters are small, i.e., in the vicinity of second-order phase boundaries between homogeneous phases or in the inhomogeneous phase close to the LP. It is possible, however, to determine the entire second-order boundary between a homogeneous and an inhomogeneous phase as well if we drop the assumption of small gradients while keeping the assumption of small amplitudes [9,14]. To this end, we start from Eq. (27), but instead of performing a gradient expansion of the mass functions, we allow the fluctuating condensates to be arbitrary periodic functions and decompose them into their Fourier modes; i.e., we write with q k belonging to the reciprocal lattice corresponding to the periodic structure. Requiring real-valued condensates in coordinate space, the Fourier components for opposite momenta are related as σ l;−q k ¼ σ Ã l;q k , etc. Using again the momentum-space representations of the propagators S f;0 as well, the integrals in Eq. (27) are readily performed. We can then decompose the thermodynamic potential as where Ω ðnÞ is of nth order in the fluctuations. The linear term Ω ð1Þ is proportional to β 1 and thus vanishes by the gap equation (35). 4 For the quadratic term, one finds with and the function and we have again used the gap equation (35) to eliminate F 1 ðM s;0 Þ from the above expressions. Note that the GL coefficients quadratic in the amplitudes could alternatively be derived as , leading to the same results as given in Sec. III A.
According to Eq. (41), the homogeneous ground state we expand about becomes unstable against developing inhomogeneities if Γ −1 l or Γ −1 s are negative in a region with nonvanishing q. Thus, as long as we are in the region where the phase diagram restricted to homogeneous condensates has a ground state with σ l ¼ 0, the second-order phase boundaries to the inhomogeneous phase correspond to the lines in the T − μ plane where one of these functions just touches the zero axis, i.e., has a vanishing minimum at some q 2 ≠ 0.

IV. NUMERICAL RESULTS
We now want to evaluate the GL coefficients and the functions Γ −1 l or Γ −1 s numerically and analyze the resulting phase structure of the model, focusing in particular on the locations of the critical and Lifshitz points. To this end, we first need to fix the model parameters. Since the NJL model is nonrenormalizable, this includes deciding upon a regularization scheme for the divergent momentum integrals. As in our previous works, e.g., Refs. [13,14,22], we choose a Pauli-Villars (PV) inspired regularization scheme, since a sharp three-momentum cutoff is known to lead to large regularization artifacts when dealing with inhomogeneous phases. More specifically, we follow Ref. [11] and regularize the vacuum contributions to the thermodynamic potential with three PV regulators, while we keep the finite medium contributions unregularized. For the integrals F n ðMÞ, we then have with E ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , the PV coefficients c j ¼ f1; −3; 3; −1g, and Fermi distribution functions n ¼ ½expððE − μÞ=TÞ þ 1 −1 andn ¼ ½expððE þ μÞ=TÞ þ 1 −1 . Similarly, we have The corresponding regulator mass Λ, the coupling constants G and K, as well as the bare quark masses m s can then be fitted to vacuum phenomenology. For a starting point, we take a parameter set [23], which has been fitted to the vacuum values of the pion decay constant and of the masses of the pion, kaon, and η 0 , while fixing the constituent light-quark mass in vacuum to a value of M ðvacÞ l ¼ 325 MeV. This yields Λ ¼ 781.2 MeV, GΛ 2 ¼ 4.90, KΛ 5 ¼ 129.8, m s ¼ 236.9 MeV, and m l ¼ 10.3 MeV. As discussed before, we work, however, in the light-flavor chiral limit; i.e., we set m l ¼ 0, instead of taking the above value. The vacuum constituent mass of the light quarks then takes a somewhat smaller value of about 310 MeV; i.e., both M vac l and Λ are roughly in the same ballpark as in similar studies in the twoflavor model. 5 Since our goals are to investigate the phase structure for both realistic and small bare strange-quark masses, we will vary m s in the following but refer to the fit value given above as the "realistic" one.
We note that the fit value of the KMT coupling in the dimensionless combination KΛ 5 is about one order of magnitude larger than typical values using a three-momentum cutoff [24][25][26]. This is consistent with Ref. [27], where a systematic comparison of various regularizations schemes has been performed and where similarly large values for KΛ 5 have been found for PV regularization. 6 Nevertheless, we will also vary the value of K in order to study its effect.
In any case, the parameters from the fit (with m l set to zero) turn out to be a reasonable starting point for our studies. In particular, when we restrict the analysis to homogeneous phases, the phase boundaries, which are displayed in Fig. 1, are in fair agreement with typical results obtained with a three-momentum cutoff. In the figure, we also show three curves associated with the roots of the GL coefficients, which, as specified in Eq. (19), determine the locations of the CP and the LP. With this parameter set, the two points turn out to be extremely close to each other, although, unlike in the two-flavor case, they do not coincide exactly (see inset), with the LP being at slightly higher temperature.
Last but not least, we show in the figure the curve where Γ −1 l ðq 2 min Þ ¼ 0, determining the onset of the instability and thus the location of the second-order phase boundary where the inhomogeneous phase meets the two-flavor chirally 5 Typical values are M vac l ¼ 300 MeV and Λ ¼ 757.0 MeV [11 , 13]. 6 A quantitative comparison with our parameters is not possible since the authors of Ref. [27] have chosen a slightly different PV scheme with two regulators and a common regulator mass. restored one. For T ¼ 0, we find that at the transition q min ∼ μ ∼ 350 MeV, and as we follow this line toward higher temperatures, its value gradually drops, until it reaches zero at the LP.
It is important to note that the instability line is located at higher chemical potential than the first-order boundary of the homogeneous analysis. Therefore, coming from the two-flavor chirally restored phase and lowering the chemical potential, the instability toward the inhomogeneous phase occurs first, and the analysis is not invalidated by a preceding homogeneous first-order transition. In turn, this means that the first-order boundary between the homogeneous phases, including the CP, cannot be trusted. In fact, we have checked by evaluating Γ −1 l that the two-flavor chirally restored solution is unstable along this first-order line, indicating that the latter is entirely inside the inhomogeneous phase.
We thus find that the phase diagram is very similar to the two-flavor one, with an inhomogeneous phase covering the entire homogeneous first-order phase boundary. Working out explicitly the left phase boundary between the homogeneous broken and the inhomogeneous phase would, however, require making an Ansatz for the shape of the modulation, which is outside the scope of this work. 7 Finally, we note that our numerical analysis finds no region of the phase diagram where Γ −1 s becomes negative; i.e., the instability is driven by the light quarks only. As discussed above, because of the flavor mixing, the dynamical strangequark mass as well as the strange condensate [see Eq. (17)] can nevertheless be inhomogeneous as well, depending on the explicit shape of the nonstrange condensates.
In order to get further insight into the nature of the splitting between the CP and the LP, we now investigate the effect of varying the model parameters associated with the third flavor, namely the KMT coupling K and the strange current mass m s on the location of the two points. A first useful observation is that the equation α 4;b ¼ 0 [see Eq. (28)] does not depend on any of these two, and therefore the location of the LP will only be affected by changes in the α 2 ¼ 0 line.
We start by showing in Fig. 2 the influence of the strange current mass on the location of the CP and the LP, varying m s from 300 MeV toward zero. As can be seen from the figure, the points move rather slowly for large and moderate values of m s but speed up when the current mass gets smaller. As expected, with decreasing m s , the location of the CP moves to lower chemical potentials, until eventually it hits the temperature axis and the phase transition between the homogeneous phases becomes first order at all chemical potentials. With our parameter set, this occurs for a strange current mass of around 10 MeV. This agrees well with the corresponding value of m s found in Ref. [18] within a three-flavor Polyakov-Nambu-Jona-Lasinio model with three-momentum cutoff, giving additional support to our parameter choice.
The LP on the other hand follows a different trajectory, significantly splitting from the CP and moving toward lower temperatures as the current strange mass decreases. In order to better understand this behavior, we show in Fig. 3 for two different values of m s the curves where the relevant GL coefficients vanish. We see that for lower m s the α 4;CP line moves toward the temperature axis (in particular around the α 2 line), thus explaining the  corresponding movement of the CP. On the other hand, recalling that α 4;LP does not change at all, the relatively mild change of α 2 suggests that the LP moves only little. Nevertheless, together with the shape of the α 4;LP ¼ 0 line, it explains why the LP goes down in temperature when m s is reduced. Again, we should, however, pay attention to the ordering of the phase transitions. Unlike for the realistic mass case in Fig. 1, we find that for small values of m s the homogeneous first-order phase boundary is to the "right" of the inhomogeneous instability line. In this case, coming from the twoflavor chirally restored phase and reducing the chemical potential, the homogeneous first-order transition takes place first and invalidates our analysis for the LP and the second-order phase transition to the inhomogeneous phase. More precisely, we find that the LPs indicated in Fig. 2 should not be trusted for m s ≤ 50 MeV, while at least part of the instability line still lies to the right of the homogeneous first-order line until m s ≃ 40 MeV. Below this value, it moves completely behind it, and it is then likely that the inhomogeneous phase disappears from the phase diagram.
In the above analysis, in order to single out the effects of varying a single parameter on the phase structure of the model, we decided not to refit completely all the vacuum observables but simply varied m s while keeping all the other parameters fixed. Changing the current strange-quark mass turns out to have a very small effect on the value of the vacuum constituent-quark mass for light quarks, 8 the value of which is known to strongly affect the phase transition in the light sector. This was, however, not the case for the sixquark coupling K, as simply decreasing it while keeping all the other parameters fixed quickly leads to unreasonably small M ðvacÞ l and consequently to the disappearance of both CP and LP from the phase diagram. In order to work around this limitation, we decided to refit the four-fermion coupling G together with K in order to give the same M ðvacÞ l as in our starting parameter set. The resulting locations of the CP and the LP for different values of K (and refitted G) are shown in Fig. 4. There, we can see that, after performing the refit, the locations of both points depend only very little on K (note the scale). Within this scale, they follow again different trajectories. As expected, the two points agree for K ¼ 0, reproducing the known two-flavor limit as the third flavor in this case decouples completely. As K increases, the two split, and the LP moves above the CP. Up to KΛ 5 ¼ 160, the situation is qualitatively similar to the one discussed in Fig. 1; i.e., the inhomogeneous phase completely covers the first-order phase boundary, invalidating the GL results for the CP. The trend reverses if we consider extremely high values of the six-quark coupling: For KΛ 5 ≳ 170, the CP reappears in the phase diagram, while the GL result for the LP can no longer be trusted. We stress again that all these variations occur in a very limited range of temperatures and chemical potentials, as can be seen by the scale in the figure.
Finally, we would like to comment on the situation at higher densities. From the two-flavor NJL model, it is known that at higher chemical potentials a second inhomogeneous phase, of which the physical nature is unclear and which we termed "inhomogeneous continent," appears [29]. The corresponding phase boundary can also be determined by the stability analysis, and we find that for the three-flavor model the inhomogeneous continent exists as well. In fact, the existence of this phase was already reported in Ref. [16] for T ¼ 0. Moreover, these authors found that, as a consequence of the KMT interaction, a third inhomogeneous phase appears between the usual inhomogeneous "island" and the continent. At the lower-μ end, this new phase is reached from the two-flavor restored phase via a second-order phase transition, whereas at the upper end, there is a first-order transition back to the restored phase. In our calculations, i.e., with our regularization and parametrization, we do not see this intermediate inhomogeneous phase, but we find a behavior which could be seen as precursor of this phase and might allow us to shed some light on its nature: Within our starting parameter set, we find that, as we increase the chemical potential beyond the second-order transition from the inhomogeneous island to the restored phase, the function Γ −1 l starts developing a new minimum for large values of q 2 , which gradually approaches zero. However, before the instability is reached, the partial chiral-symmetry restoration in the strange sector takes place, leading to a steep decrease of the dynamical strange-quark mass and pushing the minimum up away from zero. The minimum then decreases again and eventually reaches zero at the onset of the inhomogeneous continent.
At least qualitatively, this behavior and that observed in Ref. [16] can be understood from Eq. (42). Consider a point on the second-order phase boundary between the inhomogeneous phase and the two-flavor restored phase. On this point, we have Γ −1 l ¼ 0; i.e., the expression in square brackets in Eq. (42) vanishes. From this, we can see that if we now increase the value of δ, Γ −1 l becomes negative; i.e., the point moves inside the inhomogeneous phase. Hence, increasing δ increases the size of the inhomogeneous regime, moving the upper phase boundary of the inhomogeneous island upward and the lower boundary of the continent downward in μ. The situation is complicated, however, by the fact that δ depends on M s;0 and thus indirectly on T and μ. As long as M s;0 is roughly constant, we have δ ∝ K=G 2 ; i.e., increasing K (and simultaneously decreasing G by the parameter refit) enhances the value of δ and thus of the inhomogeneous region. In contrast, if we keep K fixed but increase the value of μ, the strong decrease of M s;0 related to the partial chiral-symmetry restoration in the strange sector lowers δ considerably, in our case preventing an early onset of the continent. If on the other hand the decrease of M s;0 occurs after the second inhomogeneous phase has already started, the related decrease of δ can shut this phase off again, and the continent reappears only at even higher values of μ. It is possible that this is what happened in Ref. [16], and it would be interesting to check whether this scenario is indeed what is occurring within the model regularization and parametrization used in that work.

V. CONCLUSIONS
In this work, we investigated within the NJL model how the formation of inhomogeneous chiral condensates is affected by the inclusion of strange quarks, which are coupled to the light flavors via a KMT determinant interaction, mimicking the effects of the axial anomaly. The use of a Ginzburg-Landau analysis allowed us to infer the location of the critical and the Lifshitz points in a general way without having to specify a shape for the spatial dependence of the chiral condensate. Furthermore, using a complementary stability analysis, we were able to determine the location of the second-order phase transition where inhomogeneous condensates become favored over chirally restored matter. The latter also indicates that it is favorable only for the light-quark condensates to become inhomogeneous, whereas no instability was found in the strange-quark channel. This does not exclude that an inhomogeneous strange-quark condensate is induced by the coupling to the light quarks, but this depends on the favored shape of the spatial modulations.
Our first main result is that the axial anomaly leads to a splitting between the CP of the phase transition for homogeneous matter and the LP associated with the inhomogeneous phase. This splitting turns out generally to be very small, with the exception of the limit of very small strange current quark masses, for which the CP moves toward the temperature axis and eventually disappears from the phase diagram, whereas the LP does not follow the same behavior. For a realistic parameter set fitted to vacuum phenomenology, the two points are nevertheless found to almost coincide for a wide range of values of the six-quark coupling K. It is important, however, to stress that the points do not coincide exactly, and for a wide range of parameters explored in this work, we find that the firstorder chiral phase transition which is typically present when restricting the model analysis to homogeneous matter is completely covered by the inhomogeneous phase.
A closer inspection of the behavior of the GL coefficients allowed us to better understand the behavior of the two points and the magnitude of their splitting, while our complementary stability analysis let us determine in several cases the location of the phase boundary where the inhomogeneous phase terminates. Consistently with the results of Ref. [16], we find an enhancement of the inhomogeneous phase due to the coupling with strange quarks, although, as already mentioned, the effect turns out to be quite small, and we did not witness the appearance of any additional inhomogeneous window in the phase diagram compared with the two-flavor case.
We recall that our study relies on the approximation of small amplitudes (and additionally of small gradients for the GL analysis), so it can only provide a reliable result if we are able to expand about a suitable homogeneous solution, which for simplicity we assumed to be the two-flavor restored phase. It would be tedious but straightforward to extend the method to expand about massive homogeneous solutions in the light-quark sector as well, as we have done in Ref. [14] for the two-flavor model away from the chiral limit. We note that in this case we could even study the stability of the homogeneous chirally broken phase against inhomogeneities, while this region was not accessible to our present analysis. But even then, we stay restricted to second-order phase transitions. It would therefore be desirable to perform a full numerical evaluation of the model thermodynamic potential in order to obtain the complete phase structure of the model. This would require the introduction of specific Ansätze on the spatial dependence of the chiral condensate and a full numerical diagonalization of the quark Hamiltonian, as done, e.g., in Ref. [30]. In principle, one could also perform an Ansatzfree minimization of the thermodynamic potential within lattice field theory, which has already been applied to analyze inhomogeneous phases in lower-dimensional models [31][32][33][34] but obviously becomes computationally extremely expensive in 3 þ 1 dimensions. Another way to go beyond the small amplitude constraint would be to perform an extended GL analysis, as done in Ref. [28], although the determination of higher-order GL coefficients would become significantly more involved due to the additional contributions in the model Lagrangian.
Ultimately, we are of course interested in the phase structure of QCD; i.e., one should go beyond effective models and the mean-field approximation. In this context, it is interesting to note that in Ref. [5] hints for an inhomogeneous phase have been found within a FRG study of the QCD phase diagram for 2 þ 1 flavors. It was found that in a certain regime in the vicinity and above the (homogeneous) chiral phase boundary the mesonic wave-function renormalization constant at vanishing momentum Z ϕ ð0Þ becomes negative. As the authors point out, this could indicate an instability toward a spatially modulated regime, but the signal is not fully conclusive. In fact, Z ϕ ð0Þ basically corresponds to the GL coefficient α 4;b in our language, so finding a negative value is a necessary but not a sufficient requirement for an inhomogeneous phase. This becomes clear if we recall that α 4;b is proportional to d dq 2 Γ −1 l ðq 2 Þj q 2 ¼0 . Hence, a negative value hints at a minimum of Γ −1 l at nonzero q 2 but does not tell whether its value at the minimum is positive or negative. Indeed, in a previous FRG study within the quark-meson model [35], it was found that the static pion two-point function develops a maximum (corresponding to a minimum of Γ −1 l ) at nonzero q 2 but does not change its sign, so no instability occurred.
In Ref. [5], the region with negative Z ϕ ð0Þ starts at a chemical potential well below the CP. Hence, if it was indeed related to an instability, it would mean that the CP and most likely the entire first-order phase boundary are covered by an inhomogeneous phase. Unfortunately, a fully conclusive analysis of this issue is extremely involved due to several competing order parameters. In this sense, the analysis performed in our work can act as a guide to give a qualitative picture of the phase structure of dense 2 þ 1flavor quark matter as a starting point for these more refined studies to begin with.