Inhomogeneous condensation in the Gross-Neveu model in noninteger spatial dimensions $1 \leq d<3$. II. Nonzero temperature and chemical potential

We continue previous investigations of the (inhomogeneous) phase structure of the Gross-Neveu model in a noninteger number of spatial dimensions ($1 \leq d<3$) in the limit of an infinite number of fermion species ($N \to \infty$) at (non)zero chemical potential $\mu$. In this work, we extend the analysis from zero to nonzero temperature $T$. The phase diagram of the Gross-Neveu model in $1 \leq d<3$ spatial dimensions is well known under the assumption of spatially homogeneous condensation with both a symmetry broken and a symmetric phase present for all spatial dimensions. In $d = 1$ one additionally finds an inhomogeneous phase, where the order parameter, the condensate, is varying in space. Similarly, phases of spatially varying condensates are also found in the Gross-Neveu model in $d = 2$ and $d = 3$, as long as the theory is not fully renormalized, i.e., in the presence of a regulator. For $d = 2$, one observes that the inhomogeneous phase vanishes, when the regulator is properly removed (which is not possible for $d = 3$ without introducing additional parameters). In the present work, we use the stability analysis of the symmetric phase to study the presence (for $1 \leq d<2$) and absence (for $2 \leq d<3$) of these inhomogeneous phases and the related moat regimes in the fully renormalized Gross-Neveu model in the $\mu, T$-plane. We also discuss the relation between"the number of spatial dimensions"and"studying the model with a finite regulator"as well as the possible consequences for the limit $d \to 3$.

The Gross-Neveu (GN) model is arguably one of the most simple theories that describe (self-)interacting fermions.Despite this fact and having been formulated 50 years ago [2], its chiral phase diagram in the µ, T -plane in various number of spatial dimensions d is still under investigation today.Within this work, we aim to contribute to this research by focusing on so-called inhomogeneous phases (IPs) of spatially oscillating condensates in noninteger spatial dimensions.(We refer to Ref. [3] for a review on IPs.) A. General contextualization Within the N → ∞ limit, one finds that bosonic quantum fluctuations are suppressed [2], which immensely simplifies calculations and enables mostly analytic approaches.Thus, it is not surprising that the most complete picture of the thermodynamics of the GN model is within the special N → ∞ limit, see, e.g., Refs.[4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20].Still, even with these simplifications the GN model can be seen as a prototype quantum field theory (QFT) that shares a lot of features with more realistic QFTs.It is asymptotically free and undergoes dimensional transmutation, leading to a condensation of fermion-anti-fermion pairs in the infra red (IR) in vacuum, which is similar to Quantum Chromodynamics (QCD) and QCD-like theories.However, there are also important relations between the GN model and various models from solid state physics as well as numerous extensions of the model that are used as toy model QFTs, such that studying the model within different setups remains an interesting task on its own but is also of relevance as reference work.For further reading, we refer to Refs.[21,22].

B. Recap of central results
In 1 + 1 dimensions the GN model exhibits three distinct chiral phases [9,23].At low temperatures and chemical potential the discrete chiral Z 2 symmetry is spontaneously broken and one finds the so-called homogeneously broken phase (HBP) that is characterized by a nonzero chiral condensate, which is constant in space.At moderate and high temperatures one finds the symmetric phase (SP) -a gas-like phase, which is characterized by a vanishing chiral condensate.Especially relevant for our work is the IP at low temperatures and moderate and large µ, where the chiral condensate is non-vanishing and exhibits a spatial dependence.Thus, in addition to chiral symmetry also translational invariance is spontaneously broken.This phase is associated with negative values of the bosonic two-point function for some range of spatial bosonic momenta [13,[24][25][26][27][28][29][30][31][32][33][34][35][36][37][38].Assuming a second order phase transition (PT) to the SP, a necessary condition for the IP is a negative wavefunction renormalization even though this is not a sufficient criterion [38][39][40].Nevertheless, regions of negative wave-function renormalization in the µ, T -plane are interesting on their own, since they feature a nontrivial momentum structure of the two-point function and the dispersion relation.These regions that can be larger than the actual IP were labeled as moat regimes recently and could play an important role in the hadronization process in heavy-ion collisions [40][41][42][43].As discussed in Ref. [25] the 1 + 1 GN model has a moat regime, which extends over large parts of the µ, T -plane and thus serves as a toy model for this phenomenon.
It was found that the number of spatial dimensions has a profound impact on the phase structure.In 2 + 1 di-mensions, one finds that the IP and the moat regime are only present at finite regulator and vanish in the renormalized limit [13,14,44].Still, one observes an HBP for small µ and T and an SP for large µ and T [17,18].
The situation is less clear in 3+1 dimensions, where the model is non-renormalizable (without introducing additional parameters) and the value of the regulator and the choice of the regularization scheme has a drastic impact on the phase structure [45][46][47][48][49].Most certainly one also finds an HBP and an SP, while the existence of an IP can be regarded as disputed as it heavily depends on the regularization.
In an effort to understand why the IP is absent in 2 + 1 dimensions and to what extent the existence of the phase is a regulator artifact in 3 + 1 dimensions, Ref. [1] investigated the GN model in noninteger d spatial dimensions in order to interpolate between the known integer dimensional results and to additional mimic the effect of dimensional regularization. 1 This study was conducted at T = 0, which sufficed to illuminate that (a) the IP is present for 1 ≤ d < 2 and vanishes exactly in 2 = d, and (b) one does not find the phase in 2 < d < 3 in a renormalized setup, which implies that its existence in d = 3 is caused by a finite regulator.

C. Research objective
In the present work, we extend this investigation of spatially inhomogeneous condensation in (non)integer number of spatial dimensions from zero to nonzero temperature in order to map out the d-dependence of the IP and the moat regime.This study therefore complements the previously mentioned study in Ref. [1] as well as Ref. [11].The latter study already investigated the phase diagram of the GN model in continuous dimensions 1 ≤ d < 3 under the assumption of spatially homogeneous condensation.Therefore, we aim at closing a gap in the literature about the GN model.
Furthermore, we hope that this work contributes to the general discussion about necessary conditions for the presence/absence of IPs in arbitrary models and theories.Here, especially our findings about the role of the spatial dimensionality may be essential to understand the general criteria for the formation of inhomogeneous condensates.

D. Structure
This work is structured as follows: In Section II we recapitulate some basic mathematical aspects of the GN model.We present the four-fermion action, the bosonized version of the model in the N → ∞ limit and we discuss the quantities that are relevant to map out the phase diagram.Here, we also explain our regularization and renormalization prescription as well as the stability analysis -the method to detect inhomogeneous condensation.
Afterwards, in Section III, we turn to the results.We present the evaluation of the above expressions for some points in the µ, T -plane and different dimensions d.We show sample plots for the two-point function and wavefunction renormalization.Most importantly, we present the dependence of the phase diagram and especially of the IP and moat regime on the number of spatial dimensions d.
We conclude and comment on our results in Section IV and provide a brief outlook to possible consequences of our findings and followup questions.
Our work is accompanied by a large number of appendices, where we present all relevant details of this work.We hope that the amount of technical details might help the interested reader to easily reproduce and/or build on our work.We also consider these appendices as a compilation of the most relevant formulae for the GN model in d < 3 within the N → ∞ limit.

II. THE GROSS-NEVEU MODEL IN 1 ≤ d < 3 SPATIAL DIMENSIONS IN MEDIUM
In this chapter, we introduce the GN model in d + 1 dimensions, where d is the number of spatial dimensions.We work at (non)zero temperature T and (non)zero chemical potential µ and briefly recapitulate the derivation of the grand canonical/effective potential, the bosonic two-point function as well as the bosonic wave-function renormalization.All calculations and results are in the limit N → ∞, where N is the number of fermion species.The quantities that are presented in this chapter are required for the computation of the phase diagram, the stability analysis, and the detection of a possible IP and/or moat regime, see also Refs.[13,25,44] for similar analyses and results that arise as the limiting cases for integer d.

A. The action and the potential
The microscopic action of the GN model is given by [ Here, ψ = ψ(τ, ⃗ x ) and ψ = ψ(τ, ⃗ x ) are the fermion fields, where x i ∈ (−∞, ∞), i ∈ {1, . . ., d}, are the spatial coordinates of a d-dimensional Euclidean space and τ ∈ [0, 1 T ) denotes the coordinate of the compactified temporal direction that mimics the (inverse) temperature T .The fermions have antiperiodic boundary conditions in the compact direction, come in N different species 2 and transform as spinors.We use a d γ -dimensional representation of the gamma matrices of the corresponding Clifford algebra, This generalizes to noninteger d dimensions [50] and we use the Kronecker delta as the components of the Euclidean metric.Furthermore, we use λ for the fourfermion coupling and introduce the fermion chemical potential µ in the standard way.
For a detailed discussion of the symmetries of this model in different integer dimensions, we refer for example to Refs.[22,51].
In order to study four-fermion models (especially in the N → ∞ limit) one convenient approach is to bosonize the theory in the ultra violet (UV) via a Hubbard-Stratonovich transformation [52,53].Here, the fourfermion interaction is replaced by an auxiliary real scalar bosonic field ϕ.On the level of the partition function the equivalent action is [2,4,8] where we also introduced the Yukawa coupling h to obtain canonical energy dimensions for ϕ.Now, in addition to the functional integration over fermion fields, one also has to integrate over the bosonic field ϕ.Since loop corrections for the Yukawa coupling are suppressed for N → ∞, see, e.g., Refs.[22,54], it is convenient to absorb the Yukawa coupling in the field ϕ which is afterwards of dimension energy in any number of spatial dimensions.In addition, to correctly take the limit N → ∞ one rescales the boson field with √ N .Hence, we introduce as the new bosonic degree of freedom.
2 Occasionally, these species are referred to as different colors or flavors.
In an even number of spacetime dimensions d + 1 the formation of a nonzero expectation value of σ therefore signals the breaking of the discrete Z 2 chiral symmetry of the GN action and the formation of a condensate.In an odd number of spacetime dimensions as well as for noninteger dimensions the matrix γ ch , which fulfills is either nonexisting or its definition is ambiguous (depending on d γ ).For detailed discussions of these situations, we refer to Refs.[50,51,56].Still, it is possible to study the action (3) and analyze for which µ and T one finds (non)vanishing expectation values of the scalar field ⟨σ⟩.This even holds true if one allows for spatial modulations of this condensate.Regions in the µ, T -plane with ⟨σ⟩ ̸ = 0 are denoted as phases of spatially (in)homogeneous condensation/symmetry breaking, while regions with ⟨σ⟩ = 0 are called gas like/symmetric phases.
The standard way to proceed from Eq. ( 3) is to perform the functional integration over the fermion fields and absorb the resulting fermion determinant in an effective action for the boson field σ.The resulting action S eff in the probability distribution in the thermal partition function is where we already divided by the number of species N and Det denotes a functional determinant.
Considering the limit N → ∞, this implies that the only field configurations that contribute in the partition function and the expectation values are those that minimize the effective action S eff .This is equivalent to studying the full quantum effective action (the generating functional for one-particle-irreducible vertex functions) and only taking into account the contribution of the fermionic quantum fluctuations [57][58][59].In any case, for arbitrary modulations in spatial and temporal directions, the minimization of Eq. ( 8) is a highly challenging task.It might not even be a well-posed problem for noninteger dimensions.However, assuming that the field configuration with least (effective) action is constant in space and time, e.g., σ(τ, ⃗ x ) = σ = const.the problem simplifies drastically.We define the (homogeneous) effective potential which is the effective action for homogeneous fields per species and spacetime volume.The eigenvalues of the Dirac operator for homogeneous fields σ are those of free fermions with mass m = σ.Thus, the evaluation of Ū yields, where l 0 is the Matsubara sum and momentum integral over the log of the eigenvalues.Its evaluation is presented in Appendix C. By determining the global and local minima of Eq. ( 10), we obtain the homogeneous phase diagram including the spinodal lines.We denote the homogeneous field configuration that corresponds to the global minimum of the homogeneous effective potential for a given µ and T by Σ(µ, T ).
The derivative of the homogeneous effective potential with respect to (w.r.t.) the homogeneous field σ is used to express the gap equation, which is of central importance in the renormalization of this model as discussed in Section II E. The quantity l 1 is again a Matsubara sum and spatial momentum integration.Its evaluation is discussed in Appendix D.

B. The bosonic two-point function
The bosonic two-point function at bosonic spatial momentum q = |⃗ q | for a homogeneous bosonic field σ is given by where the quantities l 1 and l 2 are fermionic Matsubara sums and loop momentum integrals that are discussed in Appendices D and E. The derivation of this quantity in the GN model is discussed in great detail in Refs.[1,13,22,25] and shall not be repeated here.Simply speaking, one obtains the bosonic two-point function, by (1.) taking two functional derivatives w.r.t. to σ of Eq. ( 8), (2.) evaluating the result at vanishing external Matsubara frequencies and for σ(τ, ⃗ x ) = σ = const..In our analysis, we evaluate the two-point function at the global minimum of the effective potential Σ(µ, T ) for the specific µ and T .

C. The bosonic wave-function renormalization
Another quantity of interest is the so-called bosonic wave-function renormalization given by which is the coefficient of the kinetic contribution 1 2 (∂ µ σ) 2 to the quantum effective action [4].After a short calculation, which is summarized in Appendix H one finds where we introduced l 3 as another Matsubara sum and integral that is evaluated in Appendix F. If the wavefunction renormalization is evaluated at the global minimum of the effective potential Σ, we denote it by Z, i.e., Z(µ, T, d) ≡ z( Σ(µ, T ), µ, T, d).

D. Regularization of vacuum contributions
Some of the previously listed quantities contain contributions with UV-divergent integrals.For d < 3, these are the vacuum parts of l 0 and l 1 , which require a UV regularization to render calculations tractable.We regularize the theory with a spatial momentum cutoff that confines the integration of spatial fermionic loop momenta to a d-dimensional sphere of radius Λ.This scheme certainly has its drawbacks in the investigation of inhomogeneous condensation as we explicitly break translational invariance.However, as is discussed in Section II E, one can renormalize the theory for d < 3. Thus, it is possible to eventually remove the regulator completely by sending Λ → ∞, allowing us to make the assumptions that Λ is large.
Regularizing the vacuum part of l 0 and l 1 with the spatial momentum cutoff and expanding the result for Λ/σ ≫ 1 results in , where the l Λ x are the regularized versions (C3) and (D3) of the respective quantities, and S d is given by Eq. (B2).These expansions are obtained by applying the formula Eq. (B23), whose origin is discussed in Appendix B 4. In the last step we also used Eq.(B4) in order to have identical prefactors of the first terms in Eqs.(15) and ( 16).

E. Renormalization
The GN model naturally experiences spontaneous symmetry breaking in the vacuum for all d for N → ∞, which is exploited in the renormalization procedure.We impose as renormalization condition that the auxiliary field σ assumes the homogeneous nonzero value σ0 in the vacuum.The coupling λ is tuned such that this condition is fulfilled and the divergences are absorbed.Within the limit of N → ∞, one finds that this is successful for spatial dimensions d < 3, see Refs.[11,19,60].This can easily be understood by looking at Eq. (15), where a divergence ∝ σ4 arises for d ≥ 3.There is no coupling in the action that can compensate this divergence and removing the cutoff is only possible by introducing more parameters.

The gap equation
Our renormalization condition can be expressed by the gap equation in vacuum Inserting Eq. ( 11) and rearranging the equation for nonzero σ0 fixes the value of the coupling (its cutoff dependence) as , where we again assumed that Λ is large, i.e., Λ/σ 0 ≫ 1 and therefore simply applied Eq. ( 16) for σ → σ0 .

Renormalization of the effective potential
Hence, for the effective potential in the presence of the UV cutoff we find which in vacuum evaluates to The last line vanishes for d < 3 in the limit of Λ → ∞, see also Ref. [1].We then obtain the renormalized homogeneous effective potential whose derivation and special limits of the parameters σ, µ, T, d are presented in Appendix G.This result is equivalent to the result presented already in Refs.[11,61].

Renormalization of the two-point function
The unrenormalized two-point function Eq. ( 12) contains the divergent contribution l 1 .By inserting the regularized value for the coupling λ and inserting the regularized expression for l 1 , one obtains ) l 2 (σ, µ, T, q, d) .
In the limit Λ → ∞, one finds that the divergent parts are absorbed by the coupling and we find for the renormalized two-point function the expression The definition of Ẽ is given in Eq. (A15), and detailed steps of the renormalization as well as special limits in the parameters σ, µ, T, q, d are presented in Appendix I.

F. The stability analysis
While the investigation of the homogeneous phase structure can be conducted via the one-dimensional minimization of the homogeneous effective potential Ū w.r.t. the variable σ, one cannot minimize the effective action for an arbitrary inhomogeneous field configuration.Typically, one either resorts to an ansatz for the field modulation, as, e.g., in Refs.[9,14,23,51,62,63], or conduct a stability analysis, which is the approach that we consider.Here, we only summarize the strategy of this technique and the final quantities.We refer to Ref. [25] for a detailed derivation and discussion at the example of the (1 + 1)-dimensional GN model and to Refs.[1,46] for some further details of the stability analysis in noninteger dimensions.Other works that relied on this or related techniques are for example Refs.[13,24,25,29,33,34,36,38,39,44,[64][65][66].
The strategy of this technique is to apply inhomogeneous perturbations to a homogeneous field configuration and expand the effective action in powers of this perturbation.When applying this expansion at the global homogeneous minimum, one finds that the first nonzero correction is the second order term quadratic in the perturbations.A negative sign of its coefficient signals that an inhomogeneous field configuration is energetically favored over the homogeneous expansion point.This coefficient is given by the bosonic two-point function Γ (2) with the external, spatial bosonic momentum corresponding to the momentum of the inhomogeneous perturbation.On a strictly formal level, we are therefore simply searching for spatial momenta q, where the two-point function (23) takes negative values, if it is evaluated at the global minimum of the potential (21).
Closely related to the IP is the so-called moat regime [41,43], which is characterized by a two-point function featuring a minimum at a finite momentum.This is also realized in an IP, where, however, the minimum of the 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -0.05 0.00 0.05 0.10 FIG. 1.The two-point function Γ (2) evaluated at the global homogeneous minimum Σ(µ, T ) at µ/σ0 = 1.2 and d = 1.8 for various temperatures T /σ0 as a function of the bosonic momentum q/σ0.two-point function is necessarily negative.A simple criterion for the detection of the moat regime is that the wave-function renormalization (13) evaluated at the homogeneous minimum of the potential ( 21) assumes a negative value. 3

III. RESULTS
The results as discussed in the following are an extension of the results presented in Ref. [1], which conducted the stability analysis of the (d+1)-dimensional GN model at T = 0. Therefore, we concentrate our presentation on the effects by nonzero T and on the phase diagram as a function of d as a whole.We refer to Ref. [1] for further results, which analyze the d-dependence of the stability analysis.Furthermore, we refer to Ref. [25] for a detailed discussion of the results for d = 1 and to Refs.[13,44] for d = 2.

A. The two-point function
We start the discussion by providing two example plots of the bosonic two-point function as a function of the 3 This criterion assumes that there are no minima in the two-point function at finite momentum that are separated from the origin by a local maximum.This situation would also correspond to a moat regime, but blind to our criterion.Such a situation can indeed occur in non-renormalizable (pseudo)scalar four-fermion models like the 3 + 1-dimensional Nambu-Jona-Lasinio model, when considering chemical potentials which are larger than the regulator [46].However, we are not aware that these models also exhibit such properties in the fully renormalized limit and, thus, the wave-function renormalization appears as an appropriate criterion.We find that the non-analytic points at q = 2µ, which are present at T = 0, see Ref. [1], are completely smoothed out already at rather low temperatures.For d < 2, the former nonanalytic point turns into a nontrivial global minimum for some temperature range, while at d > 2 no structure reminiscent of the cusp at T = 0 remains.In this way for d < 2 the instability towards an IP signaled by a negative Γ (2) at T = 0 (blue curve in Fig. 1) vanishes at some T > 0. A moat regime remains signaled by the negative curvature of Γ (2) at q = 0 (red curve in Fig. 1), which results in a nontrivial global positive minimum.For large temperatures we universally find a convex shape with no remains of the IP or a moat regime.Before we close the discussion, let us remark that the T -order of the curves in Fig. 2 is correct: The offset of the two-point function is the bosonic curvature mass, which vanishes at the second order PT.At µ/σ 0 = 1.4 one starts in the HBP at T /σ 0 = 0 and crosses the PT to the SP at approximately T /σ 0 = 0.6 by increasing T , see also Fig. 7. Of course, we could have plotted the two-point function for d = 2.5 in a region, where the global minimum of the potential is trivial, hence for some µ/σ 0 ≳ 1.5 similar to d = 1.8.This, however, would have been even less interesting and we therefore opted for a comparison of d = 1.8 and d = 2.5 at the same µ/σ 0 .

B. The wave-function renormalization
The second quantity of interest is the wave-function renormalization Z, which serves as the indicator for the moat regime.It is shown as a function of the chemical potential evaluated at the global homogeneous minimum Σ(µ, T ) at various temperatures for d = 1.8 in Fig. 3 and for d = 2.5 in Fig. 4. The two values of d are again representative for the behavior for 1 ≤ d < 2 and 2 < d respectively.The dots on the curves indicate the chemical potential that corresponds to the homogeneous PT at the given temperature.
For d < 2, one finds at T = 0 that Z is constant for small µ/σ 0 and jumps to a negative value at the homogeneous PT.The constant behavior is a consequence of the Silver blaze property as [67][68][69][70], which is no longer fulfilled for nonzero T .The jump vanishes exactly at temperatures, where one does no longer find a first order PT in µ-direction under the assumption of homogeneous condensation, i.e., at the temperature of the critical endpoint, see also Fig. 7. Above this temperature, one finds that the chemical potential beyond which the wavefunction renormalization is negative moves to higher values.This is the behavior that we would expect based on the results for d = 1 as presented in Ref. [25].However, one finds a moat regime for all temperatures for sufficiently large µ/σ 0 .
For d > 2, one finds at T = 0 that Z is constant for small µ and diverges at µ = σ0 , which is again a manifestation of Silver Blaze.This nonanalytic behavior is smoothened already for small temperatures and one finds that Z smoothly changes with µ.However, most importantly for the present investigation is, that the wave-function renormalization is always positive for arbitrary values of T and µ.This implies the total ab- sence of an IP and also the total absence of a moat regime for d > 2.
For the sake of clearness, we also provide two density plots of the wave-function renormalization in the µ, Tplane for d = 1.8 and d = 2.5 in Figs. 5 and 6. (The Figs. 3 and 4 are just sections of these density plots at constant T .)Here, it is again clearly visible that there is no moat regime and IP for d > 2, where Z is always positive.For d = 1.8, however, we find a similar structure as for d = 1 (see Ref. [25]).The moat regime is present in the region of the IP and SP below the straight line that originates from µ = T = 0 and passes through the critical point.This straight line is associated with a vanishing quartic coefficient of the effective potential.Hence, the coefficient in front of σ4 changes its sign along this line as a function of µ and T , if one expands Eq. ( 21) in σ, see Refs.[12,71,72].It can be shown [71] that this coefficient is proportional to the bosonic wave-function renormalization (14), if the wave-function renormalization z is evaluated at the trivial point σ/σ 0 = 0.This is the correct evaluation point only in the SP, which is however the important region for the moat regime.Hence, the straight line separates the regime of negative and positive Z in the SP.finds a first order PT between HBP and SP at low temperatures and a second order PT for large temperatures, if one assumes spatially homogeneous condensation in the entire µ, T -plane.On the other hand, for d ≥ 2 the PT between the HBP and SP is always of second order (except for d = 2 and T = 0 [17,18]).This result was already found in Ref. [11], where it was also observed that the critical point moves down to T = 0 and that the HBP enlarges, while going continuously from d = 1 to d = 2.However, when allowing for the formation of spatially inhomogeneous condensates and searching for these via the stability analysis, this picture is modified for 1 ≤ d < 2, while nothing happens at d ≥ 2. Already at T = 0 it was found in Ref. [1] that the stability analysis reveals an IP for 1 ≤ d < 2. This phase extends to µ/σ 0 = ∞ for d = 1 at T = 0 but shrinks and has a second order PT to the SP at some finite critical µ when 1 < d < 2. The PT between the HBP and the IP cannot be resolved correctly with this method as is discussed in detail in Refs.[25,38].Still, for d = 1 the analytic solution is well-known [9,73] and served as a test field for the stability analysis in Ref. [25].For a direct comparison, of the situation in 1 ≤ d < 2 and d > 2, we again used d = 1.8 and d = 2.5 and plotted both phase diagrams together in Fig. 7.For reference, we also included the spinodal lines (the lines that engulf the region in the phase diagram, where the effective potential has a global and local minimum/minima), plotted in blue.
However, to really observe the effect of dimensionality on the IP and the moat regime, we prepared Figs. 8 and 9, which clearly show that the phase diagram for 1 < d < 2 is similar to the phase diagram at d = 1 (except for the finite extent of the IP at T = 0) and that d = 2 is the strict upper bound for the existence of an IP.Still, it is remarkable that the IP and the moat regime vanish rather slowly as a function of d and an IP is still clearly visible for our last plotted curve at d = 1.95.Nonetheless, this behavior was actually expected, because the same slow convergence was already observed for the critical point in Ref. [11], which turns into the Lifshitz point, if one allows for spatially inhomogeneous condensation.Let us remark at this point that we did not prepare extra plots for the situation at 2 < d < 3, since we do not find any signal of spatially inhomogeneous condensation and/or a moat regime for any T and µ (this was already observed at T = 0 in Ref. [1]).Therefore, the situation of exclusively spatially homogeneous condensation is already fully covered by Ref. [11] with which our results agree.We believe that the example of d = 2.5 is absolutely sufficient within this work to understand the situation for 2 < d < 3.
We close this section by summarizing that we find moat regimes and phases, where the condensate varies in space, in the GN model for N → ∞ for 1 ≤ d < 2, while they are absent for d ≥ 2.

IV. CONCLUSIONS AND OUTLOOK
Finally, we want to summarize our results, draw several conclusions, and provide an outlook to possible follow-up projects.

A. Summary
In the present work, we investigated the phase diagram of the GN model at (non)zero fermion chemical potential µ and (non)zero temperature T in the limit of an infinite number of fermion species, N → ∞, as a function of the spatial dimension d.We focused on 1 ≤ d < 3, where the model is renormalizable and solely depends on the fixation of a single dimensionful parameter.We used the vacuum fermion mass to fix the scales and worked in the renormalized limit.The focus of this work is the detection and the dependence on the number of spatial dimensions of the IPs and a possible moat regime.We found that the well-known result for d = 1 [9,10,73,74] generalizes to 1 < d < 2 and one detects an instability of the SP that signals the presence of an IP.Furthermore, one always finds an even larger moat regime at large chemical potential.Both, the IP and the moat regime vanish, when d approaches d = 2.For d ≥ 2 we do not find any indication of an IP in terms of an instability.The same applies to the moat regime.While, we cannot exclude any kind of inhomogeneous condensation, which is not detectable via a stability analysis, such a situation is, however, highly unlikely to be present in this model.The reason is that the PT between the SP and an IP is generally expected to be of second order, which enables the use of this method (see Ref. [25] for a detailed discussion about the range of validity of this method).
Apart from these novel results, we implicitly and explicitly confirmed various existing literature results for the GN model for integer d = 1 [4,8,9,23] and d = 2 [13,14,17,18,51] as well as some of the results from Ref. [11] for continuous d.
We also provide several appendices with detailed material that may be of general use for follow-up or related projects.

B. Conclusion
Already in Ref. [1] we speculated about the relevance of the number of spatial dimensions d on the formation of spatially inhomogeneous ground states in the GN and other models.Furthermore, Refs.[13,14,45] showed that the effects of the presence of a finite regulator or an effective UV/IR cutoff in terms of a spatial lattice or a finite spatial box play an important role for the presence/absence of spatially inhomogeneous condensation.Here, we want to continue this discussion and believe that the present investigation sheds light on the previous findings.Let us therefore briefly summarize the findings for the GN model at N → ∞ up to this point: In d = 1 there is an exact solution for the phase diagram and one finds spatially inhomogeneous condensation [9,73].This feature seems to be robust at N → ∞ even in the presence of a spatial lattice etc. [28,75].For d = 2 and for N → ∞ there is also an exact solution for the phase diagram, but there are no IPs in the renormalized limit [13,14].However, in the presence of some UV or IR regulator/cutoff one recovers an IP and a phase diagram that has some similarities with the situation in d = 1 [13,14].The size of the IP and the shape of the HBP thereby strongly depends on the value of the regulator/cutoff and one finds results that are closer to d = 1 or d = 2 depending on the strength of the regularization.In d = 3 there are several models that are similar to the GN model at N → ∞ which support the presence of spatially inhomogeneous condensation, while the extent of the IP usually depends on the cutoff/model parameters (renormalization like in this work is not possible and the phase diagram usually depends on at least two parameters).If we compare these results to our findings, we come to the conclusion that spatially inhomogeneous condensation seems to be a dimensional effect.For the studies discussed previously and for the present study the situation is always the same: As soon as the (effective) dimensionality of the model is reduced to 1 + ∆d dimensions, where ∆d ∈ [0, 1), we find an IP.However, as soon as there are two or more fullfeatured spatial dimensions available, the IP vanishes.In fact, it does not play a role if the number of spatial dimensions is directly reduced via dimensional regularization or as worked out in the present study using d as a continuous parameter or even via an UV/IR cutoff as, e.g., in Refs.[13,14,24,45,46,63,76].In particular, the latter case can be viewed as a reduction of a full dimension to a fractional/part of a dimension by restricting the system to a finite spatial box (IR cutoff) or coarse lattice (UV cutoff).A similar effect is observed in a recent work Ref. [77], which analyzed the homogeneous phase structure of the 2 + 1-dimensional GN model in a finite volume.The finite volume causes the critical endpoint to be located at a finite temperature, which is a feature that is limited to d < 2 in the infinite volume.
These results also seem to be in line with the observation that one-dimensional ansatz functions for inhomogeneous condensates usually appear to be the most promising and energetically favored solutions, if an IP is present at all [36,51,63,78].
Of course, at this point the immediate question that arises is: What is the underlying nature and physical principle behind this strong relation to a single spatial dimension?So far, we were not able to come up with a conclusive answer, but we hope that this work might be an important step to start the search in the right direction.The Peirls instability, which is the origin of the IP in d = 1 [9,21] and a one-dimensional effect that cannot be directly generalized to higher dimensional systems, might be a good starting point in terms of a physics understanding and correct interpretation.

C. Outlook
Now, that we have mostly settled the situation for the GN model at N → ∞ there are basically four main directions to proceed.

One-dimensional ansatz functions
It was found in d = 1, that the stability analysis is not able to detect the portion of the IP, where it is energetically favored, but the homogeneous expansion point Σ is finite [25].As mentioned in Section III C, we expect the same situation to occur for 1 < d < 2 in our calculations with some part of the IP in the vicinity of the first order PT between HBP and SP to be missing.A way to improve on this would be to consider a one-dimensional ansatz function embedded in the d-dimensional space.Such a procedure was considered in (3 + 1)-dimensional models in Ref. [63], which treats the perpendicular space in such a general way that it can be generalized from d ⊥ = 2 to noninteger d ⊥ = d − 1 dimensions.It is expected that the additional portion of the IP is not particularly large, because it quite limited in size in d = 1 and likely shrinks even further with increasing d.Nevertheless, this step would yield the complete phase diagram of this model. 4

Finite regulator or finite volume
To solidify the concept of effective dimensionality, it would be fascinating to carry out the present investigation not in the renormalized limit, but at a finite UV regulator.In this way, one could (a) connect smoothly to the 3+1 model results by extending our analysis to d = 3 and (b) one could investigate the interplay between the explicit number of spatial dimension d and the effective dimensional reduction introduced by the UV regulator.The latter could also be investigated by considering a finite volume, which introduces an IR regularization that should also lead to an effective dimensional reduction.

Finite N
While the N → ∞ limit was essential to investigate the analytic structure of the GN model for noninteger d, this semi-classical limit is fairly different from the behavior that we would expect of a QFT.The general observation is that bosonic quantum fluctuations as they would occur for finite N weaken ordered phases in such models (see, e.g.Refs.[29,72,[79][80][81]) and therefore likely do not enable the emergence of an IP for d ≥ 2.
For d < 2, it is highly likely that the IP vanishes altogether.The most pathological aspect of the N → ∞ limit is the fact that it circumvents the Coleman-Mermin-Wagner-Hohenberg-Berezinskii theorem [28,[82][83][84][85] (or related arguments for discrete symmetries [4,22,[86][87][88]) in d = 1, which would normally forbid any type of condensation at a finite temperature.Accordingly, it was found in Refs.[22,72], that there is no symmetry breaking at finite T and N in d = 1.These effects likely suppress any IP for d < 2, which all in all suggests that these models do not exhibit an IP in any dimensions for finite N .
Nevertheless, it might be interesting to consider this model for finite N .While there is no condensation in d = 1 at finite T , one finds an HBP in d = 2 [80].Thus, considering the GN model for noninteger d might be instructive to observe how the theory evolves from a system without any symmetry breaking to the system with a broken symmetry.A functional method that admits the formulation of the theory for an arbitrary d such as the Functional Renormalization Group [89][90][91] might be the optimal framework for extending our analysis to finite N .

Consequences for higher dimensional models and QCD
Our analysis shows that four-fermion models with scalar interaction channels in the limit of N → ∞ only exhibit an IP for d < 2 or via an effective reduction of the dimensionality of the system.This suggests that an IP might exist in QCD only when the low-energy behavior of QCD is not only described by scalar-pseudoscalar four-fermion interactions, but by interactions that would admit an IP for higher dimensions in other than the scalar and pseudoscalar channel.This is most likely the case at finite chemical potential, where it was found that the relevant interaction channels are diquark interactions [92,93], which have not been systematically studied with respect to the IP.Moreover, in this regime vector interactions become relevant, which were found to mix with scalar modes at finite densities and to play an important role in the homogeneous phase transition near the critical endpoint in QCD [94].Such a mixing might even induce an instability that results in a spatial modulation of the condensates [95,96].Another mechanism that could cause the existence of an IP in QCD is that it has in fact a lower effective dimension, e.g., caused by additional strong magnetic fields.While these are aspects that are yet to be understood, they certainly imply important questions that should be answered in an effort to investigate IPs and the moat regime in QCD.However, there are also indications that a possible IP would be completely destabilized by the Goldstone bosons from chiral symmetry breaking [97]. A. K. and L. P. also like to thank M. J. Steil, because the adaptive mesh-refinement algorithm that was used to generate the data of the phase diagrams is based on his work and code.
A. K. and L. P. acknowledge support from the Helmholtz Graduate School for Hadron and Ion Research, the Giersch Foundation and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Collaborative Research Center Tran-sRegio CRC-TR 211 "Strong-interaction matter under extreme conditions" -project number 315477589 -TRR 211.
All numeric results as well as the figures in this work were obtained and designed using Mathematica [98].In this work, we use the following conventions for Fourier transformations.For the bosonic field we have while the fermion fields are Fourier-transformed according to Hereby, the corresponding Matsubara frequencies for the discretized energies are which stem from the (anti-)periodic boundary conditions in τ -direction at τ = 1 T for (fermions) bosons.

b. Fermi-Dirac distribution function
We define the Fermi-Dirac distribution function as follows [99,100] Especially for the numeric implementation we exclusively use the representation in terms of tanh.In addition, we present two useful identities, which are also part of the derivation of the explicit analytic expressions in this work as well as the numeric implementation, . (A10) "Primes" denote derivatives w.r.t.x.
For the derivation of the zero-temperature limits of some formulae of this work, we repeatedly need the following limits.lim Here, E is the energy, µ the chemical potential, and Θ is the Heaviside function.In addition, lim where δ is the Dirac-delta distribution.

Spherical symmetric integration
Most of the momentum integrals in this work are of the hyperspherical type.The integrand is usually only a function of the absolute value of the momentum such that the angular integration can be performed, Here, we introduced which is the surface of the d-dimensional sphere.

Transcendental functions
A lot of the explicit formulae for the effective potential, bosonic wave-function renormalization, and the bosonic two-point function can be expressed in terms of known functions.For the sake of a self-contained presentation, we provide these functions with links to references for further reading in this appendix.We hope that this reduces unnecessary look ups and literature searches to a minimum for the reader.

a. Gamma functions
The gamma function is given in terms of its integral representation by the following expression [101, Eq. 6.1.1], (B3) It fulfills the defining relation, For this work, we make use of the Laurent series representation [102, Eq. 5.7.1] which we use for an expansion about z = 0.The polygamma function with integer index is defined in terms of derivatives of the conventional gamma function.However, there is also an integral representation [101, Eq. 6.4.1].

b. Riemann zeta function
Some formulae of this work can be expressed in terms of the Riemann zeta function [101, Eq. 23.2.7], We also define the Dirichlet eta function via the Riemann zeta function [101,Eq. 23.2.19] and in terms of an integral, It is well-known that some integrals over Bose-Einstein or Fermi-Dirac distribution functions can be expressed in terms of (incomplete) polylogarithms.The polylogarithm is defined via the following integral [102, Eq. 25.12.10], Li which reduces to the Dirichlet eta function (B8) for z = −1, A special definition, which is used in this work is the symmetrized derivative of the polylogarithm w.r.t.its index s, Here, γ is the Euler-Mascheroni constant, ψ (n) (s) is the polygamma function (B6), and the last equality holds for 2n ≤ 0. It turned out that using the last relation is more stable and accurate when it comes to numeric evaluation [25,72].
A particular useful (linear) transformation formula is [101, Eq. 15.3.7] and is valid for |arg(−z)| < π.It can be used to expand the hypergeometric function for large |z|.

Integrals and an expansion
Next, we present two important integrals for this work as well as an expansion that is used several times.

a. First special integral
Repeatedly, we are confronted with integrals that are of the type We used the substitution as well as Eqs.(B4) and (B13).

b. Second special integral
Another integral that appears several times during our calculations is We used that E 2 = p 2 + ∆ 2 and defined The integral was evaluated using where x i are the roots of f (x).For µ 2 > ∆ 2 , thus µ > 0, is the zero of while Evaluating p a , 1 E n , and f ′ at p 0 one obtains the above results.

Expansion
At several points in this work, e.g., for sending the UV cutoff Λ to infinity or to evaluate certain expressions at the symmetric point σ → 0, we need an asymptotic expansion formula for the hypergeometric function 2 F 1 (α, β; γ; z) for large |z|.This formula is found by inserting the series representation (B12) in the linear transformation formula (B14).In particular, we find, .
Appendix C: Evaluation of l0(σ, µ, T, d) In this appendix we provide details on the l 0 (σ, µ, T, d)-Matsubara sum and integral which occurs in the expression for the effective potential (10).It is defined as follows Here, ν n denotes the fermionic Matsubara frequencies (A7) and E is the fermion energy (A13).We used contour integration to evaluate the Matsubara sum.The infinite constant term can be ignored in what follows.It corresponds to an arbitrary normalization of the effective potential.
1.For T = 0 The zero temperature limit of Eq. (C1) where we used Eq.(B1) to simplify the momentum integration with hyperspherical coordinates.Regularizing the UV divergence with a sharp UV cutoff and splitting the integral in µ-(in)-dependent parts one obtains with Eq. (B15), Here, we used the Def.(A14) to facilitate a compact notation.

For T ̸ = 0
In general, for nonzero T we can still evaluate the vacuum contribution in Eq. (C1) and find with the (µ = 0)part of Eq. (C3), Further evaluation of Eqs.(C3) and (C4) is performed after renormalization of the effective potential in Appendix G.
Appendix D: Evaluation of l1(σ, µ, T, d) In this appendix, we present explicit expressions for the l 1 (σ, µ, T, d)-Matsubara sum and integral which is part of the gap equation ( 11), the regularized effective potential (19), and the (regularized) bosonic two-point function Eqs. ( 12) and (22).It is defined and evaluated as follows Again, ν n are the fermionic Matsubara frequencies (A7) and E is the fermion energy (A13).Additionally, we introduced the Fermi-Dirac distribution (A8), while it turned out that the tanh-representation seems to be more stable and accurate for numeric computations.
1.For T = 0 In the zero-temperature limit Eq. (D1) reduces to where we made use of hyperspherical coordinates in momentum space, see Eq. (B1).Splitting µ-(in)-dependent terms, using the abbreviation (A14), and introducing the UV cutoff Λ one can make use of Eq. (B15) to arrive at .
2. For T ̸ = 0 Again, for T ̸ = 0 we solely evaluate the vacuum contribution of Eq. (D1) and again use the previous result Eq. (D3) for µ = 0.For the regularized expression one finds, Further evaluation of Eqs.(D3) and (D4) is postponed to Appendices G to I.
Appendix E: Evaluation of l2(σ, µ, T, q, d) In the expression for the bosonic two-point function ( 12) contains a ⃗ q-dependent part.Here, we show some simplifications for this contribution that reads We already used q = |⃗ q | as an argument of l 2 instead of ⃗ q.This becomes clear in the following lines.In order to get rid of the nasty vectorial ⃗ q-shift in the second propagator we use the following Feynman parameter integral (E2) Applying this to Eq. (E1) one obtains, We exchanged the order of integration, substituted ⃗ p ′ = ⃗ p + ⃗ q x, and immediately returned to the "unprimed" notation for ⃗ p.Using the Defs.(A8), (A9), (A15), (A16), and (A18), Thus, the dependence on ⃗ q is actually a dependence on its absolute value q.

For: T ̸ = 0
Of course, we also provide a simplification of Eq. (E1) for T ̸ = 0 by using the vacuum contribution of the previous result (E6), Appendix F: Evaluation of l3(σ, µ, T, d) In complete analogy to the previous appendices, we present an appendix for the partial evaluation of l 3 (σ, µ, T, d) that is part of Eq. ( 14) for the bosonic wavefunction renormalization.In terms of the (un)evaluated Matsubara sum and momentum integral it reads In this expression we used Eqs.(A7) to (A10) and (A13).

For T = 0
In the zero-temperature limit Eq. (F1) turns into where we already made use of hyperspherical coordinates via Eq.(B1).Again, it is not necessary though still useful to introduce a UV cutoff regulator Λ to use Eq.(B15) also for the vacuum contribution.Afterwards, the cutoff can be removed with the help of Eq. (B23).Splitting and evaluating the integrals in µ-(in)-dependent contributions, one finds, using integration by parts and Eqs.(B15) and (B17), Again, we used the compact notation Eq. (A14).

For T ̸ = 0
At nonzero temperature we can use the vacuum part of Eq. (F3) to simplify Eq. (F1), Appendix G: The effective potential In this appendix we turn to the detailed evaluation of the effective potential Eq. (10).To this end, we calculate the renormalized limits of Eq. ( 19), where we remove the UV cutoff by sending Λ → ∞.Explicit expressions for σ, µ, T being zero or nonzero as well as the limiting cases of d = 1 and d = 2 are provided.For the sake of the conciseness, we collected links to every special case in Table I and subdivide the appendix according to the table.
TABLE I. Quick links to the equations for explicit evaluation of the effective potential U (σ, µ, T, d).The formulae are simplified in terms of known functions as far as possible.
Considering σ ̸ = 0 there are two cases to be distinguished.
b. T ̸ = 0, σ = 0, µ = 0 It is straight forward to also set µ = 0 in the previous formulae.For general 1 ≤ d < 3 we find Next, we turn to the special cases, where T = 0.
b. T = 0, σ = 0, µ = 0 The trivial and last case is Ū (0, 0, 0, d) = 0 .(G24) (Certainly, one could always add an arbitrary constant to the potential without changing the physical observables.) Appendix H: The bosonic wave-function renormalization This appendix is dedicated to calculations as well as the presentation of detailed expressions and limiting cases for the bosonic wave-function renormalization (14).We calculate the renormalized limits, such that the final results do not contain any UV cutoff.
Step by step, we provide expressions for the cases where σ, µ, and T are zero or nonzero.Additionally, we evaluate the wavefunction renormalization for the special cases d = 1 and d = 2 and demonstrate that we reproduce known literature results.All cases are collected in Table II, which links to the explicit formulae.
However, we start by providing some useful intermediate steps for the derivation of the general formula for the bosonic wave-function renormalization Eq. ( 14).The starting point is We note that the bosonic two-point function solely depends on the absolute/square of the spatial external momentum, we can use where we defined the Matsubara sum and integral formula (F1).

T ̸ = 0
We start with the wave-function renormalization in the heat bath with T ̸ = 0. First, we study the wave-function renormalization for nontrivial background field configurations σ ̸ = 0, e.g., in the phase of symmetry breaking.
a. T ̸ = 0, σ = 0, µ ̸ = 0 We start at nonzero chemical potential.Both the vacuum and the medium contribution separately exhibit an IR divergence, which cancel each other.To account for this, we consider the vacuum part in its integral form together with the medium part and write The tricky evaluation of this expression for d = 1 is presented in Ref. [22,Eq. F.65] and one finds where we used the definition (B11).While for d = 2 integration by parts leads to b. T ̸ = 0, σ = 0, µ = 0 At vanishing chemical potential Eq. (H11) simplifies to Having discussed all cases with a heat bath, we can next turn to T = 0.
a. T = 0, σ ̸ = 0 Again, we start in the phase with a nontrivial expectation value of the bosonic field and therefore evaluate the wave-function renormalization at σ ̸ = 0.
a. T = 0, σ ̸ = 0, µ ̸ = 0 Keeping µ ̸ = 0 the general expression for continuous d reads where we used Eq.(E6) for q = 0 and Eq.(F3).For d = 1 this simplifies drastically, see also Ref. b. T = 0, σ = 0, µ = 0 Lastly, if one evaluates the wave-function renormalization in vacuum at the trivial evaluation point it is ill conditioned, when taking the respective limits from the previous results.
Appendix I: The bosonic two-point function In this appendix we calculate the bosonic two-point function (12).We use the regularized integrals (D4) and (E7) as well as (F1) and calculate the renormalized limits.Here, we present results for q ̸ = 0 and the limit q = 0.For the sake of clearness, we prepared Table III which links to the different cases with (non-)vanishing σ, µ, T , as well as the special cases with d = 1 and d = 2.The limiting cases for d = 1 are discussed in detail in Refs.[22,25], whereas the d = 2 formulae are briefly discussed in Refs.[13,44].

T = 0
Having completed the T ̸ = 0 cases, we turn to the zero-temperature limit.

FIG. 3 .
FIG. 3. The wave-function renormalization Z evaluated at the minimum of the potential for d = 1.8 for various temperatures T as a function of the chemical potential µ.The dots mark the chemical potential, where the spatially homogeneous global minimum of the potential turns (non)trivial, see Fig. 7.

FIG. 4 .
FIG. 4. The wave-function renormalization Z evaluated at the minimum of the potential for d = 2.5 for various temperatures T as a function of the chemical potential µ.The (T = 0)-curve only stays finite due to finite computational µ-resolution.The dots mark the chemical potential, where the spatially homogeneous global minimum of the potential turns (non)trivial, see Fig.7.
C. The phase diagramNext, we turn to the full phase diagram of the GN model as a function of the spatial dimension d.Indeed, already in Figs.5 and 6we basically plotted the phase structure for d = 1.8 and d = 2.5, which served as examples for 1 ≤ d < 2 and d > 2. For 1 ≤ d < 2 one always

FIG. 5 .
FIG. 5.The wave-function renormalization Z at d = 1.8 in the (µ, T )-plane.The yellow curve corresponds to the second order HBP-SP PT, the red curve to the first order HBP-SP PT, the green curve to the second order IP-SP PT, and the orange curve separates regions of positive and negative wavefunction renormalization.

FIG. 7 .FIG. 8 .
FIG. 7. The phase diagram at d = 1.8 (colored lines) and d = 2.5 (black line) in the (µ, T )-plane.The plot shows for d = 1.8: the second order PT (HBP ↔ SP) in yellow, first order phase PT (HBP ↔ SP) red, the second order PT (IP ↔ SP) in green, the spinodal lines in blue, the line of vanishing quartic coefficient of the potential (= vanishing z(σ = 0)) in orange.The black line is the second order PT (HBP ↔ SP) for d = 2.5.Black dots mark the endpoints of the lines.

FIG. 9 .
FIG. 9.The first and second order phase boundary of the HBP-SP PT for a translational invariant bosonic field and the phase boundary between the IP and SP in the (µ, T, d)space.Different colored lines correspond to different values of d.(The HBP-IP PT is not plotted and not detectable within our approach.) ACKNOWLEDGMENTSA.K. and L. P. thank J. Braun, H. Gies, G. Markó, R. D. Pisarski, D. H. Rischke, M. J. Steil, J. Stoll, M. Wagner, M. Winstel, A. Wipf, N. Zorbach for fruitful discussions about this work.A. K. and L. P. thank R. Pisarski, A. Wipf, M. Winstel and N. Zorbach for useful comments on the manuscript.A. K. and L. P. especially thank S. Floerchinger and G. Endrődi for valuable discussions and for their general support at the TPI in Jena and the faculty of physics at the University of Bielefeld, respectively.