Regulator dependence of inhomogeneous phases in the 2+1-dimensional Gross-Neveu model

The phase diagram of the Gross-Neveu model in $2+1$ space-time dimensions at non-zero temperature and chemical potential is studied in the limit of infinitely many flavors, focusing on the possible existence of inhomogeneous phases, where the order parameter $\sigma$ is non-uniform in space. To this end, we analyze the stability of the energetically favored homogeneous configuration $\sigma(\textbf{x}) = \bar\sigma = \textrm{const}$ with respect to small inhomogeneous fluctuations, employing lattice field theory with two different lattice discretizations as well as a continuum approach with Pauli-Villars regularization. Within lattice field theory, we also perform a full minimization of the effective action, allowing for arbitrary 1-dimensional modulations of the order parameter. For all methods special attention is paid to the role of cutoff effects. For one of the two lattice discretizations, no inhomogeneous phase was found. For the other lattice discretization and within the continuum approach with a finite Pauli-Villars cutoff parameter $\Lambda$, we find a region in the phase diagram where an inhomogeneous order parameter is favored. This inhomogeneous region shrinks, however, when $a$ is decreased or $\Lambda$ is increased, and finally diappears for all non-zero temperatures when the cutoff is removed completely. For vanishing temperature, we find hints for a degeneracy of homogeneous and inhomogeneous solutions, in agreement with earlier findings.


Introduction
Mapping the phase diagram of Quantum Chromodynamics (QCD) at non-zero temperature T and quark chemical potential µ is one of the major challenges in the field of strong-interaction physics [1,2]. From the theoretical side, the situation is complicated by the fact that perturbative techniques are not applicable in the regime of interest, so that non-perturbative methods must be applied. At µ = 0, precise and reliable results from lattice gauge theory with realistic quark masses are available, which have revealed that chiral symmetry, which is spontaneously broken in the vacuum, gets approximately restored in a crossover transition around T ≈ 156 MeV [3][4][5][6]. While the range of validity of present lattice QCD simulations is restricted to chemical potentials µ T (corresponding to µ B 3T with the baryon chemical potential µ B = 3µ), continuum approaches to QCD, i.e., Dyson-Schwinger equations [7] or the Functional Renormalization Group [8], predict that at higher chemical potentials there is a first-order phase transition, ending at a chiral critical endpoint (CEP) at µ CEP ≈ (1.4 . . . 2.0) T CEP . Qualitatively similar results were also obtained quite some time ago within QCD inspired models, like the Nambu-Jona-Lasinio (NJL) model [9] or the Quark-Meson (QM) model [10,11]. In the chiral limit, i.e., at vanishing bare quark mass, the crossover gets replaced by a second-order phase transition, which is joined to the first-order one at a tricritical point (TCP).
Most of these investigations have been performed assuming that the order parameter (the quark condensate ψ ψ ) is homogeneous. Several model studies have revealed, however, that, at least in the mean-field approximation, there are certain regions in the phase diagram where spatially varying condensates are favored over homogeneous ones (see Ref. [12] for a review). Such inhomogeneous phases have been analyzed in detail in models with 1 + 1 spacetime dimensions, like the 1 + 1-dimensional Gross-Neveu (GN) model [13][14][15], but are also found in 3 + 1 dimensions in the NJL and QM models [16][17][18][19] as well as in a Dyson-Schwinger approach to QCD [20]. Typically, the inhomogeneous regions cover the first-order boundary between the homogeneous phases and reach up to the CEP. In the chiral limit the TCP is then replaced by a Lifshitz point (LP), where three phases, the homogeneous symmetry-broken, the restored and the inhomogeneous phase, meet [21].
When considering inhomogeneous phases, the difficulty arises that the determination of the ground state corresponds to the functional minimization of the effective action with respect to the condensate ψ ψ (x) with arbitrary spatial shape. Obviously this is a very hard problem, which, until now, has only been solved in 1 + 1-dimensional models [13][14][15]22] but not in higher dimensions. Instead of a full minimization, most authors, therefore, use certain ansatz functions for the condensate (e.g., single plane waves [16] or embedding the known solutions from the 1 + 1-dimensional into 3 + 1 spacetime dimensions [17]) or perform Ginzburg-Landau or stability analyses [21,[23][24][25][26]. An interesting alternative is to investigate inhomogeneous phases with lattice field theory or related numerical methods, where, at least in principle, the effective action can be minimized without restricting the condensate to a specific ansatz. This has been demonstrated in Refs. [27,28], where the 1 + 1-dimensional GN model with an infinite number of fermion flavors N f was investigated numerically, reproducing the known analytical results for the phase diagram and the condensate functions within numerical precision. Recently lattice field theory was also used to simulate the 1 + 1-dimensional GN model at finite N f , showing for the first time that an inhomogeneous phase also exists in this case [29,30].
In the present article we extend these studies to investigate the GN model in 2 + 1 dimensions.
In doing so, we restrict our investigations to the limit N f → ∞, corresponding to a mean-field approximation. For homogeneous phases this model has been analyzed more than three decades ago [31,32]. It was shown that it is renormalizable in the N f expansion and it was found that the phase transition is of second order at any non-zero temperature, while at T = 0 the order parameter changes discontinuously at the phase boundary. Two of us presented a first lattice field theory study of inhomogeneous phases within that model at a recent conference [33], showing indications for the existence of such phases. Shortly afterwards, however, it was found that the inhomogeneous phase disappears in the continuum limit [34], thus, suggesting that the result of Ref. [33] was an artifact of the lattice discretization. In the present paper we investigate this more systematically, comparing the dependence of the results on the lattice spacing for two different discretization prescriptions. We also complement the lattice field theory calculations by a stability analysis in a continuum approach. In particular, guided by the idea that a non-zero lattice spacing might be similar to a finite momentum cutoff, we investigate the dependence of the inhomogeneous phase on the cutoff.
The paper is organized as follows. In Sec. 2 we introduce the GN model in 2 + 1 spacetime dimensions, discuss the symmetries of the model and address the issue of fermion representations (2-component versus 4-component spinors). In Sec. 3 we explain our continuum approach, in Sec. 4 the lattice field theory techniques. Our results are presented and discussed in Sec. 5. We draw our conclusions in Sec. 6. The Gross-Neveu (GN) model [35] is a relativistic quantum field theory describing N f fermion flavors with a four-fermion interaction. We consider the model in 2 + 1-dimensional Euclidean spacetime, where action and partition function are ψ = ψ 1 , . . . , ψ N f represents N f massless fermion fields, λ is the coupling constant and µ is the chemical potential. Representations of the γ matrices are discussed in Sec. 2.2. At non-vanishing temperature the spacetime integral is over [0, β] × V , where β = 1/T is the inverse temperature and V denotes the 2-dimensional spatial volume.
To get rid of the four-fermion interaction, we introduce an auxiliary scalar field σ and perform a Hubbard-Stratonovich transformation [36], where is the Dirac operator. One can show that the expectation value of the scalar field σ is proportional to the condensate ψ ψ , i.e., Thus, σ can be used as an order parameter for chiral symmetry breaking (for more details see Sec. 2.2).
After integrating over the fermion fields, one obtains an effective action, which only depends on the scalar field σ, Here Det denotes a functional determinant in spacetime and in spinor space, while the degenerate flavor degrees of freedom have been factored out. Strictly speaking, the dimensionful operator Q should be measured with respect to some energy scale, e.g., the temperature in thermal field theory or the inverse lattice spacing in a lattice field theory. While this scale is needed to have a dimensionless argument of the logarithm, it only leads to a constant (i.e., σ-independent) shift of S eff and does not affect our results.
In this work we restrict the dependence of σ to the spatial coordinates, i.e., σ = σ(x 1 , x 2 ). With this restriction S eff is real, which is shown in App. A.2.
As one can see from Eq. (5), the action is proportional to the number of fermion flavors N f . Since we exclusively consider the limit N f → ∞, only field configurations σ corresponding to global minima of S eff contribute to the partition function Z. Thus, instead of integrating over the scalar field σ in Eq. (5) it is sufficient to find a global minimum of S eff . Observables O(σ) are then evaluated on the minimizing field σ = σ , i.e., O(σ) = O( σ ).

Fermion representations and the discrete symmetry σ → −σ
In App. A.1 we show that the effective action (5) has a discrete symmetry i.e., S eff [σ] = S eff [−σ]. A non-vanishing σ indicates spontaneous breaking of this symmetry. Moreover, as discussed in the context of Eq. (4), σ is proportional to the condensate ψ ψ .
We also note that a suitable set of γ matrices has to fulfill the Dirac algebra in Euclidean spacetime, where 1 is the identity matrix in spinor space.

The GN model in 1 + 1 dimensions
We start with a brief discussion of the fermion representation typically used for the GN model in 1 + 1 spacetime dimensions, where the situation is less complicated than in 2 + 1 spacetime dimensions. A possible irreducible 2 × 2 representation of the Dirac algebra (7) is where τ j denote Pauli matrices. A suitable γ 5 matrix, which anti-commutes with both γ 0 and γ 1 , i.e., fulfilling {γ 5 , γ µ } = 0, can be defined according to The free fermion action is then invariant under continuous chiral transformations generated by where λ a are the generators of the U(N f ) flavor symmetry, e.g., the generalized Gell-Mann matrices and the identity, and θ a are the parameters of the transformation. The four-fermion interaction term is, however, not invariant under this chiral transformation. For example for ψ → e iθγ 5 ψ, i.e., (ψψ) 2 is invariant only for θ = nπ/2, n ∈ Z. One can show that only a discrete chiral symmetry remains, Due to Eq. (4), a non-vanishing σ implies a non-vanishing fermion condensate ψ ψ . This in turn indicates spontaneous breaking of the discrete chiral symmetry, because this symmetry implies ψ ψ → − ψ ψ .
Consequently, in 1 + 1 dimensions with fermion representation (8) σ is an order parameter for spontaneous breaking of the discrete chiral symmetry Eq. (12).

The GN model in 2 + 1 dimensions
In 2 + 1 dimensions there are two inequivalent irreducible 2 × 2 representations of the Dirac algebra (7), which can be written as Neither for the representation (14) nor the representation (15) there is an appropriate γ 5 matrix, which anticommutes with all three γ µ . Consequently, there is no discrete chiral symmetry (12) and a non-vanishing σ cannot be interpreted as indication for chiral symmetry breaking. There is, however, another discrete symmetry, changing the sign of ψ ψ , i.e., Eq. (16) implies This symmetry, which is the reflection of the x 2 coordinate, is usually referred to as parity P in 2 + 1 dimensions, since the reflection of both spatial coordinates amounts to a rotation by the angle π. Thus, a non-vanishing σ indicates spontaneous breaking of parity.
Since the GN model is often used as a toy model for chiral symmetry breaking in QCD, also reducible fermion representations with a corresponding γ 5 matrix play an important role. One possibility is to combine the representations (14) and (15) to a 4 × 4 representation, (see e.g., Refs. [37][38][39]). The three matrices are block-diagonal with the upper block corresponding to representation (14) and the lower block to representation (15). There are two linearly independent matrices that anti-commute with the three γ 0 , γ 1 and γ 2 , i.e., both fulfill the necessary properties for a suitable γ 5 matrix. 1 Chiral transformations are defined by taking both γ 4 and γ 5 into account where φ a and θ a are the parameters of the transformation. While free massless fermions are chirally symmetric according to Eq. (20), the four-fermion interaction of the GN model reduces this continuous symmetry to ψ → γ 4 ψ ,ψ → −ψγ 4 (21) or, equivalently, ψ → γ 5 ψ ,ψ → −ψγ 5 .
Note, however, that Eq. (21) and Eq. (22) are not independent, but related by a vector transformation in flavor space (see Ref. [40] for details). Thus, in the following, it is sufficient to consider the discrete chiral symmetry (22). σ is an order parameter for spontaneous breaking of this symmetry, as in the 1 + 1-dimensional case discussed in Sec. 2.2.1.

Equivalence of 2-and 4-component fermion representations
In the following we will show a simple relation between expectation values O(σ) obtained with either of the two irreducible 2-component fermion representations (14) and (15) and the 4-component fermion representation (18).
Note in particular that the phase diagram with respect to the order parameter σ is the same for all three representations. In practice this is useful, because all numerical computations can be performed with the computationally cheaper 2 × 2 fermion representation (14) (or (15)), while the corresponding results are also valid for the 4 × 4 fermion representation (18), where an interpretation in terms of chiral symmetry and its spontaneous breaking is possible. In the following we denote the dimension of the fermion representation by N d .

Continuum approach
In this section, before introducing our lattice field theory techniques in Sec. 4, we describe, how we study the phase diagram of the 2 + 1-dimensional GN model using continuum methods. As explained in Sec. 2, the ground state of the system in the limit N f → ∞ corresponds to the field configuration σ(x), x = (x 0 , x 1 , x 2 ), which minimizes the effective action S eff , given in Eq. (5). Despite being a considerable simplification compared to the situation at finite N f , finding this configuration in 2 + 1 spacetime dimensions is still extremely difficult. Since the ground state is static, we can safely ignore the (imaginary) time coordinate x 0 , but we must retain the possible dependence of the σ field on the spatial coordinates x = (x 1 , x 2 ), if we want to find the true minimum. The minimization of S eff is thus a functional minimization with respect to σ(x). This is a very hard problem, both analytically and numerically, which has not been solved so far.
While in Sec. 5.3 we will present first steps towards such a full minimization using lattice field theory, here we restrict ourselves to a simpler problem and perform a stability analysis of the lowest homogeneous state with respect to small inhomogeneous fluctuations. This method, which has already been applied in 3 + 1 dimensions to investigate inhomogeneous phases in the NJL model [16,24,25] and the QM model [26], corresponds to searching for a sufficient, although not necessary, condition for an inhomogeneous phase: If the lowest homogeneous state turns out to be unstable against small inhomogeneous fluctuations, it is clear that the ground state must be inhomogeneous. On the other hand, the lowest homogeneous state can be stable against small inhomogeneous fluctuations but still be unstable against large ones. The inhomogeneous phases found by a full minimization of the effective action could thus be larger than the unstable regions of the stability analysis, but the latter are always a part of the former. In the following we explain this method in more detail.

Stability analysis
In the first step we minimize the effective action at given chemical potential and temperature with respect to spatially constant fields, σ =σ. We will give a few technical details in Sec. 3.2, but basically this is a standard procedure, which is obviously much simpler than the functional minimization with respect to arbitrary space-dependent fields σ = σ(x). In the second step, we consider small fluctuations around these homogeneous solutions and inspect their effect on the effective action. To this end we write whereσ is a homogeneous field and δσ(x) denotes a fluctuation, which can have an arbitrary spatial shape but is assumed to have an infinitesimally small amplitude.
Decomposing the Dirac operator in the same way, Q =Q + δσ withQ = Q(σ), and noting that ln(Det(Q)) = Tr(ln(Q)), this term can straightforwardly be expanded in powers of δσ, ln(Det(Q)) = Tr(ln(Q)) + Tr ln 1 +Q −1 δσ = Tr(ln(Q)) − Accordingly, the effective action can be expanded as where S (n) eff corresponds to the contribution of the n-th order in the fluctuations. Specifically we find for the three lowest-order terms where the integrals are over the spatial coordinates x 1 and x 2 , while the factors of β originate from the x 0 integrations. Demanding thatσ minimizes S We note thatQ −1 = (γ ν ∂ ν + γ 0 µ +σ) −1 is just the propagator of a non-interacting fermion with massσ at chemical potential µ. In coordinate space it thus depends on (the difference of) two space-time variables. The functional traces above are defined as where tr denotes a trace in spinor space. These expressions are most easily evaluated using the Fourier representation, or, in the infinite-volume limit, where ω n = 2π(n − 1/2)/β are fermionic Matsubara frequencies, and is the Euclidean propagator in momentum space.
Inserting this into Eq. (35), the gap equation becomes where we have defined After performing the Matsubara sum, this takes the form with E p = p 2 +σ 2 and the Fermi functions The fluctuating field can be Fourier transformed as where the field in momentum space must obey the relation δσ(−q) = δσ * (q) to ensure that δσ(x) is a real function. For the first-order contribution to the effective action we then obtain from Eq. (33) where the second equality follows from the gap equation (40). In fact, since according to the first equality only the homogeneous (q = 0) fluctuations contribute to S (1) eff , it has to vanish if we expand about the "homogeneous ground state", i.e., the lowest homogeneous solution.
To find instabilities we therefore have to investigate the second-order contribution to the effective action. Evaluating Eq. (34) in momentum space we obtain where Carrying out the Matsubara sum, one obtains From Eq. (46) we can see that, unlike in the first-order contribution, also inhomogeneous (q = 0) fluctuations contribute to S eff . In particular, if Γ −1 (q 2 ) < 0 in some momentum region, small fluctuations in that region will lower the effective action with respect to the homogeneous ground state. A sufficient condition for an instability of the homogeneous ground state with respect to developing inhomogeneities is therefore to find Γ −1 (q 2 ) < 0 in some interval around any momentum q = 0. A second-order phase boundary between a homogeneous and an inhomogeneous phase is thus given by the values of T and µ for which Γ −1 (q) just touches the zero axis, i.e., both Γ −1 = 0 and dΓ −1 /dq 2 = 0 at some non-vanishing momentum.

Homogeneous phase diagram, tricritical points and Lifshitz points
We stress again that, in order to arrive at the above conclusions, we must expand about the homogeneous ground state, i.e., we first have to determine the constant fieldσ which minimizes S (0) eff . To this end we first solve the gap equation (40) numerically, which yields all stationary points of S (0) eff . These can be minima, maxima or saddle points, but since there is only a (small) finite number of solutions, it is easy to identify the absolute minima. Strictly speaking, since the effective action is symmetric under σ → −σ (see Eq. (6)), there is a degeneracy between S (0) eff (σ) and S (0) eff (−σ). However, this is of no further relevance and we may pick either of these solutions, e.g., the positive one. 2 Another consequence of this symmetry is the fact thatσ = 0 is always a solution of the gap equation.
Minimizing S (0) eff at a large number of points in the µ-T plane we get the homogeneous phase diagram as a by-product. The phase boundaries between the homogeneous symmetry-broken phase and the restored phase are then simply given by the lines where the minimum changes from a non-vanishing to a vanishing value ofσ and can in general be obtained by a bisection procedure to the desired accuracy. For first-order phase transitions, i.e., whenσ discontinuously drops to zero this is indeed what we do in our numerical calculations The method is less efficient for second-order phase transitions, i.e., whenσ continuously goes to zero. In this case we determine the phase boundary by a stability analysis of the restored phase against homogeneous fluctuations of the sigma field. The idea is essentially the same as in Sec. 3.1 but now we always expand around the trivial solutionσ = 0 and restrict ourselves to spatially constant fluctuations, i.e., the q = 0 mode. Defining the phase boundary is then given by the condition provided the system is stable against finite fluctuations, related to a first-order phase transition. The latter can in general be identified as described above. An exception are the regions in the close vicinity of a TCP, where a first-order phase boundary turns into a second-order one, so that their distinction becomes difficult in practice. Therefore, in order to identify these points, we perform a Ginzburg-Landau analysis. To this end we expand S (0) eff around the restored solution σ = 0 in powers ofσ, Odd powers vanish because of the symmetry (6). Assuming that the fourth-and all higherorder derivatives are positive, the phase transition is of second order and takes place at the points where the second derivative vanishes. On the other hand, if the fourth-order derivative is negative, there can be a first-order phase transition. TCPs are thus located at the points where both the second and the fourth derivative of S confirming our previous result (51) for the second-order phase boundary. For the fourth-order derivative we find i.e., the additional condition for the TCP is 2 (0)|σ =0 = 0.
Finally, let us include again the possibility of inhomogeneous phases into these considerations. Obviously, Eq. (51) only describes the location of a second-order phase boundary between a homogeneous symmetry-broken and a restored phase if the corresponding region is stable against inhomogeneities. A necessary condition for this is that the momentum dependent term in Eq. (50) does not turn Γ −1 0 negative at q 2 > 0. The LP, i.e., the point where the second-order phase boundary between the two homogeneous phases is smoothly connected to the second-order instability line with respect to inhomogeneous fluctuations is thus given by Eq. (51), together with the condition and, hence, the LP coincides with the TCP (cf. Eq. (54)). This is a well-known result from the 1+1-dimensional GN model [14] and the 3+1-dimensional NJL model [21]. Indeed, the number of spatial dimensions did not enter in an essential way into the considerations above.
It is known from the homogeneous analysis of the 2+1-dimensional GN model that the phase transition is second-order at any non-zero temperature, while at T = 0 the order parameter changes discontinuously at the phase boundary [31,32,41]. Since there is thus no TCP, we conclude that there is no LP either, and we may anticipate that we will not find an inhomogeneous phase in the model at non-zero temperature. We stress, however, that these arguments do not exclude inhomogeneous phases which are reached via first-order phase transitions or which have second-order boundaries without LP (e.g., inhomogeneous "islands" surrounded by a single homogeneous phase).

Regularization and renormalization
The integral 1 defined in Eq. (42) diverges linearly and needs to be regularized. It is well known, however, that the model is renormalizable in the large-N f expansion [32]. The renormalization is done by employing the gap equation (40) to relate the coupling constant λ to a cutoff Λ, where the scale is set by the ground-state value σ 0 ofσ at T = µ = 0, Since the model is renormalizable, all observables remain finite, even in the limit Λ → ∞. In the following sections we refer to that limit as renormalized GN model.
It is, however, instructive to study the phase diagram of the model not only in this renormalized limit but also for finite values of the cutoff parameter, in particular in view of our lattice field theory computations discussed in Secs. 4 and 5, which have to be carried out always at finite lattice spacing. Specifically, we use Pauli-Villars (PV) regularization with one regulator term, i.e., we subtract from the integrand a term with the same structure but withσ 2 replaced bȳ σ 2 + Λ 2 . We apply this prescription only to the vacuum parts, i.e., the parts which survive at T = µ = 0, but not to the medium contributions, which are related to Fermi occupation numbers in the integrand and are finite without regularization. 4 For the integral 1 we then obtain from Eq. (42) and thus, from Eq. (56), for the coupling constant For the integral 2 we get where q = |q| = q 2 . Note that 2 is finite even for Λ → ∞ but we regularize it nevertheless for consistency. In particular one can show that the unregularized expressions satisfy the relation which underlies the calculation leading to Eq. (54) and, related to this, the coincidence of a TCP and a LP. This relation is only preserved if 1 and 2 are regularized in the same way. One finds Note that this expression remains also finite in the restored phase, whereσ = 0, A more intuitive regularization scheme is to apply a sharp cutoff to the momentum integrals, Compared to PV regularization, this prescription has the disadvantage that it violates Lorentz invariance. Related to this, the result for 2 is in general not unique but changes under a shift of the loop momentum p → p + ∆p in Eq. (49). We therefore stay with the PV regularization scheme as described above. It turns out however, that for both 1 and 2 (0), which are relevant for the homogeneous second-order phase boundary as well as for the LP (see Eqs. (51) and (55)), the cutoff regularization gives identical results. 5 In this sense the PV regularization parameter Λ can be interpreted as a Lorentz invariant generalization of a momentum cutoff.
Finally, we note that one regulator term is not sufficient to render the cubically divergent effective action finite. The physically more relevant difference ∆S is however finite in this regularization scheme. In particular the determination of the homogeneous ground state as the basis of our stability analysis can always be done by calculating ∆S where the Li are polylogarithms. Inserting Eq. (59) for λ and evaluating this expression at T = µ = 0, we obtain for the vacuum solutionσ = σ 0 which is always negative, confirming the stability of the vacuum solution σ 0 with respect to the restored phase σ = 0. For Λ → ∞, the right-hand side of this equation stays finite and becomes −(N d /24π)σ 3 0 .

Analytical results for the phase diagram
We are now in the position to evaluate the relations derived in Sec. 3.2 more explicitly. In particular we can determine the second-order phase boundary between the restored and the homogeneous symmetry-broken phase by inserting Eq. (58) withσ = 0 and Eq. (59) into Eq. (51). This leads to where we have defined Here and in the following we assume the vacuum condensate σ 0 to be positive.
Solving Eq. (67) for µ = 0 yields the critical temperature which can be expanded in powers of the inverse cutoff, This is consistent with the known result T c /σ 0 = 1/2 ln 2 in the limit Λ → ∞ [31,32].
For T ≤ T c we can also solve Eq. (67) for the critical chemical potential. One finds Again this is consistent with the critical phase boundary found in Refs. [31,32] in the limit Λ → ∞, corresponding to s → σ 0 . In this case the line reaches the µ axis at µ c (T = 0) = σ 0 , which is replaced by s in the case of finite Λ.
We must keep in mind, however, that the above equations have been derived under the assumption that there is a second-order phase transition from the restored to the homogeneous symmetry-broken phase. In Ref. [32] it was shown for Λ → ∞ that this is indeed the case for all T > 0 (at least for homogeneous condensates), while at T = 0 the condensate changes discontinuously at µ c , corresponding to a first-order phase transition. For arbitrary values of Λ we can compute the corresponding TCP (which, as shown above, is equal to the LP of an inhomogeneous phase) as the simultaneous solution of Eqs. (51) and (55), i.e., as the point where the lines defined by Eq. (71) and by the condition 2 (0) = 0 cross. Using Eq. (62) and taking the limitσ → 0 one finds that the temperature at the LP is given by the equation which has the solution where W −1 is the lower branch of the Lambert W function. Inserting Eq. (72) into Eq. (71) one then obtains for the chemical potential The line of LPs in the µ-T plane for varying cutoff parameter Λ is shown in Fig. 1, where the dots indicate specific examples of Λ. For finite values of Λ there is a LP at non-zero temperature, signaling the existence of an inhomogeneous phase. 6 In the limit Λ → ∞, on the other hand, T LP → 0, and we expect a second-order phase transition between homogeneous phases which are stable against inhomogeneous fluctuations for all T > 0. The figure shows, however, that this limit is reached extremely slowly. denoted by a and the number of lattice sites is N 2 s , i.e., N s lattice sites in each of the two directions and L = N s a. Since we are interested in studying spontaneous chiral symmetry breaking, it is essential to use a chirally symmetric fermion discretization. We decided to use the naive discretization (see, e.g., the textbook [42]). Naive fermions imply fermion doubling, i.e., in our case of two spatial dimensions the number of fermion flavors N f is restricted to multiples of 4. For our work this is not a problem, because we are interested in the limit N f → ∞.
The extent of the temporal direction corresponds to the inverse temperature β = 1/T and boundary conditions are antiperiodic. In temporal direction we do not use lattice field theory, but regularize the fermion fields by a superposition of 2N 0 plane waves as discussed in detail below and in Refs. [19,28]. Since the chiral condensate does not depend on x 0 , i.e., σ = σ(x) as discussed in previous sections, plane waves allow straightforward analytical simplifications of the fermion determinant. Moreover, the chemical potential can be introduced as in the continuum by adding γ 0 µ to the Dirac operator. In particular an exponential coupling as typically used in lattice field theory is not necessary. As a consequence we expect smaller discretization errors (see Ref. [30] for a detailed discussion).

Free fermions
We define the plane-wave expansion of a fermion field representing a single flavor as where 2N 0 represents the number of modes used in temporal direction. The frequencies ω n 0 = 2π(n 0 − 1/2)/β with n 0 = −N 0 + 1, −N 0 + 2, . . . , +N 0 − 1, +N 0 imply antiperiodic boundary conditions in temporal direction. We use 1/a ≡ 1 as density of degrees of freedom 7 , i.e., 2N 0 /β = 1. Consequently, the inverse temperature and the number of modes are related according to β = 2N 0 . Inserting the plane-wave expansion into the free fermion action leads to For the two spatial directions we apply the naive lattice discretization, Note that, due to fermion doubling, χ represents four fermion flavors and, thus, cannot be interpreted as a standard fermion field ψ as, e.g., used in Eq. (76) or in Sec. 2. The relation between the components of χ and of ψ is non-trivial. N s must be even to have periodic boundary conditions in the spatial directions for all flavors. We refer to Refs. [29,43], where this is discussed in detail.
In the following we explore two suitable choices, which have the same correct continuum limit: where the corresponding Fourier transform isW 1 (k) = (cos(k) + 1)/2, and where N s is a multiple of 4 and the corresponding Fourier transform isW 1 (k) = Θ(π/2 − |k|).
Because of fermion doubling N f /4 naive fermion fields χ f are needed to represent N f fermion flavors. We stress that a specific non-diagonal structure of W 2 (x) is mandatory for a valid discretization of the GN model with naive fermions, i.e., a discretization with the correct continuum limit. The straightforward and probably more intuitive choice W 2 (x) = δ x 1 ,0 δ x 2 ,0 , which we used at an early stage of this work [33] and which was also used in Ref. [34], introduces additional four fermion interactions which are not part of the GN model, e.g., couplings between different flavors likeψ 1 ψ 2 σ,ψ 1 ψ 3 σ, etc. For a more detailed discussion we refer to Sec. 5.3.1 of this work and to appendix A of Ref. [29].
After integrating over the fermion fields in the partition function, as discussed in Sec. 2, one obtains the discretized effective action The Dirac operator Q is a matrix of size 2N 0 N 2 s N d × 2N 0 N 2 s N d with rows and columns representing spacetime and spin, This matrix is diagonal with respect to the temporal indices n 0 and n 0 . Thus, one can factorize the fermion determinant in Eq. (81) according to i.e., the problem is reduced to the computation of the determinants of 2N 0 smaller matrices of size N 2 s N d × N 2 s N d . When restricting the dependence of σ to only one of the two spatial coordinates, i.e., σ = σ(x 1 ), the numerical costs of computing S eff can be further reduced. Similar to the plane-wave expansion in temporal direction one can diagonalize the Dirac operator (82) with respect to x 2 by writing with x 2 still restricted to the sites of the spatial lattice 8 and k n 2 = 2πn 2 /L. The Dirac operator (82) then becomes Q(n 0 , x 1 , n 2 ; n 0 , x 1 , n 2 ) = = δ n 0 ,n 0 δ n 2 ,n 2 γ 0 (iω n 0 + µ)δ x 1 , where W 1 = W 1 or W 1 = W 1 . This is a diagonal matrix with respect to the temporal indices n 0 and n 0 as well as the spatial indices n 2 and n 2 . The computation of the fermion determinant in Eq. (81) is, thus, reduced to the computation of the determinants of 2N 0 N s matrices of size ln Det Q(n 0 , x 1 , n 2 ; n 0 , x 1 , n 2 ) .

Numerical evaluation of the effective action
We perform all computations with the 2 × 2 representation of γ matrices (14), i.e., N d = 2. An important part of these computations is the numerical evaluation of the effective action (81) for a given field configuration σ. Typically, this has to be repeated many times, e.g., when minimizing S eff with respect to σ, or when checking the stability of a homogeneous condensate σ =σ = const with respect to inhomogeneous perturbations. Particularly time consuming is the computation of ln(Det(Q)). To maximize efficiency, we distinguish the following three cases: The Dirac operator (82) is a block-diagonal matrix with 2N 0 blocks of size 2N 2 s × 2N 2 s . The determinant of each block is computed via a standard LU-decomposition. We use the publicly available GSL library [44].
The Dirac operator (85) is a block-diagonal matrix with 2N 0 N s blocks of size 2N s × 2N s . Again the determinant of each block is computed via a standard LU-decomposition.
• σ =σ = const: For homogeneous condensates we obtain identical results for the discretizations W 2 and W 2 , as can be seen from Eqs. (78) to (80). The eigenvalues of Q can be calculated analytically in a straightforward way. ln(Det(Q)) is then computed by summing over the eigenvalues, leading to (note that in contrast to previous equations the sum over n 0 is restricted to positive integers) with ω n 0 = 2π(n 0 − 1/2)/β and k n j = 2πn j /L.
In Sec. 4.4 we also need the second derivative of S eff with respect toσ atσ = 0, which can be calculated as Note that for σ =σ = const the effective action corresponds to the term S (0) eff of the expansion (31). Hence, Eq. (89) is the lattice field-theory version of Eq. (53). Indeed, since 2N 0 N 2 s = βV in lattice units, this is immediately obvious for the first term on the right-hand side, while the lattice version of the integral (41) is given by in consistency with the second term.

Stability analysis on the lattice
In Sec. 3.1 we have discussed in detail within the continuum approach, how to perform a stability analysis with respect to inhomogeneous perturbations of a homogeneous condensate minimizing the effective action. The same steps and calculations can also be carried out in the lattice field theory approach.

Coupling constant, lattice spacing and scale setting
As usual in a renormalizable lattice field theory, the lattice spacing a can be set by tuning the coupling constant λ, i.e., a = a(λ). 9 Moreover, in our particular regularization the number of modes in temporal direction N 0 is proportional to the inverse temperature β, i.e., N 0 = β/2a, as discussed in Sec. 4.1.1. Thus, the temperature can be adjusted by either changing N 0 or λ.
For our computations we first fix the number of modes at the critical temperature, denoted by N 0,c , where we typically use a small number, 2 ≤ N 0,c ≤ 5 (throughout this section N s = L is chosen sufficiently large, such that finite volume corrections are essentially negligible). This in turn fixes the coupling constant λ and the lattice spacing a, where the former has to be tuned in such a way that 2N 0,c a(λ) = β c . An obvious possibility is to determine λ such thatσ(λ − ) = 0 andσ(λ + ) = 0 for infinitesimal (see Fig. 2, where |σ| is plotted as a function of λ). 10 Mathematically equivalent, but more practical from a numerical point of view is to consider (89)) as a function of λ and to determine its root, which leads to Even though λ is now fixed, the temperature can still be changed in discrete steps by increasing or decreasing the number of modes, since T = 1/2N 0 . In principle, the scale could now be set via the critical temperature T c = 1/2N 0,c , i.e., we could express all dimensionful quantities in units of T c . However, we prefer to set the scale via σ 0 , which is common in the existing literature. To determine σ 0 , we compute |σ| µ=0,T for several small T by minimizing S eff with respect toσ. Using Eq. (87) this is numerically rather simple and can by done by a standard golden section search. |σ| µ=0,T quickly approaches a plateau, when decreasing T , i.e., |σ| µ=0,T is almost constant for T < ∼ T c /4, and the plateau value is identical to σ 0 (see Fig. 3, left plot). In the following we express all dimensionful quantities in units of σ 0 . For example for T c /σ 0 we obtain values rather close to the analytically known result 1/2 ln(2) [31] also at finite lattice spacing. When increasing N 0,c , which amounts to decreasing λ as well as decreasing a, and approaching the continuum limit, our results for T c /σ 0 approach 1/2 ln(2) ≈ 0.721 . . ., as can be seen in the right plot of Fig. 3. Note, however, that in contrast to our continuum approach, where this limit is reached from below when the cutoff Λ is sent to infinity (see Eq. (70)), here the data points first overshoot the limiting value and then approach the latter from above.

Results
After having introduced our different approaches, we finally turn to the discussion of our results for the phase diagram. From Refs. [31][32][33][34]45] as well as from our analytical studies in Sec. 3.4 we expect up to three phases, each characterized by a different behavior of the field σ: • A symmetric phase with σ = 0 at large µ and/or large T .
• A homogeneous symmetry-broken phase with a spatially constant (but µ and T dependent) field σ =σ = 0 at small µ and small T .
• Possibly an inhomogeneous phase, where σ(x) is a varying function of the spatial coordinates, at intermediate µ and small T [33]. This phase might only be present at a finite value of the regulator (e.g., Pauli-Villars cutoff Λ or lattice spacing a), as indicated by recent lattice field theory results reported in Ref. [34] and the behavior of the LP found in Sec. 3.4.
In the following we present and compare results, which we computed using our three regularizations, the continuum approach with Pauli-Villars cutoff Λ discussed in Sec. 3 and the two lattice discretizations with lattice spacing a discussed in Sec. 4. The advantage of the lattice approach compared to the continuum approach is that one can perform minimizations of S eff with respect to σ(x), allowing for arbitrary modulations of the condensate. There are, however, also drawbacks, namely computations are typically expensive and have to be carried out at finite spatial volume V and finite lattice spacing a, while in the continuum approach both V and the Pauli-Villars cutoff Λ can be sent to infinity. Thus, both approaches complement each other.
Our main goal is to explore and to understand the dependence of a possibly existent inhomogeneous phase in the 2 + 1-dimensional GN model on the regulators Λ and a. In particular we expect to get insights concerning the a dependence of the lattice field theory results by analogy to the Λ dependence of the continuum results, for which we have derived analytical expressions in Sec. 3.4. For a crude quantitative comparison of continuum and lattice results we relate the corresponding regulators via This equation can be obtained by equating the considered regions of spatial momenta, which are πΛ 2 (a circle in momentum space in the continuum approach; see Sec. 3.3) and (π/a) 2 (a square in momentum space in the lattice approach; see appendix A in Ref. [29]).
lattice field theory continuum approach . We also list the corresponding coupling constants λ for both the lattice field theory and the continuum approach (see Eq. (94) and Eq. (59)). All quantities are made dimensionless with the help of the vacuum value σ 0 of the scalar field σ.
The lattice spacings used in the numerical calculations presented below, together with the corresponding values of the Pauli-Villars cutoffs, are listed in Tab. 1. We also list the corresponding coupling constants λ for both cases. Despite the rather different approaches, the values of λσ 0 turn out to be quite similar, which may be considered as an additional justification of Eq. (95).

Phase diagram for homogeneous condensate σ(x) =σ = const
In this subsection we only allow spatially constant condensates, i.e., σ(x) =σ = const. The resulting homogeneous phase diagrams are presented in Fig. 4. For Λ → ∞ the phase boundary obtained in the continuum approach is of second order (with the exception of the endpoint at T = 0, where it is of first order) and given by Eq. (71) with s → σ 0 . It is identical to the known result of Ref. [31] and shown as solid blue line in both plots of Fig. 4. The two lattice discretizations with W 2 and W 2 become identical for homogeneous condensates, as discussed in Sec. 4.2. In the right plot of Fig. 4 we show results for two different lattice spacings, aσ 0 = 0.379 and aσ 0 = 0.049, but identical spatial volume (in units of the critical temperature), (LT c ) 2 = (N s /2N 0,c ) 2 = 10 2 . The phase boundary obtained with the finer lattice spacing is much closer to the result from Ref. [31] and a combined continuum and infinite-volume extrapolation indicates consistency with that result. It is interesting to note that at the coarse lattice spacing aσ 0 = 0.379 the phase transition is of second order for T /σ 0 ≥ 0.189, while it is of first order for T /σ 0 ≤ 0.165. At the fine lattice spacing aσ 0 = 0.049 we only observe second-order phase transitions, but we expect that somewhere below T /σ 0 = 0.075 (the smallest temperature we have investigated) these transitions change to first order.
Our continuum results, which are displayed in the left panel of Fig. 4, exhibit a qualitatively similar behavior. We use the finite Pauli-Villars cutoffs Λ/σ 0 = 4.68 and Λ/σ 0 = 36.17, which, according to the matching formula (95), correspond to the two lattice spacings used in the right panel of the figure (see also Tab. 1). As discussed in Sec. 3.4 and illustrated in Fig. 1, at finite Λ there always exists a TCP, separating first-and second-order phase boundaries. In Fig. 4 these points are marked by open circles. For Λ/σ 0 = 4.68 the TCP is located at µ/σ 0 = 0.881 and T /σ 0 = 0.245, which is somewhat lower in chemical potential and higher in temperature than observed in the lattice approach for aσ 0 = 0.379. A similar tendency is seen at Λ/σ 0 = 36.17, where the TCP is found at µ/σ 0 = 0.985 and T /σ 0 = 0.166 in the continuum approach, while no first-order phase transition at temperatures above T /σ 0 = 0.075 was found in the corresponding lattice calculation. We emphasize, however, that Eq. (95) is only a crude prescription to compare two very different approaches and there is no reason to expect perfect agreement. In general we find that the lattice points at low T converge faster to the Λ → ∞ line than the corresponding continuum results, while the second-order phase boundaries at low µ agree fairly well, even quantitatively.

Instabilities with respect to spatially inhomogeneous perturbations
We now relax the restriction to a homogeneous field σ and investigate the possible appearance of inhomogeneous phases in the phase diagram. To do this in a rigorous way, one has to consider arbitrary modulations of the condensate and minimize the effective action S eff with respect to σ(x). In the lattice approach this is possible, but numerically a very hard problem. It amounts to minimizing S eff in Eq. (81), which is a function in N 2 s variables σ(x). Finding the global minimum of a function in a large number of variables (in our case N 2 s = O(10 2 ) . . . O(100 2 )) is an extremely challenging task. Therefore, in a first step, we check, whether the constant σ =σ determined in Sec. 5.1 is stable with respect to spatially inhomogeneous perturbations δσ. We do this both in the continuum approach introduced in Sec. 3 as well as in the lattice approach of Sec. 4. Steps towards a rigorous minimization of the effective action with respect to arbitrarily varying fields will be discussed afterwards, in Sec. 5.3.
The starting point for the stability analysis is S (2) eff at σ =σ, which is given by Eq. (46) in the continuum approach and by Eq. (91) in the lattice approach. In both cases the perturbation δσ is expressed as a sum of plane waves. This leads to a particularly simple form for S (2) eff , where an instability with respect to inhomogeneous perturbations is indicated by a negative value of Γ −1 for at least one momentum q = 0. (For q = 0 we have Γ −1 ≥ 0 because σ =σ minimizes S eff when restricting σ to a constant.) Thus, in order to determine numerically whether there is an instability at given µ and T , one just has to find the minimum of Γ −1 . This is rather straightforward, because Γ −1 is a function with only a single argument q 2 in the continuum approach and a finite set of numbers, each corresponding to one of the discrete lattice momenta q k , in the lattice field theory approach. Instability lines separating regions of stability from regions of instability can then be computed by a simple bisection in µ direction.
In Fig. 5 we show results of such analyses obtained in lattice field theory with the discretization W 2 = W 2 (see Eq. (80)) for three different lattice spacings as well as in the continuum approach with the corresponding values of the Pauli-Villars cutoff Λ. Regions where the homogeneous ground state is unstable against inhomogeneous fluctuations and which are therefore parts of an inhomogeneous phase are shaded in green.
In the continuum approach (left column of Fig. 5) these instability regions are found in a temperature regime from T = 0 up to the LP, which, as shown in Sec. 3.2, coincides with the TCP of the corresponding homogeneous phase diagram. Accordingly, as already expected from the behavior of the LP shown in Fig. 1, the instability region shrinks when Λ is increased, and disappears completely for Λ → ∞. At finite Λ we find that instabilities only occur in the symmetric phase of the homogeneous phase diagrams, while the homogeneous symmetrybroken phase remains stable against small inhomogeneous fluctuations. In particular the "left" phase boundaries of the instability region coincide with the first-order phase boundaries of the homogeneous phase diagrams. It is quite likely, however, that the true inhomogeneous phase (at finite Λ) reaches to somewhat lower values of µ, as known, e.g., from the GN model in 1 + 1 dimensions [13][14][15] or the NJL model in 3 + 1 dimensions [16,17].
For further illustration of the continuum-approach results of Fig. 5 we show in Fig. 6 selected examples of the function Γ −1 for Λ/σ 0 = 4.68 (left plot) and Λ → ∞ (right plot). The green curve in the left plot corresponds to a point inside the instability region, as obvious from the fact that Γ −1 is negative in some momentum interval. The red curve, on the other hand, just touches the Γ −1 = 0 axis at some non-zero q 2 and therefore corresponds to a point on the boundary Figure 5: "Phase diagram" of the 2 + 1-dimensional GN model in the µ-T plane obtained by stability analyses with respect to spatially inhomogeneous perturbations. Lattice field theory results obtained with W 2 for three different lattice spacings a are shown in the right column, continuum results for the corresponding Pauli-Villars cutoff Λ = √ π/a in the left column. The blue line appearing in all six plots represents the known phase boundary from Ref. [31], which we also obtain in our continuum approach for Λ → ∞. Regions of instability, which are part of and may coincide with inhomogeneous phases, are shaded in green.
between the stable and the unstable region with respect to inhomogeneous fluctuations. The blue curve has a similar minimum as the red curve, however at q 2 = 0, thus indicating a LP (cf. Eq. (55)). Finally, the magenta curve corresponds to a point on the second-order boundary between the homogeneous symmetry-broken and the symmetric phase, where Γ −1 (0) = 0 but with a non-vanishing derivative with respect to q 2 .
The magenta curve in the right plot of Fig. 6 is qualitatively similar and corresponds to a point on the second-order boundary between the homogeneous symmetry-broken and the symmetric phase in the renormalized model, i.e., for Λ → ∞. As pointed out above, for Λ → ∞ there is no instability with respect to inhomogeneous fluctuations at any T > 0, and therefore the blue curve is characteristic for the entire phase boundary. An exceptional case is the situation at T = 0 and µ = σ 0 , indicated by the blue curve. Here Γ −1 (q 2 ) = 0 in the whole interval 0 ≤ q 2 ≤ 4σ 2 0 . This point may thus be interpreted as a point on the instability boundary, but instead of having a single unstable mode, the instability is driven by any momenta 0 ≤ |q| ≤ 2µ or superposition of them. This behavior, which clearly deserves further investigation, could be related to an observation reported in Ref. [45], where it was found that at T = 0 a specific inhomogeneous modulation based on Jacobi elliptic functions is degenerate with the favored homogeneous solutions. Going back to Fig. 5, we now turn to the diagrams shown in the right column, which have been obtained within the lattice field theory approach, using the discretization with W 2 = W 2 . As one can see by comparison with the left panels, the results are qualitatively similar to those from the continuum approach. Again, the regions of instability are only found in the symmetric phase of the homogeneous phase diagram, and they shrink when the lattice spacing is decreased. We therefore expect that in the continuum limit a → 0 the instability regions vanish (at least for T > 0), as they do for Λ → ∞ in the continuum approach.
When carefully inspecting the instability line towards the symmetric phase, one can identify small deviations from a smooth behavior. This is a finite-volume effect, similar to that observed in numerical studies of the 1 + 1-dimensional GN model [27,28]. The reason is that close to the phase boundary only non-vanishing momenta inside a small interval would lead to negative Γ −1 , but none of the numerically considered momenta, which are quantized by the finite spatial volume, is inside that interval. We also note that the LP has a higher temperature than the TCP, i.e., they do not coincide as it is the case in the continuum approach. We attribute this to the lattice discretization, which invalidates some of the steps we have performed in Sec. 3.2 to show the coincidence of the two points in the continuum approach. In particular the momentum derivative in Eq. (55) is not well defined on the lattice, where the allowed momenta are discrete. Interestingly, the positions of the LP for the different lattice spacings agree fairly well with the corresponding ones in the continuum approach, when we apply our matching prescription (95), while the agreement of the positions of the TCP is not as good (see Fig. 4).
For aσ 0 = 0.379, despite the similar positions of the LP, the instability region from the lattice approach extends to considerably higher chemical potentials than in the corresponding continuum calculation. In fact, extending the stability analysis to higher values of µ it turns out that the temperature of the boundary to the stable region rises again, leading to a much larger instability region, before it drops towards T = 0 at around µ/σ 0 = 3.5. This can be seen in Fig. 7 (upper right plot), where the results of the stability analysis are shown for the same parameters as in Fig. 5 but for a much larger µ range. In the corresponding continuum plot (upper left) we see that around µ/σ 0 = 5 a second instability region appears, which grows to very high temperatures and does not seem to be bounded at all. Such a type of inhomogeneous region is also known from the 3 + 1-dimensional NJL model, where it has been termed "inhomogeneous continent". Lowering the cutoff below Λ/σ 0 ≈ 3 (not shown in the figure), the continent merges with the finite instability region at smaller µ, which is connected to the homogeneous symmetrybroken phase (sometimes called "inhomogeneous island"), leading to a single large instability region, similar as in the lattice result at aσ 0 = 0.379 but without a decreasing boundary at large µ. On the other hand, if we increase the cutoff, the onset of the continent is pushed to higher values of µ, and it disappears completely for Λ → ∞, so that no instability region is left in the renormalized model. The behavior on the lattice is qualitatively similar, with the exception that the size of the second inhomogeneous region is finite. This difference is most likely due to the fact that the momenta of the unstable modes increase with µ, and that, in contrast to the continuum approach, the available momentum modes on the lattice are restricted.
So far we have only discussed lattice results with the discretization W 2 = W 2 . Using the lattice discretization with W 2 = W 2 (see Eq. (79)) we did not find instabilities anywhere in the µ-T plane, neither for aσ 0 = 0.379 nor for aσ 0 = 0.115. When using naive lattice fermions, momenta |k 1 | > ∼ π/2 or |k 2 | > ∼ π/2 in the Fourier expansion of σ lead to unphysical interaction terms not present in the GN model, as discussed in Sec. 4.1.2 and in more detail in Ref. [29]. Thus, these large momenta need to be suppressed, to obtain lattice actions with the correct continuum limit. This is done in different ways for the two lattice discretizations we are using. W 2 is a smooth cosine function in momentum space, while W 2 is a step function. Both discretizations yield identical results in the continuum limit, but W 2 generates a slight suppression ∝ a 2 also for small momenta, which are relevant for the formation of inhomogeneous instabilities. In contrast to that W 2 does not generate such a suppression. This could be the reason why the discretization with W 2 = W 2 leads to instability regions at finite a, while the discretization with W 2 = W 2 does not. We will come back to this issue in Sec. 5.3.

Minimization of the effective action allowing inhomogeneous modulations
As pointed out before, the true inhomogeneous phases could be larger than the instability regions found in Sec. 5.2. The latter are are fully included in the former, but the former might extend further if somewhereσ is just a local minimum of S eff with the global minimum given by an inhomogeneous σ(x). In particular the non-existence of instability regions in the renormalized model (at least for T = 0) does not generally exclude the existence of an inhomogeneous phase in this limit. In the following, in order to explore this possibility, we perform numerical minimizations of the lattice discretized effective action (81), allowing inhomogeneous modulations of the condensate. To this end, we use W 2 as well as W 2 .

Minimization using a cosine ansatz for the condensate
We start by restricting the condensate to with fixed integer n and minimize S eff with respect to α. For small α this ansatz corresponds to one of the terms in the sum of Eq. (91), i.e., it represents a particular perturbation investigated in the stability analyses of Sec. 5.2. The motivation for this subsection is to provide further support that the two lattice discretizations corresponding to W 2 = W 2 and W 2 = W 2 coincide in the continuum limit. Moreover, this particular cosine ansatz allows to make contact to and perform a cross check with existing lattice results from Ref. [34].
As an example we focus on n = 3 at (µ/σ 0 , T /σ 0 ) = (1.035, 0.110). For the lattice discretization with W 2 = W 2 and the two coarser lattice spacings aσ 0 ∈ {0.379, 0.237} the associated Γ −1 (see Eq. (92)) is negative, i.e., in these cases the ansatz leads to an inhomogeneous phase. For W 2 = W 2 and aσ 0 = 0.086 as well as for W 2 = W 2 the associated Γ −1 is positive. In Fig. 8 we plot the effective action as a function of α. The plot confirms our findings from Sec. 5.2: α = 0 is the location of a maximum of S eff for W 2 = W 2 and aσ 0 ∈ {0.379, 0.237} and of a minimum in the other cases. The qualitatively different behavior of S eff for W 2 = W 2 and W 2 = W 2 , in particular for larger a, is the reason, why the discretization with W 2 = W 2 leads to an inhomogeneous phase at finite a, while W 2 = W 2 does not. From Fig. 8 one can also see that the effective actions for W 2 = W 2 and W 2 = W 2 are approaching each other for decreasing a. This indicates that in the limit a → 0 both actions converge to the same action, the action of the renormalized 2 + 1-dimensional GN model, as theoretically expected (see the discussion in Sec. 4.1.2 and Ref. [29]).
The ansatz (96) together with a numerical minimization was already used in the previous lattice field theory study of the 2 + 1-dimensional GN model in Ref. [34]. There, however, a different discretization was used, equivalent to Eq. (78) with W 2 (x − y) = δ x,y . While this discretization does not correspond to the GN model for an arbitrary spatially varying condensate σ(x) (see Sec. 4.1.2 and the detailed discussion in Ref. [29]), it becomes identical to W 2 = W 2 if σ is restricted to a cosine-shaped modulation as in Eq. (96) with |n| < N s /4. Of particular interest is a comparison of Fig. 8 of Ref. [34] and our results from Fig. 8 for W 2 = W 2 . Since we do not use the same lattice spacings and spatial volumes as used in Ref. [34], a precise quantitative comparison is not possible. However, the two figures are qualitatively identical and indicate consistency of the results presented in Ref. [34] and our results. We show results for both lattice discretizations W 2 = W 2 and W 2 = W 2 and three lattice spacings a. S 0 is a physically irrelevant a-dependent shift chosen such that S eff | α=0 is identical for all three lattice spacings.

Minimization allowing arbitrary 1-dimensional modulations of the condensate
Finding the global minimum of the effective action for arbitrary σ = σ(x 1 , x 2 ) is time consuming, because it is quite expensive to evaluate S eff . Thus, we restrict σ to arbitrary 1-dimensional modulations, i.e., consider σ = σ(x 1 ). S eff can then be evaluated more efficiently, because the determinant of the Dirac operator factorizes into determinants of smaller matrices of size N s N d × N s N d (see Eq. (86)). Moreover, the number of variables for the minimization is reduced from N 2 s to N s . To search for the global minimum of S eff for given (µ, T ), a Fletcher-Reeves conjugate gradient algorithm is used, as implemented in the GNU Scientific Library [44]. This algorithm is suited to compute local minima of a given function. We try to find the global minimum of S eff by carrying out several local minimizations with different initial field configurations σ(x 1 ). Some of these initial configurations are proportional to a cosine as in Eq. (96), but also randomly generated initial configurations are used. For values (µ, T ), where we find local minima, the corresponding condensates are always periodic and oscillating. We note that local minimizations with randomly generated initial configurations do not lead to additional local minima, but to the same minima already found with initial configurations proportional to a cosine. We interpret this as indication that the found σ(x 1 ) with the smallest corresponding value for S eff represents the global minimum.
The shape of the condensate inside the instability region for W 2 = W 2 As discussed in Sec. 5.2, the instability regions found for the lattice discretization with W 2 = W 2 and finite lattice spacing (see Fig. 5 and Fig. 7) are regions where the condensate exhibits inhomogeneous modulations. By minimizing S eff we determine the shape of these modulations. In Fig. 9 we show the condensate σ as a function of x 1 for aσ 0 = 0.379, T /σ 0 = 0.132 and two values of the chemical potential, µ/σ 0 = 0.97 (left plot) and µ/σ 0 = 1.11 (right plot). For the smaller value of µ the number of oscillations is smaller and the amplitude is larger than for the larger value of µ. Moreover, at µ/σ 0 = 0.97 the condensate deviates from a cosine, while at µ/σ 0 = 1.11 the modulation is essentially a cosine. For even larger µ, when approaching the boundary of the instability region at µ/σ 0 = 1.134, the amplitude tends to zero, whereas the wavelength is still finite and consistent with the unstable momentum mode in the stability analysis. Note that the same behavior was observed for the 1 + 1-dimensional GN model inside the inhomogeneous phase [13,14]. Local minima of S eff in the homogeneous symmetry-broken phase Inside the σ = 0 regions shown in the right plot of Fig. 4 we find in addition to the favored homogeneous configuration σ =σ several local minima corresponding to inhomogeneous modulations. None of these minima leads to a smaller value of S eff than the constant condensate σ =σ, which seems to represent the global minimum. Only rather close to the boundary of the homogeneous symmetry-broken phase the existence and properties of the inhomogeneous minima depend on the lattice discretization, i.e., whether we use W 2 = W 2 or W 2 = W 2 , and on the lattice spacing a. Farther inside the homogeneous symmetry-broken phase they are almost independent of W 2 and of a. This indicates that these local inhomogeneous minima are also present in the renormalized 2 + 1-dimensional GN model. As an example we show in Fig. 10 the condensate σ(x 1 ) corresponding to one such local minimum at (µ/σ 0 , T /σ 0 ) = (0.60, 0.176) and lattice spacing aσ 0 = 0.237. As one might expect from the minimizations in the instability region discussed above, the wavelengths of the condensates corresponding to the local inhomogeneous minima are quite large, their amplitudes are close to σ 0 and the shape is somewhere between a cosine and a kink-antikink structure.
We also studied the differences between S eff (σ) and S eff evaluated at the local minima. These differences are positive in the homogeneous symmetry-broken phase, since σ =σ represents the global minimum of S eff . For fixed T and increasing µ the differences become smaller. For the discretization with W 2 = W 2 and T < T LP the difference to the local minimum corresponding to a kink-antikink with wavelength L approaches 0 extremely close to or exactly at the boundary to the instability region. For the discretization with W 2 = W 2 inhomogeneous minima cease to Figure 10: The condensate σ(x 1 ) corresponding to one of the local minima at (µ/σ 0 , T /σ 0 ) = (0.60, 0.176) and lattice spacing aσ 0 = 0.237. There is essentially no difference between the two lattice discretizations W 2 = W 2 and W 2 = W 2 . The dotted line represents a cosine function with the same wavelength and amplitude. exist near the boundary to the instability region. For fixed µ and decreasing T the differences between S eff (σ) and S eff evaluated at the local minima also decrease. This is consistent with Ref. [45], where the 2 + 1-dimensional GN was studied at T = 0 using a specific ansatz for σ(x 1 ) based on Jacobi elliptic functions. There it was found that such inhomogeneous modulations represent minima of S eff , which are degenerate with the minimum at σ =σ.

Phase diagram
Within the σ = 0 regions shown in the right plot of Fig. 4 and (for W 2 = W 2 ) outside the instability region we do not find any local minima corresponding to inhomogeneous modulations. Thus, the phase diagram of the 2 + 1-dimensional GN model, when allowing arbitrary 1-dimensional modulations, is identical to that already found by stability analyses (see Fig. 4, right plot for W 2 = W 2 and Fig. 5 and Fig. 7, right columns for W 2 = W 2 ). In particular we confirm and consolidate the findings from Ref. [34], where some evidence was presented that there is no inhomogeneous phase in the continuum limit.
It is, however, important to note, that we did not yet carry out minimizations for 2-dimensional modulations. Such 2-dimensional modulations of the condensate might lead to lower values of S eff and, thus, could generate larger inhomogeneous phases at finite a or inhomogeneous phases, which survive the continuum limit. We plan to study this possibility in the near future.

Conclusions
In this work we explored in detail the phase diagram of the 2 + 1-dimensional GN model in the limit of infinitely many flavors. We implemented three different regularizations, a continuum regularization with a Pauli-Villars cutoff Λ and two lattice field theory regularizations with lattice spacing a, which are based on naive fermions. Particular focus was put on studying the possible existence of inhomogeneous phases, their location in the µ-T plane and their dependence on the regularization and the corresponding regulator, Λ or a, respectively.
Our main results are the following: • For finite values of the regulator, inhomogeneous phases may exist, depending on the details of the regularization: For the continuum regularization and one of the two lattice field theory regularizations we found an inhomogeneous phase, while for the other lattice field theory regularization there is no inhomogeneous phase.
• Even if there is an inhomogeneous phase at finite values of the regulator, it seems to disappear, when the regulator is removed, i.e., in the limit Λ → ∞ or a → 0.
These results confirm existing results, e.g., continuum results for the homogeneous phase diagram [31,32] or lattice field theory results at a single lattice spacing [33] and at several lattice spacings, but with a specific cosine ansatz for the condensate [34]. Our work also substantially extends these existing results, in particular by an analytical stability analysis of homogeneous phases with respect to arbitrary inhomogeneous perturbations and by a full numerical minimization of the effective action, where arbitrary 1-dimensional modulations of the condensate are allowed, i.e., without restriction to a specific ansatz like plane waves or Jacobi elliptic functions as done in existing work [34,45].
It is important to note that our numerical minimization is currently limited to arbitrary 1dimensional modulations. Thus, an important next step is to allow arbitrary 2-dimensional modulations. This could not only lead to larger or additional inhomogeneous phases at finite values of the regulator, but also to the existence of inhomogeneous phases in the renormalized 2 + 1-dimensional GN model, i.e., in the limit Λ → ∞ or a → 0.
Our results call for a critical revision of the role of the regularization in the physically more relevant case of 3 + 1 space-time dimensions. For instance, inhomogeneous phases have also been found in the 3 + 1-dimensional NJL model. However, unlike in 2 + 1 dimensions, this model is non-renormalizable and therefore the studies have been performed using finite fixed regulators. Given that inhomogeneous phases exist in the renormalized 1 + 1-dimensional GN and NJL models, while in the 2+1-dimensional GN model they are only present at finite regulator values, one might suspect that the observed inhomogeneous phases at 3 + 1 dimensions could be regularization artifacts. The QM model, on the other hand, is renormalizable in 3 + 1 dimensions, and inhomogeneous phases have been reported to exist in that model even in the renormalized limit. Unfortunately, the model suffers from other instabilities at large cutoff values, at least in mean-field approximation. In the near future we therefore plan to extend our detailed investigations to 3 + 1-dimensional models in order to shed light on these issues.
Thus, if α j is an eigenvalue of Q[+σ], −α * j is an eigenvalue of Q[−σ]. Consequently, where we have again used that the number of eigenvalues is even. Combining Eq. (99) and Eq.