Absence of inhomogeneous chiral phases in 2+1-dimensional four-fermion and Yukawa models

We show the absence of an instability of homogeneous (chiral) condensates against spatially inhomogeneous perturbations for various 2+1-dimensional four-fermion and Yukawa models. All models are studied at non-zero baryon chemical potential, while some of them are also subjected to chiral and isospin chemical potential. The considered theories contain up to 16 Lorentz-(pseudo)scalar fermionic interaction channels. We prove the stability of homogeneous condensates by analyzing the bosonic two-point function, which can be expressed in a purely analytical form at zero temperature. Our analysis is presented in a general manner for all of the different discussed models. We argue that the absence of an inhomogeneous chiral phase (where the chiral condensate is spatially non-uniform) follows from this lack of instability. Furthermore, the existence of a moat regime, where the bosonic wave function renormalization is negative, in these models is ruled out.

In many of the above works the (chiral) order parameters are considered to be homogeneous in space. This is a reasonable choice for first investigations of the phase diagram. In strongly-interacting systems in condensed matter physics, however, crystalline-like ground states are quite common [75][76][77][78][79][80][81][82] suggesting that the assumption of homogeneous order parameters has to be reevaluated. More refined investigations of FF and Yukawa-type models have indeed revealed the presence of a chiral inhomogeneous phase (IP), where the corresponding order parameter is a function of the spatial coordinates, also in effective theories in nuclear and high energy physics, see Refs. [83][84][85], where IPs have first been found in this context. This has stimulated conjectures that inhomogeneous chiral condensates and related phenomena might also be relevant in the phase diagram of QCD. However, it has to be remarked that many of the existing model investigations have been carried out within the mean-field approximation, where the bosonic quantum fluctuations are neglected.
In 3+1 dimensions, many of the models for spontaneous chiral symmetry breaking in QCD feature an IP -typically appearing at low temperature and intermediate densities, where one would expect a first-order phase transition between homogeneous phases, see Ref. [104] for a review. However, these results are mostly obtained in the meanfield approximation and their predictive power for QCD can be questioned. In case of the FF models, we recently documented a regulator dependence of the IP in the NJL model [23], where the existence and shape of the IP depends on the chosen regularization scheme. In combination with the non-renormalizability of the NJL model this raises questions about the predictive power for QCD, especially since the chemical potentials in the region of the IP are in the order of the necessary regulator. The results in the renormalizable QM model suffer from other technical problems, such as instabilities at large field values [105,106].
Besides these difficulties of the model studies, there are still indications that inhomogeneous chiral condensates might be relevant in QCD [107,108]. A recent Functional Renormalization Group study of QCD [109] finds a socalled moat or Lifshitz regime, where the bosonic wave-function renormalization is negative and a modified dispersion relation is obtained [110]. Such a phenomenon is often related to an IP [111]. Experimental signals of this regime are discussed in Refs. [110,[112][113][114].
In 2 + 1-dimensional models, only recently the attention has moved towards the study of IPs [41,[115][116][117][118][119][120]. In the 2 + 1-dimensional GN model, IPs were observed at finite lattice spacing [115], but it turned out that these are regulator artifacts depending on the chosen regularization scheme [116,117]. After renormalization homogeneous phases are favored. This regulator dependence further highlights the problem of non-renormalizability with respect to high chemical potentials relevant in the IP of the 3 + 1-dimensional NJL model [23]. In lattice simulations of the 2 + 1dimensional GN model, oscillating, but damped correlators have been observed in Refs. [37,38], but considering the found regulator dependence of the IP [116,117] it is unclear whether these are reminiscent from the above described regulator artifacts or from underlying physical reasons, such as Fermi surface effects. We note that the existence of an IP is also not favored by the introduction of chiral imbalance [118,119] nor of a magnetic field [41].
In this paper, we rule out the instability of homogeneous condensates with respect to inhomogeneous perturbations for a large class of models with Lorentz-(pseudo)scalar isospin FF interactions. This extends the previous findings in the GN model to a variety of different FF and Yukawa models. We apply the framework of analyzing the stability of the bosonic two-point functions (as, e.g., used in Refs. [23,111,115,117,[119][120][121][122][123][124][125][126]) to these models. This analysis tests whether homogeneous field configurations are energetically unstable against inhomogeneous perturbations. Among the limitations of the method is that one can only show the existence of an IP and not the shape of energetically preferred inhomogeneous field configuration (see Ref. [111] for a detailed discussion of this method). As discussed in Ref. [111], there might still exist an IP, which is not detected by the stability analysis. However, to our knowledge no phase diagram has been found, where an IP does not at least enclose a region in parameter space with instabilities of the homogeneous minimum against inhomogeneous perturbations. Moreover, all IPs so-far observed in model studies are connected by a second-order phase transition to the chirally symmetric phase. Such a phase transition can always be detected by analyzing the stability of the symmetric minimum of the effective potential.
This work is outlined as follows. A general FF model, containing all the studied interaction channels, is introduced in Section II. Furthermore, we present an extension of these models to Yukawa models. In Section III, the stability analysis of the bosonic two-point function for finite baryon chemical potential and temperature is presented for the FF models and their extension to Yukawa models. Based on this analysis, we argue that all models with Lorentz-(pseudo)scalar interaction channels do not exhibit IPs or moat regimes. In Section IV, examples for these FF models are presented. Allowing for multiple chemical potentials, we also show the stability of homogeneous condensates for a few theories each containing a small subset of the previously discussed interaction channels. Finally, we sum up our results and conclude in Section V.

II. DEFINITION OF THE CONSIDERED MODELS
In this section, we introduce a FF model, which serves as the general prototype for the models, which will be later studied using the stability analysis introduced in Section III. We perform a bosonization with auxiliary fields and obtain the effective action. Finally, in Section II B, the Yukawa models, obtained from extending the bosonized FF models, are defined.

A. Four-fermion models
In order to set up our general analysis, we define as the most general FF action studied in this work in 2 + 1-dimensional Euclidean spacetime (x = (x, τ ) represents the spacetime coordinate). The integration over the periodic Euclidean temporal coordinate τ goes from 0 to β = 1/T , where T is the temperature, while the integration over d 2 x goes over the two-dimensional spatial plane. The vector ψ contains 2N four-component fermion fields (N identical spinors with isospin up/down respectively). We work in the chiral limit, so no bare mass term is introduced. The Dirac matrices are reducible, 4 × 4 representations of the 2 + 1-dimensional Euclidean Clifford algebra where I is the 4 × 4 identity matrix. Useful representations for computations in 2 + 1 dimensions can, e.g., be found in Refs. [117,119,127,128]. The interaction vertices c j are 8 × 8 matrices in isospin and spin space and elements of C = (c j ) j=1,...,16 = (1, iγ 4 , iγ 5 , γ 45 , ⃗ τ , i⃗ τ γ 4 , i⃗ τ γ 5 , ⃗ τ γ 45 ) , where ⃗ τ is the vector of Pauli-matrices acting on the isospin degrees of freedom. The spin matrices γ 4 and γ 5 anticommute with the γ ν , while γ 45 ≡ iγ 4 γ 5 commutes with the γ ν . All elements of C are 8 × 8 matrices, where the identity matrices in spin and isospin space are not explicitly written down whenever the matrix c j is diagonal in the corresponding space. The channels ψ (x) c j ψ(x) 2 are local, Lorentz-(pseudo)scalar FF interactions in the SU (2) isospin space and spinor space. The vector and matrix-like FF channels are not taken into account, since the analysis of these interactions differs technically from the (pseudo-)scalar ones. Therefore, these interaction terms will be left to forthcoming work. The couplings λ j of each channel have inverse energy dimension and will be set to either 0 or λ in order to study different models and allow for different symmetry groups of the action 3 . The chemical potential µ is introduced in the usual way and induces a non-vanishing baryon density ∝ψγ 3 ψ. It is well known, that the partition function of Eq. (1) is identical to a partition function with auxiliary bosonic scalar fields ⃗ ϕ, where ϕ j can be introduced via a Hubbard-Stratonovich transformation (an inverse shifted Gaussian integration) in order to get rid of the FF interaction ψ (x) c j ψ(x) 2 , up to a physically irrelevant integration constant [44]. The partially bosonized action is then given by where J is an index set containing 4 all integers 1 ≤ j ≤ 16 with λ j ̸ = 0. The auxiliary fields ϕ j have the dimension of an energy and fulfill the Ward identity 5 In order to relate to the literature about FF models as low-energy effective models (e.g. Refs. [1-3, 5, 8, 10, 12, 19, 21, 104]), we define phenomenologically motivated symbols for the ϕ j , i.e., the tuple The ordering of this tuple is used to directly map a field ϕ j to the corresponding fermion bilinearψc j ψ by reading of the index j in the tuple C in Eq. (3). As one can see from Eq. (4), non-vanishing ⟨ϕ j ⟩ give rise to dynamically generated mass terms for the fermion fields. These mass terms spontaneously break different symmetries of the action Eq. (4). The possibility of certain symmetry breaking schemes is, of course, dependent on the choice of interaction channels, i.e., on the set J (see the definition below Eq. (4)).
In Appendix A, we define the (chiral) symmetry group for 2 + 1-dimensional fermion fields. Considering no interaction terms, the fermion fields contained in ψ are invariant under transformations of the group U(4N ) composed of chiral transformation of the group U γ (2N ) and isospin transformations, which are elements of SU ⃗ τ (2). This invariance is not explicitly broken in the action (4) (or Eq. (1)) if λ j = λ for 6 j = 1, . . . , 16. When a subgroup of the interaction channels is taken into account in the partially bosonized action (4), only subgroups of the chiral symmetry transformations might be realized (depending on the choice of J). The reminiscent symmetry transformations, which are relevant in Section IV, are also defined in Appendix A.
A mass term of the form mψψ, as spontaneously generated by a non-vanishing expectation value of ϕ 1 = σ, breaks the full symmetry group U γ (2N ) down to the U I4 (N ) × U γ45 (N ) subgroup of vector transformations. Mass terms of the form im 4ψ γ 4 ψ (generated by ϕ 2 = η 4 ) and im 5ψ γ 5 ψ (generated by ϕ 3 = η 5 ) have the same symmetry breaking pattern, but, in addition also break one of the two Z 2 parity transformations 7 In addition, there is a fourth mass term m 45ψ γ 45 ψ, which does not break a continuous symmetry, but both P 4 and P 5 . This parity-odd mass can be dynamically generated by a non-vanishing expectation value of ϕ 4 = η 45 . These four different mass terms have an interpretation in condensed matter applications, such as graphene or similar systems [129][130][131], see Ref. [51] for a summary. We assign quantum numbers +, − to the auxiliary bosonic fields in Eq. (6) according to their transformation behavior under the parity transformations (7) and (8). The fields and their quantum numbers are listed in Table I. Later, we will refer to the fields by their behavior under parity, e.g., η 4 and ⃗ π 4 are fields with quantum numbers (−, +) and σ, ⃗ a 0 are (+, +) fields. After integration over the fermion fields in Eq. (4), one obtains an effective action, depending solely on the bosonic fields, where Q has been multiplied with β in order to ensure a dimensionless argument of the logarithm. This introduces only a temperature-dependent constant to the partition function. We recognize that S eff is proportional to N . Taking the limit N → ∞ for the FF models is equivalent to a mean-field approximation at finite N , which is what we consider for the rest of this work. In this approximation one takes the fermionic fluctuations fully into account through the functional Tr ln over the Dirac operator Q, while all quantum fluctuations in the bosonic degrees of freedom are neglected. This causes the global minimum ⃗ ϕ(x) = ⃗ Φ(x) of S eff to be the only relevant contribution in the partition function. Thus, expectation values of observables can be computed by evaluating them on the respective global minimum. This can become problematic when the effective action has multiple, degenerate global minima, for example in the case of a first order phase-transition or minima which are related by symmetry transformations on the level of the effective action. In the latter case, one formally has to introduce a small symmetry breaking parameter z and make the extrapolation z → 0 in order to remove ambiguities. In the mean-field approximation, however, it is common (see, e.g., Refs. [4,8,104]) and more practical to pick one of the degenerate minima whenever facing this situation. In the case of a first order phase-transition or critical end point, we will refrain from evaluating any quantities that depend on the minimum of the effective action in order to avoid ambiguities.

B. Yukawa models
In order to generalize our analysis of the FF models, as defined in Eq. (1), to corresponding Yukawa models in Section III, we introduce their action in 2 + 1-dimensional Euclidean spacetime as where ⃗ χ = (χ j ) j∈J contains scalar fields of canonical dimension energy 1/2 , h is the Yukawa coupling, κ n are the couplings of the self-interaction terms between the fields ⃗ χ and S eff is defined in Eq. (9). In the mean-field approximation, these models can be analyzed using the stability analysis in the same manner as the FF models, as is discussed in Section III B. Thus, quantum fluctuations of the fields ⃗ χ are neglected and observables are computed by evaluating them on the global minimum of S Y using the same formalism as described above for the FF models. A field χ j , that by its interaction channel in S eff corresponds to a field ϕ j in a FF model, has the same parity quantum numbers defined in Table I from ϕ j . Eq. (10) defines a QM-type of model. One can think of the fermion fields contained in S eff via the fermionic determinant as interacting through the exchange of the dynamical bosonic fields χ j , which in the QM model correspond to light mesons.
We note that the Ward identity for the expectation values of the scalar fields χ j , analogous to Eq. (5) for the auxiliary fields ϕ j , contains additional contributions from the self-interaction and kinetic terms. Thus, the expectation value of χ j does not directly yield the expectation value of a fermion bilinear (although these can nevertheless be computed). However, the spatial homogeneity of the bosonic fields still directly leads to the absence of inhomogeneous condensates. Consequently, the investigation of IPs in Yukawa models is possible by studying the stability of the fields χ j in the same way as we study the stability of ϕ j in the FF models.

III. STABILITY ANALYSIS OF MODELS WITH A BARYON CHEMICAL POTENTIAL
The stability analysis is a technique to determine whether a spatially homogeneous field configuration is unstable with respect to inhomogeneous perturbations. It follows from this instability that an inhomogeneous field configuration is energetically favored over the homogeneous expansion point. This technique has seen regular use in such investigations, e.g., in Refs. [111,117,123,126] and, thus, we recapitulate only the most relevant steps. We provide, however, a more detailed derivation in Appendix B. Note that the central result of this work, i.e., showing the stability of homogeneous condensates against inhomogeneous perturbations in all models with Lorentz-(pseudo)scalar interaction channels in a general manner, is discussed below in Section III A.
We start by considering homogeneous field configurations ⃗ ϕ and apply spatially inhomogeneous perturbations δ ⃗ ϕ to them, i.e., where the perturbations are of an arbitrary shape and assumed to be of an infinitesimal amplitude. Inserting this into the effective action (9) enables a systematic expansion of S eff in powers of δϕ. The zeroth order contribution S eff vanishes if the homogeneous expansion point ⃗ ϕ corresponds to a solution of the gap equation. The second order contribution contains the Hessian of the effective action, which is found to be diagonal in momentum space, but not necessarily diagonal in field space. The diagonalization -if possible -is done via a change of the field basis ⃗ ϕ → ⃗ φ. For the theory considered in Eq. (9), where only a baryon chemical potential is present, the Hessian is already diagonal in field space, i.e., ⃗ φ = ⃗ ϕ, if we set all non-zero λ j to λ. Nevertheless, we use the "new" basis from here on in order to be consistent with the later analysis, where a proper diagonalization is indeed needed. We then obtain as the second order contribution where δφ(q) are the Fourier coefficients of the perturbations and the magnitude of the perturbation's momentum q = |q|. The bosonic two-point function Γ (2) φj M 2 , µ, T, q 2 is the curvature of the effective action with respect to |δφ(q)|. For the curvature we find the explicit form whereν n = (ν n − iµ), the fermionic Matsubara frequencies ν n = 2π(n − 1 2 )/β and We find that a φj = a + = 4 for fields φ j with parity quantum numbers (P 4 , P 5 ) = (+, +), (−, −) and a φj = a − = 0 for fields with (P 4 , P 5 ) = (+, −), (−, +). Accordingly, we find two possible momentum dependent contributions to the two-point function. The integral ℓ 2 at T = 0 assumes the simple form Further expressions for ℓ 1 and ℓ 2 for various cases of µ, T, M 2 , q 2 can be found in Appendices B 1 and B 2. Also at finite temperature L2,−(M 2 , µ, T, q 2 = 0) = 0. However, the 0.0 contour line and its label are not drawn, because they would be obscured by the axis. In the ancillary files to this paper, we provide a script, which allows to do the numerical computation of L2,± for arbitrary values of (M 2 , µ, T, q 2 ) and that produces this figure. See Eqs. (13) and (15) for the definition of L2,±.
In order to make conclusions about whether inhomogeneous condensates are favored over homogeneous ones, we use the global homogeneous minimum as the homogeneous expansion point, i.e., the field configuration ⃗ Φ that minimizes S  φj (M 2 , µ, T, q 2 ) will signal that there is an even deeper minimum in the direction of δφ j (q). In order to provide dimensionless quantities, we express all parameters in units of the mass M 0 , which is the mass corresponding to the global homogeneous minimum in the vacuum of the theory.

A. Absence of instability
The global homogeneous minimum, that is the only relevant expansion point when searching for an IP, is stable against homogeneous perturbations, i.e., Γ (2) φj (M 2 , µ, T, q 2 = 0) ≥ 0 ∀φ j . Consequently, negative values for finite q 2 and, thus, an instability against inhomogeneous perturbations can be ruled out, if the q-dependent part L 2,φj of Γ (2) φj (M 2 , µ, T, q 2 ) is a monotonically increasing function of q 2 . This is exactly the case for L 2,± . For T = 0 the analytical form of L 2,± given by Eqs. (15) and (16) reveals its monotonically increasing behavior for all M 2 , µ, q 2 . The same is true for finite temperatures which can be verified by numerical calculations of ℓ 2 . In Fig. 1, the functional behavior of L 2,± for µ/M 0 = 1.0 and T /M 0 = 0.0, 0.05 is plotted in a color plot. We conclude that the general model defined by Eq. (9) does not experience an instability towards an IP. This conclusion holds when considering any set of interaction channels as given by J. Moreover, these models also do not feature a so-called moat regime [110]. This regime is defined by a negative wave-function renormalization Z ∝ d 2 Γ (2) /dq 2 and is often an accompanying phenomenon to an IP but can also exist independently [111]. Our computations show that L 2,± always yields Z ≥ 0, since they are monotonically increasing functions in q.
Additionally, we note that ℓ 2 ∼ Θ µ 2 − M 2 − q 2 /4 at T = 0, which causes the two-point function to be constant for momenta 0 < q 2 /4 < µ 2 − M 2 . Here not only the wave function renormalization Z vanishes, but also any higher derivative of Γ (2) with respect to q. A special point is at µ/M 0 = 1 and T = 0, where the homogeneous phase transition in the 2 + 1-dimensional GN is often regarded as being first order. However, the system rather exhibits a critical endpoint at this point [29]. This means that the effective potential (e.g. computed in Refs. [26,117]) is flat between some homogeneous field values with M ∈ [0.0, 1.0]. Such a flatness is also observed in the two-point function where the contribution 1/λ − ℓ 1 vanishes, which causes the two-point function to be constant zero for momenta 0 < q 2 /4 < µ 2 − M 2 . This vanishing curvature of the effective action is a hint for a degeneracy between homogeneous and inhomogeneous condensates. Such a behavior has already been observed in the case of the 2 + 1-dimensional GN model also at µ/M 0 = 1.0, T /M 0 = 0, which has been revealed by a study [132] with an one-dimensional ansatz for the chiral condensate. We expect this degeneracy to be restricted to (µ, T )/M 0 = (1.0, 0.0) and, especially, the homogeneous condensates to be favored against inhomogeneous ones for all non-vanishing T . The analysis in Ref. [132] suggests that also higher orders of the expansion of the effective action would vanish for a certain range of finite momenta. In Section IV A, we show that the homogeneous phase diagram of all possible models described by Eq. (9) is identical to the one of the 2 + 1-dimensional GN model. Thus, the above discussion of the flatness of the effective potential and the two-point function also applies to these models.
It should be noted that the absence of an instability does not completely rule out the existence of an IP. As discussed in Ref. [111] at the example of the 1 + 1-dimensional GN model, it is possible for the homogeneous minimum and global inhomogeneous minimum to be separated by a energy barrier. Here, the homogeneous minimum appears stable against inhomogeneous perturbations even though an IP is energetically favored. Such a phase can only be found by calculations with a guess of ansatz functions for the condensates or by explicit numerical brute-force minimizations using lattice field theory. For the 2 + 1-dimensional GN model, there is evidence in the literature that this is not realized. Lattice minimizations of this model have not found any other IP at finite lattice spacings than the ones also obtained by a stability analysis of the lattice regularized models [117,119,120]. It is important to note that these IPs vanish when taking the continuum limit [116,117], as the bosonic two-point function converges towards a momentum dependence proportional to L 2,+ . These results of the 2 + 1-dimensional GN model suggest that the stability analysis applied to our general FF model (9) also does not miss an IP. Also, to our knowledge all IPs, which are observed in model investigations, can at least in some parameter region be detected by the stability analysis. In these studies, one finds a second order phase transition between the IP and the chirally symmetric phase. At least, this second-order transition can always be detected by analyzing the stability of the symmetric minimum of the effective potential. The sum of these arguments combined with our analysis is strong evidence for the absence of inhomogeneous condensates in all of the models described by Eq. (9) (or, equivalently, Eq. (4)) independent of the considered (sub)set of interaction channels given by J.
Since this result is obtained in the mean-field approximation, one needs to consider its predictive power for the full quantum theories. It is typically found that bosonic fluctuations tend to disfavor and/or disorder (in-)homogeneous condensation 8 [41,94,95,98,99,134]. Based on these findings, we expect that the non-existence of IPs in 2 + 1dimensional FF and Yukawa models in the mean-field approximation is a clear signal that inhomogeneous ground states are not present in the corresponding full quantum theories. A scenario with an inhomogeneous ground state in a full quantum theory which is not present in the mean-field approximation has -to our knowledge -never been observed. In 2 + 1-dimensional FF models, one has so far only seen oscillating but also damped correlation functions for the auxiliary fields [38]. However, these lattice results [38] have been obtained on rather crude lattices. So, we expect these oscillations to be reminiscent of the IP, which can be found in mean-field investigations of FF models at finite lattice spacings, but vanish in the continuum limit [116,117].

B. Generalization to Yukawa models
When the models are only subjected a baryon chemical potential, it is rather straightforward to generalize the stability analysis of FF models to Yukawa models, which are defined as in Eq. (10). We only discuss the final result here and a more detailed derivation can be found in Appendix B 3.
Applying the stability analysis to the effective action as defined in Eq. (10), the second order contribution again is the first non-vanishing correction to the homogeneous action S (0) Y , when the homogeneous expansion point is a solution of the gap equation. We find the second order contribution to the effective action where ζ is the field basis that was used in the diagonalization of the Hessian 9 . The two-point function can now be identified as the curvature of the effective action in the field direction δζ j (q). Again, one finds either L 2,+ or L 2,− as the momentum dependence of the two-point function. As is obvious from Eq. (18), the additional contributions compared to the FF two-point function are either constant or monotonically increasing in q 2 . Thus, by the reasoning given in the last section, these additional terms cannot facilitate the appearance of an instability towards an IP since the corresponding FF model does not exhibit such an instability for any ⃗ Φ, µ, T, q 2 already.

IV. RESULTS OF THE ANALYSIS FOR SPECIFIC MODELS
In this section, we present examples of FF and Yukawa models where the condensates do not develop an instability towards inhomogeneous perturbations and are, thus, very unlikely to feature an IP. In Section IV A 1, FF models with only a baryon chemical potential but multiple interaction channels (see Eq. (1)) are presented. These results are obtained using the stability analysis, as explained in detail in Section III, where also the reasoning for the absence of instabilities is explained in detail. After that, we allow for multiple chemical potentials in Section IV A 2 and explain the differences of the analysis compared to only a baryon chemical potential. Again, examples for model calculations are discussed. Finally, we turn towards the extension of our findings to Yukawa models in Section IV B.

A. Four-fermion models
Before we turn towards the stability analysis of the bosonic two-point functions, we shortly present our finding for the homogeneous phase diagram of FF models described by Eqs. (1) and (4) (by considering different interactions channels in the set J).
The homogeneous phase diagram of the 2 + 1-dimensional GN model, i.e., the phase diagram when restricting the field to homogeneous field configurations, as first obtained in Ref. [26], features a phase at low temperature and chemical potential where the discrete chiral symmetry is spontaneously broken by a non-zero chiral condensate [24,26]. At finite temperature this phase is separated from the chirally symmetric phase by a second order phase transition. At The homogeneous phase diagram for all possible models, that can be obtained from Eq. (1) by setting certain couplings either to zero or λ, is identical to that of the GN model. This result is related to the present symmetries of the action, which allow to pick homogeneous minima, where only the scalar channel σ develops a non-zero expectation value (note that this is only possible when the fields are restricted to being homogeneous). Other homogeneous, global minima are connected to this one via a (chiral) symmetry transformation. Exempt from this are the fields with parity quantum numbers (−, −) namelyη 45 and ⃗ π 45 . By analyzing the homogeneous effective potential, we can nevertheless show that the homogeneous fieldsη 45 , ⃗ π 45 develop vanishing expectation values for all temperatures and baryon chemical potentials. First, we discuss this allowing for all isoscalar channels, which yields the introduction of the auxiliary fields σ, η 4 , η 5 , η 45 . Using the symmetries of Eq. (9) (or equivalently Eq. (4)) one can use the chiral transformations Eq. (A5) to obtain non-vanishing expectation values only for σ and η 45 , Then, we analyze the fermion contribution Tr ln Q, which can be represented in a block-diagonal form with 2 × 2 blocks corresponding to the irreducible spin representation in 2 + 1 dimensions. These blocks can be interpreted as GN model contributions in the irreducible fermion representation with chemical potential µ and homogeneous fieldsφ L/R ≡σ ±η 45 (see the theory sections of Refs. [117,119] for a more in-depth discussion of these irreducible blocks). The model decomposes into two GN models with the same chemical potentials, which gives usφ L =φ R and, thusη 45 = 0. An analogous analysis involving some additional isospin rotations is valid for ⃗ π 45 . Thus, we do not observe spontaneous parity breaking in the FF models for all µ and T .
The order of the homogeneous phase transition might be relevant, since the stability analysis of the bosonic twopoint function in the 1 + 1-dimensional GN revealed that a first order phase transition can cause stability of the homogeneous minimum even though there is an energetically preferred IP [111]. Here, the global energetically preferred inhomogeneous minimum was only connected to the symmetric homogeneous field configuration which was not energetically favored over the stable non-zero minimum. The fact that in 2 + 1 dimensions the critical point occurs at a single point at zero temperature increases the confidence in the presented results obtained via the stability analysis. Moreover, by inspecting all possible expansion points via varying the mass M , we can rule out a scenario as in 1 + 1 dimensions, since in 2 + 1 dimensions there is no other unstable expansion point while the physically relevant expansion point M is stable.

Models with baryon chemical potential
Table II summarizes the central findings from the stability analysis for some models that are obtained from Eq. (4) by defining J and setting λ j = λ, j ∈ J. The specific models are chosen with respect to their relevance to phenomenology and the literature. The two-point functions for models with other combinations of channels can easily be calculated as discussed in Section III. The respective symmetry transformations of the effective actions allow to obtain homogeneous expansion pointsφ j , which are vanishing except forσ. None of the two-point functions in these models exhibit a different momentum dependence than L 2,± . As discussed in Section III and illustrated by Fig. 1, L 2,± are positive, monotonically increasing functions of the square of the momentum q of the inhomogeneous perturbations. Thus, we expect the absence of IPs and moat regimes in all of these models, as discussed in detail in Section III A. Therein, we also explain that there is the possibility of a degeneracy between the homogeneous condensate and the inhomogeneous condensate, as observed at (µ, T )/M 0 = (1.0, 0.0) in the GN model [132]. By our analysis, we expect this degeneracy to be restricted to (µ, T )/M 0 = (1.0, 0.0) and, especially, the homogeneous condensates to be favored against inhomogeneous ones for all non-vanishing T .
In Table II, we discuss four FF models explicitly and give their respective symmetry groups (defined in Appendix A) and indicate whether the momentum dependence of the two-point functions are L 2,+ or L 2,− .
The first row shows the renowned GN model with one single isoscalar channel with parity quantum number (+, +) (see Eq. (7) and Eq. (8) for the definition of the parity in 2 + 1 dimensions and Table I for the quantum numbers of the fields). The two-point function exhibits the L 2,+ dependence, which was already documented in Ref. [117].
The second row shows the 2 + 1-dimensional analog of the 3 + 1-dimensional NJL model. It breaks the axial symmetry transformations Eqs. (A3) and (A4), leaving only the combined isospin and axial symmetries Eqs. (A7) and (A8) as chiral transformations. Due to the ambiguity of the γ 5 operator in 2 + 1 dimensions, it makes sense to use the generators of both transformations to construct an "NJL" Lagrangian. Similar to the case in 3 + 1 dimensions, we find that the momentum dependence of the two-point functions of the π-fields are given by L 2,− , while it is L 2,+ for the σ field [23,126].
The third row shows the 2 + 1-dimensional chiral Heisenberg-Gross-Neveu (χHGN) model [62] with an additional ψ γ 45 ψ 2 interaction term that is special to 2 + 1 dimensions 10 . Again, due to there being two axial transformations Eqs. (A3) and (A4), a Lagrangian involving both generators is an appropriate choice. A non-zero bosonic field η 45 breaks both parity symmetries Eqs. (7) and (8)    fields do not mix (as the only off-diagonal terms are ∝η 45 , compare with Eq. (B15)) and the diagonalizing field basis coincides with the auxiliary bosonic fields introduced in the bosonization. Thus, one of the two parity transformations could only be spontaneously broken if one of the fields with a negative parity quantum number (compare Table I) develops a spatial modulation. However, also for this model we found that the momentum dependent part of the two-point functions of all present bosonic fields is either L 2,+ or L 2,− and, thus, homogeneous condensates never develop an instability towards a spatially dependent condensate.
The last row lists what we call the complete Lorentz-(pseudo)scalar four-fermion (PSFF) model, since it features all interaction channels present in Section II and Eq. (1). Again, note thatη 45 = ⃗ π 45 = 0 which prevents off-diagonal second order terms from the fermionic contribution, see Eq. (B16). Even though we considered by far the largest amount of interaction channels in this model, no two-point function with a different momentum dependence than L 2,± is obtained.
As discussed in Section III A, two-point functions for fields, which have mixed parity quantum numbers ((P 4 , P 5 ) = (+, −) or (P 4 , P 5 ) = (−, +) ), have a momentum dependence proportional to L 2,− , while the others have a momentum dependence proportional to L 2,+ . This is independent of how large the symmetry group is and of the number of interaction channels considered. Even when considering all 16 Lorentz-scalar FF interactions (see Eq. (3)), we do not see a different mathematical structure of the two-point functions.
Therefore, according to the argument given in Section III A, none of the models, which can be described by Eq. (4) by defining J and setting λ j = λ, j ∈ J, exhibits an instability towards an IP for any M 2 , µ, T, q 2 . As argued in detail in Section III A, this is strong evidence for the absence of IPs in these models.

Models with multiple chemical potentials
In this section, we allow for multiple chemical potentials in FF models in addition to the baryon chemical potential, which is introduced in Eq. (1) via the usual µψγ 3 ψ term. Namely, we will study the effect of finite isospin chemical potential, introduced with the term µ Iψ γ 3 τ 3 ψ, and chiral chemical potential, that is introduced as µ 45ψ γ 3 γ 45 ψ. Although this requires an extensive analysis, we can again identify the well-known function L 2,+ in the two-point functions. Thus, also the models discussed below do not develop an instability towards an IP. Also, the existence of a moat regime with a negative wave-function renormalization can be ruled out, as Z ≥ 0 according to the discussion in Section III.
In general, the introduction of multiple chemical potentials induce that (homogeneous) expectation values of several of the auxiliary fieldsφ i will have non-zero expectation values, e.g., studying finite µ 45 can lead to a non-vanishing expectation value ofη 45 . Also, the homogeneous phase diagram has to be computed separately for each individual case and can have a more involved phase structure (compare, e.g., Ref. [119]).
This also significantly complicates the analysis of the bosonic two-point functions, although in principle the method can be applied as introduced in Section III. The main difficulty is the diagonalization of the bracketed field space matrix in Eq. (B16), which becomes quite involved depending on the studied action. For a large number of interaction channels and multiple non-vanishingφ i it can become even impossible to diagonalize this expression. Thus, a general analysis for multiple theories, as in Section III for FF models with baryon chemical potential, cannot be presented by us.
However, there are still some combinations of interaction channels and chemical potentials, where one can still obtain the momentum dependence of the two-point functions by a comparatively easy diagonalization obtained by a suitable choice of field coordinates φ j . For these models, we find the fields φ j such that the fermion propagator (B8) is block-diagonal with 2 × 2 blocks, where each corresponds to an irreducible spinor representation. Then, the contribution of the fermion-loop (as written down in, e.g., Eqs. (B12) and (B14)) also decomposes when ∆Q is written as a function of the δφ j . The vertices 11 corresponding to the interaction of the fermions with the new field basis φ j project out either one or multiple of the 2 × 2 blocks and diagonalize the bracketed field space matrix in Eq. (B16). Then one can proceed with the analysis as described in Section III and Appendix B. Table III summarizes the models, where such an analysis with respect to the stability of homogeneous condensates against inhomogeneous perturbations was performed by us. All two-point functions obtained in the models in Table III are proportional to L 2,+ . This originates in the circumstance, that in our derivation of the two-point function the field variables φ j are chosen such that one obtains a block-diagonal structure of the homogeneous fermion propagator (B8). These blocks behave as 2 + 1-dimensional GN models (see Table II for the momentum dependence of the two-point function of the σ field in the GN model) with a single effective chemical potential and a single field φ j or sums of the fields. The effective chemical potentials are linear combinations of µ, µ 45 , µ I . The momentum-dependence of the obtained two-point functions is given by linear combinations of L 2,+ with M 2 and effective chemical potentials.
The model in the first row is a GN model with an additional ψ γ 45 ψ 2 interaction channel and a chiral imbalance via finite µ 45 . A standard GN model subjected to µ 45 without the additional interaction channel has been studied in Ref. [119]. The corresponding homogeneous expectation value of the field η 45 , as induced by finite µ 45 breaks parity spontaneously. However, the new field basis ϕ L/R effectively decouples the theory into two, independent GN models with chemical potentials µ ± µ 45 . This leads to two "GN-like" two-point functions each with one independent field and chemical potential. The same procedure can be performed for the model in the second row, where instead of a ψ γ 45 ψ 2 the isospin channel ψ τ 3 ψ 2 (leading to an auxiliary bosonic field a 0,3 in the bosonic theory) is added to the GN interaction and an isospin chemical potential µ I introduces an isospin imbalance.
In the third row, we allow for both isospin and chiral imbalance introducing both corresponding chemical potentials in addition to the baryon chemical potential µ. A ψ γ 45 τ 3 ψ 2 interaction channel is studied in addition to the ψ ψ 2 interaction. The corresponding bosonic auxiliary fields π 45,3 and σ can again be combined via linear combination to obtain the diagonalizing field basis φ ± . The two-point functions Γ (2) φ± are now sums of different L 2,+ contributions, each with M 2 = Φ 2 ± , but with differing linear combinations of the chemical potentials dictating the exact form of L 2,+ . Again, we refer to Ref. [119] where an analogous competition between two chemical potentials is studied in the phase diagram and, although this is not explicitly computed therein, in the two-point function of the σ field.
In the forth row, the most involved model in Table III is considered containing both previously introduced interactions resulting in the presence of the auxiliary bosonic field η 45 and a 0,3 in addition to the σ channel simultaneously at a finite baryon, isospin and chiral density. We find the diagonalizing field basis to be ϕ L , ϕ R , a 0,3 . The momentum dependence of the two-point functions of ϕ L/R are again sums of two L 2,+ contributions each with different linear combinations of the chemical potentials and M 2 = φ L/R ±ā 0,3 2 . However, Γ a0,3 is proportional to the sum of Γ (2) ϕ L and Γ (2) ϕ R , combining four different contributions each proportional to L 2,+ . Concluding, all the studied FF models in Table III do not exhibit an instability of the homogeneous ground state when subjected to inhomogeneous perturbations and, thus, it is very unlikely that they feature an IP (c.f. Refs. [111,116,117] and Section III A). Nevertheless, it is important to state that in the investigation of FF models subjected to multiple, non-vanishing chemical potentials we restricted ourselves to a very limited set of interaction channels. We want to highlight at this point that the restriction to a few interaction channels also limits the predictive power of our models for high energy phenomenology. For example, at finite isospin chemical potential one needs to take into account for charged pion condensation (see, e.g., Ref. [135]). Attentive readers may notice that the corresponding channels are not present in the models in Table III. In order to provide an adequate description of this phenomenon one certainly has to extend our study in this direction. To study a larger set of interactions with several of the chemical   (3)) in the model. The second column gives the corresponding auxiliary bosonic fields after bosonization, that correspond to the fermion bilinear via the Ward identity (5). The third column lists the used chemical potentials. In the fourth column, the field basis φj is defined, which diagonalizes the second order correction, as described in detail in Section IV A 2. Then, the momentum dependence of the two-point functions Γ (2) φ j is given in the fifth column, right to the field definition. The last column givens an overview of the present symmetries in the model. The groups are clickable and refer to the definition of the symmetry group. potentials µ, µ 45 , µ I or even more axial imbalances 12 in models with more interaction channels, is in principle possible but requires a different, technically more involved analysis than the one done in this work. Thus, we postpone such an analysis to future works. However, our study allows us to conclude that the presence of multiple imbalances in fermion densities does not generically allow for the existence of inhomogeneous ground states.

B. Yukawa models
In this section, the results for the Yukawa model extension of the FF models discussed above are presented. The generalization of the bosonized FF models (4) to a Yukawa model is discussed in Section II B. In Eq. (10), the Yukawa action is constructed out of the FF models' effective action (9).
The homogeneous phase diagram of these models might drastically change compared to the FF models due to the introduction of additional couplings and self-interaction terms. Nevertheless, the stability analysis, as described in Section III, can be performed for all possible, homogeneous expansion points χ j =χ j such that conclusions about the stability of these homogeneous condensates against inhomogeneous perturbations can be made.
First, we discuss the generalization of the results in Section IV A 1, where the FF models are studied at finite baryon chemical potentials µ. The key observation is that all bosonic two-point functions are positive, monotonically increasing functions of q 2 and, thus, no instabilities are observed. This observation is directly obtained for the Yukawa extensions of these models, as these only affect the momentum independent offset (compare Eq. (18)) and the physical expansion point M. The momentum structure of the two-point function is still proportional to L 2,± plus an additional positive, monotonically increasing q 2 term. As discussed in Section III A and illustrated by Fig. 1, there is no expansion point for which L 2,± is not a monotonically rising function of q 2 . Thus, we can conclude that there exist no instabilities towards an IP for the Yukawa models, which are generated by extending the FF models in Table II. According to the reasoning in Section III A, the moat regime with a negative wave-function renormalization is also ruled out. Nevertheless, it might be possible, that such a Yukawa model features a first order phase transition between the homogeneous phases. Then, one might also find the existence of a first order phase transition towards an IP. In many model calculations the IP covers the homogeneous first order transition [86,104,105,125,126] and, thus, the existence of such a transition could also mean that an IP exists in these models, which is not detected by our analysis. To our knowledge, however, a model, which features no instability towards an inhomogeneous perturbation in the whole phase diagram, but which still has an IP, has never been observed before.
Next, we discuss the extension of the FF models from Table III with multiple chemical potentials to Yukawa models. In the case of multiple chemical potentials, it is not straightforward to derive the two-point functions for the Yukawa models, corresponding to the FF models by Eq. (10), starting from the FF model results. This is caused by the necessity of multiple non-vanishing homogeneous expectation valuesχ j as expansion points for the analysis when the models are subjected to multiple chemical potentials. Then, the Yukawa self-interactions of the bosonic fields cause non-vanishing second order contributions, which are off-diagonal in the field perturbations δχ j (compare Appendix B 3 and Eq. (B31)), in addition to the off-diagonal fermion contributions Eq. (B16). This is the case for all Yukawa extensions of the models in Table III.
In the analysis of such a model one needs to diagonalize the whole second-order contribution in Eq. (B31). Typically, we choose a basis ζ j for the dynamical scalar fields that corresponds to the field basis φ j , which diagonalizes the FF model part. Then, one still needs to diagonalize the off-diagonal contributions coming from the bosonic selfinteractions. The diagonalization can become very complicated, even when using symbolic diagonalization of the expressions as provided by Matlab [136]. We provide an example of this procedure in Appendix C, where we discuss the Yukawa extension of the model in the first row in Table III. In this analysis, we find a more complicated structure of the two-point functions and a diagonalizing field basis, which depends on the studied momentum q of the perturbation as well as the homogeneous expectation values of the scalar fields and the external parameters such as temperature and the chemical potentials. Nevertheless, the obtained two-point functions can be shown to be monotonically increasing functions in q and, thus, again no instabilities towards inhomogeneous perturbations are observed. For the other three models in Table III, we do not show an explicit calculation as the expressions become very lengthy and the analysis becomes very involved. However, by first calculations and due to the fact that the off-diagonal contributions coming from the scalar fields' self-interaction terms are not dependent on the momentum q of the inhomogeneous perturbations we expect that also these Yukawa models do not develop an instability towards an IP.

V. CONCLUSION AND OUTLOOK
In this work, we analyzed the stability of homogeneous condensates against inhomogeneous perturbations in a wide range of 2+1-dimensional four-fermion (FF) models and their Yukawa extensions under the influence of combinations of baryon chemical potential, isospin chemical potential and chiral chemical potential. All investigations were performed in the mean-field approximation, i.e., neglecting bosonic quantum fluctuations. The most involved model features 16 Lorentz-(pseudo)scalar FF interaction channels and the other models are subsets of this model (see Eqs. (1) and (3) for the FF model).
Our main finding is the stability of homogeneous condensates against inhomogeneous perturbations in these FF and Yukawa models at any finite baryon chemical potential and temperature. As argued in Section III A, this is strong evidence that IPs do not exist in 2 + 1-dimensional FF models with Lorentz-(pseudo)scalar interaction channels and their corresponding Yukawa extension (the correspondence is given in Section II B). Also, we can completely rule out the existence of a moat regime whose characteristic is a negative wave-function renormalization Z (see Section I for a short discussion). The momentum dependence of the obtained two-point functions only allow Z ≥ 0. We draw similar conclusions for models with a subset of the above mentioned FF interaction subjected to multiple chemical potentials. In the case of multiple chemical potentials, the extension to Yukawa models causes technical difficulties, but, in principle, we expect the result regarding stability of homogeneous phases to be the same. This is motivated in Section IV B and an example computation for an explicit model is provided in Appendix C. We note, however, that the used models are not very well suited for high energy phenomenology in the presence of multiple chemical potentials, as they lack, for example, the description of charged pion condensation for finite isospin chemical potential. Thus, one should interpret our results as a first step to generically show that multiple imbalances of fermion densities do not result in the existence of an IP.
Our results suggest that there might be a general argument, similar to Derrick's theorem [137], behind the absence of IPs in 2 + 1-dimensional models. Such a principle could possibly be found by studying the properties of a general Ginzburg-Landau free energy as done in Ref. [138]. Therefore, one would probably need to encode the fermionic determinant in our models by allowing for higher orders in the gradient expansion.
In the present study, we did not consider vector FF interactions, since their inclusion yields a technically different analysis than the one presented in Section III. The inclusion of these channels would also be relevant for low-energy effective theories of QCD, where both vector and scalar FF interaction should emerge. However, the interplay of (repulsive) vector interactions and the scalar interactions could play a crucial role in developing an instability towards an IP, especially as vector fields directly couple to fermionic momenta in the calculation 13 . In 3+1-dimensional models it was already observed that vector interactions can enlarge the IP [139,140]. The presence of such interactions causes additional contributions to appear in the two-point functions of the corresponding bosonic fields, which do not have a definite monotony as the momentum dependence of the two-point functions in the present work. Calculations for 2 + 1-dimensional FF models with vector interactions are already underway and we hope to report soon about the results.
With respect to high energy phenomenology, there are many ways to build on the present work -for example allowing for finite quark-masses, studying three-flavor versions of the FF models allowing for strangeness effects or the inclusion of bosonic quantum fluctuations using lattice field theory. With respect to the stability analysis, it is interesting to note that using finite bare quark masses disfavors an IP [126] and the inclusion of strange quarks only has marginal effects on inhomogeneous condensation [141]. However, the effect of bosonic quantum fluctuations beyond the mean-field approximation on the phase diagram is a very important aspect and also currently discussed in the literature [41,94,95,98,99,134]. The general result of these studies is that bosonic quantum fluctuations tend to disfavor ordered phases (such as condensed phases) similar to the effect of thermal fluctuations. Thus, it can be anticipated that the absence of an IP in a mean-field calculation implies that there will be no inhomogeneous ground states in the full quantum theories. To our knowledge, there is no observation of an inhomogeneous ground state being generated by bosonic quantum fluctuations when there is no IP present when these fluctuations are suppressed. An interesting aspect for 2 + 1-dimensional theories in general could be the inclusion of difermion interactions in order to allow for color-superconducting order parameters, which often are in competition with the chiral ones, see, e.g., Refs. [11,72,142,143].
Regarding the general understanding of strongly-interacting fermions in 2+1 dimensions, it would be very interesting to study the interference and competition of the different chiral chemical potentials, which can be introduced as µ 4ψ γ 3 γ 4 ψ, µ 5ψ γ 3 γ 5 ψ or µ 45ψ γ 3 γ 45 ψ (as done in the beginning of Section IV A 2) and correspond to different conserved charges. However, the first results in our work and in Ref. [119] suggest that this would not change the phase diagram with respect to the existence of IPs. In this section, we want to elaborate on the symmetries of the fermionic models, which are introduced in Section II in Eq. (1). Their partially bosonized action can be found in Eq. (4) and are invariant under certain symmetry transformations acting on the 2N four-component spinor fields (the factor of 2 comes from introducing an isospin space). The specific symmetry group, under which the model is invariant, is determined by the choice of the set J, i.e. by the choice λ j ̸ = 0, j ∈ J and, correspondingly, by which auxiliary bosonic fields ϕ j , j ∈ J have to be introduced in the bosonization. This section is intended to provide an overview of the relevant symmetry groups present in the models discussed in Section IV and the tables Table II and Table III. Taking J = {j} j=1,...,16 the N isospin up/down fermion fields in Eq. (4) will be invariant under transformation of the group U(4N ), similar to the free case. This group has (4N ) 2 generators, which we split up into three different categories. There are rotations within the internal space of the N spinors generated by N × N matrices given by the generalized Gell-Mann matrices T a with a = 1, . . . , N 2 − 1 and the identity matrix T N 2 = I N . These generators T a can be combined with any of the chiral symmetry transformations of a free fermion fields, which are generated by (I 4 , γ 4 , γ 5 , γ 45 ). Together with the internal rotations either of these can generate a U(N ) symmetry, labeled by the corresponding chiral generator, namely with real parameters α a , β a , ζ a , ι a . Thus, the whole chiral symmetry group can be defined as where U is a matrix element of U(2N ). In addition, the FF model in Eq. (1) is invariant under isospin transformations of the group U(2), which is composed of a U(1) phase factor and where ⃗ ξ are three real parameters.
Some of the models discussed in Section IV are only invariant under a subgroup of the U(4N ) transformations. We will define these symmetry transformations as a reference for Table II and Table III. As in the NJL model, one can choose J such, that one or more of the two axial transformations (A3) and (A4) are broken. Then, the action can be still invariant under a combined isospin and chiral rotation given by where the isovectors ⃗ ι ′ a and ⃗ ζ ′ a contain 3 real parameters and a = 1, . . . , N 2 . Typically in these cases, the vector symmetry Eq. (A6) is not broken.
Some of the choices of J break the continuous chiral symmetries to discrete subgroups Typically, when a model is in addition invariant under the transformation (A2), then either of (A9) and (A10) can be reproduced by a combination of (A2) and the other of the discrete symmetries. In Table III, only a remnant of the isospin symmetry transformation SU ⃗ τ (2) is present in some of the models, namely Appendix B: Derivation of the stability analysis In this section, the stability analysis for a general FF theory subjected to a baryon chemical potential is derived. This discussion is similar to model-specific discussions, as, e.g., found in Refs. [111,117,126]. The core idea of this technique is to analyze the stability of a homogeneous field configuration under inhomogeneous perturbations. Furthermore, we outline how an extension to Yukawa models can be done.
Consider an expansion of the auxiliary bosonic fields where δ ⃗ ϕ(x) is a spatially dependent inhomogeneous perturbation of the homogeneous expansion point ⃗ ϕ. These perturbations are of an arbitrary shape and assumed to be of an infinitesimal amplitude. The Dirac operator then also separates into a homogeneous and an inhomogeneous part which we can use to expand the ln DetQ as where Tr denotes a functional trace over all spaces. We insert the expansion from Eq. (B3) and Eq. (B1) into the effective action Eq. (9) to obtain the expansion where S (n) eff contains all terms of order n in the perturbations δϕ j , i.e., terms ∝ j δϕ mj j with j m j = n. The first three terms in the series are given by where the zeroth order term is proportional to the effective potential. In the position space representation, the fermion propagator depends only on the difference of two space-time variables, i.e.,Q −1 =Q −1 (x, y) ≡Q −1 (x − y). The functional traces are represented in position space as with tr denoting the trace in spinor and isospin space. In order to evaluate these traces, it is instructive to consider the Fourier representation of the homogeneous propagator whereν n = (ν n − iµ), ν n = 2π(n − 1 2 )/β are the fermionic Matsubara frequencies, and M and c ⋆ are defined in Eq. (14). Using the Fourier representation, Eq. (B6) evaluates to which we can identify to be proportional to the homogeneous gap equations where ⃗ ϕ ′ is a homogeneous field configuration that is a solution of the gap equations. Therefore, S eff vanishes when the homogeneous field configurations used as an expansion point are solutions of the gap equations. If these expansion points are used, the second order term is the first non-zero correction. We evaluate the trace in Eq. (B7) to where q = |q| and we used the Fourier representation of the spatial inhomogeneous perturbations The matrix in field space Γ F ϕj ϕ k can be cast into the form where we used that tr(γ i c j c k c l ) = 0, c 2 j = ±I and that the anti-commutator {γ i , c k } evaluates to 0 or 2c k γ i for all considered c j . Thus, we obtain for the second order correction of the effective action In order to make statements about the stability of a homogeneous field configuration one has to determine a basis φ j ( ⃗ ϕ), j ∈ J for which δ j,k λ −1 j + Γ F ϕj ϕ k (q 2 ) is diagonalized. This is not possible in general and depends on the present chemical potentials and the interactions of the model. Furthermore, we assume that all λ j are either λ or 0 according to the subset J of the field content that we consider. If this diagonalization is then indeed possible as it is the case for all channels considered in Eq. (1) at finite baryon chemical potential, one obtains the form where a ′ φj is a coefficient that is determined by the considered field φ j . In this diagonalized form, we identify Γ (2) φj (q 2 ) as the curvature of the effective action for an inhomogeneous perturbation in field direction φ j with momentum q. By writing the denominator of the integrand in Eq. (B18) in a partial fraction, we can split the integral and obtain the final form of the bosonic two-point function as given in Eq. (13), where a φj = 2(a ′ φj − 1). Note that the integral ℓ 1 is also obtained in the fermionic trace in the gap equation (B11) forφ 1 =σ (compare also the GN model gap equation in Section III of Ref. [117]). The gap equation and, correspondingly, ℓ 1 is typically used as a renormalization condition in the vacuum for the coupling constant λ j = λ.

The momentum independent part L1
We consider the integral where E 2 = p 2 + M 2 . The factor of 8 comes from the traces over the four-dimensional spinor and the two-dimensional isospin space. Performing the sum over n, we obtain the standard result where n(x) is the fermionic distribution function The vacuum part is UV-divergent, which we can regulate with a spatial momentum cutoff Λ and use to set the value of λ k via the gap-equations Eq. (B11) with the mass M 0 corresponding to the minimum of the effective action in the vacuum (defined in Section III). We are interested in the contribution of ℓ 1 to the two-point function Γ (2) , where ℓ 1 appears exclusively as L 1 = 1 λ k − ℓ 1 . Using Eqs. (B22) and (B23) we find where the medium integral over the fermionic-distribution function can be evaluated to For T = 0, this evaluates to from which we can naively take the limits µ → 0 and/or M → 0.

The momentum dependent part L2,±
In order to calculate the momentum dependent part of the two-point function, we start by carrying out the Matsubara summation in ℓ 2 and obtain ℓ 2 (M 2 , µ, T, q 2 ) = 8 where E q = M 2 + (p + q) 2 . This integral is UV-finite and, thus, we do not have to implement a regularization scheme. However, the integrand has a divergence at 2p · q = −q 2 that has to be treated with a Cauchy principal value prescription. The vacuum contribution can be calculated analytically and we obtain ℓ 2 (M 2 , µ, T, q 2 ) = 8 4πq arctan q 2|M | − q/2 0 dp p E n(β(E + µ)) + n(β(E − µ)) At any finite T , the medium contribution has to be calculated numerically. However, taking the limit T → 0 enables us to also calculate the medium contribution analytically and we find Eq. (16). From this expression, we can take either the limit q → 0 to obtain or the limit |M | → 0 to obtain While ℓ 2 is not defined for M = q = T = 0 for some values of µ, the whole momentum-dependent contribution to the two-point functions L 2,± is defined with

Generalization to Yukawa models
It is straightforward to generalize the stability analysis of FF models to Yukawa models, which are defined as in Section II B in Eq. (10). Thus, we outline only the meaningful differences with respect to the discussion of FF models in order to preserve brevity.
After an expansion of S Y in powers of inhomogeneous perturbations of the fields, one identifies again S where the second line and third line contain the additional terms resulting from the extension to Yukawa models. Note that the third line contains non-diagonal contributions from the self-interaction terms of the ⃗ χ fields. However, these are proportional toχ jχk , and, thus, vanish when either ofχ j ,χ k can be rotated to zero through a symmetry transformation, as is the case for all Yukawa-type extensions of the model Eq. (9), where only a baryon chemical potential is present. If this is not the case, one needs to take into account this off-diagonal contribution in addition to the, in principle, off-diagonal fermionic contribution. This is the case for the models discussed in Table III with additional chemical potentials. A more involved diagonalization needs to be performed, although the Yukawa contribution is not expected to change the general momentum-dependence of the two-point function, as the term is q-independent. An example for such an analysis is presented in Appendix C.
In the case of only a baryon chemical potential, we utilize symmetry transformations to obtain a homogeneous expansion point ⃗ χ = ⃗ χ such that the off-diagonal contributions in the third row of Eq. (B31) vanish. Performing similar steps as between Eq. (B15) and Eq. (13) leads to the second order correction Eq. (17) and the corresponding bosonic two-point function Eq. (18).
Appendix C: Stability analysis for a Yukawa model with more than one chemical potential In this section, we will show how the off-diagonal contribution from the Yukawa self-interactions (compare the third line of Eq. (B31)) makes the diagonalization of S (2) Y more involved. However, we will also demonstrate that this q-independent contribution does not alter the predictions coming out of the analysis.
The model that we will study is defined as the Yukawa model extension according to Eq. (10) of the FF model in the first row of Table III, i.e., it contains the σ and η 45 fields as well as a baryon chemical potential µ and a chiral chemical potential µ 45 . As documented in in Table III and [120], the FF part of the model is diagonalized by the field basis proportional to (σ ± η 45 ) . (C1) In analogy to this FF model, we study the effective action S eff [χ L , χ R ] N = − Tr ln / ∂+ γ 3 (P L µ L + P R µ R ) + hP L χ L + hP R χ R ] + (C2) where χ L/R are fields of canonical dimension and proportional to the dynamical scalar fields χ σ and χ η45 as We, again, introduced a Yukawa coupling h as well as couplings κ n for the self-interactions and define projectors and chemical potentials accordingly All terms except for the last row in Eq. (C2) either contain only 14 χ L or only χ R . Thus, the second order correction is given by |δχ j (q)| 2 h 2 Γ F χj + q 2 + n>1 κ n n 2(M 2 ) n−1 + 4(n − 1)χ 2 j (M 2 ) n−2 + (C6) + δχ L (−q)δχ R (q)4κ n n(n − 1)χ LχR (M 2 ) n−2 + L ↔ R , where M 2 =χ 2 σ +χ 2 η45 =χ 2 L +χ 2 R and Γ F χj = 1 λ − ℓ 1 + L 2,+ (h 2χ2 j , µ j , T, q 2 ) (C7) is the contribution, that also appears in the corresponding FF model (see Table III and Ref. [120]). The integrals ℓ 1 and L 2,+ are defined in Eqs. (13) and (15). The last row of Eq. (C6) contains the off-diagonal contribution of the self-interaction. This contribution is not dependent on the spatial momentum q of the perturbation, but it makes the diagonalization more complicated (compare Eq. (B31) for the form of this contribution in the more general case). In fact, we are only able to diagonalize this symbolically using Matlab [136]. Using the definitions Y j = q 2 + n>1 κ n n 2(M 2 ) n−1 + 4(n − 1)χ 2 j (M 2 ) n−2 , with j = L, R as well as 4κ n n(n − 1)χ LχR (M 2 ) n−2 (C9) and we write the diagonalization of S (2) eff in the (δχ L (q), δχ R (q))-space where B(q 2 ,χ L ,χ R , µ L , µ R , T ) is a basis changing matrix determined by Matlab, whose form is not relevant for our analysis. In this form one can determine, whether the diagonal entries of the matrix in Eq. (C11) are non-negative. For the physically relevant homogeneous expansion point M both entries are non-negative for q = 0, since otherwise the expansion point would not be a minimum when only considering homogeneous field values. Therefore, in order to prove positivity for all q = |q| it suffices again to show that the entries are monotonically increasing functions of q. We take the derivative of the entries with respect to q and require it to be non-negative where L 2,L/R = L 2,+ (h 2χ2 L/R , µ L/R , T, q 2 ) and its derivative with respect to q is non-negative, i.e., d dq L 2,L/R = L 2,L/R ′ ≥ 0, since it is a monotonically increasing function of q (compare Section III). We can rearrange Eq. (C12) and square it to obtain where obviously 0 ≤ c 2 ≤ 1 and, thus, the inequality is fulfilled for all q. Summarizing this lengthy and delicate analysis: We diagonalized the second order corrections (C6) of a Yukawa model with multiple chemical potentials given by Eq. (C2) using computer algebra systems such as Matlab [136]. Analyzing the resulting expression, we find that both eigenvalues of the relevant curvature matrix in the second order corrections are positive, monotonically increasing functions of the momentum squared q 2 of the inhomogeneous perturbation. Thus, we do not observe instabilities of the homogeneous condensates in the Yukawa model given by Eq. (C2). By the same reasoning, a negative wave-function renormalization (proportional to the second derivative of the two eigenvalues with respect to q), i.e., a so-called moat regime is not observed in the model. Similar behavior is expected for the other Yukawa models that correspond to the FF models in Table III according to Section II B and Eq. (10).