Inhomogeneous phases in the Gross-Neveu model in 1+1 dimensions at finite number of flavors

We explore the thermodynamics of the 1+1-dimensional Gross-Neveu (GN) model at finite number of fermion flavors $N_f$, finite temperature and finite chemical potential using lattice field theory. In the limit $N_f \rightarrow \infty$ the model has been solved analytically in the continuum. In this limit three phases exist: a massive phase, in which a homogeneous chiral condensate breaks chiral symmetry spontaneously, a massless symmetric phase with vanishing condensate and most interestingly an inhomogeneous phase with a condensate, which oscillates in the spatial direction. In the present work we use chiral lattice fermions (naive fermions and SLAC fermions) to simulate the GN model with 2, 8 and 16 flavors. The results obtained with both discretizations are in agreement. Similarly as for $N_f \rightarrow \infty$ we find three distinct regimes in the phase diagram, characterized by a qualitatively different behavior of the two-point function of the condensate field. For $N_f = 8$ we map out the phase diagram in detail and obtain an inhomogeneous region smaller as in the limit $N_f \rightarrow \infty$, where quantum fluctuations are suppressed. We also comment on the existence or absence of Goldstone bosons related to the breaking of translation invariance in 1+1 dimensions.


Introduction
The GN model describes Dirac fermions with N f flavors interacting via quartic interactions in 1 + 1 dimensions. It was originally introduced as a toy model that shares several fundamental features with QCD [1]: it is renormalizable, asymptotically free, exhibits dynamical symmetry breaking of the Z 2 chiral symmetry, and has a large N f limit that behaves like the 't Hooft large N c limit of QCD. The particle spectrum and thermodynamics of the theory in the N f → ∞ limit is known analytically. Similarly, the 1-flavor model is equivalent to the 1-flavor Thirring model which can be solved analytically in the massless limit [2] (it has a vanishing β-function). But for intermediate numbers of flavors 1 < N f < ∞ there is -despite many analytical and numerical studies -no complete understanding of the thermodynamics and particle spectrum.
The GN model and related four-Fermi theories in 1 + 1 dimensions have been used in particle physics, condensed matter physics and quantum information theory. For example, in condensed matter physics the GN model describes the charge-soliton-conducting to metallic phase transition in polyacetylene (CH) x as a function of a doping parameter [3]. It is equivalent to the Takayama-Lin-Liu-Maki model [4] which describes the electron-phonon interactions in CH in an effective low-energy continuum description, see Ref. [5]. Four-Fermi models are intensively studied to better understand and classify symmetry-protected topological phases of strongly interacting systems. For more details we refer to the nice summary in Ref. [6].
Recently we have seen a renewed interest in the physics of the GN model at low temperature and high baryon density, because it is the region of the QCD phase diagram which is particularly challenging for first-principles QCD approaches. Corresponding results are only available at asymptotically high densities, where the QCD coupling constant is small so that perturbation theory can be applied, at vanishing density, where lattice QCD does not suffer from the sign problem, or for unphysically large quark masses, where effective theories exist that mitigate the sign problem. At moderate densities and realistic values of the quark masses, i.e. the regime, which is probed by heavy-ion experiments and which is relevant for supernovae and compact stars, neither approach can be applied. In this regime our current picture of the QCD phase diagram is, thus, mostly based on QCD-inspired models, e.g. the GN model, the Nambu-Jona-Lasinio (NJL) model or the quark-meson model. In early calculations within these models it was assumed that the chiral condensate is homogeneous, i.e. constant with respect to the spatial coordinate(s). However, allowing for spatially varying condenates it turned out that there exist regions in the phase diagram, where inhomogeneous chiral phases are favored [7,8]. The majority of the existing calculations have been performed in the limit N f → ∞ or, equivalently, the meanfield approximation (see Ref. [9] for a review and Refs. [10][11][12][13][14][15][16][17] for examples of recent work). Using lattice field theory and related numerical methods and considering N f → ∞, the GN model has been explored in 1 + 1 and 2 + 1 dimensions, the chiral GN model in 1+1 dimensions and the NJL model in 1 + 1 and 3 + 1 dimension [18][19][20][21][22]. However, a full lattice simulation and investigation of the phase diagram of any of these models at finite number of fermion flavors, where quantum fluctuations are taken into account, is still missing. The main goal of the present work is to make a step in this direction and to explore, whether such inhomogeneous phases also exist in the 1 + 1-dimensional GN model at finite N f .
In this work we shall use naive fermions and SLAC fermions to study the multi-flavor GN model. These fermion discretizations are all chiral and no fine tuning is required to end up with a chirally symmetric continuum limit. But the theorem of Nielsen and Ninomyia [23] tells us that we have to pay a price for using (strictly) chiral fermions. And indeed, with naive fermions we can only simulate 4, 8, . . . flavors. With SLAC fermions we can simulate 1, 2, . . . flavors, but the associated Dirac operator is non-local. It has been argued elsewhere that there is no problem with SLAC fermions for lattice systems without local symmetries, see for example Ref. [24]. We shall consider GN models with 2, 8 and 16 flavors which have no sign-problem. The results for 8 and 16 flavors can be compared with the results obtained with naive fermions. We find full agreement of the results obtained with both fermion species. Note that with Wilson fermions the full chiral symmetry cannot be restored in the continuum limit with just one bare coupling. One needs to introduce a bare mass plus two bare couplings and fine tune these three parameters to arrive at a chirally symmetric continuum limit [25]. An alternative would be to use fermions which obey the Ginsparg-Wilson relation. We did not use such fermions, because we sample the full µ-T parameter space and carefully check for discretization and finite size effects. With Ginsparg-Wilson fermions this would be would be too time-consuming. This paper is organized as follows. In section 2 we summarize some known features of the GN model that are relevant for its thermodynamical properties. These include properties of the fermion determinant in the continuum and on the lattice, homogeneous and inhomogeneous phases in the N f → ∞ limit and some comments concerning the spontaneous symmetry breaking (SSB) of translation invariance. In section 3 we discuss different lattice discretizations, the scale setting and some details of the simulations. Our numerical results are presented in section 4. The main focus concerns the behavior of the two-point function of the order parameter for chiral symmetry breaking and the resulting consequences for the phase diagram in the plane spanned by the chemical potential µ and the temperature T . We shall see that the GN model with N f = 2, 8 and 16 flavors behaves qualitatively similar to the model in the N f → limit. In particular we localize three regions in µ-T parameter space, where the two-point function shows a qualitatively different dependence on the spatial separation. The model with N f = 2 is also simulated on rather large lattices with spatial extent N s up to 725 lattice points to carefully investigate the long-range behavior of the correlator. In appendix A we discuss, why the lattice GN model with naive fermions may have an incorrect continuum limit, and how to modify the interaction term to end up with an (almost) naive fermion discretization with correct continuum limit.

Theoretical basics 2.1 The Gross-Neveu model
The Gross-Neveu model (GN model) is a relativistic quantum field theory describing N f flavors of Dirac fermions with a four fermion interaction. In this work we investigate this asymptotically free model in 2 spacetime dimensions. The fermions are described by a field ψ = (ψ 1 , . . . , ψ N f ), the components of which are two-component Dirac spinors. Originally it has been studied in the 1/N f expansion. The action and partition function are where the fermion bilinears contain sums over flavor indices, e.g.ψψ = iψ i ψ i .
To be able to perform the fermion integration one follows Hubbard and Stratonovich by introducing a fluctuating auxiliary scalar field σ to linearize the operatorψψ in the interaction term, where is the Dirac operator. The four-Fermi term in (1) is recovered after eliminating σ by its equation of motion or equivalently by integrating over σ in the functional integral. In eqs. (2) and (3) we also introduced a chemical potential µ to study the system at finite fermion density. Expectation values of operators O(ψ,ψ, σ) in the grand canonical ensemble are given by Note that the integration is over fermion fields, which are anti-periodic in the Euclidean time direction, with period β = 1/T , while the auxiliary scalar field is periodic.
Integrating over the fermion fields leads to with expectation values of operators O(σ) given by Of particular interest in the present work is the chiral condensate, which distinguishes the different phases of the GN model. Translation invariance of the integral over dσ x in the (well-defined) functional integral on the lattice implies a Ward-identity which states, that the condensate is proportional to the average auxiliary field, Also of interest is the Ward-identity relating the two-point function of the condensate to the two-point function of the auxiliary field, In our analysis of the phase diagram at finite temperature and density, the two-point function of the auxiliary field on the right hand side will play a crucial role (see section 4.3).

The fermion determinant
In this section we study some relevant spectral properties of the Euclidean Dirac operator D with auxiliary field and chemical potential as defined in eq. (3).
Later on we shall discretize Euclidean spacetime on a lattice such that ∂ µ turns into a difference operator. For most discretizations one does not retain the above properties without introducing doublers -this is what the celebrated Nielsen-Ninomiya theorem tells us [23]. In the present work, however, we shall use chiral lattice fermions with the above properties, naive fermions (having doublers) and SLAC fermions. Now we investigate the spectral properties of the neither hermitian nor anti-hermitian operator D with eigenvalue equation in the continuum or on the lattice under the assumption that a) and b) hold true.
Charge conjugation: In 2 Euclidean spacetime dimensions there exists a symmetric charge conjugation matrix C with C −1 γ µ C = γ µT . Since the (Euclidean) γ-matrices are hermitian we have γ µ * = γ µT , and property b) implies It follows that all non-real eigenvalues come in complex conjugated pairs (λ, λ * ) such that det D is real. Hence there is no sign problem for an even number of flavors, since then (det D) N f is non-negative.
Chiral symmetry: The four-Fermi term breaks the U A (N f ) chiral symmetry of the kinetic term. But a discrete discrete Z 2 chiral symmetry still remains under whichψψ and σ change their signs. Under this discrete chiral symmetry the Dirac operator is conjugated with γ * = iγ 0 γ 1 : Since (on a finite lattice) the number of eigenvalues of D is even, we conclude that the determinant is an even function of the auxiliary field, Hermitean conjugation: The Dirac operator in eq. (9) is the sum of one anti-hermitean and two hermitean terms and Since the determinant is real (and the number of eigenvalues is even) it follows that det D is invariant under a simultaneous sign change of σ and µ, Together with (12) this leads to a determinant which is an even function of the chemical potential, Note that for most lattice fermions not all of the above properties hold true, for example for Wilson fermions the Z 2 chiral symmetry is explicitly broken.

Summary of existing results in large -N f limit
Before 2003 most field-theoreticians and particle physicists took for granted that in thermal equilibrium translation invariance is realized such that the chiral condensate ψ ψ is constant. Assuming translation invariance one can analytically determine the phase diagram of the GN model at finite temperature and fermion density in the large -N f limit [26]. But in the condensed matter community it has been known for a while that the Peierls instability may trigger a breaking of translation invariance. This explains, for example, the inhomogenous Fulde-Ferrell-Larkin-Ovchinnikov equilibrium state for ultracold fermions [27,28] (see [29] for a recent review).
Subsequently it was shown that for N f → ∞ the relativistic GN model exhibits an inhomogeneous condensate at low temperature and high density. In a series of interesting papers [7,8] an explicit expression for the condensate in terms of Jacobi elliptic functions has been derived.

Homogeneous phases in large N f -limit
For N f → ∞ the saddle point approximation to the functional integral with integrand exp(−N f S eff ) in eq. (5) becomes exact. This means that the condensate σ is identical to the field σ which minimizes S eff . In particular, if we assume translation invariance then we may minimize S eff [σ] on the set of constant fields. But for a constant σ the regularized action is proportional to the Euclidean spacetime volume, We renormalize the theory such that the minimum of the effective potential at zero temperature and zero chemical potential is at some σ 0 > 0. This determines the bare coupling g Λ as function of the dimensional parameter σ 0 and the momentum-cutoff Λ, The renormalized potential has the simple form with one-particle energies p = √ k 2 + σ 2 . In accordance with our previous discussion it is an even function of the auxiliary field and of the chemical potential.
The minimizing field as function of the temperature and chemical potential is depicted in Figure  1, left. Throughout the present work we use σ 0 to set the scale and thus measure the chemical potential, temperature and condensate field in units of σ 0 . The phase diagram shows a symmetric phase with vanishing condensate and a broken phase with homogeneous condensate.
The system undergoes a phase transition from the symmetric phase at high temperature or large chemical potential to the broken phase a low T and small µ [26,30]. At vanishing chemical potential the transition happens at T c = e γ /π ≈ 0.567. The second order line extends up to the Lifshitz point at (µ 0 , T ) ≈ (0.318, 0.608), where it turns into a first order line. The latter intersects the zero-temperature axis at µ c = 1/ √ 2 ≈ 0.707. At µ c and T = 0 the condensate jumps from 1 to 0.

Inhomogeneous phase in large N f -limit
At low temperature and large chemical potential the minimum of S eff does not correspond to a homogeneous but to a spatially inhomogeneous condensate σ(x) . For a time-independent but spatially varying auxiliary field σ the eigenfunctions of D have the form ψ nm (x ) = e iωnt ψ m (x) with Matsubara frequency ω n . Summing over these frequencies in log det D one arrives at the renormalized effective action where the same renormalization prescription as in the homogeneous case has been adopted. The ε m are the real eigenvalues of the hermitean Dirac-Hamiltonian which appears in the decomposition Theε m are the eigenvalues of the Dirac Hamiltonian with constant auxiliary fieldσ given bȳ The two sums over the negative one-particle energies in the first line of (19) are easily identified as difference of two divergent vacuum energies: one for the prescribed auxiliary field σ(x) and the other for the constant reference field σ defined in (22). A heat kernel regularization reveals that the first line in (19) is UV-finite if and only if the reference field is chosen as in (22). The T -and µ-dependent traces in the second line in (19) are manifestly UV-finite and represent the finite temperature and density corrections.
In the large-N f limit only auxiliary fields which minimize S eff [σ] contribute to the functional integral in (2). With the Hellman-Feynman formula for the expectation values ε m = ψ m |h σ |ψ m the variational derivative of S eff with respect to σ can be calculated and one ends up with the Gap equation This renormalized self-consistency equation is a complicated functional equation, whose solutions have been investigated at various times in the literature. Most derivations given previously derived the regularized gap equation from the regularized trace of the Green-function with bare coupling constant and cutoff parameter [31][32][33][34]. Here the point of departure is the renormalized effective action (19) with physical scale parameter σ 0 and only finite quantities enter the derivation of the gap equation.
To summarize, to calculate the chiral condensate at finite temperature and finite density in the large-N f limit one must solve the spectral problem for the σ-dependent Dirac Hamiltonian (20) and find a self-consistent solution σ(x) of the gap equation (23). At zero temperature and fermion density Dashen et al. indeed could solve the coupled system for the modes ψ n (x) and the scalar field σ(x) by using powerful inverse scattering methods [31]. They observed that a scalar field could only solve the gap equation if the solutions of the Dirac equation ψ n are not reflected. Their space-dependent solutions describe n-particle bound states with filled Dirac sea and masses Self-consistent solutions at finite temperature and fermion density have been constructed by Thies et al. [7,8] by some (nonlinear) superposition of kink-antikink solutions. They succeeded to construct periodic solutions σ(x) with associated Bloch waves ψ n (x) of the coupled system (20) and (23) in a certain region of the (T, µ) phase diagram. The Bloch waves (they are solutions of the Lamé equation) and the scalar field σ(x) are given in terms of Jacobi's elliptic functions, see Ref. [8]. The associated Dirac-Hamiltonian has one gap in the spectrum and the periodic and anti-periodic states at the band-edges are given by particular simple Jacobi functions. Thus, the property that h σ shows no reflection for baryon excitations above the vacuum is replaced by the property of having exactly one band gap in the spectrum if the system has high density.
In the large-N f limit where the saddle point approximation to the functional integral (5) becomes exact the inhomogeneous condensate σ(x) minimizes the effective action (19) and thus is given by the solution of the gap equation, i.e. by a Jacobi elliptic function. For points in the phase diagram where the inhomogeneous solution has a lower effective action as any homogeneous solution the system is in a inhomogeneous phase. The correct phase diagram in the large-N f limit is depicted in Figure 1, right. Note that the metastable phases and first order transition line (to guide the eye this line is kept as dashed line) disappear and are replaced by two second order transition lines. At low temperature and small chemical potential there is a homogeneous phase with broken chiral symmetry, at sufficiently high temperature we are in the homogeneous symmetric phase and at low temperature and large chemical potential we are in the inhomogeneous phase with an oscillating chiral condensate.
The wave length and amplitude of the condensate in the inhomogeneous phase are determined by the chemical potential or equivalently by the Fermi-momentum and by the temperature. If one moves within the inhomogeneous phase towards the symmetric phase, the amplitude of the condensate vanishes. If one moves towards the homogeneously broken phase, then the wave length of the condensate increases. In this work we mainly address the question whether there exists an inhomogeneous phase for a finite number of flavors N f or whether such a phase is an artifact of the large-N f limit.

Spontaneous breaking of a continuous symmetry in 1+1 dimensions
A well-known theorem by N. D. Mermin and H. Wagner in statistical mechanics states that a continuous symmetry cannot be spontaneously broken at finite temperature in 1-and 2dimensional statistical systems with short range interaction [35]. A similar theorem has been proven by S. Coleman for relativistic quantum field theory in d ≤ 2 dimensions [36]. Indeed, if spontaneous symmetry breaking of a continuous symmetry would occur, then as a consequence of the Goldstone theorem [37,38] one would expect to find massless Nambu-Goldstone bosons (NGBs) in the particle spectrum. But massless scalars with a relativistic dispersion relation have an infrared divergent correlation function in lower dimensions and thus should not exist.
When proving the Goldstone theorem one makes basic assumptions: the theory should be Lorentz invariant, the Hilbert space should be positive and a global symmetry group G should be broken to a subgroup H. Then there exists one massless scalar particle for each broken symmetry (or broken generator) such that there are n BG = dim(G/H) massless Goldstone bosons.
In non-relativistic systems and for spacetime symmetries the situation is more intricate, since sometimes NGBs have unusual dispersion relations or they are even redundant.
• For example, in a ferromagnet and antiferromagnet we may have the same spontaneous symmetry breaking O(3)→O (2), but in the first case we have only one NGB (the magnon) whereas in the second case there are two NGB. Since quantum field theories at high densities may be described by quasi-excitations with non-relativistic dispersion relations a similar reduction of NGB may happen in the high-density GN model.
• In addition, for a breaking of spacetime symmetries the simple counting rule does not apply. For example, crystals have phonons for spontaneously broken translations but no gapless excitations for equally spontaneously broken rotations. Again a reduction of the number of NGB may happen if we are dealing with spacetime symmetries instead of inner symmetries.
• Finally, it may happen that the NGB completely decouple from the rest of the system. Then one may evade the conclusion of Coleman's theorem about the non-existence of NGB in 2 spacetime dimensions. This seems to happen in the large N f -limit of the GN model [39], where translation invariance is definitely broken for high fermion density.
In 1976, Nielsen and Chadha [40] presented a general counting rule of NGBs valid either with or without relativistic invariance. They divided the modes into two classes, based on the behavior of their dispersion relations for small |k|: Relativistic modes are of type I and non-relativistic modes are of type II. By examining analytic properties of correlation functions they showed that where n NGB is the total number of NGBs and n I and n II are the number of type I and type II NGBs. The number of broken generators n BG = dim(G/H) agrees with the number of flat directions of fluctuations of the order parameter.
In passing we note that there exists a related, but in general slightly different division of NGBs into type A and B [41,42]. It is an algebraic classification based on the Lie algebra of symmetry generators. To each pair of non-commuting symmetry generatorsQ i ,Q j (the conserved charges belonging to the symmetry group G) a NGB of type B is associated. The NGBs of type A tend to be linearly dispersive for small |k|. There is a simple counting for these NGBs: where ρ is the Watanabe-Brauner matrix build from the conserved charges related to the sym- For more details we refer to Refs. [43,44].
Strictly speaking the above results hold for internal symmetries only. But it is believed that the NGBs originating from a spontaneously broken translation symmetry can be treated in essentially the same way as those associated with internal symmetries [43]. This leaves us with the following scenarios for the finite-temperature GN model at high density: • Only the (abelian) spatial translation symmetry is broken such that the Watanabe-Brauner matrix ρ vanishes. Then the results (27) would imply that there is just one type A NGB.
If this would be -as expected -a NGB of type I with relativistic dispersion relation then we would be confronted with infrared divergences. The way out could be that it is not of type I but of type II or that it fully decouples from the system.
• Alternatively, if the above results do not apply to the breaking of translation invariance at high density systems, then we may as well find no NGB or a NGB of type II with nonrelativistic dispersion relation ε k ∼ |k| 2 . Its correlation function is not infrared divergent and the problem with the spontaneous breaking would go away.
Since the inhomogeneous condensate appears (at least in the large-N f limit) at high density, a non-relativistic dispersion relation seems to be more likely than a relativistic one. Unfortunately, with the available ensembles on lattices of spatial extent up to N s = 725 lattice sites we cannot reliably measure the dispersion relation of the NGB (if it exists) in the model with N f = 2 flavors in the infrared and thus cannot decide whether the NGB has a non-relativistic or a relativistic dispersion relation. It may even be that for finite N f there is no SSB of translation invariance in the strict sense and that the model behaves like a simple atomic liquid, for example as liquid argon. Indeed, the correlator of the condensate on large lattices as presented in section 4 resembles the radial pair correlation function in an atomic fluid, see the reviews [45,46]. In a forthcoming accompanying publication we will further substantiate, by studying the baryon number as function of the chemical potential, that the GN model at high density is either a crystal or an extremely viscous fluid.
3 Lattice field theory techniques

Notation
The number of lattice sites in temporal and spatial direction are denoted by N t and N s , respectively. Consequently, the temperature is given by T = 1/N t a and the extent of the periodic spatial direction by L = N s a, where a is the lattice spacing.
In the following we consider spacetime averages of observables O[σ] for given field configurations We also compute ensemble averages, where the sum over σ extends over N conf field configurations σ(x ) generated by Monte Carlo sampling according to e −S eff [σ] . The "≈" sign in eq. (29) becomes the identity "=" in the limit N conf → ∞.
Moreover, we use the discrete Fourier transform either with respect to spacetimẽ or spaceF with time and space restricted to the lattice sites, i.e.
The corresponding momenta are For fermions we impose antiperiodic boundary conditions (BC) in t-direction such that the integer-spaced n 0 are half-integer valued. For bosons we impose periodic BC in t-direction such that the integer-spaced n 0 are integer valued. Both bosons and fermions are periodic in x-direction such that the integer-spaced n 1 are integer valued.

Lattice discretizations of fermions
We use two different lattice discretizations of fermions, naive fermions and SLAC fermions. Both discretizations have certain advantages but come also with subtleties, which are discussed in the following.
In section 4 we present numerical results for both discretizations and find agreement, which we consider an important cross check. Another comparison of the two discretizations can be found in section 3.2.3, where we show the N f → ∞ phase diagram with the restriction to homogeneous condensate σ.

Naive fermions
The naive discretization is at first glance the most straightforward lattice discretization of fermions (see e.g. the textbooks [47][48][49]). In contrast to other common fermion discretizations, e.g. Wilson fermions, the free massless naive fermion action is chirally symmetric, which is an essential and necessary property in the context of this work. Naive fermions, however, lead to fermion doubling according to the Nielsen-Nynomia theorem [23]. Thus, for most applications, e.g. QCD with 2, 2+1 or 2+1+1 quark flavors, naive fermions are not appropriate. In our case of 1 + 1 spacetime dimensions the number of fermion flavors is restricted to multiples of 4. This is not a severe limitation, since we are not interested in a particular number of flavors < 4, but mostly in simulating finite numbers of flavors, e.g.
Besides fermion doubling there are, however, further pitfalls, which might lead to a continuum limit different from the theory of interest. In the context of the GN model, this was first observed and discussed in Ref. [50]. In appendix A we reproduce the arguments of Ref. [50] and we derive a modification of the straightforward naive discretization of the GN model, which has the correct continuum limit. 1 This modified naive lattice action is with the well-known action for naive free fermions coupled to a chemical potential and the auxiliary field summed over neighbors with separation a and √ 2a where N f is a multiple of 4 and N t , N s are even integers such that all doublers obey the same BC (for details see appendix A). For all computations with naive fermions presented in the following sections we use the action (34).

SLAC fermions
For SLAC-fermions the non-local derivatives in the Dirac operator are easily characterized in momentum space [53], with the Fourier transform F as defined in eq. (30). We choose the discrete momenta k µ symmetric to the origin to end up with a real and antisymmetric matrix ∂ µ [54]. This means that in spatial direction (with periodic BC) the lattice has an odd number N s of lattice sites, whereas in temporal direction (with antiperiodic BC) the lattice has an even number N s of lattice sites [49], Both the naive and the SLAC derivative define chiral fermions, for which i / ∂ is hermitean and anticommutes with γ * = iγ 0 γ 1 . In contrast to naive fermions, however, there are no doublers for SLAC-fermions. Thus they can be used to study any positive integer number of fermion flavors. 2 We point out that the non-local SLAC-derivative must not be used in gauge theories, where the edge of the Brillouin zone (where k µ jumps) is probed in the functional integral, which leads to a clash with Euclidean Lorentz invariance [55]. But SLAC-fermions have been successfully applied to non-gauge theories, for example scalar field theories [56], non-linear sigma models [57], supersymmetric Yukawa models [56] and more recently to interacting fermion systems [24].
For SLAC-fermions the chemical potential µ is introduced as in the continuum theory, Note that the chemical potential µ enters linearly and not via exp(±aµ) multiplying a hopping term as e.g. for naive fermions (see e.g. eq. (35)) For some observables, for example the fermion density, this introduces an additional term ∝ µ in the continuum limit, which can be easily subtracted, since (in 2 spacetime dimensions) the term is finite and can be calculated analytically. We emphasize that this term is not a lattice artifact -it also exists in continuum theories when appropriate care is taken in manipulating divergent integrals [58]. After the subtraction is performed, results obtained with SLAC-fermions converge much faster to the continuum limit than for other fermion discretizations. 3 A similar observation has been made when using a linear µ for naive fermions in 4 spacetime dimensions [58]. All observables considered in the present work need no such subtraction.

Comparison of naive and SLAC fermions for N f → ∞ and homogeneous condensate
In Figure 2 we show the N f → ∞ phase diagram with the restriction to a homogeneous condensate σ for naive and for SLAC fermions. The corresponding computations are straightforward and computationally rather cheap, when using techniques similar to those discussed in Refs. [18][19][20]. For both discretizations we performed computations for several significantly different values of the lattice spacing, a ≈ 0.41/σ 0 and a ≈ 0.20/σ 0 (naive and SLAC) and a ≈ 0.10/σ 0 (only naive), but similar spatial extent L. When decreasing the lattice spacing, the results obtained with each of the two discretizations approach the continuum result from Ref. [26]. Note, however, that discretization errors for SLAC fermions are almost negligible, i.e. significantly smaller than discretization errors for naive fermions.

Simulation setup
We use a standard RHMC (Rational Hybrid Monte Carlo) algorithm [59] to perform numerical simulations. In detail we use the implementation described in Ref. [60], which was also used in Refs. [24,61,62].

Scale setting
We assume that at chemical potential µ = 0 and temperature T = 0 the system is in a homogeneously broken phase and use the (positive) expectation value to set the scale. In other words, we express all dimensionful quantities in units of σ 0 , e.g. we use for the chemical potential µ/σ 0 , for the temperature T /σ 0 , etc. Setting the scale via σ 0 was also done in previous analytical and numerical studies of the phase diagram of the GN model in the N f → ∞ limit (see e.g. [7,8,[18][19][20]), i.e. expressing dimensionful quantities in units of σ 0 allows a straightforward comparison of our results at finite N f to existing N f → ∞ results.
The determination of σ 0 in lattice units is technically straightforward. When increasing the number of lattice sites in temporal direction N t as well as in spatial direction N s at fixed coupling g 2 , the ensemble average |σ| | µ=0,T quickly approaches the constant σ 0 . Thus, in practice, one just has to compute |σ| | µ=0,T on a lattice with sufficiently large N t and N s , where |σ| | µ=0,T ≈ σ 0 . This is illustrated in Figure 3 for N f = 8, SLAC fermions and two different g 2 .
As e.g. in lattice simulations of 4-dimensional Yang-Mills theory or QCD, the lattice spacing a is a function of the dimensionless coupling g 2 and can be set by choosing appropriate values for g 2 . This is reflected by the two plateau values at small 1/N t in Figure 3 representing aσ 0 (the lattice spacing in units of σ 0 ), which correspond to g 2 = 0.192 (larger lattice spacing) and g 2 = 0.161 (smaller lattice spacing).

Ensembles of field configurations
To explore the µ-T phase diagram of the 1 + 1-dimensional GN model and its dependence on the number of fermion flavors N f and to exclude sizable lattice discretization and finite volume corrections, we generated a large number of ensembles of field configurations σ(x ). These ensembles are listed in Table 1.
For given coupling g 2 , i.e. for fixed lattice spacing a, we vary the temperature T = 1/N t a by changing N t , the number of lattice sites in temporal direction. Thus, at fixed g 2 the temperature T can only be changed in discrete steps. The chemical potential µ, on the other hand, is not restricted in such a way and can be set to any value.
The majority of simulations were carried out for N f = 8: • We simulated at several different spatial extents with 31 ≤ N s ≤ 128 corresponding to L = N s a to check for finite volume corrections.
• We simulated at many different values of the chemical potential, to explore the phase diagram.
• We carried out a sizable amount of these simulations using both fermion discretizations, i.e. SLAC fermions and naive fermions, to cross-check our results (the corresponding coupling constants g 2 have been tuned in such a way, that the simulated lattice spacings are almost identical).
Simulations at N f = 2 and N f = 16 were done with SLAC fermions, but not with naive fermions. For each ensemble between 300 and 10 000 configurations were generated.

Numerical results
The majority of results shown in the following (section 4.

Qualitative expectations
In the 1 + 1-dimensional GN model in the limit N f → ∞ there are three phases, a symmetric phase, a homogeneously broken phase and an inhomogeneous phase (see the discussion in section 2.3). The structure of the field configurations σ(x ) generated during our simulations by the HMC algorithm suggest that there is a similar phase structure at finite N f . At large T the field σ(x ) mostly fluctuates around zero, while at small T and small µ either σ(x ) ≈ +σ 0 or σ(x ) ≈ −σ 0 . Most interestingly, however, at small T and large µ the field σ(x ) exhibits spatial periodic oscillations similar to a cos-function, which might signal an inhomogeneous phase. 4 An example of a typical field configuration at (µ/σ 0 , T /σ 0 ) ≈ (0.450, 0.030) with such periodic oscillations is shown in Figure 4.
Thus, we expect that the field configurations σ(x ) generated by the Monte Carlo algorithm are crudely described by the following model: • Inside a symmetric phase σ(x ) = εη(x ) .
• Inside a homogeneously broken phase • Inside an inhomogeneous phase ε ≥ 0, σ 0 > 0 and A ≥ 0 are real parameters, which depend on µ and T . η(x ) are independent continuous random variables with Gaussian probability distributions p(η(x )) ∝ exp(−η(x ) 2 /2), which represent statistical fluctuations. λ = L/(q + δq) is the wavelength of σ in an inhomogeneous phase, where q ≥ 1 is an integer parameter and δq is an integer-valued discrete random variable with Gaussian probabilities p(δq) ∝ exp(−δq 2 /2∆q). ∆q q, the width of the Gaussian, is another real parameter. δx ∈ [0, L) is also a random variable, where it is a priori not clear, what kind of distribution to expect. The distribution could depend on the details of the HMC algorithm, and whether translational symmetry is spontaneously broken or not. Note, however, that the observables we are studying are constructed in such a way that they are independent of this distribution. To summarize, the model defined by eqs. (41) to (43)  in a homogeneously broken phase and around a cos-function with varying wavelength λ in an inhomogeneous phase.
With this model in mind, which is based on existing results in the N f → ∞ limit [7,8], we designed several observables, which are able to distinguish the three phases. Note that the sole purpose of this model is to provide some guidelines for the construction of observables and to develop expectations, in which way they characterize the three phases. The model is not used elsewhere in this work, in particular not for the analysis of our numerical results.

Squared spacetime average of σ(x )
A rather simple observable is the normalized ensemble average of the squared spacetime average of σ(x ). It is not suited to characterize an inhomogeneous phase, but still useful to distinguish a homogeneously broken phase from a symmetric phase and an inhomogeneous phase. 5

Within the model defined in section 4.1 one finds
• inside a symmetric phase • inside a homogeneously broken phase The red regions correspond to a homogeneously broken phase and the green regions to a symmetric and/or an inhomogenous phase. The gray lines represent the N f → ∞ phase boundary of the homogeneously broken phase [7,8].
• inside an inhomogeneous phase Thus, we expect Σ 2 ≈ 0 both inside a chirally symmetric phase and inside an inhomogeneous phase, while it should be significantly larger, Σ 2 ≈ 1, inside a homogeneously broken phase.
In Figure 5 we show Σ 2 in the µ-T plane for naive fermions and SLAC fermions. The red regions clearly indicate a homogeneously broken phase, while the green regions represents a symmetric and/or an inhomogenous phase. The two plots are very similar. The main reason for the small discrepancies are lattice discretization errors, which are expected to be significantly larger for naive fermions than for SLAC fermions (see section 3.2.3). To ease comparison with N f → ∞ results, we included the corresponding phase boundary of the homogeneously broken phase from Refs. [7,8]. It is obvious that the homogeneously broken phase at finite N f = 8 is of similar shape, but of smaller size than its analog at N f → ∞. Such a reduction in size is expected, because at finite N f there are fluctuations in σ(x ), which increase disorder and, thus, favor a symmetric phase. Note that at small T the boundary between the red and the green region starts to deviate significantly from the N f → ∞ boundary and turns towards (µ, T ) = (0, 0). Our numerical results indicate that this is caused by the finite lattice spacing (see e.g. Figure 8). A qualitatively similar behavior was observed in an N f → ∞ lattice study of the GN model [18].

The spatial correlation function of σ(x )
In the limit N f → ∞ in the inhomogeneous phase, σ(x ) is a periodic function of the spatial coordinate x. It has a kink-antikink structure with large wavelength close to the boundary to the homogeneously broken phase and is sin-like with smaller wavelength for larger µ. We expect a similar behavior also at finite N f (see also eq. (43)).
Since the action S eff is invariant under spatial translations, field configurations, which are spatially shifted relative to each other, i.e. σ(t, x) and σ(t, x + δx), contribute with the same weight e −S eff to the partition function and, thus, should be generated with the same probability by the HMC algorithm (see also the discussion on the distribution of δx in section 4.1). Consequently, simple observables like σ(x ) are not suited to detect an inhomogeneous phase in a finite system, because destructive interference should lead to σ(x ) = 0 in a finite system, even in cases, where all field configurations exhibit spatial oscillations with the same wavelength. 6 An observable, which does not suffer from destructive interference and is able to exhibit information about possibly present inhomogeneous structures, is the spatial correlation function of σ(x ), i.e.
The equality holds in thermal equilibrium and a finite box of length L with periodic boundary conditions since σ(t, x + y)σ(t, y) does neither depend on t nor on y. Actually, our HMC algorithm is able to sample all field configurations and, thus, to produce y-independent expectation values σ(t, x + y)σ(t, y) . We use the sum over t and y in Eq. (48) to decrease statistical errors in the Monte Carlo average.
The correlator C(x) is our main observable to detect and to distinguish the three expected phases, in particular an inhomogeneous phase. 5 In contrast to the typical exponential decay of correlation functions, C(x) is expected to oscillate in an inhomogeneous phase, i.e. C(x) should be positive, if x/λ is close to an integer, and negative, if x/λ is close to a half-integer, where λ denotes the wavelength of the spatial periodic structure of σ(x ). Such oscillations are also found in the N f → ∞ limit [64]. The expectation is also supported by analytical calculations within our model defined in section 4.1, where • inside a symmetric phase • inside a homogeneously broken phase • inside an inhomogeneous phase with the Jacobi ϑ function ϑ(z, τ ) = 1 + 2 ∞ n=1 e iπn 2 τ cos(2πnz) .
The cos-term in eq. (51) leads to oscillations with wave length L/q, while the factor including the ϑ function causes a damping of these oscillations for increasing separations x. This damping is due to the random fluctuations of the wave number q + δq entering eq. (43) via λ, which cause destructive interference for larger x. The damping is strong for large fluctuations, i.e. large ∆q, and not present in the limit ∆q → 0. Note that the result for C(x) inside an inhomogeneous phase, i.e. the right hand side of eq. (51), is independent of the distribution of the random variable δx introduced in section 4.1.
Of similar interest as C(x) is its Fourier transform (see eq. (31)). The expected behavior is the following: • inside a symmetric phaseC(k) is rather small and smooth without any pronounced peak.
• inside a homogeneously broken phaseC(k) has a pronounced peak at k = 0 and is rather small and smooth at k = 0.
• inside an inhomogeneous phaseC(k) has pronounced peaks at k = ±q, where q is related to the wavelength of the spatial oscillations of C(x) via λ = L/q.

This expectation is in agreement with results obtained within our model from section 4.1, where
• inside a symmetric phaseC • inside a homogeneously broken phasẽ • inside an inhomogeneous phasẽ Exemplary results for C(x) andC(k) are shown in Figure 6 inside the symmetric phase ((µ/σ 0 , T /σ 0 ) ≈ (0, 0.993)), inside the homogeneously broken phase ((µ/σ 0 , T /σ 0 ) ≈ (0, 0.083)) and inside the inhomogeneous phase ((µ/σ 0 , T /σ 0 ) ≈ (0.700, 0.083) and (µ/σ 0 , T /σ 0 ) ≈ (0.900, 0.083)). In all cases there is reasonable agreement between our results for naive fermions and for SLAC fermions. Since lattice discretization errors are expected to be significantly larger for naive fermions, as discussed in section 3.2.3, we show results obtained with naive fermions in the inhomogeneous phase for two different lattice spacings, a ≈ 0.252/σ 0 and a ≈ 0.126/σ 0 . Those corresponding to the finer lattice spacing are closer to the SLAC results, where a ≈ 0.250/σ 0 . We interpret this as indication that both discretizations agree in the continuum limit. Moreover, on a qualitative level there is agreement with our crude expectations summarized by eqs. (49) to (51) and eqs. (54) to (56), when the parameters in these equations are chosen appropriately. In particular the plots in the lower half of Figure 6 clearly indicate the existence of an inhomogeneous phase. C(x) exhibits cos-like oscillations with decreasing wavelength λ for increasing µ, as observed in the N f → ∞ limit. This is also reflected   The red regions correspond to a homogeneously broken phase, the green regions to a symmetric phase and the blue regions to an inhomogenous phase. The gray lines represent the N f → ∞ phase boundaries [7,8].
by the symmetric pair of peaks ofC(k) at the corresponding wave numbers q = L/λ. The gray curves in these plots represent the model expectations (eqs. (51) and (56)) with parameters A, q, ∆q and determined by fits to the lattice results for C(x) andC(k).
A straightforward calculation leads tõ is the Fourier transform of σ(x ) with respect to the spatial coordinate (see eq. (31)). The absolute values |σ(t, k)| are invariant under spatial translations x → x + δx, because This shows again that both C(x) andC(k) do not suffer from destructive interference, as already discussed at the beginning of this subsection. Moreover, eq. (57) shows in an explicit way that the Fourier transformed correlation functionC(k) also provides information about the absolute values of the Fourier coefficients of the field σ(x ). In particular the peaks inC(k) at nonvanishing k in the plots in the lower half of Figure 6 indicate that inside an inhomogeneous phase strong oscillations with the same wavelength are present in the majority of the generated field configurations.
From Figure 6 one can see Thus, the minimum of the correlation function C(x) is suited to plot a crude phase diagram as shown in Figure 7 both for naive and for SLAC fermions.  The red region indicates a homogeneously broken phase, the green region a symmetric phase and the blue region an inhomogeneous phase. As before, results obtained with these two different fermion discretizations are in fair agreement. Moreover, the phase diagram is qualitatively similar to the N f → ∞ phase diagram from Refs. [7,8], whose phase boundaries are also shown in Figure 7. The homogeneously broken phase and the inhomogeneous phase are, however, somewhat smaller for finite N f than for N f → ∞, presumably because quantum fluctuations at finite N f increase disorder and, thus, favor a symmetric phase. Note, however, that C min is not the expectation value of a product of local operators as for example C(x), which is the two-point function of the order parameter. In general one must be cautious, when using nonlocal quantities like C min , since they can fake non-existing phase transitions [65,66]. But in the present case transition lines have been localized by C min as well as the correlator C(x) of the local field σ.
We also checked the stability of the phase diagram with respect to variations of the lattice spacing and the spatial volume. To this end we performed simulations using SLAC fermions at three different values of the lattice spacing, a ≈ 0.410/σ 0 , 0.250/σ 0 , 0.195/σ 0 (the columns in Figure 8), and for four different numbers of lattice sites in spatial direction, N s = 31, 47, 63, 127 (the rows in Figure 8). Approaching the infinite volume limit at fixed lattice spacing corresponds to moving from the top to the bottom of the figure, while approaching the continuum limit at approximately fixed spatial volume corresponds to moving right and downwards at the same time. There is little difference in the crude phase diagrams shown in these twelve plots. We consider this as indication that our results, at our current level of accuracy, are consistent with results in the continuum and infinite spatial volume.
In the limit N f → ∞ the phase boundaries between the three phases are of second order. In the following we present selected results for finite N f and discuss, whether there are also phase transitions or rather crossovers.
• For the transition between the symmetric and the homogeneously broken phase we computed C min as function of the temperature T for vanishing chemical potential µ = 0 (see Figure 9, left plot, blue points). There is a rapid decrease of C min at around T /σ 0 = 0.25, which is qualitatively reminiscent to the N f → ∞ result and, thus, might indicate that there is also a second order phase transition at finite N f .
• For the transition between the homogeneously broken and the inhomogeneous phase we computed C min as function of the chemical potential µ for rather low temperature T /σ 0 ≈ 0.102 (see Figure 9, right plot). We observe a rapid decrease of C min at around µ/σ 0 = 0.5, which indicates a phase transition similar to the N f → ∞ case.
• As can also be seen from the phase diagram in Figure 7, the transition between the symmetric and the inhomogeneous phase is somewhat washed-out. This is also reflected by Figure 9, left plot, where the orange points represent C min as function of the temperature T for chemical potential µ/σ 0 ≈ 0.600. These results favor a weak phase transition of second or higher order or just a crossover. 7 As has been pointed out earlier, these conclusions may not be fully coherent since C min is a non-local quantity.

The long-range behavior of C(x)
We have argued in section 2.4 that SSB of a continuous (spacetime) symmetry in 1+1 dimensions is a delicate issue. To decide, whether there are NGB and if so, what type of NGB, we investigate the long-range correlations of the GN model with N f = 2 flavors. The question is, whether the long-range order (which in the present context means that C(x) oscillates with a constant nonzero amplitude for arbitrarily large |x|), which is necessary to form a crystal at N f = ∞, remains at finite N f long-range or becomes almost long-rangeà la Berezinskii, Kosterlitz and Thouless (BKT) [67,68] (in which case the amplitude decreases with distance like an inverse power). Indeed, by studying the long-range behavior of the SU(N f ) Thirring model (this is a four-Fermi theory with current-current interaction) in 1 + 1 dimensions, Witten argued that for finite N f the correlations have almost long-range order, such that only for N f → ∞ the continuous chiral symmetry is broken, i.e. ψ ψ = 0 [69]. This way the system circumvents the no-go theorems for SSB of continuous inner symmetries. In what follows we try to answer the question whether a similar mechanism is at work for translation symmetry in the GN model.
To detect SSB of translation invariance directly we could break translation symmetry explicitly in a finite box with periodic BC, for example by adding a term ε(x)σ(t, x) to the Lagrangian, perform the infinite volume limit and finally remove the source ε(x). Assuming clustering in thermal equilibrium and we conclude that C(x) is for large |x| proportional to the condensate (calculated with first L → ∞ and afterwards an adapted ε → 0), i.e.
In the inhomogeneous phase we can write where A(x) is the non-increasing amplitude function, while C periodic (x) represents the periodic oscillations. If the system forms a crystal, the amplitude function A(x) should approach a nonzero constant for sufficiently large separations |x|. In case the system has almost long-range orderà la BKT, the amplitude function decreases with |x| as an inverse power. To distinguish the two scenarios we study C(x) for small N f = 2, to detect a possible deviation of A(x) from the asymptotically constant behavior in the large-N f limit (for small N f quantum fluctuations might be strong enough to change long-range into almost long-range order as it happens in the chirally invariant SU(N f ) Thirring model for finite N f [69]). Thus, we expect one of the following amplitude functions: 1. In a BKT-like phase without SSB we expect the amplitude function A(x) to have the following behavior for large |x|: The dots indicate sub-leading terms and terms arising from the finite spatial extent of the system.
2. If there is SSB of translation invariance, C(x) oscillates with constant non-zero amplitude at large |x|, where the short-ranged contributions from excited states are suppressed. The amplitude function would then be depending on whether the excitations over the oscillating condensate are massive or massless. If the NGB decouple from the system, the amplitude function approaches the constant amplitude γ = 0 exponentially fast. If not, A(x) will approach γ = 0 with an inverse power of |x|.
Before discussing the long-range behavior of C(x) we present the phase diagram from C min for N f = 2 in Figure 10. Again we recognize the same phases as in the large-N f limit: a homogeneously broken phase which (as expected) is significantly smaller than for N f = 8 and N f → ∞, a symmetric phase for sufficiently large temperature and a region, where C min is clearly negative. One unexpected and striking feature of the latter phase is that the temperature range with negative C min grows with increasing µ, which is qualitatively different from the situation at large N f . This already happens for N f = 8 in a small region of parameter space (see Figure 8) but is so pronounced at N f = 2 that up to µ/σ 0 = 1.4 we found no evidence that the transition line separating the inhomogeneous and the restored phase will bend down to the µ axis. The red region corresponds to a homogeneously broken phase, the green region to a symmetric phase and the blue region to an inhomogenous phase. The grey lines represent the N f → ∞ phase boundaries [7,8].
Now we investigate the long-range behavior of the amplitude function A(x) defined in eq. (64) for N f = 2 on lattices with a rather large number of sites in spatial direction. Figure 11 shows the correlator C(x) and its Fourier transformC(k) for (µ/σ 0 , T /σ) = (0.500, 0.030) and various N s up to 725 corresponding to spatial extents up to L ≈ 297.3/σ 0 . We observe 32 periods of statistically significant oscillations in C(x) over the whole range of separations and a pronounced peak ofC(k) at the corresponding wave number. The position of the peak is essentially the same for all L demonstrating once more that the wave length is independent of the spatial extent.
The amplitude function A(x) is extracted from the peaks of the correlation function C(x) (which we identified using the scipy.signal.find peaks method [70] with prominence=0.01) for all N s , as exemplified for N s = 725 in the left plot of Figure 12. The peaks of C(x) for various N s are depicted in the right plot of Figure 12. There is a rapid drop for small separations x that flattens out for asymptotically large x. We performed χ 2 -minimizing fits of symmetrized    versions of the expectations for the amplitude function A(x) (eqs. (65) and (66)) to the extracted peaks for N s = 725 (see Figure 12, left plot). We only used data points with x ≥ x min in the fitting procedure, because all three fit functions are expected to model the amplitude function for large |x|. Since χ 2 red as a function of x min is almost constant for large |x| (see Figure 13, upper plot), we chose (independently for each of the three fits) the minimal value x min , where χ 2 red is consistent with the asymptotic constant, i.e. that value, where the constant behavior sets in. In the lower plot of Figure 13 we show the extracted parameters γ appearing in the two fit functions in eq. (65) as functions of x min . It is reassuring that both parameters γ are essentially independent of x min as long as the corresponding χ 2 red is small. Since there are massive excitations in the SSB model and massless excitations in the SSB' and BKT models, one should not expect to find the same x min for the three models, but a smaller x min for the SSB model, where the exciations are short-range. This expectation is confirmed by our numerical results (see column "x min " in Table 2).
The fit results for the parameters of the amplitude functions A(x) from eqs. (65) and (66) for N s = 725 are collected in Table 2. There are several points to note: • The SSB model admits the most stable fits with resulting parameters almost independent of x min and the initial values used in the fitting algorithm. γ is cleary different from zero.
• Fits for the SSB' model are less stable, which is reflected by the large uncertainties obtained for α and β. γ, however, can be determined in a reliable way and the result is again different from zero. Moreover, it is in excellent agreement with the corresponding result for the SSB model. It is also interesting to note that the SSB' model, which differs from the BKT model by the additive constant γ, leads to significantly smaller χ 2 red (and γ = 0) than the BKT model, when the same x min is used.
• The BKT model is only able to describe the amplitude function for rather large |x|, i.e.
the corresponding x min is significantly larger than for the SSB model and SSB' model. The fit result for the exponent is β = 0.66 ± 0.13, which is similar to the corresponding analytically known value β = 1/N f = 1/2 of the SU(N f ) Thirring model.
To summarize, it seems that the excitations are probably massive, because the SSB model is able to describe the extracted peaks for significantly smaller separations x, than it is possible with the SSB' model or the BKT model. When allowing for a constant γ (as in the SSB model and the SSB' model), the fits lead to stable and clearly non-zero results. Thus, our current data is best described by the SSB scenario, where the NGB completely decouples from the system. However, this does not imply that this is the physical situation in the thermodynamic limit, since we saw that even N s = 725 lattice sites in spatial direction are not sufficient, to rule out the 1/|x| β almost long-range behavior of the BKT scenario. In other words, a constant behavior ∝ γ and an inverse power ∝ 1/|x| β are very similar for large |x|, in particular in a periodic spatial volume, such that even larger lattices or higher accuracy is needed to clearly distinguish between these scenarios.
Despite the fact that we could not fully reveal the nature of the inhomogeneous phase, we can still argue that there exists a phase transition between the inhomogeneous low-temperature phase and the symmetric high temperature phase. This can be seen from Figure 14, where we show the spatial correlation function C(x) together with cosh-fits for three different (µ, T ) in the symmetric phase. It is evident that the cosh-functions perfectly fit the data points, which indicates that in the symmetric phase the excitations are massive. In contrast to that, at low temperature in the inhomogeneous phase C(x) is long-range or almost long-range (see Figure 11). Since C(x) behaves qualitatively different in the inhomogeneous phase and in the symmetric phase, i.e. long-range or almost long-range versus exponentially decaying, we expect a phase transition, either a symmetry-restoring transition or a BKT-like transition from the low-temperature inhomomgeneous phase to the high-temperature symmetric phase.

4.5
Approaching the N f → ∞ results with computations at finite N f In sections 4.2 to 4.4 we have presented results at N f = 2 and N f = 8, which are similar to the analytically known N f → ∞ results [7,8], e.g. for the phase diagram. To check and to confirm that results at finite N f approach for increasing N f the N f → ∞ results, we also performed simulations at N f = 16. An exemplary plot is shown in Figure 15, where Σ 2 is shown as a function of the temperature T for vanishing chemical potential µ = 0 and N f = 2, 8, 16, ∞.
While results for N f = 2 agree with the N f → ∞ result only for rather small T /σ 0 , there is agreement also for larger T /σ 0 , when N f is increased, indicating that one can approach the analytically known N f → ∞ results with computations at finite N f .
In the present work we could localize three regimes in the space of thermodynamic control parameters T and µ, in which the two-point function of the order parameter shows qualitatively different behaviors. We spotted a homogeneously broken phase, a symmetric phase and a region with oscillating correlation function C(x) defined in eq. (48). The results of our Monte Carlo simulations with two different types of chiral fermions for systems with 2, 8 and 16 flavors on lattices with sizes up to N t = 80 and N s = 725 were presented, analyzed and discussed in the main body of the text. Although we could not answer the question, whether in GN models with a finite number of flavors translation invariance is spontaneously broken at low T and large µ, or whether the system is in a Berezinskii-Kosterlitz-Thouless like phase, we clearly spotted a low-temperature and high-density region, where the model exhibits oscillating spatial correlators. The wave-length of the spatial oscillation is determined by the chemical potential and temperature and not by the lattice spacing and the spatial extent of the lattice, i.e. spatial oscillations are neither a lattice artifact nor a finite size effect. We argued that there is a transition between the inhomogeneous phase and the symmetric phase which could be an infinite order transition (according to the Ehrenfest classification). In an accompanying paper we shall demonstrate that the ratio of the system size and the dominant wave length of the condensate oscillations is equal to the number of baryons in the systems. This further substantiates the physical picture that the GN model in equilibrium at low temperature and large fermion density either forms a crystal of baryons or a viscous fluid of baryons. In this work we also showed that the amplitude of oscillations stays constant or decays very slowly as suggested by a related result [69]. The first behavior is expected for a baryonic crystal the second behavior for a viscous baryonic fluid. To better understand, how the long range behavior at low temperature and high density does not clash with the absence of Nambu-Goldstone bosons in 1 + 1 dimensions, needs further high-precision results on the two-point function of the order parameter. If the dispersion relation is non-relativistic or if the massless modes fully decouple from the system then there should be no problem.
Independent of whether the oscillating correlator C(x) points to a baryonic crystal or a baryonic liquid for finite N f we have seen that mean-field/large-N f approximations may keep more information on the physics at finite N f than one might expect. This is reassuring since in particle physics and even more so in solid state physics we often rely on mean-field type approximations. An important question is, of course, whether our results have any relevance for QCD at finite baryon density. On the one hand, we established that the interpretation as baryonic matter is not spoiled by taking quantum fluctuations into account. On the other hand, although recent numerical investigations of four-Fermi theories in 2 + 1 dimensions and for N f → ∞ spotted inhomogenous condensates [21], the spatial modulation is related to the cutoff scale and seems to disappear in the continuum limit [22]. Clearly, if this happens then we cannot expect a breaking of translation invariance for a finite number of flavors. Thus, extending our lattice studies to higher dimensions is of relevance for QCD. Simulations of interacting Fermi theories in 2 + 1 dimensions are under way and we hope to report on our findings soon.

A Lattice discretization of the GN model with naive fermions A.1 Free naive fermions
The action of free naive fermions with chemical potential µ is given in eq. (35). The Fourier representations of the fermion fields are where the discrete 2-momenta k = (k 0 , k 1 ) are in the first Brillouin zone, −π ≤ k µ a ≤ π, and are chosen such that BC in t direction are antiperiodic and in x direction are periodic (see section 3.1, in particular eq. (30)). Inserting these Fourier representations into eq. (35) leads to where we abbreviated k 0 = cosh(µa) sin(k 0 a) a − i sinh(µa) cos(k 0 a) a ,k 1 = sin(k 1 a) a .
In the limit a → 0 the sums over k 0 and k 1 can be restricted to the "soft modes", where both |k 0 a| 1 and |k 1 a| 1. There are four regions of soft modes in the first Brillouin zone, and they are denoted by R uv with u, v ∈ {0, 1}. The momenta of the soft modes in region R uv are in the neighborhood of the four momenta at which (for µ = 0) the lattice momentak 0 andk 1 vanish. For k ∈ R uv we have Now we define the soft modes in the four regions according toχ uv (k ) =χ(k uv +k ) (and analogous forχ) with small |k µ a|. Neglecting the O(a 2 ) corrections in (71) we can approximate the free lattice action (68) by This short calculation exhibits the well-known fermion flavor doubling for each spacetime dimension. It also shows that both N t and N s must be even to obey anti-periodic boundary conditions in t direction and periodic boundary conditions in x direction for each of the four fermion flavors.
It is important to note that the action (72) differs in a couple of minus signs in front of the γ matrices for flavors (u, v) = (0, 0) from the corresponding continuum expression for four free fermion flavors. These minus signs can be eliminated by changing field coordinates viã Then eq. (72) becomes This shows that the lattice action (35) corresponds in the continuum limit to four massless non-interacting fermion flavors.

A.2 Naive fermions and the GN model
Discretizing the GN model (2) with N f = 4 flavors in a straightforward way using the fields χ andχ via actually results in a theory different from the GN model. To show this, we insert again the Fourier representations of the fermionic fields (67) as well as of the real scalar field σ, whereσ(−k ) =σ * (k ) and the discrete momenta k are chosen such that σ is periodic in x 0 and x 1 direction. The action (75) becomes S σ [χ,χ,σ] = S free [χ,χ] + i √ N t N s k k χ (k )σ(k − k )χ(k ) + N f 2g 2 k σ(k ) 2 , N f = 4 .
(77) In the limit a → 0 only the soft fermion modes contribute, as discussed in appendix A.1. Note, however, that there is no kinetic term for the field σ and, thus, no corresponding suppression of σ modes. Consequently, for a → 0 the interaction term in eq. (77) can be written as with symmetric kernel in momentum spacẽ In terms of the usual field coordinatesψ uv , related toχ uv via eq. (73), the interaction term is Now it is obvious that the action (75) is not a discretization of the GN model with N f = 4 fermion flavors. While the four terms with uv = u v in eq. (80) represent the correct GN interaction for the four fermion flavors, there are twelve additional terms not present in the GN model, where the field σ couples two different fermion flavors, i √ N t N s k ,k ψ 10 (k )γ 0σ (π/a + k 0 − k 0 , k 1 − k 1 )ψ 00 (k ) + eleven more terms , as was already pointed out in Ref. [50]. Including these twelve terms in a numerical simulation, i.e. using the action (75), corresponds to studying a different theory and leads to results significantly different from those obtained with a correct discretization of the GN model (examples are shown at the end of this section). Now we derive a proper lattice discretization of the GN model. To this end, we note that only the soft fermion modes contribute in the limit a → 0 and that the four correct interaction terms are proportional toσ uv,uv , while the twelve spurious interaction terms are proportional toσ uv,u v with uv = u v . Thus, one can eliminate the spurious terms by replacingσ(k ) in the interaction term in eq. (77) byW (k )σ(k ). HereW is a weight-function with •W (k ) → 1 for k ≈ k 00 = (0, 0) (i.e. in region R 00 ), •W (k ) → 0 for k ≈ k uv with (u, v) = (0, 0) (i.e. in the other regions R uv ).
Expressing this modified action in terms of χ(x),χ(x) and σ(x) is straightforward, where W is the Fourier transform ofW and given in eq. (36). All calculations and arguments presented in this section apply to N f = 8, 12, 16, . . . flavors in a trivial way. The generalization of eq. (84) is eq. (34).
In Figure 16 we show numerical evidence that using the straightforward naive discretization of the GN model (75) leads to incorrect results, i.e. results not corresponding to the GN model. We plot Σ 2 as a function of the temperature T for chemical potential µ = 0. The blue and orange curves correspond to the SLAC discretization (see section 3.2.2) and the correct naive discretization (34) (or equivalently (84)). These curves are rather similar and get closer, when decreasing the lattice spacing, indicating that they have the same continuum limit. The green curve, on the other hand, corresponding to the straightforward naive discretization (75) is quite different and does not approach the blue and orange curves, when decreasing the lattice spacing. We obtained similar results also for non-vanishing chemical potenial.