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

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

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

I. INTRODUCTION
Monte Carlo simulations based on path-integral quantization are a powerful and widely used tool for the study of strongly coupled quantum systems. They are applied in many different areas of physics, ranging from highenergy physics, where they are employed e.g. to study the phase diagram and particle spectrum of Quantum Chromodynamics (QCD), to condensed matter physics, where they are used to study strongly correlated electron systems. Quite often, they are the only way to obtain reliable information from first principles. Unfortunately, their applicability is restricted to a very special class of systems, namely those where the path integral exhibits a positive-definite measure which can be interpreted as a probability density. This excludes most fermionic systems away from half filling (unless the complex parts of the measure cancel exactly due to an anti-unitary symmetry) as well as quantum systems which evolve in real (as opposed to Euclidean) time. This restriction is known as the sign problem and is a long-standing problem of numerical physics.
In principle, brute-force reweighting techniques can be used to extract information about systems with charge imbalance from simulations of a corresponding system at charge neutrality. These are plagued by a severe signal-to-noise ratio problem however, originating from a loss of overlap between the ensembles at zero and nonzero charge density when the thermodynamic limit is approached, and typically fail well below µ/T ≈ 1 except on very small systems. While a theorem by Troyer and Wiese states that the sign problem is a NP-hard problem in a generic spin-glass system [1] (making the discovery of a complete universal solution to all sign problems unlikely), many attempts have been made to construct specialized solutions for specific systems (typically by introducing a set of dual variables), or improved general approaches which outperform reweighting, with some success. Most notably, in QCD, Taylor expansions of the partition function with respect to µ/T have now been pushed to µ/T ≈ 2 or 3 [2].
A promising, quite different, idea to deal both with ergodicity problems in Monte Carlo simulations of sys-tems close to a first order phase transition, and overlap problems resulting from a non-positive probability measure, is to use non-Markovian random walk simulations, which do not rely on importance sampling with respect to a positive Gibbs factor. A particularly interesting class of algorithms use the inverse density-of-states as a measure for updating configurations. This measure is positive (semi-)definite by definition, and produces a random walk which efficiently samples configuration space even in "deprived" regions with very low probabilistic measure. Prominent examples are the multicanonical algorithm by Berg and Neuhaus [3] and the Wang-Landau approach [4], which both were designed for theories with discrete degrees of freedom.
The Linear Logarithmic Relaxation ("LLR") method, first described in Ref. [5], generalizes this idea to quantum systems with continuous degrees of freedom. The goal of LLR is to estimate the slope a(X) = d dX ln(ρ(X)), where X is some observable and ρ(X) is the "generalized density of states" (gDOS). Once a(X) is obtained, ρ(X) can be reconstructed up to a multiplicative factor by numerical integration and can be used to compute thermodynamic observables. A crucial property of LLR is that it features exponential error suppression: The estimate for a(X), and by extension of ln(ρ(X)), can be obtained with roughly the same statistical error, independent of the exact value of X, even if X is from a region of overall low weight.
In recent years, LLR was successfully applied to obtain ρ(E) in SU(2) and U(1) [5] as well as SU(3) [6] gauge theories and was shown to be effective at dealing with ergodicity issues arising at first-order phase transitions in U(1) gauge theory [7] and the q = 20 state Potts model [8]. In Ref. [9] LLR was applied to obtain the Polyakov loop distribution for two-color QCD (which has no sign problem) with heavy quarks at finite densities. To deal with the sign problem, one needs to compute the gDOS for the imaginary part of the Euclidean action ρ(S I ), or some related observable. This was achieved using LLR for a Z 3 spin model at finite charge density [10] and for QCD in the heavy-dense limit [11]. To date, LLR has never been applied to a sign problem of a system with fully dynamical fermions however.
In this work, we apply LLR to the fermionic Hubbard model on the honeycomb lattice away from half filling within a Hybrid Monte Carlo framework. Despite its simplicity, the Hubbard model, which describes fermionic quasi particles with contact interactions, continues to be of profound interest, as it remains the quintessential example of an interacting fermion system, and can qualitatively describe many non-perturbative phenomena such as dynamical mass-gap generation or superconductivity. On the honeycomb lattice, this model exhibits a second order phase transition in the Gross-Neveu universality class [12][13][14][15] and extended versions, which include longrange interactions, realistically describe the physics of both mono-and bilayer graphene [16,17]. Using LLR, we compute the gDOS for the average of a real-valued auxiliary field, which is introduced in Hybrid Monte Carlo simulations to transmit inter-electron interactions. We demonstrate that this result can be used to reconstruct the particle density as a function of chemical potential. We show that, in its present form, LLR can probe much further into the finite density regime than standard reweighting, that the relative advantage of LLR grows as the interaction strength is increased, and argue that future improvements might put the van Hove singularity in the single-particle bands within reach. This paper is structured as follows: We start in Sec. II by introducing the basic lattice setup and illustrating the sign problem away from half filling. Subsequently, we introduce the generalized density of states of the average Hubbard field ρ(s) in Sec. III. In Sec. IV, we discuss the reconstruction of the particle density n from ρ(s). We describe two different reconstruction schemes, whereby n(µ) is obtained from both the canonical and grand-canonical partition functions. As a benchmark, we apply both schemes to test data obtained for the noninteracting tight-binding theory. Full LLR calculations of the interacting theory, including additional numerical details, are then presented in Sec. V. We summarize and conclude in Sec. VI.

II. LATTICE SETUP AND THE SIGN PROBLEM
We consider the repulsive Hubbard model on the hexagonal (honeycomb) lattice with fermionic creation operatorsĉ † x ≡ (ĉ † x,↑ ,ĉ † x,↓ ) for two spin components at site x, which is defined by the effective Hamiltonian for the grand canonical ensemble: Here κ is the hopping parameter, the sum over x, y sums all independent pairs of nearest-neighbor sites, m s is the sublattice-staggered mass term (with alternating sign on the two triangular sublattices) for explicit sublattice-symmetry breaking with spin-density-wave order, µ is the charge chemical potential andρ x =ĉ † xĉx − 1 the charge operator. The constant U controls the interaction strength, which is positive in the repulsive Hubbard model. The creation and annihilation operators satisfy the fermionic anticommutation relations {ĉ x ,ĉ † y } = δ x,y 1. Lattice simulations of (1) using Hybrid Monte-Carlo by now have a long history already [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33], we thus summarize only the essential steps here. 1 To derive the functional integral representation of the partition function at inverse temperature β = 1/T , we first write the exponential in terms of N t identical factors and split the Hamiltonian into the free tight-binding part plus interactions,Ĥ =Ĥ 0 +Ĥ int . A symmetric Suzuki-Trotter decomposition of each of the factors then yields This introduces a finite step size δ = β/N t in Euclidean time and a discretization error of O(δ 2 ). As we will see shortly, it is convenient to include the chemical-potential term in the definition ofĤ int here, The four-fermion interaction inĤ int is then converted to bilinears by Hubbard-Stratonovich transformation, whereby the auxiliary ("Hubbard-Coulomb") field φ x,t is introduced. Finally, the trace over the fermionic operators is performed by integrating the fermionic coherent states [29], which yields Different versions of the fermion matrix M (φ) have been used in the past which are either equivalent or at least yield the same continuum limit. In this work we use h xx = δ xx m s − κ n δ x ,x+ n . 1 In particular we omit the discussion of fermionic coherent states and the partial particle-hole transformation that is applied. These and other details can be found e.g. in Refs. [21,29,33]. for Ns = Nt = 6, β = 2.7κ −1 , U = κ/2 at different µ. The results are modelled with Gaussian (for µ = 0.0148κ and 0.0667κ) and uniform (for µ = 0.2κ) distributions respectively. The inlay shows the adjusted R 2 for a constant fit to data at different µ. For µ 0.15κ the numerical data is well described by a uniform distribution, indicating a hard sign problem. An analogous figure as indication of a sign problem is obtained for graphene with long-range interactions [27].
in which the free tight-binding hopping contributions of the form e −δh are linearized, in order to be able to work with sparse matrices, but the diagonal couplings to Hubbard field and chemical potential are not.
It is clear that Eq. (5) is sign-problem free at half filling, i.e. for µ = 0, as det(M M † ) = | det M | 2 . This is no longer true for µ = 0. By writing we can consider the complex ratio of determinants with unlike-sign chemical potentials as an observable in the "phase-quenched" theory (defined by the modulus of the fermion determinant) with partition function Z pq and obtain This ratio is unity for µ → 0 and is a measure of the severity of the sign problem, as it quantifies the effectiveness of brute-force reweighting. Fig. 1 shows histograms of the phase 2 of the ratio det M (φ, µ)/ det M (φ, −µ) for different non-zero values 2 The modulus also deviates from unity at µ = 0 but need not be considered here.
of µ, obtained on a lattice of N s × N s unit cells, with 2 sites per unit cell, and N s = N t = 6, at β = 2.7κ −1 and U = κ/2, together with fit-model curves. The adjusted R 2 of a constant fit (corresponding to a uniform distribution) shows a rather rapid crossover and approaches values close to 1 at µ ≈ 0.15κ. This indicates that the signal is lost in the noise rather quickly already on small lattices (signalling a hard sign problem), and for rather high temperatures and weak interaction strengths. This effect is enhanced with larger lattice sizes, lower temperatures and larger couplings. To present a quantitative comparison of brute-force reweighting and the LLR method for different sysem sizes and interaction strengths is one of the main objectives of this work.

III. GENERALIZED DENSITY OF THE HUBBARD FIELD
Assume we have a quantum system with action βS(φ). Defining the density of states ρ(E) as we can then express the partition function as and the vacuum expectation value of an observable O(E) becomes If ρ(E) is known, then O can be obtained by (numerically or analytically) integrating Eq. (11). Here we have assumed that we know how to express O in terms of E, which is in general not the case however. Moreover, Eq. (9) is ill-defined if S(φ) is not strictly real. To compute different observables in a generic setting, the concept of the density of states can thus be generalized to quantities other than the action. The basic idea of LLR is to obtain ρ(X) (where X is some observable) by carrying out a sequence of "microcanonical" Monte-Carlo simulations, in which X is forced to assume a set of different (sufficiently dense) values X i . Obtaining the partition function or thermodynamic expectation values then essentially amounts to computing a Fourier or Laplace transform of ρ(X). To alleviate the sign problem, ρ(X) must reflect the phase fluctuations of the path-integral measure. To this end, we consider ρ(s = Φ) in this work, where Φ is the spacetime-volume average of the auxiliary field, see below. First, we apply the transformation to Eq. (5), which then leads to In this formulation, the complex part of the action has been shifted completely to the bosonic sector. Eq. (13) is now rewritten as where we have introduced the average Hubbard field where N c denotes the number of unit cells with 2 sites each. Finally, we introduce and rewrite the partition function as where ρ(s) is the generalized density of states of the average Hubbard field Φ and will be the target of our LLR calculation.

IV. RECONSTRUCTING THE PARTICLE DENSITY
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 fact that LLR obtains ρ(s) only for a discrete and finite set of points and with finite numerical precision. Our ultimate goal is to obtain the particle density n(µ). We present two reconstruction schemes which achieve this in the following, which both operate in the in the frequency domain and avoid the instabilities of Eq. (17).
We note that ρ(s) has a periodicity of 2π/β = 2πT and can thus be expanded in Fourier series. For later convenience we will first introduce a dimensionless variable and density via If we furthermore introduce an imaginary chemical potential via we observe that up to a Gaussian smearing with variance U/2N c T the generalized density of states is in fact essentially the same as the partition function at imaginary chemical potential, We will obtainρ(x) only at a discrete set of points x n = 2πn/K, where n = {0, . . . , K − 1} andρ n ≡ρ(x n ). We must hence truncate the Fourier series, naturally leading to a discrete Fourier transform which can be used for interpolation viã On the other hand, inserting (21) into (17), we obtain (22) and the exact result is recovered for K → ∞. In fact, in this limit, Eq. (22) becomes the fugacity expansion and we can identify for k = N , as the corresponding canonical partition function with particle number N . In the infinite volume limit N c → ∞ for fixed N , or equally so for T U , we may therefore neglect the exponential factor and essentially identify the Fourier series coefficientsρ k of our generalized DOS with the canonical partition functions at N = k. At the same time, it is also evident from Eq. (22) that the generalized DOS 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 verifies that the truncated coefficientsρ k at finite K, obtained from the discrete Fourier transform in (21), then yield pseudo-canonical partition functions,ρ k ∝ Z K c (T, k), which represent ensembles with particle number N = k mod K. Likewise, the discrete sampling ofρ(θ) provides us with an interpolation of Z I (θ) which agrees with the exact result for imaginary chemical potential at the discrete values θ n = 2πn/K. The general relation between ρ(s) and the partition function at imaginary chemical potential of course also follows from Eq. (22), with µ = is (and K → ∞), In a finite volume, on the other hand, i.e. at any finite number N c of unit cells, the particle numbers N are restricted to values between ±2N c , with N = 0 at half filling, corresponding to an average of one of the maximally possible two electrons on each of the 2N c sites. We then obtain the exact canonical partition functions from Eq. (23) already for and with particle-hole symmetry at half filling, one actually only needs K max = 2N c + 1.
In principle, the particel number N (µ) can be directly obtained from Eq. (22), which is free of oscillating terms, by taking the deritvative with respect to µ, Computing theρ k from ρ(s n ) can be done with high numerical precision using modern FFT libraries.
Alternatively, we can also compute the chemical potential from the canonical partition functions, as the free energy difference of ensembles with subsequent particle numbers. From Eq. (22) we then obtain and obtain the density in form of the number of particles per unit cell, n(µ) ≡ N (µ)/N c , by inversion. The exact calculation would again require K max = 2N c + 1 Fourier coefficients. This is then similar in spirit to Refs. [34,35], which carried out canonical calculations of QCD at finite charge density, or Ref. [36] which followed essentially the same strategy for finite isospin density from the lowest states in multi-pion correlators. With truncating at K < K max , we strictly speaking obtain canonical ensembles at particle number N modulo K as discussed above. The term ∝ U/2N c T in Eq. (26) represents an explicit finite volume effect which, as we will discuss below, only contributes in trivial way and can be dropped. In tight-binding or mean-field calculations, there is no such term in the first place, and the generalized DOS can be calculated analytically. The result is of the form where ε ≥ 0 is the single-particle energy with spectral density ρ ε (ε) for which an analytic expression is known in the infinite system [37]. In a finite system with periodic (Born-von Kármán) boundary conditions we use the dispersion relation ε = ε(k) instead, and simply sum over the corresponding discrete set of points k n within the first Brillouin zone, with energies ε n = ε(k n ). The same can be done to compute the exact density in the finite system with N c unit cells which then yields for the number of particles per unit cell, We have carried out a set of benchmark calculations in which we compared the canonical and grand-canonical reconstruction schemes. Thereby, a discete set of values l n = ln ρ(s n ) for s n = 2πT n/K, N = {0, . . . K − 1}, was produced as mock data from the tight-binding calculation, which can efficiently be done with arbitrary numerical precision. High-precision calculations are especially important in the reconstrucion of the density because we need with high precision the discrete Fourier transform of {ρ n = e ln } 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 the canonical approach (26). We have then compared both results with the exact calculation of the density based on the tight-binding formula (28). This was done for different setups, whereby the production of {l n = ln ρ(s n )} was done with different levels of floating point precision. The application of the reconstruction scheme was done with a 1024 digit accuracy in each case to avoid additional errors. We find that both methods yield comparable results, with the canonical procedure having a very slight advantage for a given precision of ln ρ(s). We thus choose to use this procedure exclusively in the following sections to process our LLR results. Fig. 2 shows an example calculation of n(µ) for two different temperatures on a lattice with N c = 36 unit cells, where ln ρ(s) was produced for U = 0 with 1024 digit accuracy and processed using (26). This illustrates that our method can in principle cover the entire width of the valence band, from the empty valence band at half fillg up to saturation when it is completely filled. The van Hove singularity will emerge at µ = κ in the infinite volume limit which can here be anticipated already by the rapid increase in the number density at the lower temperature around µ = κ.
In practice, the leading source of errors is of course the precision with which the Fourier coefficientsρ k can be obtained, which in turn is highly sensitive to statistical errors of ln ρ(s) in our LLR calculations. In order to get to saturation, with a completely filled lattice, we would obviously need the maximum number K max = 2N c + 1 of coefficients. So the double challenge here will be to compute as many of them as accurately as possible.

A. The algorithm
The goal of LLR is to calculate derivatives of ln ρ(s) at a sufficiently dense set of supporting points with high precision and to reconstruct ρ(s) by integration. We divide the domain of support of ρ(s) into K intervals of size δ s . At the center of each of these intervals, the slope a k = d ds ln ρ(s)| s=s k can be calculated from a stochastic non-linear equation [5]. A key element of this equation is the restricted and reweighted expectation value 3 Here Z LLR is a normalization constant, Φ was introduced in Eq. (15), a is an external parameter and θ [s k ,δs] is a window function which restricts Φ to an interval of size δ s around s k .
With the choice W (Φ) = Φ − s k , the coefficients a k are solutions of This equation can be solved through Robbins-Monro iteration [38]: The sequence converges to the correct result for any choice of α n that fulfills This is true, even if W (Φ)(·) k is approximated by an estimator, as we do in Monte-Carlo calculations. Moreover, if the iteration is terminated at some finite number N and repeated many times, the final values a (N ) k are Gaussian distributed around the true value a k and can be processed by a standard bootstrap analysis.
The window function can be chosen in different ways. The straight-forward choice is a step function, but for HMC a Gaussian window function is more appropriate, as its derivative can be taken, which implies that its effect can be reproduced by a molecular-dynamics force term. In this work, we choose where The full procedure is then summarized as follows: 1) For a given s k , initialize a k with some random value a  In practice, we start with several repetitions of steps 3−6 with fixed α = 1 and switch to underrelaxed interations with α n+1 = α n /(1 + n) after some time. Also, the whole procedure is terminated after some finite iteration number N and repeated several times for each fixed s k , to produce a sample of final 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 fixed set of external parameters. We choose N s = N t = 6, β = 2.7 κ −1 , m s = 0.185 κ, U = 1.0 κ, s = 1.33 κ for illustration. 4 For each set of parameters considered in this work, we first obtain such a sample of a k values. We then obtain ln ρ(s k ), and by extension ρ(s k ),ρ k and n(µ) together with errorbars, by feeding bootstrap averages of the final a computing the Fourier transform of ρ(s k ) and applying the canonical reconstruction scheme described in Sec. IV.

B. Nt dependence
We begin by studying the effect of the timediscretization δ. To this end, we carry out LLR calculations at N s = 6, β = 2.7 κ −1 , U = 1.0 κ, m s = 0.185 κ for different values of N t . Fig. 4 shows the results for 4 Note that the phenomenological value of the hopping parameter in the tight-binding model for graphene typically is κ ≈ 2.7 eV, so this would correspond to a temperature of T ≈ 1 eV in graphene.  Our first observation is that the dependence on N t is very mild for our choice of parameters. It is practically invisible in a(s) and n(µ). A very small difference between different N t can be seen in ln ρ(s k ) and lnρ k , which is of a similar magnitude as the statistical uncertainty however. On the other hand, our results clearly demonstrate exponential error suppression, whereby the relative error of ln ρ(s) is roughly the same across several orders of magnitude. We find that lnρ k is extremely sensitive to this small error however, to a degree that only the first few Fourier modes lnρ k can be computed accurately. This can be traced back to the fact that ρ(s) enters into the Fourier transform and not ln ρ(s). It is also reflected in our computation of n(µ), which exhibits a loss of signal at µ ≈ 0.5 κ, indicating the onset of a hard sign problem.
We note that for U = 1.0 κ which is well inside the weak-coupling phase of the model, and the temperature considered here, n(µ) basically fully agrees with the infinite-volume limit in the non-interacting theory when the linear term in Eq. (26) is dropped. We take this as an indication that this extra term represents the dominant finite-volume effect at finite U which however is a rather trivial one to correct. Further confirmation of this is provided by a comparison between results from N s = 6 and N s = 12 lattices, which also reveals a faster convergence to the thermodynamic limit without this term. We thus drop this term for all results presented in the following. We expect that deviations from the non-interacting limit will become visible at stronger couplings, of course. This is investigated further, and ultimately confirmed, below.

C. ms dependence
We now turn to studying the dependence on the explicit sublattice and spin-staggered mass term m s . Given that such a term already opens an explicit gap in the energy spectrum, we carry out this study at the comparatively weak coupling strength of U = 0.1 κ. We find that, again, the number density n(µ) coincides with the noninteracting theory and shows no significant dependence on m s . The linear term in Eq. (26) has a negligible effect here, due to the small value of U . Fig. 6 shows the results for a(s) and ln ρ(s) with N s = N t = 6, β = 2.7 κ −1 and three different choices of m s . We refrain from showing any additional figures for lnρ k and n(µ), as these fully agree (within statistical errors) with the results shown in the lowest panels of Figs. 4 and 5.
An interesting observation here is that m s has a quite strong effect on both a(s) and ln ρ(s), which turns out not to carry over to n(µ) at all. The underlying reason is that this dependence is only present in regions where ρ(s) is strongly suppressed. It is only visible due to the logarithmic scale, and thus has no significant effect on the computation of the Fourier modes.

D. U dependence
Having validated our numerical procedure at weak coupling, we now turn to a more detailed study of the dependence on the interaction strength U . This represents the central part of this work to which the bulk of our computing resources were dedicated. We thereby computed a(s), ln ρ(s), lnρ k and n(µ) again with β = 2.7 κ −1 , m s = 0.185 κ for several different choices of U . To have control over finite volume and time-discretization effects we have studied two different lattice sizes, N s = N t = 6 and 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.7 κ −1 , m s = 0.185 κ. For ln ρ(s) we include the tight-binding result to illustrate the approach to the non-interacting limit. The first observation is that a(s) gets suppressed when U is increased, which ultimately makes simulations more expensive at strong coupling. On the other hand, we clearly see a deviation from the non-interacing limit in the Fourier modes lnρ k for the strongest interaction strength U = 2.0 κ. To underscore that this deviation is absent for all weaker interactions, we show a seperate plot in Fig. 8 which directly compares lnρ k for U ≤ 1.0 κ to the tight-binding theory. Our N s = N t = 6 results for n(µ) are shown in Fig. 9. They clearly show a corresponding drop of the number density at fixed µ for the strongest coupling. N s = 12 results are shown in Fig. 10 for a(s), ln ρ(s) and lnρ k and Fig. 11 for n(µ). These confirm the qualitative changes at U = 2.0 κ. Furthermore, a direct comparison with N s = 6 suggests that finite volume effects on n(µ) are rather mild.
We point out here that the sign problem sets in at much smaller µ for the larger system (as expected). While we are able to reliably compute n(µ) up to µ ≈ 0.35 κ for N s = 6 with U = 1.0 κ, we only reach µ ≈ 0.1 κ for N s = 12. On the other hand, in both cases LLR drastically outperforms brute-force reweighting: With comparable numerical resources we obtain a signal for the determinant ratio (8) up to µ ≈ 0.14κ on N s = 6 and µ ≈ 0.075κ on N s = 12 using the brute-force method. While the relative advantage of LLR becomes smaller on the larger lattice, we can reach much larger values of µ for U = 2.0 κ (µ ≈ 0.5 κ on N s = 6 and µ ≈ 0.2 κ on N s = 12). In contrast, the µ range of reweighting is drastically diminished at stronger coupling (cf. Fig. 13 in Sec. VI). It is this last feature which ultimately makes LLR in its present form a promising method and deserving of further attention.

E. Compressed sensing
Lastly, we report on our attempts to improve our results by using fit functions for ln ρ(s), a procedure re- ferred to as compressed sensing in the LLR literature. The basic idea is to, instead of processing the raw data for ln ρ(s) pointwise at the supporting points s k , fit the entire data set with a series expansion in some complete set of functions, and use the model curve to compute observables instead. The hope is that an appropriate set of functions, which reflects the true (but a priori unknown) physics of the theory, will both suppress noise in the numerical data for ln ρ(s) and effectively generate an interpolation to a much denser set of supporting points. This in turn should allow for the computation of lnρ k at larger k and hence the number density at larger µ. Fig. 12 shows two such attempts, where ln ρ(s) was fit with a Fourier series and a series of Chebyshev polynomials of the first kind respectively. The fit function was subsequently evaluated at a much denser set of points than the original s k and used to compute lnρ k and subsequently n(µ). In each case, higher order terms were added to the expansion until the final result stabilized. We do not obtain any errorbars. Results from the direct calculation are included for comparison and represented by dashed lines.
We find that these attempts do not improve the calculation of n(µ) significantly. At best, one or two additional points (at higher densities) can be computed before the results scatter in an uncontrolled fashion. The Fourier series thereby seems to work only slightly better than the Chebyshev polynomials. We take this as an indication that additional qualitative information about the system, which places additional contraints on the choice of functions to use, is a necessary requirement for compressed sensing to be effective here.

VI. CONCLUSION AND OUTLOOK
In this work, we have applied the Linear Logarithmic Relaxation method to the repulsive fermionic Hubbard model on the honeycomb lattice, in order to assess its utility for alleviating the hard sign problem of an unbalanced dynamical fermion system. A central problem thereby is the proper choice of a target observable, which adequately reflects the complex part of the action and yields a generalized density of states which is suitable for further processing. We used the average value Φ of the auxiliary (Hubbard) field to this end, which appeared as the natural choice, as it allows for the shifting of the complex part of fermion determinant into the bosonic sector and provides a simple integral expression (Eq. (17)) for obtaining the partition function and hence the particle density. To deal with an oscillating contribution to this integral, we chose to work in the frequency domain and divised two methods to extract the particle density from the Fourier modes of the gDOS of Φ which essentially yields the partition function at imaginary chemical potential. Due to a slightly better performance in benchmark calculations, of these we chose a method based on the canonical ensembles to further process our LLR results.
We have carried out LLR calculations for a fixed temperature of β = 2.7 κ −1 , two different lattice sizes (6 3 and 12 3 ) and different interaction strengths in the weak and intermediate coupling regime, and obtained the particle density as a function of chemical potential. We thereby observed significant deviations from the non-interacting theory for the largest interaction strength considered, possibly signalling the onset of spontaneous mass-gap for- mation. We found that using LLR in its present form, on the smaller 6 3 lattice we are able to probe at least twice as far into the finite-density regime as with bruteforce reweighting. While the relative advantage of LLR is smaller on the 12 3 lattice, we find that LLR performs much better when the interaction strength is increased. Fig. 13 shows a quantitative comparison of the effective µ-ranges for the different parameter sets considered in this work. Attempts to reach into higher-density regions were made using different forms of compressed sensing, i.e. by fitting ln ρ(s) with Fourier series and Chebyshev polynomials and using the model curves for interpolation. While this allows us to reach slightly higher densities, we suspect that this procedure introduces an uncontrolled systematic error, as the physics at higher densities is strongly sensitive to the high-frequency modes of ρ(s), which such interpolations cannot account for.
The results presented here should be taken as a proof of principle. There are several different directions for future improvements. First and foremost, the computational resources spent for the final calculations (i.e. the sum of all results shown in Figs. 9 and 11) were around two months of runtime on a total of 18 GTX 980 Ti GPUs, which leaves much space for larger-scale projects. We estimate that, using the most modern hardware and libraries for sparse linear algebra, the precision for ln ρ(s) can be increased by at least an order of magnitude. More advanced techniques for compressed sensing could also be applied, such as Gaussian and telegraphic approximations or an advanced moments approach, which were proposed in Ref. [39]. Quite possibly, introducing a complex instead of a real auxiliary field has advantages, and in fact, it was shown that an optimal mixing factor between real and imaginary Hubbard fields exists, for which the sign problem is the mildest [40]. LLR might also be more effective with a discrete Hubbard field, which is used in BSS Quantum Monte-Carlo calculations [41,42]. In addition, an alternative time-discretization with a symmetry of time reversal times sublattice exchange, which was proposed already in Ref. [18] and recently used in a grand canonical HMC simulation [43], was shown to have strongly suppressed discretization effects in Ref. [26] and might positively impact the performance of LLR as well. And finally, there has been much recent progress regarding the Lefschetz thimble method [40,44,45], and constructing a hybrid approach, which combines the advantages of both methods, might be feasible. Specifically, one could attempt to apply the Lefschetz thimble decomposition directly to Eq. (17), in order to avoid the use of reconstruction schemes altogether and obtain a cleaner signal for n(µ).
Taken together, we find it not unreasonable to expect that future developments might put the van Hove singularity (VHS) of the single-particle spectrum within reach, which is of great interest in the context of superconducting phases. A crucial point thereby is the apparent stability of the LLR technique against increases of the coupling strength U . Experiments on charge-doped graphene systems have revealed a strong bandwidth renormalization (narrowing of the width of the π-bands) due to interactions [46], which suggests that the VHS can be probed at smaller µ for larger U . Furthermore, a HMC study of graphene at finite spin density revealed that the electronic Lifshitz transition at the VHS can become a true thermodynamic phase transition in the presence of interactions, with a critical temperature which increases with the coupling strength [27]. A study of an analogous transition at finite charge-carrier density thus might be feasible at large U , in particular as the sign problem becomes milder at higher temperatures.
There are many possibilities to linearize the quartic fermionic interaction using auxiliary fields. The choice in this paper is inspired by the observation that an explicit analytic continuation (i.e. Eq. (12)) was sufficient to split off the complex part of the fermion determinant. The formulation is elegant: The calculation of one (real) density of states, ρ(s), is sufficient to relay the calculation of the partiton function to one integration for each given value of the chemical potential. Note, however, that the use of a Hubbard field φ with a compact domain of support implies that the domain of the density of states is also compact. The calculation of such an "intensive" density of states to sufficient precision is difficult [39] and the most successful LLR calculations for theories with a sign problem are based upon non-compact densities [10]. The use of a non-compact formulation is left to future work.
Lastly, we should mention that extending our work to the QCD sign problem remains an open conceptual challenge. The system considered here was special since we succeeded to remove the complex part the fermion determinant by a simple analytic continuation. For gauge theories no such simple transformation exists, and measuring a proper extensive phase is much more involved. It may well be that this step is the most computationally demanding and contains the central computational complexity of the QCD sign problem.