The QCD sign problem as a total derivative
aa r X i v : . [ h e p - l a t ] S e p The QCD sign problem as a total derivative
Jeff Greensite , Joyce C. Myers and K. Splittorff Physics and Astronomy Dept., San Francisco State University, San Francisco, CA 94132, USA Discovery Center, The Niels Bohr Institute, University of Copenhagen,Blegdamsvej 17, DK-2100, Copenhagen Ø, Denmark (Dated: September 27, 2018)We consider the distribution of the complex phase of the fermion determinant in QCD at nonzerochemical potential and examine the physical conditions under which the distribution takes a Gaus-sian form. We then calculate the baryon number as a function of the complex phase of the fermiondeterminant and show that the exponential cancellations produced by the sign problem take theform of total derivatives that the full baryon number is orthogonal to this noise. These insightsallow us to define a self-consistency requirement for measurements of the baryon number in latticesimulations. Introduction:
Perhaps the simplest way to identifythe asymmetry between matter and anti-matter in QCDis to measure the free energy of a static quark and com-pare it to that of an anti-quark. These two free energiesare linked to the expectation values of an observable, thePolyakov loop, and its complex conjugate. In a baryonasymmetric world the chemical potential, µ , favors prop-agation of quarks and as a result the expectation valueof the Polyakov loop will be different from that of itscomplex conjugate. If the action was real the two ex-pectation values would be equal. Hence the asymmetrybetween matter and anti-matter can only occur if theQCD action takes complex values. This fact is known as the QCD sign problem , see [1, 2] for excellent reviews.As the physical origin suggests, sign problems likewiseoccur naturally for systems with strongly correlated elec-trons [3]. In both cases the complex nature of the actioninvalidates the Vafa-Witten theorem [4]. Vector symme-tries such as baryon number and isospin are therefore notprotected. In this way the sign problem allows for a muchricher phase structure.The QCD sign problem enters through the fermion de-terminant, det( / D + µγ + m ) = re iθ . (1)For µ = 0, the determinant is real due to the γ Her-miticity of the Dirac operator (cf. [5]). The chemicalpotential term breaks γ Hermiticity and for this reasonthe determinant is in general complex (apart from certaingauge groups such as SU (2) and G with real or pseudo-real representations). While the physical nature and theexplicit origin of the sign problem is clear we are onlybeginning to understand the consequences. The reasonis that vacuum expectation values in such non statisticalsystems have entirely new ways to manifest themselves.The chiral condensate, for example, results from extremecancellations [6, 7].The cancellations due to the sign problem are expo-nential in the volume and thus the use of lattice QCDsimulations is severely limited. Despite this, remarkablenumerical results obtained by working at fixed phase of the fermion determinant, and subsequently integratingover the distribution of the phase have been presented[8–10]. This approach goes by several names: the factor-ization method [8], the density of states method [9] andthe histogram method [10]. Here we derive the distribu-tion of the complex phase of the fermion determinant atnonzero chemical potential, which is at the core of thisapproach. In particular, we examine the conditions underwhich the distribution takes the Gaussian form assumedin [10], see [11] for a review.The results presented here follow directly from a gen-eral assumption on the analytic structure of the free en-ergies involved. This assumption is identical in nature tothat of the replica method, see e.g. [12], and hence weexpect that the results applies equally generally.The central challenge to lattice QCD in this contextis to compute the asymmetry between matter and anti-matter, that is, the average baryon number h n B i N f with N f dynamical flavors. In the approach discussed hereone first computes the distribution h n B δ ( θ − θ ′ ) i N f ofthe baryon number as a function of the complex phaseof the fermion determinant. (Notation: θ is the value ofthe phase which we keep fixed and θ ′ denotes the valueof the phase of the determinant for the gauge field con-figurations in the average.) Then the baryon number isthe integral of this distribution over θ , h n B i N f = Z d θ h n B δ ( θ − θ ′ ) i N f . (2)In this paper we determine the analytic form of the distri-bution h n B δ ( θ − θ ′ ) i N f of the baryon number as a func-tion of the complex phase of the fermion determinant.The sign problem manifest itself in the complex natureof the distribution h n B δ ( θ − θ ′ ) i N f . In particular, as wewill demonstrate below, terms which are sub-leading inthe volume at fixed θ contribute at leading order afterintegration over θ . This makes the θ -integral hard tocontrol numerically. The analytic form of the distribu-tion of the baryon number as a function of the complexphase, however, reveals two important properties: first ,the cancellations produced by the sign problem takes theform of total derivatives wrt. θ , second , the full baryonnumber (the signal) is orthogonal to this noise. Basedon this we formulate a self-consistency requirement formeasurements of the baryon number.The general assumption and results obtained are ex-emplified within the hadron resonance gas model, andin the combined strong coupling and hopping parameterexpansion. Distribution of the phase:
We first consider thedistribution, h δ ( θ − θ ′ ) i N f , of the complex phase. Usingthe basic definitions we have [13] h δ ( θ − θ ′ ) i N f = 1 Z N f Z D A δ ( θ − θ ′ )det N f ( / D + γ µ + m ) e − S Y M = 1 Z N f Z D A δ ( θ − θ ′ ) | det( / D + γ µ + m ) | N f e N f iθ ′ ( A ) × e − S Y M = e N f iθ Z | N f | Z N f h δ ( θ − θ ′ ) i | N f | . (3)The subscript | N f | refers to the phase quenched theorywhere the fermion determinant is replaced by its absolutevalue. The phase quenched distribution, h δ ( θ − θ ′ ) i | N f | ,is of course real and positive. The distribution in the fulltheory, however, has an explicit phase factor exp( N f iθ ).This explicit phase must give rise to exponential cancel-lations in order to fulfill the normalization Z d θ h δ ( θ − θ ′ ) i N f = 1 (4)because Z | N f | /Z N f is exponentially large in the volume.To compute the distribution of the complex phase an-alytically we express it through a Fourier transform [14] h δ ( θ − θ ′ ) i N f = 12 π ∞ X p = −∞ e − ipθ h e ipθ ′ i N f , (5)where the moments are h e ipθ ′ i N f = 1 Z N f (cid:28) det N f + p/ ( / D + γ µ + m )det p/ ( / D − γ µ + m ) (cid:29) , (6)and the bracket without a subscript is the averagew.r.t. the Yang-Mills action. When p is even, the mo-ments are simply partition functions with p/ p/ p isodd.It is natural to assume that the associated freeenergies have a series expansion in p , ie. that log h e ipθ ′ i is a polynomial in p . The form of this polynomial is constrained by the following requirements: Since partition functions are even in µ the moments,Eq. (6), are invariant under p/ → − p/ − N f If we take p/ → p/ − N f / h e i ( p − N f ) θ ′ i N f = Z | N f | Z N f h e ipθ ′ i | N f | . (7)To ensure these requirements the unquenched momentstake the form h e ipθ ′ i N f (8)= e − p/ p/ N f ) X − ( p/ p/ N f )) X − ( p/ p/ N f )) X − ... . The values of X j determine the distribution of the com-plex phase through Eq. (5). A few properties of the X j can be deduced from first principles: First of all, sincethe moments are partition functions the coefficients X j are extensive. Moreover, as noted above we have h e ipθ ′ i | N f | (9)= Z N f Z | N f | e − ( p − N f ) X / − ( p − N f ) X / − ( p − N f ) X / ... , and in particular for p = 01 = Z N f Z | N f | e ( N f / X − ( N f / X +( N f / X ... . (10)The sum in the exponent is therefore the free energy dif-ference between the full and the phase quenched theoryand hence must be positive since Z N f ≤ Z | N f | . Gaussian distribution:
In [10] it was found numeri-cally that the distribution of the phase to good approxi-mation assumes a Gaussian form. To obtain the Gaussianform of h δ ( θ ′ − θ ) i | N f | requires that only X is nonzero,ie. that X j = 0 for all j >
1. In that case we get h δ ( θ − θ ′ ) i N f = e iN f θ √ πX e ( N f / X ∞ X k = −∞ e − ( θ +2 πk ) /X . (11)The sum is a compactified Gaussian with a width thatgrows as the square root of the volume (recall that X isextensive). Note that this fits perfectly with the generalrelation, Eq. (3), using Eq. (10). The Gaussian form isthe key approximation in the approach of [10, 11]. There-fore it is desirable to understand under which conditionswe can expect a Gaussian distribution of the phase.To examine if conditions exist such that X j = 0 for all j >
1, analytically from first principles is at present notfeasible since it is clearly as hard as to evaluate the fullQCD partition function. However, from the form of themoments (9) it is clear that any Feynman diagram withmore than three quark lines which contributes to the freeenergy of the moments will lead to a non-Gaussian form.Let us exemplify this general analytic understanding ofthe conditions under which the Gaussian follows by twoexamples.
First , suppose we evaluate the moments Eq. (6) us-ing the hadron resonance gas model. This involves a1-loop computation of all possible mesons and baryonswhich can be formed from the p + N f quark flavors(In practice these states correspond to the decompo-sition of n ⊗ ¯ n and n ⊗ n ⊗ n in SU (2( p + N f )) to SU ( p + N f ) flavor × SU (2) spin , followed by collection ofweights with each spin degeneracy. All details will beprovided in [15].). Since all possible meson and baryonstates have either 2 or 3 internal quark lines the combina-torics at 1-loop is unable to generate terms proportionalto ( p/ p/ N f )) or higher. Hence, the summationupon p in Eq. (5) necessarily leads to the Gaussian formin the hadron resonance gas as long as we work at 1-looporder. Second , if we instead compute the moments Eq. (6)in the combined strong coupling and hopping parameterexpansion, then to third order in the hopping parameter(see [16] for the form of the hopping parameter expan-sion to third order) it is again impossible to generate( p/ p/ N f )) and higher order terms, simply becausethe number of winding loops is insufficient. Explicitly inthe confined phase, we find order by order in the hoppingparameter h (see [15] for details), X = 4 a h sinh ( µ/T ) N s and X j> = 0 (12)to second order, X = − a h sinh ( µ/T )18 λ N s and X j> = 0 (13)to 4th order and finally X = − a h sinh ( µ/T )(5 + 450 λ ) N s (14)and X j> = 0to 6th order. N s is the lattice volume and thus the X j are extensive as expected. The coupling λ is the leadingcoupling in an effective Polyakov line action derived fromthe Wilson action via a strong coupling expansion. Itsprecise relation to the gauge coupling is given in ref. [16].The constant a , which is O (1), depends on the choice oflattice action.We see that, starting at fourth order in the hoppingparameter expansion terms of order ( p/ p/ N f )) appear naturally. Furthermore, from Eqs. (12)-(15) weobserve that the corrections to the Gaussian form set inat increasing order in µ/T . This is completely general: Ifwe Taylor expand the free energy of the moments Eq. (6)in µ/T , terms of order ( p/ p/ N f )) j are present atorder ( µ/T ) j and higher. To order ( µ/T ) j the terms of order ( p/ p/ N f )) j are proportional to the 2 j thorder cumulant of θ ′ considered in [17].We conclude that one generically finds non-Gaussiancorrections to the distribution of the phase in QCD. How-ever, as we have seen, the Gaussian approximation has anatural interpretation in terms of combinatorical factorsresulting from quark lines. If the combinatorics in eval-uating the moments of the phase factor involve at mostthree quark lines then the Gaussian form results.As we shall now see, the distribution of the baryonnumber over θ helps us to understand what it takesto make the Gaussian approximation self-consistent formeasurement of the baryon number. Distribution of the baryon number:
To accessthe baryon number distribution we again proceed via theFourier transform h n B δ ( θ − θ ′ ) i N f = 12 π ∞ X p = −∞ e − ipθ h n B e ipθ ′ i N f , (15)where the moments are given by h n B e ipθ ′ i N f = 1 Z N f (16) × lim ˜ µ → µ ∂∂ ˜ µ (cid:28) det p/ ( / D + γ µ + m )det p/ ( / D − γ µ + m ) det N f ( / D + γ ˜ µ + m ) (cid:29) . The p -dependence of the log of these moments will ofcourse reduce to a polynomial in p/ p/ N f ) at ˜ µ = µ ,however, before we take this limit we must differentiatewrt ˜ µ . For ˜ µ = µ the p -dependence of free energy in themoments can be a more general polynomial in p Z N f (cid:28) det p/ ( / D + γ µ + m )det p/ ( / D − γ µ + m ) det N f ( / D + γ ˜ µ + m ) (cid:29) = exp( k + k p + k p + . . . ) , (17)for suitable functions k i . This directly leads to h n B e ipθ ′ i N f (18)= ( c + c p + c p + . . . ) e − p/ p/ N f ) X + ... , where c j ≡ lim ˜ µ → µ ∂ ˜ µ k j . We conclude that the distribu-tion of the baryon number takes the form h n B δ ( θ − θ ′ ) i N f = 12 π ∞ X p = −∞ e − ipθ (19) × (cid:0) c + c p + c p + . . . (cid:1) e − p/ p/ N f ) X + ... . This expression looks rather complicated, however, as weshall now demonstrate the entire contribution to the av-erage baryon number comes from the c term. The re-mainder of the terms are noise. Total derivative:
To see how the integration over θ singles out the c term we now rewrite the distributionof the baryon number in terms of total derivatives of thedistribution of the phase itself. One easily verifies fromEq. (19) that h n B δ ( θ − θ ′ ) i N f (20)= 12 π ∞ X p = −∞ (cid:18) c + c − i ∂∂θ + c ( − i ) ∂ ∂θ + . . . (cid:19) × e − ipθ e − p/ p/ N f ) X + ... = (cid:18) c + c − i ∂∂θ + c ( − i ) ∂ ∂θ + . . . (cid:19) h δ ( θ − θ ′ ) i N f . This amazing relation is completely general, our only as-sumption is that the free energy of the moments are poly-nomials in p . The integration over θ needed to obtain thefull baryon number is now trivial h n B i N f = Z dθ h n B δ ( θ − θ ′ ) i N f = c . (21)The terms proportional to c i with i ≥ c is to be measured. As we will demonstrate belowit is essential to exclude the noise in order to obtain theaverage baryon number. Baryon number in the Gaussian approximation:
To see that it is essential to exclude the noise let us againconsider the special case where X j = 0 for j > h n B δ ( θ − θ ′ ) i N f = 1 √ πX e N f X / N f iθ (22) × ∞ X k = −∞ e − ( θ +2 πk ) /X (cid:2) c − c ( N f + i θ + 2 πk ) X )+ c (cid:0) ( N f + i θ + 2 πk ) X ) + 2 X (cid:1) + . . . (cid:3) . The total derivatives of Eq. (20) include terms of differentpowers of volume (recall that X is extensive). This is thecore of the sign problem for measurements of the baryonnumber : Terms which are down with powers of volume inthe distribution contribute at leading order to the baryonnumber. If the contributions from (22) that go like 1 /X ,1 /X , ... are dropped the result is c − c N f + c N f + . . . rather than c as required, cf. Eq. (21). Self-consistency requirement:
As we have justseen it is essential to control the noise terms in orderto measure the baryon number. In other words we mustmeasure the distribution of the phase and check that wecan control the error on the derivatives of the distribu-tion. Clearly if we would use brute force to measureall derivatives the numerical problem would be exponen-tially hard. However, let us again turn to the analytic cases of the1-loop hadron resonance gas and the strong coupling ex-pansion to third order in the hopping parameter, wherethe distribution of the phase automatically becomes aGaussian: in both cases the polynomial in the exponentof Eq. (17) is at most third order in p and no higher thana second order derivative can appear in Eq. (20). There-fore, to make the Gaussian approximation self-consistentfor the measurement of the baryon number it is sufficientto measure the distribution of θ to an accuracy which isquadratic in the inverse volume. While this is a demand-ing numerical task it is clearly advantageous comparedto exponential accuracy.Similarly, it is possible to consistently include correc-tions to the Gaussian coming from X j with j > j = J if we are able to control the errorson the derivatives up to order 2 J . Orthogonality of signal and noise:
Such higherorder terms in θ vary increasingly rapidly with θ and thusare harder to capture. The Gaussian example, however,suggests a possible way out of this. Let us rewrite thedistribution of the baryon number in Eq. (22) as h n B δ ( θ − θ ′ ) i N f = 1 √ πX ∞ X k = −∞ (cid:20) c H (cid:16) iX N f / − θ − πk √ X (cid:17) + ic √ X H (cid:16) iX N f / − θ − πk √ X (cid:17) + . . . (cid:21) × e − ( iX N f / − θ − πk ) /X , (23)where H n are the Hermite polynomials (It is of coursenot accidental that the Hermite polynomials appear, itfollows because the Gaussian is the generating functionand we take derivatives, cf. Eq. (20).). After a shift of thecontour by iX N f / θ in (21) followsdirectly from the orthogonality of the Hermite polynomi-als Z d x H m ( x )H n ( x ) e − x = 2 n n ! √ πδ mn , (24)and the fact that H = 1. Note that the shift along theimaginary axis is proportional to the volume. Moreover,due to the orthogonality of the Hermite polynomials onthe Gaussian weight, the signal c H is in this sense or-thogonal to the noise produced by P i ≥ c i H i .The orthogonality is also present in general where thenoise is orthogonal to the signal in the following sense Z dθ h n B δ ( θ − θ ′ ) i c h δ ( θ − θ ′ ) i N f h n B δ ( θ − θ ′ ) i c j > h δ ( θ − θ ′ ) i N f h δ ( θ − θ ′ ) i N f = 0 . (25)This follows directly from the general relation, Eq. (20),where the signal is given by the c term and the noiseis the remainder, ie. the terms with j >
0. Note thatsince the pions do not carry baryon number they onlycontribute to the c j with j >
0, see [14, 15] for details.As demonstrated for unitary fermions [18] analytic in-sights on the form of the noise can be turned into a prac-tical numerical tool. If the orthogonality found abovecan be implemented numerically it has the potential tocancel the entire pion noise from the measurement of thebaryon number.
Conclusions:
The analysis of the moments has givenus a physical understanding of the Gaussian approxima-tion which is central to the histogram method: if thecombinatorics in evaluating the moments of the phasefactor involve at most three quark lines then the Gaus-sian form results. The Gaussian approximation thereforepotentially captures the dynamics of the 1-loop hadronresonance gas model as well as that of the strong couplingexpansion to third order in the hopping parameter. How-ever, we also note that non-Gaussian terms are generic,and increase in importance as the hopping parameter,chemical potential µ/T and lattice coupling β increase.These terms may therefore be quite significant for realQCD, which of course entails light quark masses and the β → ∞ limit.The results of this paper follow from the assumptionthat the free energy of the moments can be written as apolynomial in p . Since this is equivalent in nature to thecentral assumption of the replica method, see e.g. [12],we expect that it applies equally generally. Finally, when µ > m π where m π is the mass of the pion the momentsare dominated by Bose-Einstein condensation of pions.Such a condensate gives rise to a linear term in p , see[14]. It would be most interesting to extend the resultsof this paper also to this region. Acknowledgments:
We would like to thank JanRosseel, Simon Hands, Tim Hollowood, Gert Aarts, aswell as organizers and participants of the INT programon ’Quantum Noise’ for discussions. This work was sup-ported by U.S. DOE Grant No. DEFG03-92ER40711(JG) and the
Sapere Aude program of The Danish Coun-cil for Independent Research (JCM and KS). [1] P. de Forcrand, PoS LAT , 010 (2009)[arXiv:1005.0539 [hep-lat]]. [2] G. Aarts, PoS LATTICE , 017 (2012)[arXiv:1302.3028 [hep-lat]].[3] S. Chandrasekharan and U. -J. Wiese, arXiv:1108.2461[cond-mat.str-el].[4] C. Vafa and E. Witten, Phys. Rev. Lett. , 535 (1984).[5] C. Gattringer and C. B. Lang, Lect. Notes Phys. , 1(2010).[6] J. C. Osborn, K. Splittorff and J. J. M. Verbaarschot,Phys. Rev. Lett. , 202001 (2005) [hep-th/0501210]. ;J. C. Osborn, K. Splittorff and J. J. M. Verbaarschot,Phys. Rev. D , 065029 (2008) [arXiv:0805.1303 [hep-th]].[7] J. R. Ipsen and K. Splittorff, Phys. Rev. D , 014508(2012) [arXiv:1205.3093 [hep-lat]].[8] J. Ambjorn, K. N. Anagnostopoulos, J. Nishimuraand J. J. M. Verbaarschot, JHEP , 062 (2002)[hep-lat/0208025].[9] Z. Fodor, S. D. Katz and C. Schmidt, JHEP , 121(2007) [hep-lat/0701022].[10] S. Ejiri, Phys. Rev. D , 014508 (2008) [arXiv:0706.3549[hep-lat]]; S. Ejiri et al. [WHOT-QCD Collabo-ration], Central Eur. J. Phys. , 1322 (2012)[arXiv:1203.3793 [hep-lat]]; PoS LATTICE , 089(2012) [arXiv:1212.0762 [hep-lat]].[11] S. Ejiri, arXiv:1306.0295 [hep-lat].[12] P. H. Damgaard and K. Splittorff, Phys. Rev. D ,054509 (2000) [hep-lat/0003017].[13] K. Splittorff and J. J. M. Verbaarschot, Phys. Rev. D ,014514 (2008) [arXiv:0709.2218 [hep-lat]].[14] M. P. Lombardo, K. Splittorff and J. J. M. Verbaarschot,Phys. Rev. D , 054509 (2009) [arXiv:0904.2122 [hep-lat]].[15] J. Greensite, J. C. Myers and K. Splittorff, in prepara-tion .[16] J. Langelage and O. Philipsen, JHEP , 055 (2010)[arXiv:1002.1507 [hep-lat]].[17] S. Ejiri et al. [WHOT-QCD Collaboration], Phys. Rev.D , 014508 (2010) [arXiv:0909.2121 [hep-lat]].[18] M. G. Endres, D. B. Kaplan, J. -W. Lee and A. N. Nichol-son, PoS LATTICE , 017 (2011) [arXiv:1112.4023[hep-lat]].[19] In [14] the distribution of the phase and of the baryonnumber as a function of the phase was computed within1-loop chiral perturbation theory. Since there are onlytwo quark lines in the loop the computation necessar-ily resulted in a Gaussian distribution. Moreover, sincepions are baryon number neutral only the c1