Thirring model at finite density in 0+1 dimensions with stochastic quantization: Crosscheck with an exact solution
aa r X i v : . [ h e p - l a t ] M a y Thirring model at finite density in dimensionswith stochastic quantization: Crosscheck with an exact solution
Jan M. Pawlowski
1, 2 and Christian Zielinski ∗ Institut für Theoretische Physik, Universität Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany ExtreMe Matter Institute EMMI, GSI, Planckstraße 1, D-64291 Darmstadt, Germany (Dated: July 8, 2018)We consider a generalized Thirring model in dimensions at finite density. In order todeal with the resulting sign problem we employ stochastic quantization, i.e., a complex Langevinevolution. We investigate the convergence properties of this approach and check in which parameterregions complex Langevin evolutions are applicable in this setting. To this end we derive numerousanalytical results and compare directly with numerical results. In addition we employ indirectindicators to check for correctness. Finally we interpret and discuss our findings.
PACS numbers: 05.50.+q, 71.10.FdKeywords: Thirring model, finite density field theory, complex Langevin evolution, stochastic quantization
I. INTRODUCTION
Despite all efforts, one of the outstanding problems oflattice field theory until this day is the sign problem. Theintroduction of a finite chemical potential µ > rendersthe path integral measure complex and rapidly oscillatingin many theories of interest, like quantum chromodynam-ics (QCD) in dimensions. The oscillatory behaviorsignificantly increases the numerical costs, in particularin the continuum limit. A hard sign problem exists fortheories where the costs grow more than polynomiallywith the volume. This obstacle hinders numerical ab ini-tio studies of strongly interacting matter under extremeconditions and the understanding of the phase diagram ofQCD. There is no satisfactory solution known to the signproblem, despite the large number of proposed solutions.Among the proposed solutions we can find reweightingtechniques, Taylor expansions about µ = 0 , extrapola-tions from imaginary chemical potential, the introductionof dual variables and a canonical ensemble approach. Forrecent reviews see e.g. [1, 2]. However, due to the over-lap problem reweighting techniques are computationallyexpensive and can only be used for small µ , while thenumerical determination of Taylor coefficients is noisyand the expansion converges slowly [3]. Also the contin-uation from imaginary chemical potential is a nontrivialtask [4]. The application of dual variables and the canon-ical ensemble approach is still under active research, seee.g. [5, 6] for dual observables and [7, 8] for simulationswith canonical ensembles.In this paper we employ a different approach. Parisiproposed already in 1983 that stochastic quantization[9]—for a review see e.g. [10]—could circumvent the signproblem in terms of a complex Langevin evolution [11].However, it is well known that the Langevin evolutionmay converge towards unphysical fixed points. It has ∗ Permanent address: Division of Mathematical Sciences, NanyangTechnological University, Singapore 637371 been successfully applied to the SU(3) spin model [12, 13],to an effective theory of QCD in the strong-coupling limit[14], simple models of quantum chromodynamics [15, 16]and to the relativistic Bose gas [17, 18] at finite den-sity. Furthermore it has been applied to quantum fieldsin Minkowski time [19, 20], also in nonequilibrium [21].Counterexamples are given by the three-dimensional XYmodel at finite chemical potential for small β [22] and incases of gauge theories with static charges [23]. Early in-vestigations of complex Langevin evolutions can be foundin [24–26], while for reviews see e.g. [27, 28]. Recently,a set of consistency conditions indicating correct conver-gence could be derived [29–31]. When truncating thisinfinite tower of identities one obtains necessary condi-tions for correctness.In this work we apply a complex Langevin evolutionto a generalized Thirring model at finite density. Here itserves us as a model theory to check for the applicabil-ity of this method. Our results extend the studies car-ried out in [32], which led to ambiguous results. In thispaper we restrict ourselves to the case of dimen-sions and deal with the question of whether a complexLangevin evolution can enable finite density calculationsin this setting. Further investigations of this approachin the -dimensional generalized Thirring model arepresented in [33]. The -dimensional model appearsfor example in effective theories of high temperature su-perconductors and graphene, see e.g. [34] and referencesgiven therein. It is also worth mentioning that in thecase of the three-dimensional massless Thirring model, afermion bag approach was successfully applied in [35].We organized the paper as follows: In Sec. II we intro-duce a generalized Thirring model and its formulation onthe lattice. We discuss the Langevin equation and its nu-merical implementation. In Sec. III we present a closedexpression for the partition function of the lattice theoryand derive some observables of interest. We also discussadditional indicators to evaluate the convergence prop-erties of the complex Langevin evolution. In Sec. IV wediscuss the results of the numerical part of this work andaim to answer the question in which parameter regimeresults are reliable. We end this paper with concludingremarks in Sec. V. II. THE GENERALIZED THIRRING MODELA. Continuum formulation
We consider a generalization of the Thirring model.The historical model was introduced in 1958 by WalterE. Thirring and is one of the rare examples of an exactlysolvable quantum field theory [36]. While the originalmodel describes self-interaction fermions in dimen-sions, we consider N f fermion flavors at finite density.We begin with a generalization to d dimensions andthen later specialize to the case of dimensions. TheEuclidean Lagrangian in the continuum reads L Ψ = N f X i =1 Ψ i (cid:0) /∂ + m i + µ i γ (cid:1) Ψ i + g N f N f X i =1 Ψ i γ ν Ψ i . (1)The index i = 1 , . . . , N f enumerates fermion flavors, m i and µ i denote the bare mass and bare fermion chemi-cal potential of the respective flavor and g is the barecoupling strength. The γ matrices satisfy the Cliffordalgebra { γ µ , γ ν } = 2 δ µν .The model shows breaking of chiral symmetry at µ = 0 in dimensions [37]. For the -dimensionalThirring model, the equivalence to the sine-Gordonmodel can be shown [38, 39].The four-point interaction can be resolved with theintroduction of an auxiliary field A ν . This formulationreads L = X i Ψ i (cid:0) /∂ + i /A + m i + µ i γ (cid:1) Ψ i + N f βA ν . (2)Here we introduced the inverse coupling β = 1 / (cid:0) g (cid:1) .When integrating A ν out, we recover (1). Although theauxiliary field A ν is not a gauge field, the model can beinterpreted as a more general gauge theory after gaugefixing, see e.g. [40]. After integrating out the fermionicdegrees of freedom we find Z = ˆ D A Y i det K i ! e − S A = ˆ D A e − S eff ,S A = N f β /T ˆ d t ˆ d d − x A ν . (3)Here we introduced the temperature T and K i = /∂ + i /A + m i + µ i γ . Including the fermion determinant in the exponential term yields S eff = S A − X i Tr log K i . (4)For the fermion determinant the relation det K i ( µ ) = [det K i ( − µ ⋆ )] ⋆ (5)holds, thus rendering the path integral measure complexfor µ > . At vanishing or purely imaginary chemicalpotential, the determinant is real and the theory is freeof a sign problem. If the fermion determinant is replacedby its modulus, we refer to this as the phase-quenchedcase. Physically this corresponds to the introduction ofan isospin chemical potential.Like quantum chromodynamics, the Thirring modelexhibits Silver Blaze behavior [41, 42]. It implies thatat vanishing temperature there is a threshold µ c , so thatobservables are independent of the chemical potential µ for µ < µ c . While in the full theory the onset is given bythe physical fermion mass m phys , in the phase-quenchedtheory we have µ c = m π / , where m π is the physicalpion mass. B. Lattice formulation
We consider the case of dimensions—corresponding to a quantum mechanical system—withlattice spacing a and N t lattice points. We employ stag-gered fermions [43–46] and denote the number of latticeflavors, i.e., the number of staggered fermion fields, by N .Furthermore we assume that N t is even, as otherwise theformulation of staggered fermions is conceptually prob-lematic and (5) is violated. In order to introduce a finitechemical potential µ , we use the prescription by Hasen-fratz and Karsch [47]. For notational ease we refer tothe one-component auxiliary field as A t = A ( x = t ) and all dimensionful quantities are scaled dimensionlessby appropriate powers of a . Furthermore we introducethe hopping parameter κ = 1 / (2 m ) . The temperaturecorresponds to the inverse temporal extension T = N − t .Using this formulation the lattice partition functionreads Z = ∞ ˆ −∞ N t Y t =1 d A t Y i det K i ! e − S A (6)with S A = N β P t A t and flavor index i = 1 , . . . , N .The fermion matrix takes the form K i ( t, τ ) = 12 (1 + i A t ) e µ i δ t +1 ,τ −
12 (1 − i A τ ) e − µ i δ t − ,τ + m i δ tτ , (7)where we impose antiperiodic boundary conditions,cf. [32, 48]. In our analysis we focus on a few observ-ables, namely the fermion density and condensate, theenergy density and the phase factor of the fermion de-terminant. In the following, sums over the flavor index i are not implied. The fermion density of a given flavor isgiven by h n i i = 1 N t (cid:18) ∂ log Z∂µ i (cid:19) T = 1 N t (cid:28) Tr (cid:18) ∂K i ∂µ i K − i (cid:19)(cid:29) . (8)The fermion condensate follows from h χ i χ i i = 1 N t (cid:18) ∂ log Z∂m i (cid:19) T,µ i = 1 N t (cid:10) Tr K − i (cid:11) (9)and the energy density reads h ε i i = − (cid:18) ∂ log Z∂N t (cid:19) µ i + µ i h n i i , (10)which we normalize to h ε i i ( µ = 0) = 0 .The phase factor of the determinant is defined by exp ( i φ ) = det K/ | det K | . It can be expressed in termsof the partition function Z N = ∞ ˆ −∞ Y t d A t (det K ) N e − S A , (11)for N degenerated flavors and Z pq N for the phase-quenched case, where the fermion determinant in (11)is replaced by its modulus. The expectation value of exp ( i N φ ) follows in the N flavor phase-quenched theory[49, 50] as (cid:10) e i N φ (cid:11) pq N = Z N Z pq N ∈ [0 , . (12)A value close to zero indicates a rapidly oscillating pathintegral measure with a severe sign problem. C. Complex Langevin evolution
The idea of stochastic quantization is that observablesin a Euclidean quantum field theory can be obtained asthe equilibrium values of a statistical system coupled toa heat bath [10]. The problem of quantizing a field the-ory is then reduced to finding the static solutions of anassociated Langevin equation. If the action is real andbounded from below, correctness of this approach can beensured. We can also formally generalize to the case of acomplex action [11]. This situation naturally arises whenconsidering field theories at finite density. Until this daythere is a lack of rigor mathematical understanding re-garding the validity of this procedure. However, in caseswhere it is converging correctly one has a very elegantsolution for the sign problem at hand.We aim to check for the applicability of complexLangevin evolutions to the Thirring model. To this end we have to find the static solution of the Langevin equa-tion ∂∂ Θ A t (Θ) = − δS eff [ A ] δA t (Θ) + √ η t (Θ) , (13)where Θ denotes a fictitious time. The noise term η t (Θ) follows a Gaussian distribution with h η t (Θ) i = 0 , h η t (Θ) η t ′ (Θ ′ ) i = δ ( t − t ′ ) δ (Θ − Θ ′ ) . (14)A simple approach to solve the Langevin equation nu-merically is a first order integration scheme with fixedstepsize ǫ L . Higher order integration schemes of O ( ǫ / L ) have been employed in the literature too [13, 51]. How-ever, in some models fixed stepsize integration schemesfail due to the occurrence of run-away trajectories, whichcan be avoided by the use of an adaptive stepsize [52, 53].Although a constant stepsize proved here to be sufficient[32], we employ an adaptive stepsize algorithm due tobetter convergence properties. For N degenerated fla-vors our discretization of (13) reads A t (Θ + ǫ L ) = A t (Θ) + ǫ L D t (Θ) + √ ǫ L η t (Θ) (15)with drift term D t (Θ) = −N βA t (Θ)+ N i (cid:2) K − ( t + 1 , t ) e µ + K − ( t, t + 1) e − µ (cid:3) . (16)After each integration step the stepsize ǫ L will be up-dated according to ǫ L ≡ ǫ L (Θ) = δ max t | D t (Θ) | (17)with stepsize parameter δ = 10 − (compare to [32]).It is possible to generalize the real noise term in (15)to an imaginary one [32] via the replacement η t (Θ) → √I + 1 Re η t (Θ) + i √I Im η t (Θ) (18)with I ≥ . The noise correlators then read h Re η t (Θ) Re η t ′ (Θ ′ ) i = h Im η t (Θ) Im η t ′ (Θ ′ ) i = δ ( t − t ′ ) δ (Θ − Θ ′ ) (19)and h Re η t (Θ) Im η t ′ (Θ ′ ) i = 0 . Assuming correctness ofthe complex Langevin evolution and numerical stability,we expect expectation values to be independent of I . III. ANALYTICAL RESULTSA. Exact partition function
We begin with the partition function for one staggeredfermion field, i.e., N = 1 . We incorporate antiperiodicboundary conditions and for brevity we introduce B ± = 12 (2 κ ) N t (cid:16) ± p B c /β (cid:17) N t , (a) Plot of the fermion density h n i .(b) Plot of the fermion condensate h χχ i . Figure 1: In the phase structure in d = 0 + 1 we find acondensed phase for large µ . B c = β + 4 ( β + 1) κ . (20)Then the partition function (6) reads Z = 2 (cid:18) π β (cid:19) N t / [ B + + B − + cosh ( N t µ )] . (21)This can be shown for example by systematic saturationof the Grassmann integral or the help of the determinant identities in [54]. For the fermion density we find h n i = sinh ( µ/T ) B + + B − + cosh ( µ/T ) , (22)while the fermion condensate is given by h χχ i = 2 κ p β/ B c ( B + − B − ) B + + B − + cosh ( µ/T ) . (23)The expression for the energy density is rather lengthyand we will not quote it here explicitly. Figure 1 showsthe dependence of these observables on both β and µ . Forlarge µ we find a condensed phase, which is well separatedfor large N t .As it turns out, we can take the continuum limit ofthe density h n i analytically. To this end we recover thephysical units of all dimensionful quantities by reintro-ducing the lattice spacing a . We fix the dimensionfultemperature T − = aN t , express the lattice spacing a asa function of the number of lattice points N t and takethe limit N t → ∞ . We obtain h n i cont = sinh (cid:0) µT (cid:1) exp (cid:16) κ − β T βκ (cid:17) (cid:2) (cid:0) T κ (cid:1)(cid:3) + cosh (cid:0) µT (cid:1) , (24)where all units are explicitly dimensionful. In the zerotemperature limit T → we find h n i cont = Θ ( µ − m phys ) with physical fermion mass m phys = m + g and bare mass m . B. Several flavors
We also considered the case of more than one latticeflavor, but only quote here the simplest case of (11) for N t = 2 lattice points and N = 2 staggered fermion fields Z = π β κ (cid:2) β + 4 βκ + 8 β κ + 5 κ + 12 βκ + 12 β κ + 8 β κ cosh (2 µ )+ κ cosh (4 µ ) + 8 βκ cosh (2 µ ) − βκ cosh (4 µ )+16 β κ cosh (2 µ ) + 4 β κ cosh (4 µ ) (cid:3) . (25)We observe that for two flavors the density, conden-sate and energy density have plateaus in the range of β ≈ . – . , see also Fig. 9. They appear between theonset to the condensed phase and saturation of the den-sity. The height of the plateau depends on the massesof the flavors, and in the case of degenerated flavors itis exactly at half of the maximum value of h n i or h χχ i ,respectively. If β is very small or large, the plateauseventually disappear. For large β this can be understoodfrom the weak coupling limit, see Sec. III C. The exis-tence of these plateaus on the lattice is further confirmedby a heavy dense limit [33] and the Monte-Carlo studiesin Sec. IV D. In the general case of N flavors, we can findup to N − intermediate plateaus in these observables.We give a natural explanation of these structures in theAppendix. C. Weak coupling limit
In the limit of β ≫ the path integral measure has astrong peak at the origin and can be approximated by aDirac δ function. For N flavors we find Z weak = (cid:18) π N β (cid:19) Nt (cid:18) N t (cid:19) N × Y i (cid:2) B + i + B − i + cosh ( N t µ ) (cid:3) (26)for the partition function in the weak coupling limit.Here we introduced B ± i = 12 (2 κ i ) N t (cid:18) ± q κ i (cid:19) N t . (27)This can be shown again using the identities in [54]. Notethat in this limit the contribution from different flavorsfactorize and the phase-quenched case equals the full the-ory. For N degenerated flavors, i.e., B ± i = B ± , the totalfermion density is given by h n i = N sinh ( µ/T ) B + + B − + cosh ( µ/T ) (28)and the fermion condensate by h χχ i = 2 κ N ( B + − B − ) √ κ ( B + + B − + cosh ( µ/T )) . (29)While approaching the noninteracting limit, the plateaustructures observed for N > disappear. D. Analyticity in µ An observable O which is even in µ can be interpretedas a function of µ . Assuming analyticity in µ , we cananalytically continue O to purely imaginary chemical po-tential, i.e., µ ≤ . In this case the fermion determinantis real and free of a sign problem due to (5) and we canemploy a real Langevin evolution. For µ > we employa complex Langevin evolution.Assuming the correctness of the complex Langevin evo-lution, O should be analytic at µ = 0 . Any nonanalytic-ity would be an indicator for incorrect convergence. Thiscriterion was previously employed in [13] and in this workwe apply it to the condensate h χχ i . E. Consistency conditions
In [31] the authors derived a set of identities indicatingcorrect convergence of expectation values obtained by acomplex Langevin evolution. These consistency condi-tions state that for all entire holomorphic observables O the expectation value h L Oi = 0 has to vanish. Here L t = (cid:18) dd A t + D t (cid:19) dd A t (30)
0 1 2 3 4 5 6 P r opaga t o r | G ( ∆ ) | Distance ∆ Full theory - β = 1, m = 1, N t = 4 µ = 0.0 µ = 0.1 µ = 0.2 Figure 2: The propagator at different values of µ .denotes the Langevin operator and D t = − d S eff / d A t thedrift term.As this defines an uncountable number of identities,we restrict the analysis to a finite subset. We follow [31]and consider here observables O ( t, k ) = exp ( i kA t ) . Theresulting conditions take the form h L t O ( t, k ) i = (cid:10) i k [ i k + D t ] e i kA t (cid:11) = 0 , (31)which have to hold for ∀ t and ∀ k . Without loss of gen-erality we set t = 1 . Besides O ( t, k ) , we also check theconsistency conditions for the fermion density and thecondensate. F. Propagator at finite density
We define the propagator at finite temperature T = N − t and chemical potential µ by h χ ( t ) χ ( t ) i = 1 Z (cid:10) K − ( t , t ) (cid:11) (32)with partition function Z given by (21). It is helpful tointroduce the notation G (∆) = D χ (1) χ (cid:16) (cid:17)E (33)with ˆ∆ = ∆ mod N t . For small lattices, the inversion ofthe fermion matrix and the calculation of the expectationvalue can be done analytically. A typical example canbe found in Fig. 2. For small distances the propagator G (∆) falls off exponentially, but due to the finite latticeextension eventually rises again. The introduction of afinite chemical potential µ shifts the minimum and resultsin G (∆) = G ( − ∆) . -0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 F e r m i on den s i t y < n > Chemical potential µ Full theory - β = 3, m = 1, N t = 4TheoryCLEMC Figure 3: Typical errors estimated by the standarddeviation.
IV. NUMERICAL RESULTSA. Implementation
In the previous section we derived numerous analyt-ical results, which allow us to benchmark the complexLangevin evolution and check for its correctness. Forthe numerical simulations we implemented an adaptivenumerical integration scheme of the associated Langevinequation as described in Sec. II C. Using a computer alge-bra system, we can determine the inverse of the fermionmatrix in (7) and find an analytical expression of thedrift term in (16) for small lattices in order to minimizenumerical errors.The evaluation of the observables for a given set ofparameters begins with a hot start, i.e., the random ini-tialization of the auxiliary field, followed by typically steps for thermalization. The field configuration is thensampled about O (cid:0) (cid:1) times and used to evaluate theobservables. To reduce potential autocorrelation effects,two subsequent samples are separated by ten dummysteps. In all following plots, numerical values obtained bya complex Langevin evolution will be denoted by “CLE,”while analytical results are denoted by “Theory.”A simple method to estimate the error is to take thestandard deviation of the different samples of a givenobservable in a particular run. However, the resultingerrors are much larger than the empirical observed sta-tistical fluctuations between different runs and are over-estimated, see Fig. 3 for a typical example. In order toobtain more reliable error bounds, we employ a bootstrapanalysis [55]. Besides the statistical error, we also have todeal with systematic errors induced by a finite stepsize.To this end we checked that the stepsize parameter δ in(17) is sufficiently small.As an additional test we used an adaptive quasi-MonteCarlo strategy [56] for the direct evaluation of expecta-tion values, which works well for sufficiently small N t and µ . On larger lattices the sign problem is severeand the algorithm fails. As the path integral measurefalls off rapidly for large field configurations, we replacethe numerically difficult noncompact integration over theauxiliary field by a compact domain [ − Λ , Λ] N t . Weparametrize this ad hoc cutoff Λ by Λ = r σ N β . (34)If the geometric mean of the field configuration is Λ , thenthe path integral measure is dampened by a factor of exp ( − σ ) . If we choose σ too small, we introduce a largecutoff effect. If σ is very large, we still have the problemthat the integrand is close to zero in most of the inte-gration region. As a compromise here we use σ = 10 .In figures, we refer to numerical results obtained by thismethod as “MC.” B. Comparison of results
We can see the results of a typical run for two differentsets of parameters in Figs. 4 and 5. Because of the highsample size, the error bars tend to be very small. We seethat the results obtained with an adaptive Monte-Carlomethod are in good agreement with analytical results.The numerical results obtained with a complexLangevin evolution show some deviations for intermedi-ate µ . We systematically observe that the gap widens fordecreasing values of β , i.e., stronger couplings. In addi-tion the numerical evaluation becomes more noisy as therelative magnitude of the noise term η t increases. Forlarge β the agreement becomes better and the numericsare very stable.From the phase factor of the fermion determinant, wesee how the sign problem gets more severe for increasinglattice sizes. However, our approach is not affected bythis. It is interesting to note that the sign problem seemsto be less pronounced when using complex Langevin evo-lutions.In Fig. 5 we can already see how in the limit T → the Silver Blaze behavior becomes apparent. In the limit N t → ∞ the observables will eventually become indepen-dent of the chemical potential µ up to some threshold µ c .It seems that the position of this onset differs slightlyfrom analytical predictions. However, for large β theyare in good agreement.That the observed deviations are not just an effect ofa finite stepsize can be seen in Fig. 6. For a typical setof parameters we calculate the fermion density and con-densate for varying stepsize factors δ . The extrapolationshows that also in the limit δ → there is a discrepancybetween theory and complex Langevin evolutions. -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 F e r m i on den s i t y < n > Chemical potential µ Full theory - β = 3, m = 1, N t = 4TheoryCLEMC (a) Fermion density h n i . F e r m i on c onden s a t e < χ - χ > Chemical potential µ Full theory - β = 3, m = 1, N t = 4TheoryCLEMC (b) Fermion condensate h χχ i . P ha s e f a c t o r R e < e i φ > Chemical potential µ Full theory - β = 3, m = 1, N t = 4CLEMC (c) Phase factor of determinant (cid:10) e i φ (cid:11) . Figure 4: Benchmarking results for β = 3 . -0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 F e r m i on den s i t y < n > Chemical potential µ Full theory - β = 1, m = 1, N t = 8TheoryCLEMC (a) Fermion density h n i . -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.5 1 1.5 2 F e r m i on c onden s a t e < χ - χ > Chemical potential µ Full theory - β = 1, m = 1, N t = 8TheoryCLEMC (b) Fermion condensate h χχ i . -0.2 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 P ha s e f a c t o r R e < e i φ > Chemical potential µ Full theory - β = 1, m = 1, N t = 8CLEMC (c) Phase factor of determinant (cid:10) e i φ (cid:11) . Figure 5: Benchmarking results for β = 1 . F e r m i on den s i t y < n > Stepsize δ Full theory - β = 1, m = 1, µ = 1, N t = 4 δ = 0, TheoryCLELinear Fit (a) Fermion density h n i . F e r m i on c onden s a t e < χ - χ > Stepsize δ Full theory - β = 1, m = 1, µ = 1, N t = 4 δ = 0, TheoryCLELinear Fit (b) Fermion condensate h χχ i . Figure 6: Extrapolation to vanishing stepsize δ → . F e r m i on den s i t y < n > Chemical potential µ Full theory - β = 100, m = 1, N t = 4Weak Coupling LimitCLEMC Figure 7: Density in the weak coupling regime. F e r m i on den s i t y < n > Inverse Coupling β Full theory - m = 1, µ = 1, N t = 4TheoryCLEMC (a) Fermion density h n i . -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.5 1 1.5 2 2.5 3 F e r m i on c onden s a t e < χ - χ > Inverse Coupling β Full theory - m = 1, µ = 1, N t = 4TheoryCLEMC (b) Fermion condensate h χχ i . Figure 8: Observables as function of β . C. Coupling parameter dependence
In the weak coupling limit of Sec. III C, we observevery good agreement of all numerical and analytical re-sults. For increasing β the numerical results are asymp-totically approaching analytical results. In Figure 7 wesee that for β & O (100) the method is then able to de-liver reliable results. In the XY model at finite density[57] a similar behavior was observed. For small β theauthors observed incorrect convergence, while for large β they found agreement. However, in the XY model thereis a clear separation between both regimes.As already pointed out, the applicability of complexLangevin evolutions depend on the magnitude of β . InFigure 8 we find the density and condensate as a functionof β . We observe that for weak couplings, i.e., in thelimit β → ∞ , it slowly converges towards the analyticalresults. For β = O (100) and above we then find goodagreement. Contrariwise for strong couplings we observea large discrepancy. Moreover for β . the evaluationof observables is very noisy, resulting in large error bars. -0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 1 2 3 4 5 F e r m i on den s i t y < n > Chemical potential µ Full theory - N f = 2, β = 0.6, m = 3, N t = 4TheoryCLEMC (a) Fermion density h n i . -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 1 2 3 4 5 F e r m i on c onden s a t e < χ - χ > Chemical potential µ Full theory - N f = 2, β = 0.6, m = 3, N t = 4TheoryCLEMC (b) Fermion condensate h χχ i . Figure 9: Observables for N = 2 degenerated flavors. D. Several flavors
In Sec. III B we found that in a certain range of β ,the density, the condensate and the energy density showup to N − intermediate plateaus when considering N flavors.We give an example for N = 2 flavors at a couplingof β = 0 . in Fig. 9. For small lattice sizes and small µ Monte Carlo studies confirm these predictions. Becauseof the sign problem we are restricted to µ . with theseparticular parameters, before the numerical evaluationbreaks down.When trying to reproduce the plateaus with a com-plex Langevin evolution, we notice that for every valueof β the result qualitatively resembles the N = 1 case.Directly after the onset the fermion density rises untilit reaches saturation. Only in the limit of large β theplateaus disappear in the analytical solution and there isagreement with numerical results. F e r m i on c onden s a t e < χ - χ > Chemical potential µ Full theory - m = 1, N t = 4 β = 0.5, Theory β = 1, Theory β = 3, Theory β = 10, Theory β = 0.5, Langevin β = 1, Langevin β = 3, Langevin β = 10, Langevin Figure 10: Logarithmic plot of h χχ i over µ . E. Analyticity at µ = 0 As described in Sec. III D, we check for possible non-analytic behavior of the condensate h χχ i at µ = 0 . InFig. 10 we see that for µ ≤ the real Langevin evolu-tion is, as expected, in agreement with theory. Becauseof the absence of a sign problem numerical results arereliable in this regime. Also for small µ & there is nostatistically significant deviation. Hence the condensate h χχ i is analytic within our numerical accuracy.If we increase µ further, we observe that at some pointa disagreement becomes apparent. The resulting gap ismore pronounced for small β . It does not appear sud-denly, but develops smoothly and is visible for all non-large β . F. Consistency conditions
We checked the consistency conditions of Sec. III E forthe density h n i , the condensate h χχ i , and O ( t = 1 , k ) for different sets of parameters at δ = 10 − . Before webegin with the interpretation of the consistency condi-tions, we point out again the problem of using the stan-dard deviation to estimate error bounds. The result-ing error is much larger than the actual observed sta-tistical fluctuations and typically we obtain results likee.g. Re h L O (1 , i = 0 . ± . for N t = 4 , I = 0 , N = β = m = µ = 1 . Here N t denotes the temporalextension of the lattice, N the number of degeneratedflavors, I the imaginary noise, β the inverse couplingconstant, m the mass and µ the chemical potential. Theoverestimation of the error makes a meaningful interpre-tation of the conditions impossible. Hence we estimatethe error with a bootstrap analysis.We begin with the consistency conditions of the densityand the condensate. The evaluation is extremely noisyand only becomes slightly more stable for large valuesof β . Without loss of generality we restrict ourselves to0
100 1000 10000 -3 -2 -1 0 1 2 3 C oun t s Re LnFull theory - β = 1, m = 1, µ = 1, N t = 4Samples (a) Histogram of h Ln i .
100 1000 10000 100000 -3 -2 -1 0 1 2 3 C oun t s Re L χ - χ Full theory - β = 1, m = 1, µ = 1, N t = 4Samples (b) Histogram of h Lχχ i .
10 100 1000 10000 -3 -2 -1 0 1 2 3 C oun t s Re Le iA Full theory - β = 1, m = 1, µ = 1, N t = 4Samples (c) Histogram of h L O (1 , i . Figure 11: Histograms for different consistencyconditions.
100 1000 10000 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 C oun t s Im A-Full theory - β = 1, m = 1, µ = 1, N t = 4Samples Figure 12: Distribution of (cid:10) Im A (cid:11) .the case of t = 1 . As an example we quote Re h Ln i =(143 ± × and Re h Lχχ i = ( − ± × for N t = 4 , I = 0 , N = β = m = µ = k = 1 . Theresulting expressions seem to be numerically ill-definedand despite large sample sizes, we are unable to draw anyconclusions about the validity of the conditions. In ourcase both observables proved to be unsuitable to checkfor the correctness of the complex Langevin evolution.The evaluation of the conditions for the observable hO ( t, k ) i on the other hand is stable for all checked pa-rameter sets. If we take the error bounds seriously, wehave to conclude that some conditions are not compat-ible with a vanishing value and are violated. We quotehere Re h L O ( t, k ) i = − . ± . for N t = 4 , I = 0 , N = β = m = k = 1 , µ = 3 as a typical example of aviolated condition and Re h L O ( t, k ) i = 0 . ± . for N t = 4 , I = 0 , N = β = m = µ = k = 1 as one whichis compatible with a vanishing value. We also checkedthe conditions for k = 2 , , which turned out to be morenoisy compared to the k = 1 case. We interpret the vio-lated conditions as an indicator for incorrect convergenceof the complex Langevin evolution.When considering the distributions of the evaluatedconsistency conditions, we can gain further insights. InFig. 11 we find histograms of Ln , Lχχ , and L O (1 , .In general the distributions are asymmetrical and non-Gaussian. In the case of the density and condensate wesaw that the distribution falls off extremely slowly and weobserved values with an absolute value of up to O (cid:0) (cid:1) ,resulting in a very noisy evaluation. In contrast the dis-tribution of L O (1 , falls off rapidly.It is also interesting to check the distribution of theimaginary part of the auxiliary field itself. In thederivation of the consistency condition, one assumes thatboundary terms of the field would vanish. It is then im-portant to check that the imaginary part falls off rapidlyenough. In Fig. 12 we find a typical histogram of theaverage imaginary part Im A = N − t P t Im A t with samples. We see that the resulting distribution appears1Figure 13: Eigenvalues of the fermion matrix.to be compatible with our hypothesis. G. Eigenvalues of the fermion matrix
In Figure 13 we find a scatter plot of the eigenvalues ofthe fermion matrix K defined in (7). We sampled approx-imately × eigenvalues during a complex Langevinevolution. Eigenvalues come in point-reflected pairs atpoint ( m, , where m denotes the mass of the fermion.For a vanishing chemical potential µ = 0 all eigenvalueslie on the line Re λ = m . H. Imaginary noise
Assuming correct convergence of the complex Langevinevolution, observables should turn out independent of theimaginary noise term controlled by I . However, in Fig. 14we actually observe such a dependence. This was ob-served previously in other systems too [30]. Attemptsto fine-tune I so that the complex Langevin evolutionis deformed and correctly reproduces analytical resultsdid not succeed. In almost all cases an imaginary noise I > caused a more severe disagreement between nu-merical and analytical results. I. Periodic continuation of the imaginary part
A necessary condition for the correctness of the com-plex Langevin evolution is that the distributions of the F e r m i on c onden s a t e < χ - χ > Chemical potential µ Full theory - β = 3, m = 1, N t = 4TheoryI = 0.0, CLEI = 0.1, CLEI = 0.3, CLEI = 1.0, CLE Figure 14: Condensate h χχ i with imaginary noise I . F e r m i on c onden s a t e < χ - χ > Chemical potential µ Full theory - β = 3, m = 1, N t = 4Theory Λ = ∞ , CLE Λ = 1.0, CLE Λ = 0.2, CLE Λ = 0.1, CLE Figure 15: Condensate h χχ i with periodic cutoff Λ .imaginary parts of the fields have to fall off rapidlyenough. Although in the generalized Thirring model thisseems to be the case, we can restrict the imaginary partof the field to the interval [ − Λ , +Λ] and make a periodiccontinuation. A similar approach was employed in [30](see also [29]), where a fine-tuning of the cutoff enabledsimple toy models to reproduce correct results. However,in Fig. 15 we see that every finite value of Λ widens thegap to the correct results. Hence in this simple form thisapproach cannot be applied to the generalized Thirringmodel. V. CONCLUSIONS
In this paper we applied a complex Langevin evolutionto a generalized Thirring model in dimensions inorder to deal with the resulting sign problem for µ > . For intermediate values of the chemical potential wefound a gap between analytical and numerical results,which size depends on the inverse coupling β . While for2small β the discrepancy is large, we observe agreementfor large β . In particular, for β & O (100) we are usuallyable to reproduce analytical results with high accuracy.Furthermore for small µ & we did not observe anysignificant deviations to theoretical predictions and thefermion condensate is analytic at µ = 0 .However, in the case of more than one flavor we ob-served a qualitative disagreement. Our approach seemsto be unable to reproduce the plateaus we found forcertain ranges of β . Another interesting observation isthe violation of several consistency conditions, indicatingthat the complex Langevin evolution might not convergecorrectly in general. Attempts to force correct conver-gence by an ad hoc fine-tuning of a periodic cutoff Λ oran imaginary noise term I did not succeed.In a subsequent paper, we will present our findings forthe generalized Thirring model in dimensions [33].Further investigations have to deal with the question ofhow to address the aforementioned problems. In par-ticular, coordinate transformations as suggested in [58]and gauge cooling procedures like the one employed in[59] might allow a stabilization of the complex Langevinevolution. ACKNOWLEDGMENTS
We thank I.-O. Stamatescu for uncountable discussionsand useful remarks on the manuscript. We also thankC. Gattringer for a proof of (21) by systematic satura-tion of the Grassmann integral and V. Kasper for a proofusing the determinant identities in [54]. Furthermore weacknowledge E. Seiler’s derivation of the plateau struc-tures as presented in the Appendix. Finally, we thankG. Aarts and D. Sexty for discussions. This work is sup-ported by the Helmholtz Alliance HA216/EMMI and byERC-AdG-290623. C. Z. thanks the German NationalAcademic Foundation for financial support.
APPENDIX: PLATEAUS
As previously discussed, in the case of N > flavors weobserve intermediate plateaus in the considered observ-ables. In the nonlinear O(2) sigma model, the authors of[60] found similar structures. They explained this finite size behavior with the crossing of energy levels at finitechemical potential. However, here we reproduce themwith the help of a continuum model, which correspondsto the dimensional generalized Thirring model. Tothis end we consider the Hamiltonian H = N f X i =1 m i ( n i + n i ) + g N f X i =1 Q i (35)with flavor index i = 1 , . . . , N f and bare masses m i . Weintroduced Q i ≡ n i − n i , the particle number operator F e r m i on den s i t y < n > Chemical potential µ Toy model - β = 10, g = 1, m = 1.0, m = 0.9
A22 , 457 (2007),arXiv:hep-lat/0509180 [hep-lat].[3] F. Karsch, B.-J. Schaefer, M. Wagner, and J. Wambach,Phys.Lett.
B698 , 256 (2011), arXiv:1009.5211 [hep-ph]. [4] M. P. Lombardo, PoS
CPOD2006 , 003 (2006),arXiv:hep-lat/0612017 [hep-lat].[5] A. Schmidt, Y. D. Mercado, and C. Gattringer, PoS
LATTICE2012 , 098 (2012), arXiv:1211.1573 [hep-lat].[6] Y. D. Mercado, C. Gattringer, and A. Schmidt, Com-put.Phys.Commun. , 1535 (2013), arXiv:1211.3436 [hep-lat].[7] A. Alexandru, M. Faber, I. Horvath, and K.-F. Liu,Phys.Rev. D72 , 114513 (2005), arXiv:hep-lat/0507020[hep-lat].[8] A. Alexandru and U. Wenger, Phys.Rev.
D83 , 034502(2011), arXiv:1009.2197 [hep-lat].[9] G. Parisi and Y. Wu, Sci.Sin. , 483 (1981).[10] P. H. Damgaard and H. Huffel, Phys.Rept. , 227(1987).[11] G. Parisi, Phys.Lett. B131 , 393 (1983).[12] F. Karsch and H. W. Wyld, Phys.Rev.Lett. , 2242(1985).[13] G. Aarts and F. A. James, JHEP , 118 (2012),arXiv:1112.4655 [hep-lat].[14] N. Bilic, H. Gausterer, and S. Sanielevici, Phys.Rev. D37 , 3684 (1988).[15] G. Aarts and I.-O. Stamatescu, JHEP , 018 (2008),arXiv:0807.1597 [hep-lat].[16] G. Aarts and K. Splittorff, JHEP , 017 (2010),arXiv:1006.0332 [hep-lat].[17] G. Aarts, Phys.Rev.Lett. , 131601 (2009),arXiv:0810.2089 [hep-lat].[18] G. Aarts, JHEP , 052 (2009), arXiv:0902.4686 [hep--lat].[19] J. Berges, S. Borsanyi, D. Sexty, and I.-O. Stamatescu,Phys.Rev.
D75 , 045007 (2007), arXiv:hep-lat/0609058[hep-lat].[20] J. Berges and D. Sexty, Nucl.Phys.
B799 , 306 (2008),arXiv:0708.0779 [hep-lat].[21] J. Berges and I.-O. Stamatescu, Phys.Rev.Lett. ,202003 (2005), arXiv:hep-lat/0508030 [hep-lat].[22] G. Aarts and F. A. James, JHEP , 020 (2010),arXiv:1005.3468 [hep-lat].[23] J. Ambjorn, M. Flensburg, and C. Peterson, Nucl.Phys. B275 , 375 (1986).[24] H. W. Hamber and H. Ren, Phys.Lett.
B159 , 330 (1985).[25] J. Flower, S. W. Otto, and S. Callahan, Phys.Rev.
D34 ,598 (1986).[26] E.-M. Ilgenfritz, Phys.Lett.
B181 , 327 (1986).[27] G. Aarts, PoS
LATTICE2012 , 017 (2012),arXiv:1302.3028 [hep-lat].[28] G. Aarts, PoS
LAT2009 , 024 (2009), arXiv:0910.3772[hep-lat].[29] G. Aarts, E. Seiler, and I.-O. Stamatescu, Phys.Rev.
D81 , 054508 (2010), arXiv:0912.3360 [hep-lat].[30] G. Aarts, F. A. James, E. Seiler, and I.-O. Stamatescu,Eur.Phys.J.
C71 , 1756 (2011), arXiv:1101.3270 [hep-lat].[31] G. Aarts, F. A. James, E. Seiler, and I.-O. Stamatescu,PoS
LATTICE2011 , 197 (2011), arXiv:1110.5749 [hep--lat].[32] D. Spielmann,
Aspects of confinement in QCD fromlattice simulations , Ph.D. thesis, Ruprecht–Karls–Universität Heidelberg (2010).[33] J. M. Pawlowski and C. Zielinski, Phys.Rev.
D87 , 094509(2013), arXiv:1302.2249 [hep-lat]. [34] H. Gies and L. Janssen, Phys.Rev.
D82 , 085018 (2010),arXiv:1006.3747 [hep-th].[35] S. Chandrasekharan and A. Li, Phys.Rev.Lett. ,140404 (2012), arXiv:1111.7204 [hep-lat].[36] W. E. Thirring, Annals Phys. , 91 (1958).[37] S. Christofi, S. Hands, and C. Strouthos, Phys.Rev. D75 , 101701 (2007), arXiv:hep-lat/0701016 [hep-lat].[38] S. R. Coleman, Phys.Rev.
D11 , 2088 (1975).[39] G. Benfatto, P. Falco, and V. Mastropietro, Com-mun.Math.Phys. , 713 (2009), arXiv:0711.5010 [hep--th].[40] T. Itoh, Y. Kim, M. Sugiura, and K. Yamawaki,Prog.Theor.Phys. , 417 (1995), arXiv:hep-th/9411201[hep-th].[41] T. D. Cohen, Phys.Rev.Lett. , 222001 (2003),arXiv:hep-ph/0307089 [hep-ph].[42] K. Splittorff and J. J. M. Verbaarschot, Phys.Rev.Lett. , 031601 (2007), arXiv:hep-lat/0609076 [hep-lat].[43] J. B. Kogut and L. Susskind, Phys.Rev. D11 , 395 (1975).[44] T. Banks, L. Susskind, and J. B. Kogut, Phys.Rev.
D13 ,1043 (1976).[45] T. Banks et al. (Cornell-Oxford-Tel Aviv-Yeshiva Collab-oration), Phys.Rev.
D15 , 1111 (1977).[46] L. Susskind, Phys.Rev.
D16 , 3031 (1977).[47] P. Hasenfratz and F. Karsch, Phys.Lett.
B125 , 308(1983).[48] L. Del Debbio and S. Hands, Phys.Lett.
B373 , 171(1996), arXiv:hep-lat/9512013 [hep-lat].[49] J. Han and M. A. Stephanov, Phys.Rev.
D78 , 054507(2008), arXiv:0805.1939 [hep-lat].[50] J. O. Andersen, L. T. Kyllingstad, and K. Splittorff,JHEP , 055 (2010), arXiv:0909.2771 [hep-lat].[51] C. C. Chang, Math. Comp. , 523 (1987).[52] J. Ambjorn and S. K. Yang, Phys.Lett. B165 , 140(1985).[53] G. Aarts, F. A. James, E. Seiler, and I.-O. Stamatescu,Phys.Lett.
B687 , 154 (2010), arXiv:0912.0617 [hep-lat].[54] L. G. Molinari, Linear Algebra and its Applications ,2221 (2008).[55] T. DeGrand and C. E. Detar,
Lattice methods for quan-tum chromodynamics (World Scientific, 2006).[56] W. J. Morokoff and R. E. Caflisch, Journal of Computa-tional Physics , 218 (1995).[57] G. Aarts and F. A. James, PoS
LATTICE2010 , 321(2010), arXiv:1009.5838 [hep-lat].[58] G. Aarts, F. A. James, J. M. Pawlowski, E. Seiler,D. Sexty, et al. , JHEP , 073 (2013), arXiv:1212.5231[hep-lat].[59] E. Seiler, D. Sexty, and I.-O. Stamatescu, (2012),arXiv:1211.3709 [hep-lat].[60] D. Banerjee and S. Chandrasekharan, Phys.Rev.