Inhomogeneities in the $2$-Flavor Chiral Gross-Neveu Model

We investigate the finite-temperature and -density chiral Gross-Neveu model with an axial U$_A$(1) symmetry in $1+1$ dimensions on the lattice. In the limit where the number of flavors $N_\mathrm{f}$ tends to infinity the continuum model has been solved analytically and shows two phases: a symmetric high-temperature phase with a vanishing condensate and a low-temperature phase in which the complex condensate forms a chiral spiral which breaks translation invariance. In the lattice simulations we employ chiral SLAC fermions with exact axial symmetry. Similarly to $N_\mathrm{f}\to\infty$, we find for $8$ flavors, where quantum and thermal fluctuations are suppressed, two distinct regimes in the $(T,\mu)$ phase diagram, characterized by qualitatively different behavior of the two-point functions of the condensate fields. More surprisingly, at $N_\mathrm{f}=2$, where fluctuations are no longer suppressed, the model still behaves similarly to the $N_\mathrm{f}\to\infty$ model and we conclude that the chiral spiral leaves its footprints even on systems with a small number of flavors. For example, at low temperature the two-point functions are still dominated by chiral spirals with pitches proportional to the inverse chemical potential, although in contrast to large-$N_\mathrm{f}$ their amplitudes decrease with distance. We argue that these results should not be interpreted as the spontaneous breaking of a continuous symmetry, which is forbidden in two dimensions. Finally, using Dyson-Schwinger equations we calculate the decay of the U$_A$(1)-invariant fermion four-point function in search for a BKT phase at zero temperature.


I. INTRODUCTION
A surprising amount of physical phenomena in particle-and condensed-matter physics are well described by four-Fermi theories. For instance, they are employed to model low-energy chiral properties of Quantum Chromodynamics (QCD). The effective four-Fermi theory describing the dynamics of nucleons and mesons goes back to Nambu and Jona-Lasinio (NJL) [1] and is built upon interacting Dirac fermions with chiral symmetry, paralleling the construction of Cooper pairs from electrons in the BCS theory of superconductivity.
In fact, most of our knowledge about QCD at intermediate baryon densities stems from the study of NJLtype effective theories, since in this regime one needs nonperturbative methods but cannot use lattice field theory techniques due to the complex-action problem. In a similar spirit, a four-Fermi current-current interaction among leptons (and quarks) was proven to give an accurate phenomenological description of the weak interaction at low energy p 2 m 2 W . In the pioneering work by E. Fermi the currents are made up from the proton, neutron, electron and neutrino fields [2]. In four spacetime dimensions interacting Fermi theories, such as the NJL model or Fermi theory, are non-renormalizable and thus can only serve as effective (low-energy) approximations which need to be UV completed. For the two examples given these completions are of course known.
The dynamical creation of a condensate from strong fermion interactions as seen in NJL-type models inspired many theories of the breaking of electroweak symmetry, such as technicolor (see the review [3]) and the top-quark condensate [4].
Four-Fermi theories in two spacetime dimensions are renormalizable and asymptotically free (some are integrable or even soluble) and share certain features with their cousins in four dimensions. The most prominent examples are the Thirring model with a current-current interaction [5], which is S-dual to the sine-Gordon model, and the Gross-Neveu (GN) model with a scalar-scalar interaction [6], which serves as a toy model for the theory of strong interactions.
With the discovery of novel materials (like Dirac and Weyl semimetals in two and three spatial dimensions) and the development of experimental techniques (for example optical lattices to trap atoms) we have witnessed a steadily increasing interest in models describing interacting fermions. Such models in lower dimensions describe one-dimensional and planar systems, such as polymers [7][8][9][10][11], graphene [12,13] or high-T c superconductors [14,15], to name some prominent examples.
Interacting Fermi theories at finite temperature and density were mainly investigated in the limit of a large number of fermion flavors N f . For N f → ∞ the saddlepoint approximation becomes exact and one can solve the corresponding gap equation analytically on the set of homogeneous condensates. But, for the (1 + 1)-dimensional GN model at low temperature and large chemical potential the relevant solutions of the gap equation are actually inhomogeneous in space. They have been constructed in [16] for the GN model with discrete and in [17,18] for the chiral GN model with continuous chiral symmetry. These remarkable analytic results for N f → ∞ prove the existence of inhomogeneous phases, which are regions in parameter space where the chiral condensate acquires a spatial dependence, indicating the spontaneous breakdown of not only chiral symmetry alone but in a combination with spacetime symmetries (see [19] for a review).
Are these inhomogeneous phases at large densities an artifact of the large-N f limit as suggested by various no-go theorems in two spacetime dimensions? To address this question, a better understanding of interacting Fermi systems at finite N f with regards to inhomogeneous phases is required. But the spontaneous breaking of translation invariance is not merely of academic interest: Systems where an inhomogeneous state develops spontaneously have been extensively discussed in the condensed-matter literature. A prominent example is the inhomogeneous pairing inside a superconductor in a magnetic field, predicted by Larkin, Ovchinnikov, Fulde and Ferrell (LOFF phase) [20,21]. Similar types of pairings can occur in many other physical systems, ranging from supersolids to ultracold atomic gases (see the reviews [22,23]). The UV cutoff which is inherent in all condensed-matter systems inhibits a direct translation of these findings to quantum field theory and particle physics where one removes the cutoff during the process of renormalization.
A first attempt to investigate the fate of inhomogeneous phases at finite N f has been made in recent lattice studies [24,25], where the existence of spatially varying chiral condensates in the (1 + 1)-dimensional GN model with 2, 8, and 16 flavors was confirmed. The present work serves as a follow-up, providing a similar analysis of the chiral Gross-Neveu (cGN) model with a continuous axial symmetry, characterized by the Lagrangian where g 2 denotes a dimensionless coupling constant and the two-dimensional matrix γ * = iγ 0 γ 1 is the analog of γ 5 in two spacetime dimensions. The summation over N f flavors of fermions is implied in the fermion bilinears entering Eq. (1). Below we shall see that the results of our simulations with chiral SLAC fermions resemble the analytical findings of the large-N f limit [17,18]. The analysis of the GN model in [24] has already given clear evidence that the chiral and doubler-free SLAC fermions and naive fermions yield comparable results in the continuum limit, with the former converging considerably faster. 1 Using SLAC fermions has the additional advantages that the lattice cGN model is invariant under axial U A (1) transformations and that we can study the system with N f = 2 without encountering a sign problem. With naive fermions the GN and cGN models have no sign problems only for N f a multiple of 8. In the present work, however, we want to investigate how much the models at finite flavor number differ from the analytic solutions at infinite N f , for which N f = 8 might be too large, see [24]. We 1 The same observation applies to supersymmetric Yukawa models [26,27].
do not use Wilson fermions since we are mainly interested in the chiral properties of cGN models. Staggered fermions, on the other hand, may lead to wrong results for interacting Fermi systems, as has been demonstrated in [28][29][30]. Our work is organized as follows. In Sec. II we summarize relevant facts about the finite-temperature and -density cGN model with Lagrangian (1) in the continuum, which will be used in the subsequent sections. In Sec. III the lattice cGN model with chiral SLAC fermions is presented, relevant observables are introduced and the lattice setup is discussed. Section IV contains our simulation results on the inhomogeneous condensation of the scalar and pseudo-scalar bilinears and their interrelation. We calculate the phase diagram in the (T, µ) plane for various lattice sizes and lattice constants in order to study the thermodynamic and continuum limits. We shall see that even for the smallest accessible value N f = 2 the results resemble those for the exact solution of the system with N f → ∞. Towards the end we exploit Dyson-Schwinger equations to study the U A (1)-invariant fermion four-point function in the infrared.

A. Symmetries and reformulations
The chiral GN model with Lagrangian (1) most prominently features a global axial U A (1) symmetry, with a continuous parameter α ∈ R. In this work we denote spacetime coordinates by bold letters, for example The continuous axial symmetry is to be compared with the discrete Z 2 symmetry of the model considered in [24,25]. Further symmetries of the model include a flavor-vector symmetry that ensures the factorization of the fermion determinant, parity and charge conjugation symmetry responsible for the absence of the sign problem for even N f (see [24] for details) and, of course, Euclidean spacetime symmetry.
As is usually done we introduce the complex auxiliary field ∆ in order to bring the Lagrangian (1) to the equivalent form where P ± = 1 2 (1 ± γ * ) are the chiral projectors. This Lagrangian is invariant under the axial transformations (2) supplemented by One can show the equivalence of Lagrangians (4) and (1) by using the equations of motion for the auxiliary field ∆. This equivalence persists on the quantum level because the ∆ integration in the path integral is Gaussian and can be done analytically, leading back to Eq. (1).
It is no more difficult to obtain the following Dyson-Schwinger (DS) equations relating the expectation values of the auxiliary fields to the symmetry-breaking chiral condensates: 2 For later use we introduce two further parametrizations of ∆ in terms of its real and imaginary parts σ and π and in terms of its absolute value ρ and phase θ: In order to study finite baryon densities we also introduce a chemical potential µ for the fermion number densitȳ ψγ 0 ψ, such that the Lagrangian takes the form where the Dirac operator D is defined as It is understood that this operator acts on all flavors in the same way, such that in the multi-flavor case we may use the same symbol as for one flavor. While there is no gauge invariance in this model, one can still trade the compact field θ for an imaginary vector potential in the following sense: where the covariant derivative D µ is defined as Since the main focus in our work is on homogeneous and inhomogeneous phases of the finite-temperature and finite-density cGN model we impose that ψ,ψ are antiperiodic and ∆, ∆ * are periodic in Euclidean time with period β, where β is the inverse temperature. We furthermore impose that all fields are periodic in the spatial direction with period L.
Integrating out the fermions in the partition function yields an effective bosonic theory in which the auxiliary bosons become dynamical via fermion loops, with the effective action We used that the fermion determinant of the multi-flavor model is (det D) N f with the one-flavor operator D appearing in Eq. (14). A convenient (and widely adopted) way of renormalizing this formal expression is a choice of the bare coupling g 2 such that S eff for T = 0 and µ = 0 takes its global minimum at some prescribed positive value ρ(t, x) = ρ 0 . The corresponding gap equation in the thermodynamic limit, yields the cutoff dependence of the bare coupling.

B. Large-N f results
In the large-N f limit the saddle-point approximation to the path integral (13) becomes exact and the grand potential Ω proportional to the minimum of the effective action (14) on the space of auxiliary fields, This means that in the large-N f limit the path integral is localized at the minimizing configuration ∆ min . It follows, for example, that the expectation value of ∆ is equal to ∆ min . The condition of a (local) minimum, maximum or saddle point is expressed by the gap equation and for all µ the cGN model (in the large-N f limit) is in a symmetric phase with a vanishing condensate field [31]. Here γ is the Euler-Mascheroni constant. More surprising is the fact that below T c and for all µ = 0 there are no homogeneous solutions of the gap equation which minimize S eff . Instead, the minimizing configurations are helixes with pitch π/µ, so-called chiral spirals, with a temperature-dependent amplitudeρ(T ) and k(µ) = −µ in the large-N f limit. For vanishing chemical potential the chiral spiral degenerates to a homogeneous configuration, which relates to the large-N f solution of the Z 2 GN model at µ = 0. We conclude that the profile functionρ(T ) is just the condensate of the GN model at µ = 0, which decreases monotonically in T until it vanishes at T c . The large-N f phase diagram in the (T, µ) plane is depicted in Fig. 1. The phase diagram of the cGN model in the large-N f limit (see [18]). One critical temperature Tc for all µ marks the transition from chiral spirals at T < Tc to the symmetric phase at T > Tc. Units are set by the condensate ρ 0 at zero temperature and zero chemical potential.

C. Spontaneous symmetry breaking in low dimensions
Under rather natural assumptions the existence of perfect long-range order (as opposed to quasi -long-range order) in lower dimensions is excluded by the celebrated Coleman-Hohenberg-Mermin-Wagner theorem [32][33][34][35]. This theorem states that continuous symmetries cannot be spontaneously broken at finite T in low-dimensional systems with short-range interactions. In particular, for zero-temperature systems the theorem says: The continuous symmetries cannot be spontaneously broken in (1 + 1)-dimensional quantum systems. If a continuous symmetry were spontaneously broken, then the system would contain Goldstone bosons, which is impossible in two spacetime dimensions because massless scalar fields have an IR-divergent behavior [35]. Discrete symmetries, on the other hand, can still be spontaneously broken in two dimensions.
There is a domain-wall proof of the theorem, of which the basic intuition is to rotate the field or, in a spin-model language, the values of spins in a finite region with an arbitrarily small energy cost. This is achieved by creating a domain wall of finite thickness interpolating between the regions with rotated and unrotated spins. If the symmetry group were discrete, there would be no smooth interpolation and hence a finite cost for creating domain walls.
The increasing strength of fluctuations (thermal and quantum) in the IR with decreasing dimension d is known from the (Euclidean) free scalar field with propagators The interpretation of the IR divergence in d = 1 and 2 is that the field fluctuations cannot stay centered around a mean. It implies that far away from a given spacetime point the field takes completely different values than at the given point. This happens in one and two dimensions where the fluctuations move the field arbitrarily far from an initial value such that it has no well-defined average. This reasoning should apply to translation invariance as well: If the distance between two neighboring particles on a wire fluctuates by δx, then the nth particle's separation fluctuates as √ n δx and thus diverges for large n.
These large fluctuations destroy any long-range order in the position of the particles and R. Peierls concluded that a one-dimensional equally spaced chain with one electron per ion is unstable [36]. In higher dimensions (d ≥ 3) the fluctuation-induced correlations fall off at large distances and are not strong enough to destroy long-range order. Furthermore, based on the powerful energy-entropy argument it has been argued that any spontaneous symmetry breaking (SSB) should be disallowed in 1 + 1 dimensions at finite temperature [32]. In the argument one considers a small number N of local perturbations of an ordered state (e.g. aligned spins in the Ising model). The entropic contribution of these perturbations to the free energy is ∝ N ln N while the energy penalty is only ∝ N . Thus, the entropic contribution can overcome the energy barrier and destroy the order. This perspective is directly applied to the discrete GN model in [37].
Hence, the breaking of translation invariance in the (1 + 1)-dimensional GN model seems to be excluded on general grounds. On the other hand, the no-go theorems do not apply in the large-N f limit where the analytical solution shows that the finite-temperature and finite-density equilibrium state is not translation invariant. What may happen at finite N f is a subtle issue and has been discussed (including the underlying assumptions of certain no-go theorems) in [24].
Besides the questions raised in [24] there are more points to be clarified with regards to the applicability of the no-go theorems: It is not obvious whether the effective action S eff [∆] containing the non-local fermion determinant is short ranged enough to ensure the convergence of certain integrals, which is assumed in [33]. Although [34] treats fermions as well, the result is based on sufficient convergence (in form of f sum rules) and gives itself an example of violation.
We emphasize that the no-go theorems allow for a BKT phase with quasi-long-range order expressed by slowly decaying correlations ∝ 1/|x | α and a BKT transition to a massive phase with short-range correlations ∝ e −m|x | [38,39]. There is no symmetry breaking and no order parameter involved in the strict sense, but the slowly decaying correlations of a BKT phase allow for large regions of one distinguished local state.

D. Perturbations of chiral spirals
How are the inhomogeneous phases of the GN and cGN models in the large-N f limit compatible with the no-go theorems discussed above? In a way the large parameter N f takes over the role of an extra spatial dimension. For example, in the domain-wall argument given above the energy penalty is multiplied by the large number N f and in the limit N f → ∞ may overcome the entropy gain.
This and further heuristic arguments can be substantiated by a systematic expansion in the small parameter 1/N f , whereby one assumes that for finite N f the continuous U A (1) axial symmetry is not spontaneously broken. Under an axial rotation the radial field ρ is left invariant and θ is shifted by a constant. This means that an invariant effective action is a functional of the form [40] This effective action is used to calculate expectation values of functions of ∆ = ρe iθ and its complex conjugate field ∆ * . However, in the continuum model a condensate ∆ cannot form (it would break the axial symmetry) and with chiral SLAC fermions and the ergodic rHMC algorithm it averages out in lattice simulations, see Sec. III B. Thus, following our previous studies [24,25], the correlator will be of particular interest to us. For N f → ∞ the path integral is localized at the chiral spiral (19) and we find Clearly, for finite N f we must admit small deviations from the chiral spiral, and expand the effective action in powers of the fluctuation fields δρ and δθ. An explicit calculation at zero temperature and in an infinite volume shows that the term linear in the fluctuation fields vanishes if the bare GN coupling depends onρ according to The first relation is recognized as the gap equation of the Z 2 GN model. For large volumes the wave number k becomes continuous and the second relation can be fulfilled for all µ. Since the effective action only depends on k via k + µ, this relation implies that S eff is independent of both k and µ. In a finite box with quantized k, however, the background fieldρ and the effective action will generically depend on k + µ. The contribution quadratic in the fluctuation fields is rather lengthy and has divergent terms which all cancel when one uses the gap equation (25). If in addition the wave number of the chiral spiral obeys k + µ = 0, then one obtains where the dots indicate higher-order terms, the integrals extend over the spacetime volume and we introduced the (non-local) operator K ∆ containing the Laplace operator, In a low-energy approximation we may perform the gradient expansion, which yields the simple expression containing the standard kinetic terms plus higher derivative terms. The first term under the first integral is just the second-order term in the expansion of U eff (ρ + δρ) in powers of δρ. Thus, up to second order the effective action for ρ =ρ + δρ and δθ at low energies has the form where we inserted the explicit form of the effective potential at zero temperature and the dots indicate cubic and higher-order terms and higher derivative terms. We see explicitly that ρ describes a massive field and δθ a massless would-be Nambu-Goldstone mode. At large N f the latter decouples from the system while at finite N f it destroys perfect long-range order.
To study long-range correlations we can safely neglect contributions from the massive field and obtain, for large but finite N f , the valid approximation It holds information about the dominant wave numbers of typical configurations in an ensemble. Due to the logarithmic divergence in the correlator of the massless scalar field one finds for x → ∞ such that in a BKT phase with quasi-long-range order the amplitude of the oscillating correlator decays fairly slowly, following a power law. At finite temperature the correlation length, given by where K 1 denotes a modified Bessel function of the second kind, is finite. The coefficient α increases monotonically with the inverse temperature β from α = 0 to α = 1. This means that the correlation length diverges in the large-N f limit or for T → 0.

A. Objectives and observables
The previous discussion makes clear that we should not expect to see SSB in the cGN model with U A (1) symmetry. Indeed, there are stronger arguments against perfect long-range order in this model than in the GN model with Z 2 symmetry. However, the difference between a spontaneously broken and a BKT phase at zero temperature most likely appears on exponentially large length scales that cannot be reached in our lattice simulations, see, for instance, [41]. It could very well happen that on physically relevant length scales one can hardly distinguish between quasi-long-range and perfect long-range order. Furthermore, we shall see that even the distinction between a massive symmetric phase and a BKT phase at low temperatures is non-trivial if one allows for contributions of the first excited state.
Either way we will find striking similarities between the cGN model with only two flavors and the model with N f → ∞, which, for µ = 0, shows SSB of translation invariance. If similar observations apply to more realistic models in higher dimensions then this could be relevant for the physics of compact neutron stars, heavy-ion collisions or condensed matter in small systems.
We shall see that for 8 flavors the correlation function C(x) in (22) has the form (30) predicted by the effective low-energy Lagrangian and can be hardly distinguished from the large-N f result (23). For example, at low temperature its discrete Fourier transform F[C](k) is peaked at the dominant wave number which for large N f is given by the chemical potential, while for N f = 2 we find k max < µ. Notice that we have included a factor of 1/2 in Eq. (33), in line with the introduction of k in Eq. (19) as half the wave number. We will use this convention for k and k max in the remainder of this work. The spatial correlation function C(x) also encodes the distinction between the massive symmetric and BKT phases in its decay properties, For a comparison we included the asymptotic behavior in a symmetry-broken phase. The temperature-dependent correlation length ξ β was defined in Eq. (32).

B. Lattice setup
We discretize two-dimensional Euclidean spacetime to a finite lattice with N s and N t lattice sites in the spatial and temporal directions respectively, such that L = N s a is the spatial extent, T = 1/N t a is the temperature and a denotes the lattice constant.
In our simulations we employ chiral SLAC fermions [42,43], which discretize the dispersion relation in momentum space, leading to a non-local kinetic term in position space. They have proven advantageous over other fermion discretizations for low-dimensional fermionic theories, see e.g. [24]. The use of SLAC fermions restricts N s to be odd and N t to be even. For further details we refer to sections 2.1 and 4.1 of [26]. Note that our lattice setup is the same as in [24], with the only difference that besides a scalar field σ we now have an additional pseudo-scalar field π and both fields enter the complex condensate field ∆ via Eq. (7).
For an easy comparison with the analytic results we express physical quantities in units of the expectation value ρ at T = µ = 0, denoted by ρ 0 . This is analogous to the scale σ 0 in our previous studies [24,25]. One should stress that this neither assumes any form of symmetry breaking nor is in conflict with any no-go theorem because a non-vanishing expectation value of ρ does not break any symmetry. Fig. 2 shows histograms of x ∆(x ) in the  complex plane for µ = 0 and 12 different temperatures. For these histograms we used ensembles with O(10 4 ) configurations each. We clearly observe that the distribution of ∆ is angle independent or U A (1) invariant. At low temperature it is ring shaped with its maximum at ρ > 0, while at high temperature it turns into a Gaussian-like distribution and the maximum moves to ρ = 0.
In the actual simulations, however, the quantity ρ 0 is surprisingly hard to determine. App. A sheds some light on the details of this procedure. In summary, we used 3 with τ = 1, . . . , N MC enumerating the Monte Carlo (MC) configurations. This yields a good signal at low temperatures where ρ is required. For most of our simulations, we used one of three different spatial extents N s = 63, 127, 255 and lattice constants aρ 0 ≈ 0.46, 0.19, 0.08 in order to study both the continuum limit and the infinite-volume limit. We vary the temperature by changing the number of lattice points in the temporal direction, N t , at fixed a and we vary a by changing the coupling 1/g 2 in Eq. (8). For these lattices we map out phase diagrams in the (T, µ) plane. More details as well as a table summarizing all parameter settings are given in App. D.
Experience with interacting fermion models teaches us [24] that systems exhibiting (quasi-)long-range inhomogeneous structures can have rather long thermalization times when running simulations with randomly generated initial configurations, e.g., using a Gaussian distribution. As a way to counteract this problem, we employ a different approach for the majority of results presented in this work and perform a systematic "freezing out" in the following way: Starting at high temperatures with N t N s , where thermalization times are not an issue, we generate at least 1000 configurations to ensure proper thermalization. We then map the last of these configurations to a lower-temperature lattice with N t > N t by simply repeating the data in the temporal direction and use it as a seed configuration on the larger lattice. This reduces the thermalization period (where no measurements are performed) if the temperature step is small. In our simulations we systematically approach lower and lower temperatures using this "freeze-out" procedure. This way we experienced significantly less "getting stuck" in some far from typical configurations, although it could still not be completely prevented from happening.
A cross-check with thermalized results using standard Gaussian-distributed initial configurations yields equivalent results, with the "freezing" method having noticeably better thermalization properties and thus overall smoother results. As an additional cross-check we also performed the inverse procedure, i.e., a "heating", for a handful of parameters in order to exclude any hysteresis effects caused by the "freezing" method.
As can be seen from Fig. 3, where we show the Fourier transform of C σσ (x) (to be defined in Eq. (37)) computed via each of the three methods, the "frozen" and "heated" results agree very well, indicating that hysteresis effects are negligible. The fact that the "independent" results, i.e. the ones obtained by using Gaussian-distributed initial configurations, show some deviation is likely to be attributed to their lower statistics and worse thermalization properties. The temperature is the second lowest considered, i.e. we compare several "freezing" steps with a single "heating" step. The vertical lines indicate the maxima of the "frozen" results at the lowest temperature.
The vertical lines in Fig. 3 show the peak positions, which were estimated with the "freezing" method, for the lowest temperature considered. We see that for the highest density (µ/ρ 0 ≈ 1.31 in the figure) the peak positions of the two lowest temperatures differ. This dependence on temperature is not seen in the large-N f limit and is caused by bosonic fluctuations.
For small µ on our smallest lattice, where homogeneous configurations dominate, we furthermore compared with a "cold start", which amounts to starting the simulation from ∆ (0) (t, x) = 1 + i. Again we found matching results except for the lowest temperatures where the cold start is expected to suffer from severe autocorrelation problems. A detailed analysis of autocorrelation effects can be found in App. B.

C. Lattice estimators
We have argued that spatial correlation functions are useful tools to probe for inhomogeneous phases since they avoid the destructive interference one would encounter when directly calculating chiral condensates on the lattice. We consider the two spatial correlators where the sums extend over all lattice sites and · denotes the Monte Carlo average. If these correlators show an oscillating behavior, one can infer the existence of inhomogeneities. The unbroken U A (1) symmetry (5) implies that for any temperature and chemical potential Also note that the fermion determinant is invariant when σ and µ both change their signs, such that from which we conclude that We see that additional correlators that arise from interchanging σ ↔ π in Eq. (37) are not independent and we refrain from using them in subsequent equations to save some space. In the measurements, however, we do not implement the symmetries (38) and instead use all four correlators C σσ , C σπ , C πσ and C ππ in order to reduce statistical correlations. From Eq. (22) one obtains and the property (40) means that C is real for vanishing µ.
In [24] we introduced the minimal value to map out the entire phase diagram of the (discrete) GN model. This parameter is negative if there is (quasi-) long-range order with oscillating C σσ (x) and is also useful for discussing the physics of the chiral GN model. For the chiral model the choice of C σσ might seem arbitrary but because of (38) any quadratic correlator of a linear combination of σ and π would serve the same purpose.
It is important to note that taking the minimum is a global operation that disqualifies this quantity as a local observable. Furthermore, this minimum might (and actually commonly will) be taken for small spatial separations x. In such cases, C min does not probe the longrange behavior of the system.
We estimate the dominant wave number k max as given by Eq. (33), but calculated from C σσ instead of C. Sometimes we quote results in terms of the integer-valued dominant winding number ν max , related to k max via From analytical studies [40,44] it is expected that the U A (1)-invariant fermionic four-point function of the GN model, at zero temperature and zero fermion density should have a power-law behavior in the limit of large separations, where c is some constant. Similarly to the spatial correlation functions (37) for the condensate fields we introduce the spatial correlation function for the N f fermionic lattice fields, The asymptotic form (45) would imply a power-law decay Dyson-Schwinger equations (see App. C) relate C 4F (x) to the spatial correlation functions of the condensate fields, since the contact term in Eq. (C5) does not contribute for large x. Since C σσ and C σπ are easily accessible in lattice simulations we shall use this relation to study the infrared properties of C 4F . For µ = 0 the latter is real, see (40). From the effective low-energy approximation outlined in Sec. II D we expect that the phase of the complex condensate field, θ = arg(∆), holds important information about the existence of inhomogeneous structures. We thus studied the space dependence of its expectation value, defined in the following way: where the bar indicates time averaging, i.e., We chose this (unusual) order of time-and MC averaging to suppress statistical uncertainties. Although the two averages in (49)

IV. NUMERICAL RESULTS
In previous studies of the discrete GN model [24,25] To be more precise, if the BKT scenario were correct, then, for instance, in order to obtain a decay to half the amplitude a crude estimate using C(x) ∼ |x| −1/8 yields at the very least. Thus, in order to make any meaningful statements about such an amplitude decay we would require around O(10 3 ) lattice points at sufficiently small temperature (large temporal extent). This does not take into account severe autocorrelation problems, finite-size effects and contributions from excited states that might all spoil the signal. This crude estimate motivated us to study the long-range behavior for N f = 2 in [24], for which the same estimate yields feasible 40 lattice points.
Although our focus is on 2 flavors we performed one parameter scan in (T, µ) for N f = 8, N s = 63 and aρ 0 ≈ 0.41 in order to compare with results for the discrete GN model. Some of our results are depicted in Fig. 5. Fig. 5a shows the phase diagram extracted from C min (see Eq. (42)), which is to be compared with Fig. 1 for infinite flavor number and is also the equivalent of Fig. 7 in [24]. We see that the phase diagram agrees well with the large-N f prediction for small chemical potential (µ < 0.5ρ 0 ) and at least shows the predicted structure at larger µ. At vanishing chemical potential C min is positive for small temperatures, indicating predominantly homogeneous configurations with non-vanishing amplitudes. They relate to the homogeneously symmetry-broken phase at large N f . In Fig. 5b we see that in this regime C σσ (x) is a positive and monotonically decaying function (blue curve) and C σπ (x) ≈ 0 in agreement with (40). Raising the temperature we find a small temperature regime around T ∼ 0.3ρ 0 where we observe a sudden drop of the amplitude such that the µ = 0 data mimic a second-order phase transition. In the hightemperature regime the non-zero correlator C σσ falls off even more rapidly. This (would-be) transition temperature at µ = 0 is approximately equal to the one found in the discrete GN model in [24]. This was to be expected since in the large-N f limit the GN and cGN models at vanishing chemical potential have the same critical temperature. It is also not surprising that for N f = 8 the transition temperature is significantly lower than in the large-N f limit (cf. Eq. (18)), where quantum fluctuations are suppressed. The symmetric high-temperature regime at µ = 0 extends to non-vanishing chemical potential (orange curve in Fig. 5b).
At low temperature and non-vanishing fermion density we can clearly confirm that the dominant contributions to the path integral come from chiral-spiral-like configurations. An example of this is shown in Fig. 5b (green curve). Such configurations are the cause of the large region of negative values in Fig. 5a. The transition line to the region where oscillations are no longer dominant is roughly a line of constant temperature for small chemical potential (µ < 0.5ρ 0 ), as expected from the large-N f solution. For large chemical potential it tilts upwards unexpectedly, thereby enlarging the regime where inhomogeneities are found. This effect was also observed in [24] for N f = 2 in the discrete GN model and is related to short-and intermediate-range phenomena that will be discussed later. Nevertheless, the fact that we encounter it already for N f = 8 strengthens the point that quantum fluctuations are much stronger in the chiral model compared to the discrete one. For N f = 8 the winding numbers (43) of the inhomogeneous configurations match the large-N f expectation very well if one accounts for the discretization of the wave number due to the finite box size, as can be seen in Fig. 5c (note that ν max is integer valued by definition). As in [25] there is a tendency for the lattice data to lie slightly below the N f = ∞ expectation. The linear fit through the origin yields a slope of roughly 7.91, which is lower than the large-N f value L/π ≈ 8.27, but well within the expected accuracy of the large-N f expansion We remark that autocorrelations appear to be under control. However, due to limited statistics we cannot rule out the existence of another, larger, autocorrelation scale at low temperatures, see App. B for details.  To monitor the fluctuations in the system at finite temperature and density, we measure the dominant wave number (33) of the equilibrium ensemble. It characterizes important configurations for the given set of control parameters and tells us which chiral spiral is favored in the rough landscape defined by the effective action with its many local minima. This analysis presupposes that chiral spirals are the dominant configurations even for N f = 2 or that the non-dominant winding numbers are suppressed. We shall see that this is a valid assumption at small temperatures. Fig. 6 shows such histograms for 3 values of T and 3 values of µ. As expected, the data show three distinct peaks, one for each value of µ. At the lowest temperature and µ = 0 the peaks are pronounced with over 80% of the configurations sharing the same dominant wave number. Increasing the temperature then broadens the peaks. Concerning the question of spontaneous symmetry breaking, one should stress three features: 1. While the peaks flatten significantly for rising temperature, they do not vanish completely. At temperatures as high as ∼ 0.5ρ 0 we could still make out small bumps (not shown in Fig. 6). Thus, even at these high temperatures the system knows about the inhomogeneity arising from its finite density.
2. There is no qualitative (or even sudden) change in these distributions that could characterize a second-order phase transition. Instead the flattening of the peaks is a rather smooth process.
3. Even at low temperatures (e.g. T ≈ 0.05ρ 0 ), well inside the would-be symmetry-broken regime, the contributions from concurrent frequencies are significant (around 10 − 20% in the example). The interference of these contributions is the mechanism which, crucially, prevents a breaking of symmetry.
The features 1 and 2 discussed above are clearly visible in the spatial correlators depicted in Fig. 7 T and non-vanishing µ we clearly observe remnants of a chiral spiral, see Fig. 7a. From Fig. 7b we see that both correlators are oscillating with a phase shift of π/2. This is the parameter region where there are sharp peaks in Fig. 6. At higher temperature the peaks flatten and the correlators show damped oscillations as shown in Fig. 7c. Notice, however, that even in this regime we still find C min < 0, i.e., this is classified as a region of spatial inhomogeneities according to our definition. Here we observe a clear deviation from the large-N f solution. Since the oscillations in Fig. 7c are only seen on short scales we must be cautious when interpreting a negative C min as a signal for inhomogeneities. As already stressed before, a negative C min ensures that there are oscillations on some length scale but this scale can be -and certainly is for large parts of the blue region of the phase diagram -a short or intermediate one. deviate considerably from those in the large-N f limit and those for N f = 8, cf. Fig. 5c. One might conjecture that the winding numbers decrease with decreasing N f .

C. Phase diagram for N f = 2
One could expect a qualitatively different "phase diagram" for the cGN model with 2 flavors as compared to the large-N f diagram depicted in Fig. 1. In order to test this expectation we calculated C min defined in Eq. (42) on a grid in the space of control parameters µ and T on lattices with N s = 63, 127 and 255 lattice points in the spatial direction. We studied both an infinite-volume extrapolation at (approximately) fixed lattice constant aρ 0 ≈ 0.46 and a continuum extrapolation at (approximately) fixed physical volume.
The diagrams for systems with constant lattice spacing in Fig. 9 show that the infinite-volume limit significantly shrinks the red (C min > 0, i.e. predominant homogeneous contributions) region without affecting the blue and green region of predominant inhomogeneous resp. symmetric configurations. There are three rather different phenomena at work here: 1. The simplest one is just geometrical: When we enlarge the spatial volume, we can fit larger wavelengths into the finite box. For small µ the pitch of the chiral spiral would exceed the box size, which means that the chiral spiral does not fit into the box. Such a suppression of chiral spirals with large pitches disappears for larger volumes. Hence, the region of predominant homogeneous configurations must shrink in the direction of non-vanishing µ.
2. Finding a shrinking of the C min > 0 region in the temperature direction is clear evidence against spontaneous symmetry breaking. In fact, the (qualitative) behavior of the apparent transition temperature and the condensate is well understood by a comparison of the analytically known correlation length Eq. (32) with the box size. We can thereby clearly identify the remnant condensates that were measured as finite-size effects.
3. The transition from the blue (C min < 0, i.e. predominant inhomogeneities) to the green (C min ≈ 0, i.e. predominantly symmetric) regime can be easily understood as the following short-range effect: At finite temperature there are two length scales in the system (besides the finite box size), i.e. the temperature-induced finite correlation length ξ β from Eq. (32) and the predominant wavelength inversely proportional to µ (up to discretization due to the finite box size). Obviously, for C min to be negative, the amplitude, which decays due to ξ β , must not have dropped to (almost) vanishing values at separations where the first minimum of the oscillations occurs. Since the latter is given by µ (up to a constant factor), the transition line from blue to green signals that the temperature scale takes over as the shortest relevant scale from the chemical potential. This is not really a qualitative change. As this is independent of the much larger box size, it is not affected by the infinite-volume limit.
An interesting, but unfortunately hard to quantify observation is the following. While on smaller lattices (e.g. N s = 63) the data tend to fluctuate around only one background configuration, like a chiral spiral with a fixed winding number, larger lattices admit changing the winding number more often as opposed to less often. An example is depicted in Fig. 11, which shows a time series of the modulus ρ (τ ) of the spacetime average of∆ For most of the MC time ρ (τ ) fluctuates about a constant value. During this time the real part C σσ of C defined in (22) is almost constant and the imaginary part C σπ is small (recall that C σπ = 0). But at several MC times, e.g. τ ≈ 1100 or 3860, the field ρ (τ ) drops and the real and imaginary parts of C(x) show the profiles typical for a chiral spiral. While the lower left plot in Fig. 11 is representative for most of the configurations, the sudden drops in the time series are strongly correlated with the appearance of inhomogeneous configurations as seen in the lower right panel. That this is only seen on large lattices is counterintuitive at first since autocorrelation times usually increase with the system size and it is also the opposite of what was observed for the Z 2 GN model during the work on [24,25]. However, the fact that considerable phase fluctuations on large scales are allowed is the analytically predicted mechanism for avoiding spontaneous symmetry breaking, see Sec. II. From that per-spective, it supports the analytical claims. Obviously, Fig. 11 showcases a large autocorrelation time τ , which, however, is under control due to good statistics of the order of 20τ .
The phase diagrams for systems with approximately constant physical volume and successively smaller lattice spacing are shown in Fig. 10. As can be seen, we find inhomogeneities 5 for all our lattice spacings and the results are consistent with their existence in the continuum limit. Unfortunately, setting the scale in a partially conformal system is a very subtle issue as the dominant fluctuations have no scale at all (at zero temperature). Since the details of this scale setting procedure are highly non-trivial (see App. A), we must leave a more detailed analysis to a future publication.
As is discussed in detail in App. B, our simulations suffer from large autocorrelations. For a large region in parameter space on all geometries these autocorrelations are under control due to sufficient statistics. However, close to the critical line at T = 0 autocorrelation times diverge and the shown data should only be regarded as qualitative in the sense that they surely capture the important phenomena found in large but finite regions of space but might be off quantitatively due to autocorrelation-related suppression of subdominant local minima. However, the discussion of App. B makes it clear that these will not affect our conclusions.
We conclude that we observe inhomogeneous structures in the cGN model with only 2 flavors -similarly as in the large-N f model. The notable difference is that -in accordance with pertinent no-go theorems -these are incoherent on sufficiently large scales thereby hindering spontaneous symmetry breaking. A comparable study of the Z 2 GN model with 8 flavors in [24] led to a similar conclusion: inhomogeneous structures persist in the infinite volume limit. We cannot say for certain whether this remarkable feature survives the continuum limit of the cGN lattice models.

D. Decay properties of C4F
We analyzed the decay properties of C 4F (x) as given by Eq. (48) on a 72 × 63 lattice with aρ 0 ≈ 0.46, i.e. at a temperature T ≈ 0.03ρ 0 . In order to study its infrared behavior we computed the connected correlation function. Motivated by the asymptotic forms (31) predicted by the low-energy effective action we fit the data points via a (symmetrized) algebraic function as well as a double-cosh ansatz, and show the results in Fig. 12.  (48)) with algebraic and exponential fits whose fit parameters are given in (54) and (55).
The fit parameters for a power-law decay are These results confirm similar findings obtained for the Z 2 GN model, namely that it is very difficult to distinguish between power-law and exponential decays on the lattices with N s = 63, which was also to be expected following the previous discussion and [41]. However, from the perspective of our analytical discussion, where we predicted a massive phase for any T > 0 with the mass vanishing in the limit T → 0, there is a very well-motivated explanation. Eq. (32) predicts which is reasonably close to the fitted value (remember that we expanded in O(1/N f ) ∼ 50% for N f = 2) and explains its seemingly unnaturally small magnitude.
On the other hand, we find that the value β ≈ 0.52 is only marginally larger than the theoretical zerotemperature prediction of β = 0.5 for two flavors in Eq. (35). This result -although not precisely a proof -is in astonishing agreement with the analytical prediction coming from an expansion around N f = ∞ 2 and furthermore beautifully reveals how the massive phase more and more approximates the conformal behavior at zero temperature by the unexpectedly large mass ratio m 2 /m 1 ≈ 10.

E. The phase field θ
In this section we analyze θ(x) . This discussion should be regarded as complementary to the previous analysis of the correlators in the sense that we now use a quantity directly related to the fields. It will further substantiate our previous findings.
We show the x dependence of the average θ(x) , as defined in Eq. (49), on a 72 × 63 lattice for aρ 0 ≈ 0.46 and for three values of the chemical potential in Fig. 13. For vanishing µ the argument of the averaged complex condensate field ∆ is constant, which means that the latter does not wind. For the intermediate value µ ≈ 0.88ρ 0 the average angle is an almost linear function of x and the complex condensate winds 6 times when one moves along the spatial direction. When one further increases the chemical potential to µ ≈ 1.31ρ 0 , the slope of the (almost) linear mapping x → θ(x) increases and the condensate winds 9 times. We see that the winding number of the chiral spiral around the spatial direction increases with increasing µ. In fact, the winding number extracted from the averaged field θ(x) perfectly agrees with the dominant winding number defined in Eq. (43) and depicted in Fig. 8. To summarize: at low temperature θ(x) calculated from (49) is almost a linear function of x with a slope proportional to µ. In agreement with the analysis based on the dominant wave number we detect a chiral spiral with a winding number proportional to µ.

V. CONCLUSIONS
In the present work we studied the (1 + 1)-dimensional chiral Gross-Neveu model with chiral SLAC fermions and exact axial U A (1) symmetry on the lattice. Our two main results are summarized in the following.
First, we have strong and multiple evidence that the analytical prediction from an expansion in 1/N f well describes the qualitative features of the cGN model with 2 flavors. As expected we see no spontaneous symmetry breaking with long-range order in the strict sense, and our results suggest that at T = 0 there is a BKT phase with quasi-long-range order. For example, the low-temperature regime, where we have signals of (in)homogeneous ordering over the whole lattice, is well explained by the analytically predicted correlation length exceeding the finite box size and it shrinks consistently in the thermodynamic limit. Additionally, the decay properties of pertinent correlation functions are well fitted by the analytically predicted ansätze with reasonable parameter values.
The latter suggests that for T = 0 fluctuations of the phase field θ on all scales exist and are responsible for the restoration of the axial symmetry. We demonstrated that θ is uniformly distributed on unit circles in complex field space and that large system sizes allow for long-range phase fluctuations strong enough to change the winding number. This behavior is predicted by the effective lowenergy theory for θ which has been taken from [44] and extended to µ = 0.
Despite this, our second finding is that, rather unexpectedly, our simulations at finite temperature and density reveal that the cGN model with only N f = 2 flavors resembles the analytic large-N f solution in many ways. The chiral spirals are still seen in the dominant configurations and their winding numbers increase linearly with the chemical potential. The only qualitative difference at low temperatures is that these structures are only coherent in finite but -depending on the temperature -potentially very large regions of space. Instead of a temperature-driven phase transition at intermediate temperatures, we found a competition of the two important scales in the system, viz. the temperature-induced finite correlation length and the density-induced wavelength. So, the question whether or not oscillating behavior was observed (on potentially short scales) can be answered only from comparison of the wavelength with the correlation length. Or, put differently, it is very likely that oscillating behavior can be found for any temperature and non-vanishing chemical potential as long as the wavelength is shorter than the correlation length of the system. This is qualitatively different from the large-N f behavior where there is a strict critical temperature above which no oscillation can be observed.
We have verified these results mainly via the correlator C in (22) and by analyzing the phase of the averaged field ∆, defined in (49). We generated many ensembles for the control parameters T and µ on grids with up to 192 points. To quantify finite-size and discretization effects the simulations were repeated on lattices with 63, 127 and 255 points in the spatial direction. While we have good signals for the behavior in the thermodynamic limit, whether the inhomogeneities remain after the continuum limit has been taken is less clear. With the chosen scale setting, which is a subtle issue in a theory with quasi-long-range order, we observe that inhomogeneities remain in the limit aρ 0 → 0. We hope to gain a more thorough understanding of this limit in the future.
Although we found strong evidence that consistently supports the analytical predictions, our method of MC simulations will never be able to prove this in a rigorous sense. Therefore, it would be interesting to compare our findings with results from other methods, for example the functional renormalization group. It would be valuable to continue the study of the (1 + 1)-dimensional Gross-Neveu-Yukawa model in [45] to related systems in finite volumes and inhomogeneous background fields.
The mechanism of how the cGN-model realizes the U A (1) symmetry is similar to the flattening of the constraint effective potential for a spacetime-averaged order parameter∆ [46]. For example, in the Ising model at low temperature, if we impose that the spatially averaged spin vanishes in the sum over spin configurations, then in a typical configuration we observe large regions with spin up and large regions with spin down. Despite the surface energy stored in the walls separating the "up" and "down" regions, this is the energetically preferred way of fulfilling the external constraint.
Models with a continuous symmetry react differently to the constraints. For example, in the 3-dimensional O(2) model with a Mexican hat potential for a complex scalar field ∆ the constraint |∆| < |∆| is met by inhomogeneous spin-wave-like configurations with |∆(x )| ≈ |∆| [47]. These configurations resemble the chiral spiral in the cGN model, for which the modulus of∆ can be much smaller than |∆| . In the 2-dimensional cGN model the constraint∆ ≈ 0 is not imposed by hand but by general theorems which ensure that ∆ = 0. In a typical configuration the modulus of ∆(x ) is near the minimumρ of the effective potential -in order to minimize the bulk energy -but the real and imaginary parts σ and π have vanishing expectation values caused by large phase fluctuations about the relevant chiral spiral. The main difference between the 3-dimensional O(2) model and the 2-dimensional cGN model is that in the former model the wavelength of the inhomogeneity is given by the box size [47] and in the latter by the inverse chemical potential.
In [48] it has been emphasized that the occurrence of correlation functions exhibiting damped oscillations in the spatial directions is directly related to particular features of the dispersion relations. The associated quantum spin liquid behavior, which we also spotted in the 2-flavor cGN model, may thus be observed in a larger class of field theories.
After publishing the initial draft of our manuscript a similar study of the (1+1)-dimensional cGN model using the naive fermion discretization was published in [49]. Its results are in qualitative agreement with ours.

ACKNOWLEDGMENTS
We thank Björn Wellegehausen for providing the code base used in the present work and for fruitful discussions. In addition, we thank Laurin Pannullo, Marc Wagner and Marc Winstel for many discussions on four-Fermi theories and the collaboration during a previous work on the Z 2 GN models. This work has been funded by the Deutsche Forschungsgemeinschaft ( For an easy comparison of our results with the analytic large-N f solution we use ρ 0 = ρ T =µ=0 to set the scale. Unfortunately, it is difficult to obtain an accurate estimate for ρ in our simulations. In this appendix, we first explain the (statistical) problems with direct approaches to measure ρ and afterwards present our solution.
From a field theory perspective, the direct lattice estimator for ρ would be ρ t,x for any (fixed) point (t, x) on the lattice. Now, ρ t,x should be homogeneous up to fluctuations and, hence, one can improve the statistics by combining the data from all estimators ρ t,x for all lattice points. Example data for this estimator can be seen in Fig. 14a.
The final estimate for ρ would then read where τ is the MC time, ∆ t,x the field value at site (t, x) of the τ 'th configuration and mean # means averaging with respect to the respective subscript. In order to actually show the distribution from which the final estimates are calculated, we present (here and in the following) the histograms one obtains by stripping the means after the absolute values have been taken.
The histogram of the straightforward estimator shown in Fig. 14a is dominated by its broad variance (as is expected for a local estimator). More importantly, since the field θ is quasi-long-range, it requires many sweeps through the lattice to obtain a θ-independent distribution of ∆ like the ones depicted in Fig. 2. In fact, a typical configuration in the simulations is not distributed symmetrically around the origin but rather around some finite value ∆ 0 . The center of the configurations moves slowly (in Monte Carlo time) around the origin in field space. For this reason, taking the modulus right in the beginning leads to a significant bias towards larger values in the estimator (A1).
The broad variance mentioned above is a known statistical phenomenon in MC simulations and is usually cured by averaging over the spacetime lattice before taking the absolute value, schematically This sharpens the distribution but is less well motivated from a field theory perspective. The choice (A2) can be justified if there is spontaneous symmetry breaking and a small trigger is sufficient to align the values of the field on the lattice sites. In this case the absolute value does not change the result if we take the limits in the correct order, i.e. the spatial volume to infinity before removing the trigger. In the symmetric phase, on the other hand, already the spatial average should vanish in the thermodynamic limit and again taking the absolute value does not make a difference. Example distributions of this estimator are shown in Fig. 14c. Note the different scales on the x axes. It may come as a surprise that there is a second peak visible that distorts the mean of this distribution. This is due to the fact that at any non-zero temperature there are contributions from inhomogeneous configurations, which average out over the lattice to a very good approximation, see also Fig. 11. While for these data the distortion might be mild, we are not willing to take the risk of severely underestimating the observable for scale setting.
In the present work, what is even more problematic is that long-range (quasi-periodic) inhomogeneities must not be averaged over the spatial direction before taking absolute values. But, since we have to improve statistics as much as possible we will compromise by using where we, similarly as in the spatial correlation functions (37), first average over time.
As Fig. 14b indicates, this yields acceptable statistics while only using the assumption of temporal homogeneity which is a feature of all large-N f results we know of and was checked to be valid in our MC data, see, for example, Fig. 4. One should note that this procedure does not work in the high-temperature regime as the distribution in this case approaches that of Eq. (A1).
In future works other scale settings could be used and the corresponding results should be compared with those obtained in the present work. For example, the mass of the field ρ(t, x) may serve as an energy scale. The drawback of choosing a scale different from the minimum of the effective potential U eff (ρ) (at zero temperature and density) is that it is less straightforward to relate to the analytic results for large N f . In the large-N f limit the field ρ becomes infinitely heavy.

Appendix B: Autocorrelation analysis
During our simulations, we had to tackle severe autocorrelations similar to those described in [25]. In this appendix we summarize our extensive analysis of autocorrelation functions (ACFs) of various lattice estimators and provide details on how we arrive at the conclusion that our qualitative statements are robust despite large autocorrelations for certain parameter regions.

Identifying autocorrelation scales in an example
To facilitate such a discussion it is useful to visualize the topography of the effective action of the theory. In the infinite volume and -less important for this argument -continuum case, [31] found the general form of the saddle points of the effective action. The spatial profiles of the order parameter ∆ are given as a continuous family with four parameters related to overall scale, amplitude and phase variations, and a phase offset which is tightly related to translations. The finite volume we work in subjects these self-consistent solutions to the imposed boundary conditions such that for some of these parameters only a discrete subset of allowed solutions yields saddle points in finite volume. This entails a ragged landscape of valleys with local minima of the effective action that are separated by ridges stemming from the finitevolume effects and melting in the infinite-volume limit. From the analytical results, we expect chiral-spiral-like local minima (including the degenerate one, i.e. the constant order parameter) to be most important and our simulations confirm this expectation.
The above discussion immediately suggests that there are three qualitatively different kinds of autocorrelations: Sampled configurations will typically tend to fluctuate around one local minimum correlated within this valley on a (MC-time) scale τ fluct . During this process the reference chiral spiral will rotate the overall phase offset on a time scale τ U(1) which, for non-degenerate chiral spi-rals, is equivalent to translating this spiral. Eventually, the algorithm will climb (or tunnel through) a ridge and arrive in another valley on a time scale τ kmax .
From these three time scales, τ U(1) is of minor importance to us because we carefully crafted all of our observables to respect the U(1) (and closely related translational) symmetry. From the notable exceptions Fig. 2 and Fig. 13, however, we learned that it is quite sizable but clearly under control as the almost-perfect circles of Fig. 2 illustrate.
The other autocorrelation scales can be clearly distinguished in Fig. 15. For one exemplary parameter set, the figure shows ACFs of C σσ (x) for some randomly chosen lattice points x as well as the average and the (local in MC-time separation) maximal autocorrelations obtained over all lattice points. The latter rather unconventional quantity can be considered as a worst-case scenario for autocorrelations in C σσ . All of these are well described by an ansatz where b, τ 1 and τ 2 are free parameters. While the detailed numbers obviously depend on the data set chosen for fitting, the orders of magnitudes are consistent (cf. Tab. I).  Fig. 15). b τ1 τ2 x/a = 5 0.65 0.3 111.5 x/a = 19 0.86 0.8 101.5 x/a = 27 0.94 0. In order to associate the two numerical values τ 1 and τ 2 with the mechanisms described above, we consider representative time evolutions on each time scale in Fig. 16. Fig. 16a shows the Fourier spectrum of C σσ over 20 × τ 1 MC configurations where we conservatively use τ 1 ≈ 1. It is clearly seen that there is a constant peak at ν max = 3, while the MC evolution produces small fluctuations around this reference configuration. We conclude that the small time scale for this parameter set is generated from fluctuations around one local minimum, i.e. τ fluct = τ 1 .
To probe the larger time scale, we show a MC-time window of 20 × τ 2 configurations in Fig. 16b, where we conservatively estimated τ 2 ≈ 126. For visual clarity, we averaged blocks of τ 2 MC configurations which should be thought of as a coarse graining integrating out the high frequencies similar to an RG transformation. On this scale, the correlator spectra are smooth (due to the coarse graining) and sharp peaks and the MC evolution produces jumps in the dominant frequency k max . This finding relates the long time scale to τ kmax that was suggestively named after its effect of jumping between valleys changing k max . One should also stress that it is a non-trivial statement that 126 configurations tend to be rather coherent -again strengthening the association of τ 2 with τ kmax .

Analysis and reasoning about the rest of parameter space
The previous example indicates two important facts: Firstly, our choice of algorithmic parameters rendered τ fluct negligible while τ kmax is of considerable size. Secondly, besides τ kmax being large it is still under control in the sense that we have a statistically significant number of independent configurations N MC /τ kmax ≈ 380 even in the worst case discussed above. As the chosen parameter set is well inside the region of intermediate-scale inhomogeneous order we conclude from a statistically robust ensemble that our claims of intermediate-scale inhomogeneities without spontaneous symmetry breaking of any kind are robust with respect to autocorrelation effects. We checked for similar examples on all lattice sizes and lattice spacings. However, the above example was taken from a moderate temperature region. As we confirmed in this study, at the T = 0 line the system is critical which implies diverging correlation lengths -also in the MC-time evolution as is well known around practitioners [50]. We therefore expected and a posteriori verified huge autocorrelations for temperatures close to zero. One should stress that this is a physical effect; it can likely be circumvented by an appropriately adapted algorithm but still bares physically relevant information. Still, for a small region of very low temperatures τ kmax can easily exceed our greatest efforts of up to 8 · 10 4 configurations generated for some parameter sets. We therefore suggest to view the very-low-temperature results with a grain of salt quantitatively: They surely give a good impression of what phenomena to expect in the exceedingly large regions of space that are correlated for these temperatures but they might be quantitatively off due to autocorrelations suppressing the interference from subdominant local minima.
We remark that τ kmax has a clear tendency to decrease in the infinite-volume limit. This is the opposite behavior of what is typically found in symmetry-breaking systems and considered further evidence in support of the existence of a BKT phase and against spontaneous symmetry breaking, as was also mentioned in the main text.
The effect of larger flavor numbers is to reduce quantum fluctuations or, in the pictorial language from above, to deepen the valleys and grow the ridges. This effect is responsible for ultimately obtaining actual spontaneous symmetry breaking in the limit of infinite number of flavors. It also greatly enhances autocorrelations, particularly in suppressing jumps between different valleys. For N f = 8 the temperatures necessary to clearly observe inhomogeneities on average and jumps in the dominant wave number during the MC evolution occur concurrently are higher than in the 2-flavor case. This strengthens our finding in this and previous papers [24,25] that the convergence to mean-field behavior is quite rapid with the number of flavors. While technically this casts some doubt on the quantitative accuracy of the N f = 8 data, we want to point out again that this is driven by physical properties of the system more than technical difficulties.
Finally, we want to share some thoughts on how to improve on the situation: Due to extensive analytical results about the effective action, we were able to obtain a clear picture of the cause of autocorrelations in this model. One can easily imagine algorithmic improvements leveraging this knowledge. As the local minima can be enumerated by their dominant wave number k max , local updates in Fourier space might suffice, e.g. local Metropolislike updates or swapping of Fourier components. As the current updating scheme is very efficient to reduce some part of the autocorrelations, it would probably be advan- tageous to combine both kinds of updates into a single update step. Similar ideas are already used, e.g. [51]. Another approach could be to constrain the simulations to a single sector and a posteriori combining the results into a weighted sum. However, these approaches would be very specific and hardly generalizable to related problems, e.g. topological freezing in lattice QCD [52]. A modern approach that is agnostic of analytical knowledge, which is usually not as easily employed in more realistic theories, could be independence samplers from generative models, i.e. independently drawing the configurations from an efficient approximation of the probability distribution. Promising results in this direction were presented in [53], where they overcame topological freezing in 1+1D U(1) gauge theory.

Appendix D: Parameters
N f Ns = L/a 1/g 2 aρ0 Nt In order to calculate the various phase diagrams we generated many ensembles characterized by the control parameters (N f , T, L, µ) or (N f , N t , N s , ρ 0 µ), plus the four-Fermi coupling g 2 tuned to the required lattice spacing measured in units of ρ 0 = ρ T =µ=0 . We summarize the lattice spacings corresponding to the different values of N f , N s and g 2 in Tab. II.
As explained in the main text, we used different initial conditions for the fields to deal with thermalization problems: We performed scans with Gaussian-distributed seeds with mean zero, a freeze-out from high temperatures to reduce thermalization times and a heat-up procedure from the lowest temperature to exclude any hysteresis effects from the freeze-out. We also used a homogeneous cold start, in the sense of setting the initial configuration to ∆(x ) = 1 + i for all x , at small µ, where inhomogeneous configurations are suppressed. In Tab. III we collect the control parameters N t and µ for which we generated ensembles in equilibrium for each of these methods. Notice that we use the same lattice spacings as in Tab. II, which were determined via the freeze-out procedure, irrespective of the initial conditions. TABLE III: Parameter sets used in the simulations. Note that the uncertainty of aρ 0 (from Tab. II) propagates to the values of µ/ρ 0 , although we did not make this explicit for the sake of readability.