# A density of states approach to the hexagonal Hubbard model at finite density

Michael Körner, Kurt Langfeld, Dominik Smith, Lorenz von Smekal

AA density of states approach to the hexagonal Hubbard model at ﬁnite density

Michael K¨orner, Kurt Langfeld, Dominik Smith, and Lorenz von Smekal Institut f¨ur Theoretische Physik, Justus-Liebig-Universit¨at, Heinrich-Buﬀ-Ring 16, 35392 Giessen, Germany School of Mathematics, University of Leeds, Leeds, LS2 9JT, UK

We apply the Linear Logarithmic Relaxation (LLR) method, which generalizes the Wang-Landaualgorithm to quantum systems with continuous degrees of freedom, to the fermionic Hubbard modelwith repulsive interactions on the honeycomb lattice. We compute the generalized density of statesof the average Hubbard ﬁeld and divise two reconstruction schemes to extract physical observablesfrom this result. By computing the particle density as a function of chemical potential we assess theutility of LLR in dealing with the sign problem of this model, which arises away from half ﬁlling.We show that the relative advantage over brute-force reweighting grows as the interaction strengthis increased and discuss possible future improvements.

I. INTRODUCTION

Monte Carlo simulations based on path-integral quan-tization are a powerful and widely used tool for the studyof strongly coupled quantum systems. They are appliedin many diﬀerent areas of physics, ranging from high-energy physics, where they are employed e.g. to studythe phase diagram and particle spectrum of QuantumChromodynamics (QCD), to condensed matter physics,where they are used to study strongly correlated electronsystems. Quite often, they are the only way to obtain re-liable information from ﬁrst principles. Unfortunately,their applicability is restricted to a very special class ofsystems, namely those where the path integral exhibitsa positive-deﬁnite measure which can be interpreted asa probability density. This excludes most fermionic sys-tems away from half ﬁlling (unless the complex parts ofthe measure cancel exactly due to an anti-unitary sym-metry) as well as quantum systems which evolve in real(as opposed to Euclidean) time. This restriction is knownas the sign problem and is a long-standing problem of nu-merical physics.In principle, brute-force reweighting techniques can beused to extract information about systems with chargeimbalance from simulations of a corresponding systemat charge neutrality. These are plagued by a severesignal-to-noise ratio problem however, originating froma loss of overlap between the ensembles at zero and non-zero charge density when the thermodynamic limit is ap-proached, and typically fail well below µ/T ≈ µ/T have now beenpushed to µ/T ≈ Linear Logarithmic Relaxation (“LLR”) method,ﬁrst described in Ref. [5], generalizes this idea to quan-tum systems with continuous degrees of freedom. Thegoal of LLR is to estimate the slope a ( X ) = ddX ln( ρ ( X )),where X is some observable and ρ ( X ) is the “generalizeddensity of states” (gDOS). Once a ( X ) is obtained, ρ ( X )can be reconstructed up to a multiplicative factor by nu-merical integration and can be used to compute thermo-dynamic observables. A crucial property of LLR is thatit features exponential error suppression: The estimatefor a ( X ), and by extension of ln( ρ ( X )), can be obtainedwith roughly the same statistical error, independent ofthe exact value of X , even if X is from a region of overalllow weight.In recent years, LLR was successfully applied to obtain ρ ( E ) in SU(2) and U(1) [5] as well as SU(3) [6] gaugetheories and was shown to be eﬀective at dealing withergodicity issues arising at ﬁrst-order phase transitions inU(1) gauge theory [7] and the q = 20 state Potts model[8]. In Ref. [9] LLR was applied to obtain the Polyakovloop distribution for two-color QCD (which has no signproblem) with heavy quarks at ﬁnite densities. To dealwith the sign problem, one needs to compute the gDOSfor the imaginary part of the Euclidean action ρ ( S I ), orsome related observable. This was achieved using LLRfor a Z spin model at ﬁnite charge density [10] and forQCD in the heavy-dense limit [11]. To date, LLR hasnever been applied to a sign problem of a system withfully dynamical fermions however. a r X i v : . [ h e p - l a t ] A ug In this work, we apply LLR to the fermionic Hubbardmodel on the honeycomb lattice away from half ﬁllingwithin a Hybrid Monte Carlo framework. Despite itssimplicity, the Hubbard model, which describes fermionicquasi particles with contact interactions, continues to beof profound interest, as it remains the quintessential ex-ample of an interacting fermion system, and can qualita-tively describe many non-perturbative phenomena suchas dynamical mass-gap generation or superconductivity.On the honeycomb lattice, this model exhibits a secondorder phase transition in the universality class of the chi-ral Gross-Neveu model in 2+1 dimensions [12–15]. Withits relativistic dispersion relation for the low-energy ex-citations in the Dirac-cone region it therefore also pro-vides a convenient lattice regularization, with minimaldoubling, of relativistic theories for chiral fermions withlocal four-fermi interactions such as the Gross-Neveu orNambu-Jona-Lasinio models which are of continued in-terest in searches for inhomogeneous phases [16] as pre-dicted also for the QCD phase diagram, mainly frommean-ﬁeld studies of the NJL model [17]. Extended ver-sions of the hexagonal Hubard model, which include long-range interactions, are also used to realistically describethe physics of both mono- and bilayer graphene [18, 19].Using LLR, here we compute the gDOS for the averageof a real-valued auxiliary ﬁeld, which is introduced in Hy-brid Monte Carlo simulations to transmit inter-electroninteractions. We demonstrate that this result can be usedto reconstruct the fermion density as a function of chem-ical potential. We show that, in its present form, LLRcan probe much further into the ﬁnite density regimethan standard reweighting, that the relative advantage ofLLR grows as the interaction strength is increased, andargue that future improvements might put the van Hovesingularity in the single-particle bands within reach.This paper is structured as follows: We start in Sec. IIby introducing the basic lattice setup and illustrating thesign problem away from half ﬁlling. Subsequently, weintroduce the generalized density of states of the aver-age Hubbard ﬁeld ρ ( s ) in Sec. III. In Sec. IV, we dis-cuss the reconstruction of the particle density n from ρ ( s ). We describe two diﬀerent reconstruction schemes,whereby n ( µ ) is obtained from both the canonical andgrand-canonical partition functions. As a benchmark, weapply both schemes to test data obtained for the non-interacting tight-binding theory. Full LLR calculationsof the interacting theory, including additional numericaldetails, are then presented in Sec. V. We summarize andconclude in Sec. VI. II. LATTICE SETUP AND THE SIGNPROBLEM

We consider the repulsive Hubbard model on thehexagonal (honeycomb) lattice with fermionic creationoperators ˆ c † x ≡ (ˆ c † x, ↑ , ˆ c † x, ↓ ) for two spin components atsite x , which is deﬁned by the eﬀective Hamiltonian for the grand canonical ensemble:ˆ H = − κ (cid:88) (cid:104) x,y (cid:105) (ˆ c † x ˆ c y + h.c.)+ (cid:88) x (cid:16) m s ˆ c † x σ ˆ c x + U ρ x − µ ˆ ρ x (cid:17) . (1)Here κ is the hopping parameter, which quantiﬁes theenergy cost of fermionic quasi-particles propagating be-tween nearest-neighbor sites. Its phenomenological valuein the tight-binding description of the electronic prop-erties of graphene on a substrate is κ ≈ . κ , which eﬀectively corresponds to setting κ ≡

1. The sum over (cid:104) x, y (cid:105) sums all independent pairsof nearest-neighbors, m s is the sublattice-staggered massterm (with alternating sign on the two triangular sub-lattices) for explicit sublattice-symmetry breaking withspin-density-wave order. The chemical potential µ cou-ples to the charge operator ˆ ρ x = ˆ c † x ˆ c x − U controlsthe interaction strength, which is positive in the repul-sive Hubbard model. The creation and annihilation op-erators satisfy the fermionic anticommutation relations { ˆ c x , ˆ c † y } = δ x,y . Lattice simulations of (1) using HybridMonte-Carlo by now have a long history already [22–37],we thus summarize only the essential steps here. To derive the functional integral representation of thepartition function at inverse temperature β = 1 /T , weﬁrst write the exponential in terms of N t identical factorsand split the Hamiltonian into the free tight-binding partplus interactions, ˆ H = ˆ H + ˆ H int . A symmetric Suzuki-Trotter decomposition of each of the factors then yields Z = Tr (cid:16) e − β ˆ H (cid:17) = Tr (cid:16) e − δ ˆ H e − δ ˆ H int e − δ ˆ H . . . (cid:17) + O ( δ ) . (2)This introduces a ﬁnite step size δ = β/N t in Euclideantime and a discretization error of O ( δ ).As we will see shortly, it is convenient to include thechemical-potential term in the deﬁnition of ˆ H int here,i.e. deﬁning ˆ H int ≡ (cid:88) x (cid:16) U ρ x − µ ˆ ρ x (cid:17) . (3)The four-fermion interaction in ˆ H int is then converted tobilinears by Hubbard-Stratonovich transformation, e − δ ˆ H int ∼ = (cid:90) D φ e − δ U (cid:80) x φ x e − iδ (cid:80) x ( φ x + iµ )ˆ ρ x , (4) In particular we omit the discussion of fermionic coherent statesand the partial particle-hole transformation that is applied.These and other details can be found e.g. in Refs. [25, 33, 37]. whereby the auxiliary (“Hubbard-Coulomb”) ﬁeld φ x,t isintroduced. Finally, the trace over the fermionic opera-tors is performed by integrating the fermionic coherentstates [33], which yields Z = (cid:90) D φ det (cid:2) M ( φ, µ ) M † ( φ, − µ ) (cid:3) exp (cid:40) − δ U (cid:88) x,t φ x,t (cid:41) . (5)Diﬀerent versions of the fermion matrix M ( φ ) have beenused in the past which are either equivalent or at leastyield the same continuum limit. In this work we use M ( φ, µ ) ( x,t ) , ( x (cid:48) ,t (cid:48) ) = δ xx (cid:48) exp { iδ ( φ x,t + iµ ) } δ tt (cid:48) − (cid:16) δ xx (cid:48) − δh xx (cid:48) (cid:17) δ t +1 ,t (cid:48) , (6) h xx (cid:48) = δ xx (cid:48) m s − κ (cid:88) (cid:126)n δ x (cid:48) ,x + (cid:126)n . in which the free tight-binding hopping contributions ofthe form e − δh are linearized, in order to be able to workwith sparse matrices, but the diagonal couplings to Hub-bard ﬁeld and chemical potential are not.It is clear that Eq. (5) is sign-problem free at half ﬁll-ing, i.e. for µ = 0, as det( M M † ) = | det M | . This is nolonger true for µ (cid:54) = 0. By writing Z = (cid:90) D φ (cid:12)(cid:12) det M ( φ, µ ) (cid:12)(cid:12) det M ( φ, µ )det M ( φ, − µ ) × exp (cid:8) − δ U N t − (cid:88) t =0 (cid:88) x,y φ x,t (cid:9) , (7)we can consider the complex ratio of determinants withunlike-sign chemical potentials as an observable in the“phase-quenched” theory (deﬁned by the modulus of thefermion determinant) with partition function Z pq andobtain Z ( µ ) Z pq ( µ ) = (cid:68) det M ( φ, µ )det M ( φ, − µ ) (cid:69) pq . (8)This ratio is unity for µ → of the ratiodet M ( φ, µ ) / det M ( φ, − µ ) for diﬀerent non-zero valuesof µ , obtained on a lattice of N s × N s unit cells, with 2sites per unit cell, and N s = N t = 6, at β = 2 . κ − and U = κ/

2, together with ﬁt-model curves. The adjusted R of a constant ﬁt (corresponding to a uniform distribu-tion) shows a rather rapid crossover and approaches val-ues close to 1 at µ ≈ . κ . This indicates that the signal The modulus also deviates from unity at µ (cid:54) = 0 but need not beconsidered here. μ = κμ = κμ = κ - π - π - π - π π π π π - π - π - π - π π π π π phase p r o p a b ili t y d i s t r i b u t i o n ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● Adjusted R² μ in κ FIG. 1. Histograms of the phase of det M ( φ, µ ) / det M ( φ, − µ )for N s = N t = 6, β = 2 . κ − , U = κ/ µ . Theresults are modelled with Gaussian (for µ = 0 . κ and0 . κ ) and uniform (for µ = 0 . κ ) distributions respec-tively. The inlay shows the adjusted R for a constant ﬁt todata at diﬀerent µ . For µ (cid:38) . κ the numerical data is welldescribed by a uniform distribution, indicating a hard signproblem. An analogous ﬁgure as indication of a sign problemis obtained for graphene with long-range interactions [31]. is lost in the noise rather quickly already on small lat-tices (signalling a hard sign problem), and for rather hightemperatures and weak interaction strengths. This eﬀectis enhanced with larger lattice sizes, lower temperaturesand larger couplings. To present a quantitative compar-ison of brute-force reweighting and the LLR method fordiﬀerent sysem sizes and interaction strengths is one ofthe main objectives of this work. III. GENERALIZED DENSITY OF THEHUBBARD FIELD

Assume we have a quantum system with action βS ( φ ).Deﬁning the density of states ρ ( E ) as ρ ( E ) = (cid:90) D φ δ ( S ( φ ) − E ) , (9)we can then express the partition function as Z = (cid:90) dE ρ ( E ) e βE , (10)and the vacuum expectation value of an observable O ( E )becomes (cid:104) O (cid:105) = 1 Z (cid:90) dE O ( E ) ρ ( E ) e βE . (11)If ρ ( E ) is known, then (cid:104) O (cid:105) can be obtained by (numeri-cally or analytically) integrating Eq. (11). Here we haveassumed that we know how to express O in terms of E , which is in general not the case however. Moreover,Eq. (9) is ill-deﬁned if S ( φ ) is not strictly real. To com-pute diﬀerent observables in a generic setting, the con-cept of the density of states can thus be generalized toquantities other than the action.The basic idea of LLR is to obtain ρ ( X ) (where X issome observable) by carrying out a sequence of “micro-canonical” Monte-Carlo simulations, in which X is forcedto assume a set of diﬀerent (suﬃciently dense) values X i .Obtaining the partition function or thermodynamic ex-pectation values then essentially amounts to computinga Fourier or Laplace transform of ρ ( X ). To alleviate thesign problem, ρ ( X ) must reﬂect the phase ﬂuctuationsof the path-integral measure. To this end, we consider ρ ( s = Φ) in this work, where Φ is the spacetime-volumeaverage of the auxiliary ﬁeld, see below. First, we applythe transformation φ x,t → φ x,t − iµ (12)to Eq. (5), which then leads to Z ( µ ) == (cid:90) D φ det [ M ( φ, µ ) M † ( φ, − µ )] exp (cid:40) − δ U (cid:88) x,t φ x,t (cid:41) = (cid:90) D φ | det M ( φ, | exp (cid:40) − δ U (cid:88) x,t ( φ x,t − iµ ) (cid:41) . (13)In this formulation, the complex part of the action hasbeen shifted completely to the bosonic sector. Eq. (13)is now rewritten as Z ( µ ) = (cid:90) D φ | det M ( φ, | × exp (cid:40) − δ U (cid:88) x,t ( φ x,t − Φ) − δV U (Φ − iµ ) (cid:41) , (14)where we have introduced the average Hubbard ﬁeldΦ = 1 V (cid:88) x,t φ x,t , V = 2 N c N t , (15)where N c denotes the number of unit cells with 2 siteseach. Finally, we introduce ρ ( s ) = (cid:90) D φ | det M ( φ, | δ (Φ − s ) × exp (cid:40) − δ U (cid:88) x,t ( φ x,t − s ) (cid:41) , (16)and rewrite the partition function as Z ( µ ) = (cid:90) ds ρ ( s ) exp (cid:26) − N c U T ( s − iµ ) (cid:27) , (17)where ρ ( s ) is the generalized density of states of the av-erage Hubbard ﬁeld Φ and will be the target of our LLRcalculation. IV. RECONSTRUCTING THE PARTICLEDENSITY

Assume we have obtained ρ ( s ) using some method.Due to the oscillating contribution of the exponential,it is clear that Eq. (17) will be hard, if not impossible,to evaluate numerically. This is exacerbated by the factthat LLR obtains ρ ( s ) only for a discrete and ﬁnite setof points and with ﬁnite numerical precision. Our ul-timate goal is to obtain the particle density n ( µ ). Wepresent two reconstruction schemes which achieve this inthe following, which both operate in the in the frequencydomain and avoid the instabilities of Eq. (17).We note that ρ ( s ) has a periodicity of 2 π/β = 2 πT andcan thus be expanded in Fourier series. For later conve-nience we will ﬁrst introduce a dimensionless variable anddensity via x = sT , and ¯ ρ ( x ) = T ρ ( xT ) . (18)If we furthermore introduce an imaginary chemical po-tential via µ = iθT , and Z I ( θ ) ≡ Z ( iθT ) , (19)we observe that up to a Gaussian smearing with variance U/ N c T the generalized density of states is in fact es-sentially the same as the partition function at imaginarychemical potential, Z I ( θ ) = (cid:90) dx ¯ ρ ( x ) exp (cid:26) − N c TU ( x − θ ) (cid:27) . (20)We will obtain ¯ ρ ( x ) only at a discrete set of points x n =2 πn/K , where n = { , . . . , K − } and ¯ ρ n ≡ ¯ ρ ( x n ). Wemust hence truncate the Fourier series, naturally leadingto a discrete Fourier transform which can be used forinterpolation via˜ ρ k = 1 K K − (cid:88) n =0 ¯ ρ n e πi nk/K , ¯ ρ ( x ) ≈ K − (cid:88) k =0 ˜ ρ k e − ikx . (21)On the other hand, inserting (21) into (17), we obtain Z ( µ ) ≈ (cid:90) dx (cid:32) K − (cid:88) k =0 ˜ ρ k e − ikx (cid:33) exp (cid:26) − N c TU (cid:16) x − i µT (cid:17) (cid:27) . = (cid:114) πUN c T K − (cid:88) k =0 ˜ ρ k exp (cid:26) − U N c T k − µT k (cid:27) , (22)and the exact result is recovered for K → ∞ . In fact, inthis limit, Eq. (22) becomes the fugacity expansion andwe can identify for k = N , Z c ( T, N ) = (cid:114) πUN c T ˜ ρ N exp (cid:26) − U N c T N (cid:27) (23)as the corresponding canonical partition function withparticle number N . In the inﬁnite volume limit N c → ∞ for ﬁxed N , or equally so for T (cid:29) U , we may thereforeneglect the exponential factor and essentially identify theFourier series coeﬃcients ˜ ρ k of our generalized DOS withthe canonical partition functions at N = k . At the sametime, it is also evident from Eq. (22) that the generalizedDOS itself then becomes equal, up to a constant factor,to the partition function at imaginary chemical potential,i.e. ˜ ρ k ∝ Z c ( T, k ), and ¯ ρ ( θ ) ∝ Z I ( θ ) or ρ ( s ) ∝ Z ( is ).Moreover, one easily veriﬁes that the truncated coeﬃ-cients ˜ ρ k at ﬁnite K , obtained from the discrete Fouriertransform in (21), then yield pseudo-canonical partitionfunctions, ˜ ρ k ∝ Z Kc ( T, k ), which represent ensembleswith particle number N = k mod K . Likewise, the dis-crete sampling of ¯ ρ ( θ ) provides us with an interpolationof Z I ( θ ) which agrees with the exact result for imaginarychemical potential at the discrete values θ n = 2 πn/K .The general relation between ρ ( s ) and the partitionfunction at imaginary chemical potential of course alsofollows from Eq. (22), with µ = is (and K → ∞ ), Z ( is ) = (cid:114) πUN c T ∞ (cid:88) k = −∞ ˜ ρ k exp (cid:26) − U N c T k (cid:27) e − isk/T → (cid:114) πU TN c ρ ( s ) , UN c T → . (24)In a ﬁnite volume, on the other hand, i.e. at any ﬁ-nite number N c of unit cells, the particle numbers N arerestricted to values between ± N c , with N = 0 at halfﬁlling, corresponding to an average of one of the maxi-mally possible two electrons on each of the 2 N c sites. Wethen obtain the exact canonical partition functions fromEq. (23) already for K = K max = 4 N c + 1 , and with particle-hole symmetry at half ﬁlling, one actu-ally only needs K max = 2 N c + 1.In principle, the particle number N ( µ ) can be directlyobtained from Eq. (22), which is free of oscillating terms,by taking the deritvative with respect to µ , N ( µ ) = T ddµ ln Z ( µ ) (25)= − (cid:80) N c k =0 k ˜ ρ k exp (cid:110) − U N c T k − µT k (cid:111)(cid:80) N c k =0 ˜ ρ k exp (cid:110) − U N c T k − µT k (cid:111) . Computing the ˜ ρ k from ρ ( s n ) can be done with high nu-merical precision using modern FFT libraries.Alternatively, we can also compute the chemical po-tential from the canonical partition functions, as the freeenergy diﬀerence of ensembles with subsequent particlenumbers. From Eq. (22) we then obtain µ ( N + 1 /

2) = − T (cid:0) ln Z c ( T, N + 1) − ln Z c ( T, N ) (cid:1) (26) ≈ T (cid:20) ln ˜ ρ N − ln ˜ ρ N +1 + U N c T (cid:0) N + 1 / (cid:1)(cid:21) , and obtain the density in form of the number of particlesper unit cell, n ( µ ) ≡ N ( µ ) /N c , by inversion. The exactcalculation would again require K max = 2 N c + 1 Fouriercoeﬃcients. This is then similar in spirit to Refs. [38–40], which carried out canonical calculations of QCD atﬁnite charge density, or Ref. [41] which followed essen-tially the same strategy for ﬁnite isospin density fromthe lowest states in multi-pion correlators. With truncat-ing at K < K max , we strictly speaking obtain canonicalensembles at particle number N modulo K as discussedabove. The term ∝ U/ N c T in Eq. (26) represents an ex-plicit ﬁnite volume eﬀect which, as we will discuss below,only contributes in trivial way and can be dropped.In tight-binding or mean-ﬁeld calculations, there is nosuch term in the ﬁrst place, and the generalized DOS canbe calculated analytically. The result is of the formln ρ ( s ) = 2 N c (cid:90) dε ρ ε ( ε ) ln (cid:16) cosh ε T − sin s T (cid:17) , (27)where ε ≥ ρ ε ( ε ) for which an analytic expression is knownin the inﬁnite system [42]. In a ﬁnite system with pe-riodic (Born-von K´arm´an) boundary conditions we usethe dispersion relation ε = ε ( k ) instead, and simply sumover the corresponding discrete set of points k n withinthe ﬁrst Brillouin zone, with energies ε n = ε ( k n ). Thesame can be done to compute the exact density in theﬁnite system with N c unit cells which then yields for thenumber of particles per unit cell, n ( µ ) = 1 N c (cid:88) n (cid:16) tanh ε n + µ T − tanh ε n − µ T (cid:17) . (28)We have carried out a set of benchmark calculationsin which we compared the canonical and grand-canonicalreconstruction schemes. Thereby, a discete set of values l n = ln ρ ( s n ) for s n = 2 πT n/K , N = { , . . . K − } , wasproduced as mock data from the tight-binding calcula-tion, which can eﬃciently be done with arbitrary numer-ical precision. High-precision calculations are especiallyimportant in the reconstrucion of the density because weneed with high precision the discrete Fourier transform of { ρ n = e l n } rather than that of { l n } . The number density n ( µ ) was subsequently computed from the FFT result { ˜ ρ k } , using both the fugacity expansion via (25) and thecanonical approach (26). We have then compared bothresults with the exact calculation of the density based onthe tight-binding formula (28). This was done for dif-ferent setups, whereby the production of { l n = ln ρ ( s n ) } was done with diﬀerent levels of ﬂoating point precision.The application of the reconstruction scheme was donewith a 1024 digit accuracy in each case to avoid additionalerrors. We ﬁnd that both methods yield comparable re-sults, with the canonical procedure having a very slightadvantage for a given precision of ln ρ ( s ). We thus chooseto use this procedure exclusively in the following sectionsto process our LLR results.Fig. 2 shows an example calculation of n ( µ ) for twodiﬀerent temperatures on a lattice with N c = 36 unit β = κ - , m s = κ , Volume = β = κ - , m s = κ , Volume = μ in κ n ( μ ) FIG. 2. The number of particles n ( µ ) per unit cell on a latticewith N c = 36 unit cells for two diﬀerent temperatures at U = 0 from the tight-binding calculation (solid lines) andfrom the canonical reconstruction procedure based on (26)using input data of 1024 digit accuracy (discrete points). cells, where ln ρ ( s ) was produced for U = 0 with 1024digit accuracy and processed using (26). This illustratesthat our method can in principle cover the entire widthof the valence band, from the empty valence band at halfﬁllg up to saturation when it is completely ﬁlled. Thevan Hove singularity will emerge at µ = κ in the inﬁnitevolume limit which can here be anticipated already bythe rapid increase in the number density at the lowertemperature around µ = κ .In practice, the leading source of errors is of course theprecision with which the Fourier coeﬃcients ˜ ρ k can beobtained, which in turn is highly sensitive to statisticalerrors of ln ρ ( s ) in our LLR calculations. In order to getto saturation, with a completely ﬁlled lattice, we wouldobviously need the maximum number K max = 2 N c + 1of coeﬃcients. So the double challenge here will be tocompute as many of them as accurately as possible. V. LLR RESULTSA. The algorithm

The goal of LLR is to calculate derivatives of ln ρ ( s )at a suﬃciently dense set of supporting points with highprecision and to reconstruct ρ ( s ) by integration. We di-vide the domain of support of ρ ( s ) into K intervals ofsize δ s . At the center of each of these intervals, the slope a k = dds ln ρ ( s ) | s = s k can be calculated from a stochasticnon-linear equation [5]. A key element of this equation is the restricted and reweighted expectation value (cid:104)(cid:104) W (Φ)( a ) (cid:105)(cid:105) k = 1 Z LLR (cid:90) D φ θ [ s k ,δ s ] (Φ) | det M ( φ ) | × W (Φ) e − βS ( φ ) e − a Φ . (29)Here Z LLR is a normalization constant, Φ was introducedin Eq. (15), a is an external parameter and θ [ s k ,δ s ] is awindow function which restricts Φ to an interval of size δ s around s k .With the choice W (Φ) = Φ − s k , the coeﬃcients a k aresolutions of (cid:104)(cid:104) W (Φ)( a k ) (cid:105)(cid:105) k = 0 . (30)This equation can be solved through Robbins-Monro it-eration [43]: The sequence a ( n +1) k = a ( n ) k + α n δ s (cid:104)(cid:104) W (Φ)( a ( n ) k ) (cid:105)(cid:105) k (31)converges to the correct result for any choice of α n thatfulﬁlls ∞ (cid:88) n =0 α n = ∞ , ∞ (cid:88) n =0 α n < ∞ . (32)This is true, even if (cid:104)(cid:104) W (Φ)( · ) (cid:105)(cid:105) k is approximated by anestimator, as we do in Monte-Carlo calculations. More-over, if the iteration is terminated at some ﬁnite number N and repeated many times, the ﬁnal values a ( N ) k areGaussian distributed around the true value a k and canbe processed by a standard bootstrap analysis.The window function can be chosen in diﬀerent ways.The straight-forward choice is a step function, but forHMC a Gaussian window function is more appropriate,as its derivative can be taken, which implies that its eﬀectcan be reproduced by a molecular-dynamics force term.In this work, we choose (cid:104)(cid:104) Φ − s (cid:105)(cid:105) ( a ) = 1 Z LLR (cid:90) D φ det M ( φ ) det M † ( φ ) (Φ − s ) × exp (cid:40) − δ U (cid:88) x,t ( φ x,t − s ) − δ s ( s − Φ) − a Φ (cid:41) , (33)where Z LLR ( a ) = (cid:90) D φ det M ( φ ) det M † ( φ ) × exp (cid:40) − δ U (cid:88) x,t ( φ x,t − s ) − δ s ( s − Φ) − a Φ (cid:41) . (34)The full procedure is then summarized as follows: The double-bracket notation is customary in the LLR literatureand should be understood as deﬁned by Eq. (29). It is not impliedhere that an expectation value is taken twice. - n a n Volume 6x6x6, β = κ - , m s = κ , U = κ , s = κ FIG. 3. Illustration of stochastic Robbins-Monro iteration.A set of 20 starting values a (0) k are generated and are eachupdated according to Eq. (31). Underrelaxation is switchedon at n = 15. The procedure is terminated at n = 105 toobtain the ﬁnal values used for bootstrapping.

1) For a given s k , initialize a k with some random value a (0) k not too far from zero.2) Initialize Hubbard ﬁeld (e.g with a value whichminimizes the window function).3) With ﬁxed a k , thermalize Hubbard ﬁeld with HMCtrajectories according to Eq. (34), i.e. Z LLR ( a k ).4) With additional HMC trajectories, compute an es-timate of (cid:104)(cid:104) Φ − s (cid:105)(cid:105) ( a k ).5) Update a k using Eq. (31).6) Continue from step 3.In practice, we start with several repetitions of steps 3 − α = 1 and switch to underrelaxed interationswith α n +1 = α n / (1+ n ) after some time. Also, the wholeprocedure is terminated after some ﬁnite iteration num-ber N and repeated several times for each ﬁxed s k , toproduce a sample of ﬁnal a ( N ) k values.Fig. 3 shows one example of a stochastic Robbins-Monro iteration, taken from our actual production runs,where the procedure described above is applied for a ﬁxedset of external parameters. We choose N s = N t = 6, β = 2 . κ − , m s = 0 . κ , U = 1 . κ , s = 1 . κ forillustration. For each set of parameters considered inthis work, we ﬁrst obtain such a sample of a k values. Wethen obtain ln ρ ( s k ), and by extension ρ ( s k ), ˜ ρ k and n ( µ ) Note that the phenomenological value of the hopping parameterin the tight-binding model for graphene typically is κ ≈ . T ≈ together with errorbars, by feeding bootstrap averages ofthe ﬁnal a ( N ) k intoln ρ ( s k ) = k − (cid:88) i =0 a i δ s + 12 a k δ s , (35)computing the Fourier transform of ρ ( s k ) and applyingthe canonical reconstruction scheme described in Sec. IV. B. N t dependence We begin by studying the eﬀect of the time-discretization δ . To this end, we carry out LLR calcula-tions at N s = 6, β = 2 . κ − , U = 1 . κ , m s = 0 . κ for diﬀerent values of N t . Fig. 4 shows the results for a ( s ), ln ρ ( s k ) and ln ˜ ρ k , while the ﬁnal results for n ( µ )are shown in Fig. 5. The latter ﬁgure includes twosubﬁgures, whereby the linear ( ∼ U ) contribution toEq. (26) is included or neglected respectively. Fig. 5 alsoshows a corresponding calculation of n ( µ ) in the non-interacting tight-binding theory. All errorbars were ob-tained through bootstrap analysis.Our ﬁrst observation is that the dependence on N t is very mild for our choice of parameters. It is practi-cally invisible in a ( s ) and n ( µ ). A very small diﬀerencebetween diﬀerent N t can be seen in ln ρ ( s k ) and ln ˜ ρ k ,which is of a similar magnitude as the statistical uncer-tainty however. On the other hand, our results clearlydemonstrate exponential error suppression, whereby therelative error of ln ρ ( s ) is roughly the same across sev-eral orders of magnitude. We ﬁnd that ln ˜ ρ k is extremelysensitive to this small error however, to a degree thatonly the ﬁrst few Fourier modes ln ˜ ρ k can be computedaccurately. This can be traced back to the fact that ρ ( s )enters into the Fourier transform and not ln ρ ( s ). It isalso reﬂected in our computation of n ( µ ), which exhibitsa loss of signal at µ ≈ . κ , indicating the onset of ahard sign problem.We note that for U = 1 . κ which is well inside theweak-coupling phase of the model, and the tempera-ture considered here, n ( µ ) basically fully agrees with theinﬁnite-volume limit in the non-interacting theory whenthe linear term in Eq. (26) is dropped. We take this as anindication that this extra term represents the dominantﬁnite-volume eﬀect at ﬁnite U which however is a rathertrivial one to correct. Further conﬁrmation of this is pro-vided by a comparison between results from N s = 6 and N s = 12 lattices, which also reveals a faster convergenceto the thermodynamic limit without this term. We thusdrop this term for all results presented in the following.We expect that deviations from the non-interacting limitwill become visible at stronger couplings, of course. Thisis investigated further, and ultimately conﬁrmed, below. N t = N t = N t = π / π π / π - - κ a ( s ) Volume 6x6xN t , β = κ - , U = κ , m s = κ N t = N t = N t = π / π π / π - - -

50 s in κ l n ρ ( s ) Volume 6x6xN t , β = κ - , U = κ , m s = κ N t = N t = N t = - - - -

50 k l n ρ k Volume 6x6xN t , β = κ - , U = κ , m s = κ FIG. 4. LLR result: N t dependence of a ( s ), ln ρ ( s ), ln ˜ ρ k for N s = 6, β = 2 . κ − , U = 1 . κ , m s = 0 . κ . Individualbootstrap averages are shown for ln ˜ ρ k to illustrate loss ofsignal for the higher modes. C. m s dependence We now turn to studying the dependence on the ex-plicit sublattice and spin-staggered mass term m s . Giventhat such a term already opens an explicit gap in the en-ergy spectrum, we carry out this study at the compara-tively weak coupling strength of U = 0 . κ . We ﬁnd that, N t = N t = N t = - - - - - μ in κ n ( μ ) Volume 6x6xN t , β = κ - , U = κ , m s = κ N t = N t = N t = - - - - - μ in κ n ( μ ) Volume 6x6xN t , β = κ - , U = κ , m s = κ FIG. 5. LLR result: N t dependence of n ( µ ) for N s = 6, β = 2 . κ − , U = 1 . κ , m s = 0 . κ . Top ﬁgure includesthe linear ( ∼ U ) term in Eq. (26), bottom ﬁgure does not.Errorbars computed by boostrapping. Solid line shows thenon-interacting tight-binding theory. again, the number density n ( µ ) coincides with the non-interacting theory and shows no signiﬁcant dependenceon m s . The linear term in Eq. (26) has a negligible eﬀecthere, due to the small value of U . Fig. 6 shows the resultsfor a ( s ) and ln ρ ( s ) with N s = N t = 6, β = 2 . κ − andthree diﬀerent choices of m s . We refrain from showingany additional ﬁgures for ln ˜ ρ k and n ( µ ), as these fullyagree (within statistical errors) with the results shown inthe lowest panels of Figs. 4 and 5.An interesting observation here is that m s has a quitestrong eﬀect on both a ( s ) and ln ρ ( s ), which turns outnot to carry over to n ( µ ) at all. The underlying reasonis that this dependence is only present in regions where ρ ( s ) is strongly suppressed. It is only visible due to thelogarithmic scale, and thus has no signiﬁcant eﬀect onthe computation of the Fourier modes. D. U dependence Having validated our numerical procedure at weak cou-pling, we now turn to a more detailed study of the depen- m s = κ m s = κ m s = κ π / π π / π - - - κ a ( s ) Volume 6x6x6, β = κ - , U = κ m s = κ m s = κ m s = κ π / π π / π - - - - - -

50 s in κ l n ρ ( s ) Volume 6x6x6, β = κ - , U = κ FIG. 6. LLR result: m s dependence of a ( s ) and ln ρ ( s ) for N s = N t = 6, β = 2 . κ − , U = 0 . κ . dence on the interaction strength U . This represents thecentral part of this work to which the bulk of our com-puting resources were dedicated. We thereby computed a ( s ), ln ρ ( s ), ln ˜ ρ k and n ( µ ) again with β = 2 . κ − , m s = 0 . κ for several diﬀerent choices of U . To havecontrol over ﬁnite volume and time-discretization eﬀectswe have studied two diﬀerent lattice sizes, N s = N t = 6and N s = N t = 12, respectively.Fig. 7 shows the U dependence of a ( s ), ln ρ ( s ) and ln ˜ ρ k for N s = N t = 6, β = 2 . κ − , m s = 0 . κ . For ln ρ ( s )we include the tight-binding result to illustrate the ap-proach to the non-interacting limit. The ﬁrst observationis that a ( s ) gets suppressed when U is increased, whichultimately makes simulations more expensive at strongcoupling. On the other hand, we clearly see a devia-tion from the non-interacing limit in the Fourier modesln ˜ ρ k for the strongest interaction strength U = 2 . κ . Tounderscore that this deviation is absent for all weakerinteractions, we show a seperate plot in Fig. 8 which di-rectly compares ln ˜ ρ k for U ≤ . κ to the tight-bindingtheory. Our N s = N t = 6 results for n ( µ ) are shown inFig. 9. They clearly show a corresponding drop of thenumber density at ﬁxed µ for the strongest coupling. U = κ U = κ U = κ U = κ U = κ π / π π / π - - κ a ( s ) Volume 6x6x6, β = κ - , m s = κ U = κ U = κ U = κ U = κ U = κ TB π / π π / π - - - -

50 s in κ l n ρ ( s ) Volume 6x6x6, β = κ - , m s = κ U = κ U = κ U = κ U = κ U = κ - - - - -

50 k l n ρ k Volume 6x6x6, β = κ - , m s = κ FIG. 7. LLR result: U dependence of a ( s ), ln ρ ( s ), ln ˜ ρ k for N s = N t = 6, β = 2 . κ − , m s = 0 . κ . Individual boot-strap averages are shown for ln ˜ ρ k . Result for non-interactingtight-binding theory is included for ln ρ ( s ). N s = 12 results are shown in Fig. 10 for a ( s ), ln ρ ( s )and ln ˜ ρ k and Fig. 11 for n ( µ ). These conﬁrm the quali-tative changes at U = 2 . κ . Furthermore, a direct com-parison with N s = 6 suggests that ﬁnite volume eﬀectson n ( µ ) are rather mild.We point out here that the sign problem sets in at muchsmaller µ for the larger system (as expected). While weare able to reliably compute n ( µ ) up to µ ≈ . κ for N s = 6 with U = 1 . κ , we only reach µ ≈ . κ for0 Tight - BindingU = κ U = κ U = κ U = κ - - - - - - -

50 k l n ρ k Volume 6x6x6, β = κ - , m s = κ FIG. 8. LLR result: ln ˜ ρ k for N s = N t = 6, β = 2 . κ − , m s = 0 . κ and diﬀerent U , compared with non-interactingtight-binding theory. U = κ U = κ U = κ U = κ U = κ - - - - - μ in κ n ( μ ) Volume 6x6x6, β = κ - , m s = κ FIG. 9. LLR result: U dependence of n ( µ ) for N s = N t = 6, β = 2 . κ − , m s = 0 . κ . Errorbars computed by boost-rapping. Solid line shows the non-interacting tight-bindingtheory. N s = 12. On the other hand, in both cases LLR dras-tically outperforms brute-force reweighting: With com-parable numerical resources we obtain a signal for thedeterminant ratio (8) up to µ ≈ . κ on N s = 6 and µ ≈ . κ on N s = 12 using the brute-force method.While the relative advantage of LLR becomes smaller onthe larger lattice, we can reach much larger values of µ for U = 2 . κ ( µ ≈ . κ on N s = 6 and µ ≈ . κ on N s = 12). In contrast, the µ range of reweighting isdrastically diminished at stronger coupling (cf. Fig. 13 inSec. VI). It is this last feature which ultimately makesLLR in its present form a promising method and deserv-ing of further attention. E. Compressed sensing

Lastly, we report on our attempts to improve our re-sults by using ﬁt functions for ln ρ ( s ), a procedure re- U = κ U = κ π / π π / π - - κ a ( s ) Volume 12x12x12, β = κ - , m s = κ U = κ U = κ π / π π / π - - - - -

100 s in κ l n ρ ( s ) Volume 12x12x12, β = κ - , m s = κ U = κ U = κ - - -

50 k l n ρ k Volume 12x12x12, β = κ - , m s = κ FIG. 10. LLR result: U dependence of a ( s ), ln ρ ( s ), ln ˜ ρ k for N s = N t = 12, β = 2 . κ − , m s = 0 . κ . Individualbootstrap averages are shown for ln ˜ ρ k . ferred to as compressed sensing in the LLR literature.The basic idea is to, instead of processing the raw datafor ln ρ ( s ) pointwise at the supporting points s k , ﬁt theentire data set with a series expansion in some completeset of functions, and use the model curve to computeobservables instead. The hope is that an appropriateset of functions, which reﬂects the true (but a priori un-known) physics of the theory, will both suppress noise inthe numerical data for ln ρ ( s ) and eﬀectively generate aninterpolation to a much denser set of supporting points.1 U = κ U = κ - - - - - μ in κ n ( μ ) Volume 12x12x12, β = κ - , m s = κ FIG. 11. LLR result: U dependence of n ( µ ) for N s = N t = 12, β = 2 . κ − , m s = 0 . κ . Errorbars computed by boost-rapping. Solid line shows the non-interacting tight-bindingtheory. This in turn should allow for the computation of ln ˜ ρ k atlarger k and hence the number density at larger µ .Fig. 12 shows two such attempts, where ln ρ ( s ) was ﬁtwith a Fourier series and a series of Chebyshev polyno-mials of the ﬁrst kind respectively. The ﬁt function wassubsequently evaluated at a much denser set of pointsthan the original s k and used to compute ln ˜ ρ k and sub-sequently n ( µ ). In each case, higher order terms wereadded to the expansion until the ﬁnal result stabilized.We do not obtain any errorbars. Results from the directcalculation are included for comparison and representedby dashed lines.We ﬁnd that these attempts do not improve the calcu-lation of n ( µ ) signiﬁcantly. At best, one or two additionalpoints (at higher densities) can be computed before theresults scatter in an uncontrolled fashion. The Fourierseries thereby seems to work only slightly better thanthe Chebyshev polynomials. We take this as an indi-cation that additional qualitative information about thesystem, which places additional contraints on the choiceof functions to use, is a necessary requirement for com-pressed sensing to be eﬀective here. VI. CONCLUSION AND OUTLOOK

In this work, we have applied the Linear LogarithmicRelaxation method to the repulsive fermionic Hubbardmodel on the honeycomb lattice, in order to assess itsutility for alleviating the hard sign problem of an un-balanced dynamical fermion system. A central problemthereby is the proper choice of a target observable, whichadequately reﬂects the complex part of the action andyields a generalized density of states which is suitable forfurther processing. We used the average value Φ of theauxiliary (Hubbard) ﬁeld to this end, which appeared asthe natural choice, as it allows for the shifting of the com-plex part of fermion determinant into the bosonic sector U = κ U = κ U = κ U = κ - - - - μ in κ n ( μ ) Volume 6x6x6, β = κ - , m s = κ U = κ U = κ U = κ U = κ - - - - - μ in κ n ( μ ) Volume 6x6x6, β = κ - , m s = κ FIG. 12. LLR result: µ -dependence of particle density for N s = N t = 6, β = 2 . κ − , m s = 0 . κ and diﬀerent U (results include the linear ∼ U term in Eq. (26)). Dashedlines were obtained directly from ln ˜ ρ k , while dots employed compressed sensing : ln ρ ( s ) was ﬁt with a Fourier series (topﬁgure) and Chebyschev polynomials (bottom ﬁgure) respec-tively. Solid line shows non-interacting tight-binding theory. and provides a simple integral expression (Eq. (17)) forobtaining the partition function and hence the particledensity. To deal with an oscillating contribution to thisintegral, we chose to work in the frequency domain anddivised two methods to extract the particle density fromthe Fourier modes of the gDOS of Φ which essentiallyyields the partition function at imaginary chemical po-tential. Due to a slightly better performance in bench-mark calculations, of these we chose a method based onthe canonical ensembles to further process our LLR re-sults.We have carried out LLR calculations for a ﬁxed tem-perature of β = 2 . κ − , two diﬀerent lattice sizes (6 and12 ) and diﬀerent interaction strengths in the weak andintermediate coupling regime, and obtained the particledensity as a function of chemical potential. We therebyobserved signiﬁcant deviations from the non-interactingtheory for the largest interaction strength considered, U/κ = 2 .

0, signalling strong correlations which mighteventually lead to spontaneous mass-gap formation which2

X XX X brute force N C = N C = N C = N C = μ in κ U i n κ FIG. 13. Comparison of eﬀective µ -range of brute-forcereweighting and LLR (results shown for 6 and 12 lattices at β = 2 . κ − ). For each value of U the phase distribution ofdet (cid:102) M ( φ, µ ) / det (cid:102) M ( φ, − µ ) was measured at diﬀerent µ untilthe signal was lost. A roughly equal amount of computer timewas spent for corresponding LLR calculations. is known to occur at around U/κ ≈ . lattice we are able to probe at least twiceas far into the ﬁnite-density regime as with brute-forcereweighting. While the relative advantage of LLR issmaller on the 12 lattice, we ﬁnd that LLR performsmuch better when the interaction strength is increased.Fig. 13 shows a quantitative comparison of the eﬀective µ -ranges for the diﬀerent parameter sets considered inthis work.Attempts to reach into higher-density regions weremade using diﬀerent forms of compressed sensing, i.e. byﬁtting ln ρ ( s ) with Fourier series and Chebyshev poly-nomials and using the model curves for interpolation.While this allows us to reach slightly higher densities, wesuspect that this procedure introduces an uncontrolledsystematic error, as the physics at higher densities isstrongly sensitive to the high-frequency modes of ρ ( s ),which such interpolations cannot account for.The results presented here should be taken as a proofof principle. There are several diﬀerent directions forfuture improvements. First and foremost, the computa-tional resources spent for the ﬁnal calculations (i.e. thesum of all results shown in Figs. 9 and 11) were aroundtwo months of runtime on a total of 18 GTX 980 TiGPUs, which leaves much space for larger-scale projects.We estimate that, using the most modern hardware andlibraries for sparse linear algebra, the precision for ln ρ ( s )can be increased by at least an order of magnitude. Moreadvanced techniques for compressed sensing could alsobe applied, such as Gaussian and telegraphic approxima-tions or an advanced moments approach, which were pro-posed in Ref. [44]. Quite possibly, introducing a complex instead of a real auxiliary ﬁeld has advantages, and infact, it was shown that an optimal mixing factor betweenreal and imaginary Hubbard ﬁelds exists, for which thesign problem is the mildest [45]. LLR might also be moreeﬀective with a discrete Hubbard ﬁeld, which is usedin BSS Quantum Monte-Carlo calculations [46, 47]. Inaddition, an alternative time-discretization with a sym-metry of time reversal times sublattice exchange, whichwas proposed already in Ref. [22] and recently used ina grand canonical HMC simulation [48], was shown tohave strongly suppressed discretization eﬀects in Ref. [30]and might positively impact the performance of LLR aswell. And ﬁnally, there has been much recent progressregarding the Lefschetz thimble method [45, 49, 50], andconstructing a hybrid approach, which combines the ad-vantages of both methods, might be feasible. Speciﬁcally,one could attempt to apply the Lefschetz thimble decom-position directly to Eq. (17), in order to avoid the use ofreconstruction schemes altogether and obtain a cleanersignal for n ( µ ).Taken together, we ﬁnd it not unreasonable to expectthat future developments might put the van Hove singu-larity (VHS) of the single-particle spectrum within reach,which is of great interest in the context of superconduct-ing phases. A crucial point thereby is the apparent stabil-ity of the LLR technique against increases of the couplingstrength U . Experiments on charge-doped graphene sys-tems have revealed a strong bandwidth renormalization(narrowing of the width of the π -bands) due to interac-tions [20], which suggests that the VHS can be probedat smaller µ for larger U . Furthermore, a HMC studyof graphene at ﬁnite spin density revealed that the elec-tronic Lifshitz transition at the VHS can become a truethermodynamic phase transition in the presence of in-teractions, with a critical temperature which increaseswith the coupling strength [31]. A study of an analogoustransition at ﬁnite charge-carrier density thus might befeasible at large U , in particular as the sign problem be-comes milder at higher temperatures.There are many possibilities to linearize the quarticfermionic interaction using auxiliary ﬁelds. The choicein this paper is inspired by the observation that an ex-plicit analytic continuation (i.e. Eq. (12)) was suﬃcientto split oﬀ the complex part of the fermion determinant.The formulation is elegant: The calculation of one (real)density of states, ρ ( s ), is suﬃcient to relay the calculationof the partiton function to one integration for each givenvalue of the chemical potential. Note, however, that theuse of a Hubbard ﬁeld φ with a compact domain of sup-port implies that the domain of the density of states isalso compact. The calculation of such an “intensive” den-sity of states to suﬃcient precision is diﬃcult [44] and themost successful LLR calculations for theories with a signproblem are based upon non-compact densities [10]. Theuse of a non-compact formulation is left to future work.Lastly, we should mention that extending our work tothe QCD sign problem remains an open conceptual chal-lenge. The system considered here was special since we3succeeded to remove the complex part the fermion de-terminant by a simple analytic continuation. For gaugetheories no such simple transformation exists, and mea-suring a proper extensive phase is much more involved.It may well be that this step is the most computation-ally demanding and contains the central computationalcomplexity of the QCD sign problem. ACKNOWLEDGMENTS

This work was supported by the Helmholtz Interna-tional Center known as HIC for FAIR and its successor, the Helmholtz Research Academy Hessen for FAIR.We are gratefull to Pavel Buividovich, John Graceyand Maksim Ulybyshev for helpful discussions and com-ments. [1] M. Troyer and U.-J. Wiese,

Computational complexityand fundamental limitations to fermionic quantumMonte Carlo simulations , Phys. Rev. Lett. (2005)170201, [ cond-mat/0408370 ].[2] A. Bazavov et al., The QCD Equation of State to O ( µ B ) from Lattice QCD , Phys. Rev.

D95 (2017)054504, [ ].[3] B. A. Berg and T. Neuhaus,

Multicanonical ensemble: ANew approach to simulate ﬁrst order phase transitions , Phys. Rev. Lett. (1992) 9–12, [ hep-lat/9202004 ].[4] F. Wang and D. P. Landau, Eﬃcient, multiple-rangerandom walk algorithm to calculate the density of states , Phys. Rev. Lett. (Mar, 2001) 2050–2053,[ cond-mat/0011174 ].[5] K. Langfeld, B. Lucini and A. Rago, The density ofstates in gauge theories , Phys. Rev. Lett. (2012)111601, [ ].[6] C. Gattringer and K. Langfeld,

Approaches to the signproblem in lattice ﬁeld theory , Int. J. Mod. Phys.

A31 (2016) 1643007, [ ].[7] K. Langfeld, B. Lucini, R. Pellegrini and A. Rago,

Aneﬃcient algorithm for numerical computations ofcontinuous densities of states , Eur. Phys. J.

C76 (2016)306, [ ].[8] K. Langfeld,

Density-of-states , PoS Lattice2016 (2017)010, [ ].[9] K. Langfeld and J. M. Pawlowski,

Two-color QCD withheavy quarks at ﬁnite densities , Phys. Rev.

D88 (2013)071502, [ ].[10] K. Langfeld and B. Lucini,

Density of states approachto dense quantum systems , Phys. Rev.

D90 (2014)094502, [ ].[11] N. Garron and K. Langfeld,

Anatomy of thesign-problem in heavy-dense QCD , Eur. Phys. J.

C76 (2016) 569, [ ].[12] F. F. Assaad and I. F. Herbut,

Pinning the order: thenature of quantum criticality in the Hubbard model onhoneycomb lattice , Phys. Rev. X3 (2013) 031010,[ ].[13] Y. Otsuka, S. Yunoki and S. Sorella, Universal quantumcriticality in the metal-insulator transition oftwo-dimensional interacting dirac electrons , Phys. Rev.X (Mar, 2016) 011029.[14] F. Parisen Toldin, M. Hohenadler, F. F. Assaad andI. F. Herbut, Fermionic quantum criticality in honeycomb and π -ﬂux hubbard models: Finite-sizescaling of renormalization-group-invariant observablesfrom quantum monte carlo , Phys. Rev. B (Apr,2015) 165108.[15] M. Hohenadler, F. Parisen Toldin, I. F. Herbut andF. F. Assaad, Phase diagram of the kane-mele-coulombmodel , Phys. Rev. B (Aug, 2014) 085146.[16] J. Lenz, L. Pannullo, M. Wagner, B. Wellegehausen andA. Wipf, Inhomogeneous phases in the Gross-Neveumodel in 1+1 dimensions at ﬁnite number of ﬂavors , Phys. Rev. D (2020) 094512, [ ].[17] M. Buballa and S. Carignano,

Inhomogeneous chiralcondensates , Prog. Part. Nucl. Phys. (2015) 39–96,[ ].[18] T. O. Wehling, E. S¸a¸sıo˘glu, C. Friedrich, A. I.Lichtenstein, M. I. Katsnelson and S. Bl¨ugel, Strengthof eﬀective coulomb interactions in graphene andgraphite , Phys. Rev. Lett. (Jun, 2011) 236805.[19] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi,K. Kuroki and L. Fu,

Maximally localized wannierorbitals and the extended hubbard model for twistedbilayer graphene , Phys. Rev. X (Sep, 2018) 031087.[20] S. Ulstrup, M. Sch¨uler, M. Bianchi, F. Fromm,C. Raidel, T. Seyller et al., Manifestation of nonlocalelectron-electron interaction in graphene , Phys. Rev. B (Aug, 2016) 081403.[21] D. K. Efetov and P. Kim, Controlling electron-phononinteractions in graphene at ultrahigh carrier densities , Phys. Rev. Lett. (Dec, 2010) 256805.[22] R. C. Brower, C. Rebbi and D. Schaich,

Hybrid MonteCarlo simulation on the graphene hexagonal lattice , PoSLattice2011 (2011) 056, [ ].[23] P. V. Buividovich and M. I. Polikarpov,

Monte-Carlostudy of the electron transport properties of monolayergraphene within the tight-binding model , Phys. Rev.

B86 (2012) 245117, [ ].[24] M. V. Ulybyshev, P. V. Buividovich, M. I. Katsnelsonand M. I. Polikarpov,

Monte-Carlo study of thesemimetal-insulator phase transition in monolayergraphene with realistic inter-electron interactionpotential , Phys. Rev. Lett. (2013) 056801,[ ].[25] D. Smith and L. von Smekal,

Monte-Carlo simulation ofthe tight-binding model of graphene with partiallyscreened Coulomb interactions , Phys. Rev. B (2014) ].[26] D. Smith and L. von Smekal, Hybrid Monte-Carlosimulation of interacting tight-binding model ofgraphene , PoS Lattice2013 (2013) 048, [ ].[27] D. Smith, M. Koerner and L. von Smekal,

On thesemimetal-insulator transition and Lifshitz transition insimulations of mono-layer graphene , PoS Lattice2014 (2014) 055, [ ].[28] P. Buividovich, D. Smith, M. Ulybyshev and L. vonSmekal,

Interelectron interactions and the rkky potentialbetween h adatoms in graphene , Phys. Rev. B (Oct,2017) 165411, [ ].[29] P. Buividovich, D. Smith, M. Ulybyshev and L. vonSmekal, Competing order in the fermionic Hubbardmodel on the hexagonal graphene lattice , PoSLattice2016 (2016) 244, [ ].[30] T. Luu and T. A. L¨ahde,

Quantum monte carlocalculations for carbon nanotubes , Phys. Rev. B (Apr, 2016) 155106.[31] M. Koerner, D. Smith, P. Buividovich, M. Ulybyshevand L. von Smekal, Hybrid Monte Carlo study ofmonolayer graphene with partially screened Coulombinteractions at ﬁnite spin density , Phys. Rev. B (2017) 195408, [ ].[32] S. Beyl, F. Goth and F. F. Assaad, Revisiting theHybrid Quantum Monte Carlo Method for Hubbard andElectron-Phonon Models , Phys. Rev.

B97 (2018)085144, [ ].[33] P. Buividovich, D. Smith, M. Ulybyshev and L. vonSmekal,

Hybrid-Monte-Carlo study of competing orderin the extended fermionic Hubbard model on thehexagonal lattice , Phys. Rev. B (Dec, 2018) 235129,[ ].[34] P. Buividovich, D. Smith, M. Ulybyshev and L. vonSmekal, Numerical evidence of conformal phasetransition in graphene with long-range interactions , Phys. Rev.

B99 (2019) 205434, [ ].[35] J.-L. Wynen, E. Berkowitz, C. K¨orber, T. A. L¨ahde andT. Luu,

Avoiding Ergodicity Problems in LatticeDiscretizations of the Hubbard Model , Phys. Rev.

B100 (2019) 075141, [ ].[36] S. Krieg, T. Luu, J. Ostmeyer, P. Papaphilippou andC. Urbach,

Accelerating Hybrid Monte Carlo simulationsof the Hubbard model on the hexagonal lattice , Comput.Phys. Commun. (2019) 15–25, [ ].[37] D. Smith, P. Buividovich, M. Koerner, M. Ulybyshevand L. von Smekal,

Quantum phase transitions on thehexagonal lattice , .[38] A. Alexandru, M. Faber, I. Horvath and K.-F. Liu, Lattice QCD at ﬁnite density via a new canonicalapproach , Phys. Rev. D (2005) 114513,[ hep-lat/0507020 ].[39] P. de Forcrand and S. Kratochvila, Finite density QCDwith a canonical approach , Nucl. Phys. B Proc. Suppl. (2006) 62–67, [ hep-lat/0602024 ].[40] A. Nakamura, S. Oka and Y. Taniguchi,

QCD phasetransition at real chemical potential with canonicalapproach , JHEP (2016) 054, [ ].[41] W. Detmold, K. Orginos and Z. Shi, Lattice qcd atnonzero isospin chemical potential , Phys. Rev. D (Sep, 2012) 054507.[42] J. P. Hobson and W. A. Nierenberg, The statistics of atwo-dimensional, hexagonal net , Phys. Rev. (Feb,1953) 662–662. [43] H. Robbins and S. Monro, A stochastic approximationmethod , Ann. Math. Stat. (1951) 400.[44] N. Garron and K. Langfeld, Controlling the SignProblem in Finite Density Quantum Field Theory , Eur.Phys. J.

C77 (2017) 470, [ ].[45] M. V. Ulybyshev and S. N. Valgushev,

Path integralrepresentation for the Hubbard model with reducednumber of Lefschetz thimbles , .[46] R. Blankenbecler, D. J. Scalapino and R. L. Sugar, Monte carlo calculations of coupled boson-fermionsystems. i , Phys. Rev. D (Oct, 1981) 2278–2286.[47] R. T. Scalettar, D. J. Scalapino and R. L. Sugar, Newalgorithm for the numerical simulation of fermions , Phys. Rev. B (Dec, 1986) 7911–7917.[48] J. Ostmeyer, E. Berkowitz, S. Krieg, T. A. Laehde,T. Luu and C. Urbach, The Semimetal-Mott InsulatorQuantum Phase Transition of the Hubbard Model on theHoneycomb Lattice , .[49] M. Ulybyshev, C. Winterowd and S. Zafeiropoulos, Lefschetz thimbles decomposition for the Hubbard modelon the hexagonal lattice , .[50] M. Ulybyshev, C. Winterowd and S. Zafeiropoulos, Taming the sign problem of the ﬁnite density Hubbardmodel via Lefschetz thimbles ,1906.02726