Inhomogeneous phases in the chirally imbalanced $2+1$-dimensional Gross-Neveu model and their absence in the continuum limit

We study the $\mu$-$\mu_{45}$-$T$ phase diagram of the $2+1$-dimensional Gross-Neveu model, where $\mu$ denotes the ordinary chemical potential, $\mu_{45}$ the chiral chemical potential and $T$ the temperature. We use the mean-field approximation and two different lattice regularizations with naive chiral fermions. An inhomogeneous phase at finite lattice spacing is found for one of the two regularizations. Our results suggest that there is no inhomogeneous phase in the continuum limit. We show that a chiral chemical potential is equivalent to an isospin chemical potential. Thus, all results presented in this work can also be interpreted in the context of isospin imbalance.


Introduction
The Gross-Neveu (GN) model describes a theory of N f fermion flavors with a quartic interaction. It is a rather simple model commonly used to explore and describe the spontaneous breaking of chiral symmetry [1] in the µ-T plane, where µ denotes the chemical potential and T the temperature. In the limit N f → ∞ (corresponding to the meanfield approximation or, equivalently, the neglect of bosonic quantum fluctuations) the 1 + 1-dimensional GN model exhibits three phases: a symmetric phase (with a vanishing chiral condensate), a homogeneous symmetry-broken phase (with a non-zero, but spatially constant condensate) and an inhomogeneous phase, where the chiral condensate is an oscillating function of space [2][3][4]. The phase diagrams of the GN model and related theories were also investigated at finite N f , i.e., with bosonic quantum fluctuations included, using lattice Monte-Carlo simulations [5][6][7][8][9][10][11] and the functional renormalization group (FRG) [12].
Inhomogeneous phases are not limited to the GN model, but were found in several models in the mean-field approximation in 1 + 1 dimensions [13][14][15][16] and in 3 + 1 dimensions [17][18][19][20][21][22][23]. For a review, we refer to Ref. [24]. In recent works [25][26][27][28] it has been discussed that inhomogeneous phases might be related to so-called moat regimes, where the bosonic wave function renormalization Z is negative. Such a regime has been found in an FRG study of the phase diagram of quantum chromodynamics (QCD) [29]. A similar regime has been observed in the 1 + 1-dimensional GN model (see Figure 9 of Ref. [30] , where the wave function renormalization Z is plotted in the µ-T plane). One finds that a negative Z accompanies the instability of a homogeneous condensate with respect to inhomogeneous perturbations as a necessary, but not a sufficient condition. The possible existence of a moat regime in the QCD phase diagram encourages us to improve our understanding of inhomogeneous condensation and related phenomena in QCD-inspired models.
Recently, the existence of inhomogeneous phases was also explored in the 2 + 1dimensional GN model in the mean-field approximation [31][32][33]. Such 2 + 1-dimensional four-fermion theories are of interest both in high energy physics [34][35][36][37][38] and in condensed matter physics [39][40][41][42][43][44][45][46], but also to study conceptual questions, e.g., renormalizability in the 1/N expansion or in a perturbative approach [47][48][49][50]. Hence, confirming the existence of an inhomogeneous phase in such a model could have a significant impact. Early seminal studies of the µ-T phase diagram of the 2 + 1-dimensional GN model [51,52] reported a second-order phase transition between the symmetric and the homogeneous symmetrybroken phase at finite T and µ and a first-order phase transition at T = 0. However, in these studies only a homogeneous order parameter was considered. In our recent publication Ref. [33] we studied the existence of an inhomogeneous phase in the 2 + 1dimensional GN model within the mean-field approximation. Our main findings were that an inhomogeneous phase is present at finite regulator and for certain regularization schemes (a Pauli-Villars cutoff and a specific lattice discretization), but it disappears when the regulator is removed, as previously observed in Ref. [38].
In this work we continue our investigations from Ref. [33] by extending the 2 + 1dimensional GN model with a chiral chemical potential. We studied its phase diagram, where our main focus was on possibly existing inhomogeneous phases. While the GN model might be too simple to realistically describe the effect of chiral imbalance on QCD, it might still improve our conceptual understanding of inhomogeneous condensation in the presence of chiral imbalance, which is an important problem. A difference in the densities of left-and right-handed quarks is relevant in physical systems such as heavy-ion collisions [53,54] or compact stars [55,56]. The impact of chiral imbalance on chiral symmetry breaking has been studied (e.g., in 1 + 1-dimensional models) [14], where it had no influence on the existence of the inhomogeneous chiral spiral, and in 3 + 1-dimensional models [57][58][59], where only homogeneous order parameters were considered. A chirally imbalanced 2 + 1-dimensional GN model, extended by a quartic difermion interaction, was explored in Refs. [60,61] with the aim to clarify the competition of homogeneous fermion-fermion condensation and homogeneous chiral condensation. In recent two-color and three-color QCD studies [62,63] a chiral chemical potential was investigated and found to increase the chiral transition temperature. This result is supported by a Nambu-Jona-Lasinio (NJL) model study [64].
This paper is structured as follows. We start in Section 2 by discussing the theoretical basics of the GN model in 2 + 1 dimensions including details on the underlying chiral symmetry. We also add a chiral chemical potential to the model and show the equivalence of chiral imbalance and isospin imbalance. In Section 3 we discretize the effective action of the model using lattice field theory. Section 4 is the main part of our paper, where numerical results are presented and discussed. Finally, we conclude in Section 5. Preliminary results from this project were presented at a recent conference [65].

Action and Partition Function
The action of the GN model in 2 + 1 dimensions with N f fermion flavors is where ψ n represents N f massless fermion fields; µ is the chemical potential; and g 2 is the coupling of the four-fermion interaction.
x with d 2 x = dx 1 dx 2 and T denoting the temperature given by the inverse extent of the periodic temporal direction of Euclidean space-time.
The action (1) is equivalent to where σ is a scalar boson field; λ = N f g 2 is the rescaled coupling; and is the Dirac operator. Integration over the fermion fields leads to the so-called effective action and the corresponding partition function One can show that σ(x) is related to the condensate ψ n (x)ψ n (x) according to As in previous studies of chiral inhomogeneous phases, we restricted the dependence of σ to the spatial coordinates, i.e., σ = σ(x). With this restriction Det Q is real, which is shown in the Appendix of Ref. [33]. For even N f , the effective action S eff [σ] is then real. (Our numerical calculations of the determinant showed that Det Q is exclusively positive, i.e., the effective action S eff [σ] is real for all values of N f ).
Since S eff [σ] ∝ N f , the limit N f → ∞ reduces the relevant configurations in the partition function (4) to the global minima of S eff [σ]. Thus, the computation of a path integral is reduced to an optimization problem. In the case of degenerate global minima, spontaneous symmetry breaking selects one of these minima. Consequently, an expectation value O(σ) is identical to the value of O evaluated at the corresponding global minimum, i.e., O(σ) → O(σ). In particular, σ → σ. For the remainder of this paper we exclusively consider the limit N f → ∞.

Representation of the Dirac Matrices and Chiral Symmetry
Typically one uses either an irreducible 2 × 2 representation or a reducible 4 × 4 representation of the Dirac algebra for the γ matrices appearing in the Dirac operator (3) (for details see, e.g., Refs. [33,36,37,66]). In the case of an irreducible 2 × 2 representation, there is no symmetry, which can be interpreted as chiral symmetry, because it is impossible to define a matrix γ 5 , which anticommutes with γ 0 , γ 1 and γ 2 . Therefore, a reducible 4 × 4 representation is more appropriate in our context, e.g.: where τ j denote the Pauli matrices. The three matrices +τ 1 , +τ 2 and +τ 3 as well as the three matrices −τ 1 , −τ 2 and −τ 3 form irreducible 2 × 2 representations, which are inequivalent. The corresponding upper two and lower two entries of the fermion fields ψ n can be interpreted as left-handed and right-handed components, respectively.
The Dirac operator (3) is then block-diagonal, where represent Dirac operators for left-handed and right-handed fermion fields ψ L/R n (see also Equations (14) and (15) in Ref. [33]). One can show that Det Q (2) [µ, σ] and DetQ (2) [µ, σ] are invariant under both µ → −µ and σ → −σ. Using the latter one can show The action with Q = Q (4) [µ, σ] is invariant under the discrete chiral transformations. (For free fermions one can define continuous axial chiral symmetry transformations with both γ 4 and γ 5 . The four-fermion interaction in Equation (1) breaks the corresponding symmetries explicitly, see, e.g., Ref. [33]). with Both γ 4 and γ 5 anticommute with γ 0 , γ 1 and γ 2 , thus fulfilling the necessary properties for generating an axial chiral transformation. The symmetries (11) and (12) are also present for the action (2), where the corresponding transformation of σ is in both cases σ → −σ. Thus, σ is an order parameter for chiral symmetry breaking. (When using an irreducible 2 × 2 fermion representation, there is no chiral symmetry; however, σ can still be interpreted as an order parameter for parity breaking).
Thus, there is only one independent Z 2 symmetry, i.e., the structure of chiral symmetry A chiral chemical potential µ 45 can be introduced in a straightforward way by extending and replacing the Dirac operator in (3) or equivalently (7) according to µ 45 contributes to the chemical potentials of the left-handed (upper two) components and the right-handed (lower two) components with opposite sign, thus causing chiral imbalance.
We note that there are other possibilities to define chirality and chiral imbalance (see, e.g., Refs. [60,61] and Section 5) differing from our definition, where left-and right-handed fermion fields ψ L/R n are projected from the fermion fields as with P L/R denoting the corresponding projectors.
As done for µ 45 = 0 in appendix A of Ref. [33], one can show that Det Q (4) Det Q (4) [µ, µ 45 , σ] is also invariant under the exchange of the ordinary and the chiral chemical potential, µ ↔ µ 45 . Clearly, S eff [σ] as well as the phase diagram share this invariance. In Section 4 we use this property to cross-check our numerical results. We note that the effective action can be written as the sum of a left-handed and a right-handed part, with µ L = µ + µ 45 and µ R = µ − µ 45 . Of course, the two parts are not independent but coupled via σ. Moreover, both parts are equivalent to the chirally balanced effective action (see Section II C. of Ref. [33]), i.e., This property will be useful when we discuss our numerical results in Section 4.

Equivalence of Isospin and Chiral Imbalance
In this subsection we consider an even number of fermion flavors N f , again in the 4-component reducible representation, and assign half of them (the "u flavors") a chemical potential µ + µ I , the other half (the "d flavors") a chemical potential µ − µ I . Clearly, µ I generates an imbalance between the u and the d flavors and, thus, can be interpreted as isospin chemical potential.
The corresponding effective action is where the Dirac operator is an 8 × 8 matrix in spin and isospin space, Using Equation (20) one can show Consequently, the effective action for the GN model with isospin imbalance, Equation (23), is identical to the the effective action for the GN model with chiral imbalance, Equation (21), when identifying µ I = µ 45 . Thus, all numerical results presented in Equation (4) can either be interpreted in the context of chiral imbalance or of isospin imbalance. We note that this equivalence of isospin and chiral imbalance is specific to the GN model in 2 + 1 dimensions.

Lattice Discretization
We used a lattice discretization of the effective action (21), which was similar to the discretization discussed in Section IV of our previous work [33] . The key difference is that we use the naive fermion discretization also in temporal and not only in the spatial directions.
We considered a 3-dimensional space-time volume βV, where β = 1/T was the inverse temperature and V = L 2 the quadratic spatial volume. The boundary conditions were periodic in the 2 spatial directions and periodic and antiperiodic in temporal direction for the fields σ and ψ n ,ψ n , respectively. We used a cubic lattice with N t × N 2 s lattice sites and lattice spacing a, i.e., β = aN t and L = aN s . In the following, all dimensionful quantities are expressed in units of the lattice spacing, e.g., a ≡ 1, β ≡ β/a. Because of the finite space-time volume, the 3-dimensional momenta were quantized, with and η = 0, 1/2, corresponding to periodic and antiperiodic boundary conditions in the temporal direction. In our numerical implementation, the effective action and fields were treated in momentum space, are the Fourier coefficients of the field σ(x).
is the Dirac operator in momentum space. On the lattice, this operator is a matrix with columns and rows labeled by the momenta p and q, respectively. The sin for ν = 0 contains the matrix γ 45 but can be simplified according to The Dirac operatorQ (4) p,q has eight regions of soft modes, where the dispersion relation is approximately linear, in the first Brioullin zone and where each region describes a fermion flavor. In the continuum limit, these fermion flavors do not interact with each other but with the scalar field σ in the same manner (for details see Ref. [5] and the Appendix of Ref. [7]). Thus, in order to study N f fermion fields on the lattice, where N f is restricted to a multiple of 8, one has to use N f /8 naive fermions in the discretization of the fermion bilinear in Equation (2) resulting in the factor 1/8 in Equation (27) (compare Section IV A. in Ref. [33]). An appropriately chosen weight functionW 2 (p) was necessary to ensure the correct continuum limit (see Refs. [5,7,33] for details). We investigated and compared two possible choices, with Θ denoting the Heaviside function.
Because of the restriction of σ to the spatial coordinates, i.e., σ = σ(x), the Dirac operator (29) was block-diagonal with respect to p 0 and q 0 . This simplified the computation of DetQ (4) [µ, µ 45 , σ] to the computation of N t determinants of smaller matrices of size 4N 2 s × 4N 2 s .

Numerical Results
Using lattice field theory the phase diagram of the 2 + 1-dimensional GN model with µ 45 = 0 was extensively explored in Refs. [31][32][33]. There is a symmetric phase with σ = 0 at large µ or large T and a homogeneous symmetry-broken phase with a constant σ =σ at small µ and small T. Moreover, at finite lattice spacing and for certain discretizations (e.g., W 2 =W 2 ) there is additionally an inhomogeneous phase, where σ(x) is a varying function of the spatial coordinates. However, this inhomogeneous phase shrinks, when decreasing the lattice spacing, and seems to vanish in the continuum limit.
The main focus of this paper is to investigate in particular the phase structure for µ 45 = 0 to clarify whether inhomogeneous phases exist. At first, we recalled that the effective action S eff [σ] can be written as the sum of a left-handed part S L eff [σ] and a right-handed part S R eff [σ] with chemical potentials µ L and µ R , respectively (see Equation (21)). Moreover, each of the two parts was equivalent to the action of the chirally balanced GN model, which was investigated in detail in our previous work [33]. Thus, for |µ L | > µ c (T) and |µ R | > µ c (T) both S L eff [σ] and S R eff [σ] had their respective minima at σ = 0 (µ c (T), denoting the location of the phase boundary of the symmetric phase at µ 45 = 0 and temperature T).
(We ignored the existence of an "inhomogeneous island" or "inhomogeneous continent" at large chemical potentials where cutoff effects were particularly strong [33].). Consequently, the minimum of S eff [σ] also corresponded to σ = 0. In other words, from numerical results obtained in Ref. [33] at µ 45 = 0 we could conclude that chiral symmetry was restored in the chirally imbalanced GN model for |µ L | > µ c (T) and |µ R | > µ c (T). In the remaining regions of the (µ, µ 45 , T) space, S L eff [σ] and S R eff [σ] compete and the behavior of the condensate needed to be investigated numerically. We started doing that in Section 4.1 by restricting our computations to a homogeneous condensate. After that, in Section 4.2, we carried out a stability analysis of the favored value of the homogeneous condensate with respect to inhomogeneous perturbations. Finally, in Section 4.3, we performed numerical minimizations of the effective action, allowing arbitrary inhomogeneous modulations.
The lattice spacing a was a function of the coupling λ. As explained in our previous work [33] we tuned λ such that the temporal extent N t,c a corresponded to the inverse critical temperature β c = 1/T c , which separated at µ = µ 45 = 0, the symmetric and homogeneous symmetry-broken phase. Then, at fixed λ, the temperature T = 1/N t a could be changed in discrete steps by increasing or decreasing N t . A summary of the lattice parameters used to generate all following numerical results is given in Table 1. We note that throughout this section dimensionful quantities are expressed in units of the vacuum expectation value of σ, σ 0 = σ| µ=0,µ 45 =0,T=0 .

Restriction to a Homogeneous Condensate
In the case of a homogeneous condensate, i.e., σ =σ or equivalentlyσ(p) =σδ p,0 , the two lattice discretizations withW 2 andW 2 (Equations (31) and (32)) were identical. The Dirac operator corresponded to a block-diagonal matrix with N t N 2 s blocks of size 4 × 4. Thus, the ln DetQ (4) [µ, µ 45 , σ] term in the effective action (27) could be computed quite efficiently by summing over N t N 2 s determinants of 4 × 4 matrices. Moreover, the effective action at given µ, µ 45 and T is a function of just a single variableσ; hence, it could be minimized numerically in a straightforward and rather cheap way to obtain the physically preferred value of the homogeneous condensate. Figure 1 shows the phase diagram in (µ, µ 45 , T) space for aσ 0 = 0.2327 and Lσ 0 = 120 aσ 0 = 27.92. For µ 45 = 0.0 the phase boundary is quite similar to the analytically obtained continuum result [51] with slight deviations due to discretization and finite volume effects. At high temperature T/σ 0 > ∼ 0.4 the phase boundary exhibited an approximate rotational symmetry in the µ-µ 45 plane, i.e., it was crudely described by µ 2 + µ 2 45 ≈ (µ c (T)) 2 . In contrast to that, at low temperature, the phase boundary approached a square-like shape in the µ-µ 45 plane.
The left plot of Figure 2 shows sectional views of the phase diagram at fixed lattice spacing aσ 0 = 0.2327 for two different spatial extents, Lσ 0 = 60 aσ 0 = 13.95 and Lσ 0 = 120 aσ 0 = 27.92. At low temperature the phase boundary exhibits an oscillatory behavior, which was more pronounced for the smaller lattice volume. We expected that the oscillations would disappear in the infinite volume limit. The right plot of Figure 2 shows sectional views of the phase diagram for two different lattice spacings, aσ 0 = 0.3649 and aσ 0 = 0.2327, at fixed ratio N s /N t,c = 20 implying similar spatial extents Lσ 0 = 80 aσ 0 = 29.19 and Lσ 0 = 120 aσ 0 = 27.92. There were visible discrepancies due to discretization effects, in particular when both T is small and µ ≈ µ 45 . Continuum results at µ 45 = 0 from Ref. [51] as well as our lattice results at various small temperatures, lattice spacings and spatial volumes point towards T = 0 phase boundaries at µ/σ 0 = 1 for 0 ≤ µ 45 /σ 0 < 1 and at µ 45 /σ 0 = 1 for 0 ≤ µ/σ 0 < 1 in the continuum limit. We also studied the behavior ofσ at small temperature T/σ 0 = 0.0716, aσ 0 = 0.2327 and Lσ 0 = 120 aσ 0 = 27.92.σ is shown as function of µ and µ 45 in the left plot of Figure 3 and as function of µ L = µ + µ 45 and µ R = µ − µ 45 in the right plot of Figure 3. To explain these results, we noted that the effective action (21) was the sum of a left-handed part S L eff [σ] and a right-handed part S R eff [σ] with chemical potentials µ L and µ R , respectively. At the beginning of Section 4, we had already concluded thatσ = 0, if |µ L | > µ c (T) and |µ R | > µ c (T). Similarly, we argue now that both parts favorσ ≈ σ 0 , if |µ L | < µ c (T) and |µ R | < µ c (T). The numerical results from Figure 3 are consistent with that expectation. In particular, the yellow regions, whereσ ≈ σ 0 , corresponds to |µ L | < µ c (T) and |µ R | < µ c (T). In the remaining regions of the µ-µ 45   We note that the lattice data shown in this subsection also represents a non-trivial cross-check of our implementation: all numerical results were consistent with the symmetry µ ↔ µ 45 within machine precision.

Stability of a Homogeneous Condensate
Now, we relaxed the constraint that σ was a homogeneous condensate. To determine the preferred modulation of the condensate in a possibly existing inhomogeneous phase, we had to allow arbitrary spatial modulations of σ, i.e., consider σ = σ(x), and minimize the effective action with respect to these modulations. In lattice field theory this is possible but numerically very challenging. As a first step, therefore, we explored whether the homogeneous minima σ =σ, which were determined in Section 4.1 for many different (µ, µ 45 , T), were stable or unstable with respect to spatially inhomogeneous perturbations δσ(x). Boundaries between stable and unstable regions in (µ, µ 45 , T) space were identical to phase boundaries if the amplitude of the inhomogeneity became infinitesimal when approaching the boundary. However, a stability analysis failed to detect inhomogeneous condensates in regions of the phase diagram where the homogeneous minimum (found, e.g., as described in Section 4.1) corresponded to a local, but not global, minimum of S eff [σ(x)]/N f . This was, e.g., the case in the 1 + 1-dimensional GN model [30,67].
A detailed derivation of the formalism to probe the stability of a homogeneous condensate σ =σ with respect to arbitrary spatial perturbations δσ(x) can be found in Ref. [33] for a continuum approach. This formalism can be transferred to lattice discretizations in a straightforward way, which is discussed in the same reference. A quantity of central importance is where ∑ p runs over all 3-dimensional lattice momenta (26), q = (0, q k ), the trace refers to spinor space andQ p [µ, µ 45 ,σ] is defined viaQ (4) p,q [µ, µ 45 ,σ] = δ p,qQp [µ, µ 45 ,σ] and Eq. (29), i.e.,Q Negative values of Γ −1 (q k )/N f with q k = 0 indicate instability of the condensate σ =σ with respect to harmonic perturbations with momentum q k . Such perturbations decrease S eff [σ]; consequently, an inhomogeneous condensate was preferred. By evaluating Γ −1 (q k )/N f for suitably chosen parameters (µ, µ 45 , T) we could identify regions that were part of an inhomogeneous phase.
We searched extensively for regions, whereσ was unstable, using both discretizations (31) and (32). ForW 2 =W 2 such regions did not seem to exist. ForW 2 =W 2 and finite lattice spacing there was a region of instability at small T consistent with the findings at µ 45 = 0 reported in Ref. [33]. Figure 4 shows the lattice with the finer lattice spacing, aσ 0 = 0.2327, and spatial extent Lσ 0 = 100 aσ 0 = 23.27. The region of instability was located within the tetrahedral shape. At smaller temperature it had a larger extent in the µ-µ 45 plane. A somewhat unexpected result was the large extent of the region of instability in µ 45 direction (e.g., for T/σ 0 = 0.076 and µ/σ 0 ≈ 1.0 up to µ 45 /σ 0 ≈ 0.5). Its boundary is plotted in Figure 5 in the µ R -µ L plane. The plot shows that the instability region extended up to µ L = µ + µ 45 ≈ 1.5 and at the same time down to µ R = µ − µ 45 ≈ 0.5. A symmetric phase was preferred by S L eff [σ] with chemical potential µ L ≈ 1.5, while S R eff [σ] with chemical potential µ R ≈ 0.5 prefered a homogeneous symmetry-broken phase. Thus, neither of the two parts of the effective action (21) favored an inhomogeneous phase, but in combination they did. This highlighted the non-trivial interplay of S L eff [σ] and S R eff [σ] gave rise to a rather large inhomogeneous phase at finite lattice spacing for certain discretizations. Note that Figure 5 also reveals that the homogeneous phase boundary was engulfed by the region of instability, which was not the case in our previous study at µ 45 = 0 [33], where a different lattice regularization was used.  In Figure 6 we show sectional views of the region of instability for the discretizatioñ W 2 =W 2 and various temperatures. The upper row corresponds to the larger lattice spacing aσ 0 = 0.3649 and the lower row to the smaller lattice spacing, aσ 0 = 0.2327, while the left column corresponds to smaller and the right column to larger spatial extent Lσ 0 . From comparing the upper and the lower row, it is obvious that the instability region shrank when the lattice spacing decreased, most prominently in the µ direction. This is particularly evident in the right column, where the boundaries are significantly less distorted by finite volume effects. In the plots in the left column, however, there are pronounced oscillations that seem to be caused by small spatial volume. These oscillations are reminiscent of those observed in the µ-T plane in lattice studies of the chirally balanced GN model in 1 + 1 and 2 + 1 dimensions [31,33,67].
In summary, we found no region of instability forW 2 =W 2 but a shrinking region of instability for decreasing lattice spacing forW 2 =W 2 . This strongly suggests that there is no region of instability in the chirally imbalanced 2 + 1-dimensional GN model in the continuum limit.

Arbitrary Spatial Modulations of the Condensate
Now, we discuss the minimization of the effective action (27) with respect to the condensate allowing arbitrary spatial modulations, i.e., arbitrary Fourier coefficientsσ(p). We did this for selected parameters (µ, µ 45 , T) by carrying out several conjugate gradient minimizations, which differed in starting values forσ(p). For each (µ, µ 45 , T) we found only a small number of local minima although a significantly larger number of different starting values forσ(p) were provided to the minimization algorithm. This might have indicated that for all considered (µ, µ 45 , T) the corresponding global minimum was among the found local minima. We note that in our previous work [33] only 1-dimensional modulations were studied, i.e.,σ =σ(p 1 ). In this work we relaxed that constraint to allow arbitrary 2-dimensional modulations, i.e.,σ =σ(p). Since this was a numerically difficult and computer time-intensive task, we used a rather small lattice with coarse lattice spacing aσ 0 = 0.3649, temperature T/σ 0 = 0.114, Lσ 0 = 28 aσ 0 = 10.22 and discretization corresponding toW 2 =W 2 . First, we studied the chirally balanced model, i.e., µ 45 = 0. As an example, the upper left plot in Figure 7 shows a configuration σ(x) corresponding to one of the global minima of the effective action at µ/σ 0 = 1.041. An inhomogeneous condensate was favored, as we learned from the stability analysis discussed in Section 4.2. Even though the minimization algorithm allows arbitrary 2-dimensional modulations, the resulting global minimum was just a plane wave with wave vector q k = 2π(1, 2)/L. Similarly, the upper-right plot in Figure 7 shows a minimizing condensate at a larger chemical potential µ/σ 0 = 1.083. Again, we found a plane wave, but this time with wave vector q k = 2π(2, 2)/L, i.e., with a smaller wavelength. Even though the found plane waves were 1-dimensional structures, by allowing arbitrary 2-dimensional modulations it reduced finite volume corrections. This was so because, in contrast to our previous work [33], the direction of the wave vector was no longer restricted to being parallel to one of the coordinate axes; thus, its magnitude could be changed in finer steps.
The minimization algorithm also found 2-dimensional modulations. These, however, corresponded exclusively to local minima of the effective action (27). An example is shown in the center of Figure 7.
We also investigated, how chiral imbalance, i.e., µ 45 = 0, affected the preferred modulation of the condensate. The plots in the lower row of Figure 7 show the global minima of the effective action for (µ/σ 0 , µ 45 /σ 0 ) = (1.041, 0.05) and (µ/σ 0 , µ 45 /σ 0 ) = (1.041, 0.25). At fixed µ/σ 0 = 1.041 the frequency was almost independent of µ 45 (cf. the upper left plot, the lower left plot and the lower right plot of Figure 7). The amplitude, however, decreased when increasing |µ 45 |, suppored a second-order phase transition, also at µ 45 = 0. We searched extensively for inhomogeneous condensates outside the instability region explored and discussed in Section 4.2, but did not find any. Thus, we concluded that the boundaries of the instability region were identical to phase boundaries and that the corresponding phase transitions were of the second order.

Conclusions
In this work we studied the phase diagram of the 2 + 1-dimensional GN model with chiral imbalance introduced via a chiral chemical potential µ 45 using the mean-field approximation. Our lattice field theory results indicated that an inhomogeneous phase exists at finite lattice spacing a, when using a specific lattice discretization (W 2 =W 2 ). Non-vanishing µ 45 , however, seems to disfavor inhomogeneous modulations. Moreover, the inhomogeneous phase shrank for decreasing a and was expected to disappear in the continuum limit. These findings are consistent with our previous work [33], which was restricted to µ 45 = 0. (For completeness we note that inhomogeneous condensates at T = 0 and µ 45 = 0, which are energetically degenerate to the homogenous condensate in the homogeneous symmetry-broken phase, were found in Ref. [68]. Our numerical lattice studies were, however, limited to T > 0). Investigations of 1 + 1-dimensional GN-type models [5][6][7][8][9][10][11][12] at finite N f showed that bosonic quantum fluctuations can influence the extent and existence of (in-)homogeneous phases significantly. Such fluctuations are likely to affect the extent of the homogeneous symmetry-broken phase in the present 2 + 1dimensional GN model. The inhomogeneous phase was expected to remain absent since bosonic quantum fluctuations tend to disfavor ordered phases even more.
Moreover, for our chirally imbalanced 2 + 1-dimensional GN model, we showed that an isospin chemical potential µ I is equivalent to the chiral chemical potential µ 45 . Thus, all results presented can either be interpreted in the context of chiral imbalance or of isospin imbalance. In particular the µ-µ 45 -T phase diagram was identical to the µ-µ I -T phase diagram. Interestingly, a recent study of the 3 + 1-dimensional NJL model in the large-N c limit [57] conjectured a similar approximate duality of the phase diagram.
Color-superconductivity might play an important role in such studies [69,70], especially at finite isospin chemical potential [71,72]. In the considered 2 + 1-dimensional GN model this could, however, not be investigated, because the necessary difermion interaction was not present. Thus, we are not yet in a position to compare our results to up-to-date lattice QCD simulations at finite µ I (see, e.g., Refs. [73][74][75][76][77]), where a phase with Bose-Einstein condensation of charged pions was observed. As a next step, it might therefore be interesting to establish contact with Refs. [60,61], where a color-superconducting channel was added to the chirally imbalanced 2 + 1-dimensional GN model. It should be noted, however, that these references introduced the chiral chemical potential in a conceptually different way using spin matrices γ 0 γ 4 as well as γ 0 γ 5 , instead of γ 0 γ 45 , as in Equation (18).
Author Contributions: All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript. produced with the matplotlib package [82]. Calculations on the GOETHE-HLR and on the FUCHS-CSC high-performance computer of the Frankfurt University were conducted for this research. We would like to thank HPC-Hessen, funded by the State Ministry of Higher Education, Research and the Arts, for programming advice.

Conflicts of Interest:
The authors declare no conflict of interest.