Divergences in the quark number susceptibility: The origin and a cure

Contents lists available atScienceDirect

Physics Letters B

www.elsevier.com/locate/physletb

Divergences in the quark number susceptibility:

The origin and a cure

Rajiv V. Gavai

^a

, Sayantan Sharma

^b,∗

aDepartmentofTheoreticalPhysics,TataInstituteofFundamentalResearch,HomiBhabhaRoad,Mumbai400005,India bFakultätfürPhysik,UniversitätBielefeld,D-33615Bielefeld,Germany

a r t i c l e i n f o a b s t ra c t

Articlehistory:

Received10April2015

Receivedinrevisedform4July2015 Accepted16July2015

Availableonline21July2015 Editor:J.-P.Blaizot

Quarknumbersusceptibilityonthelattice,obtainedbymerelyaddingaμN termwithμasthechemical potential and N astheconservedquarknumber,hasaquadraticdivergenceinthecut-offa.Weshow thatsuchadivergencealreadyexists forfreefermionswithacut-offregulator.Whileonecaneliminate it in the free latticetheory by suitably modifyingthe action, as is popularly done, it can simply be subtractedoffaswell.Computationsofhigherordersusceptibilities,neededforestimatingthelocation oftheQCDcriticalpoint,thenneedalotfewernumberofquarkpropagatorsatanyorder.Weshowthat thismethodofdivergenceremovalworksintheinteractingtheory.

©2015TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).FundedbySCOAP³.

1. Introduction

Thephasediagramofthestronglyinteractingmatterdescribed byQuantumChromodynamics(QCD)hasbeenasubjectofintense research in the recent years. Usual weak coupling perturbative approach may work for sufficiently high temperatures. However, the gauge interactions are likely to be strong enough for tem- peratures close to QCD, the typical scale of QCD, necessitating strong coupling techniques. Lattice QCD is the most successful non-perturbativetechnique whichhasprovideduswithsome in- terestingresultspertaining to thephase diagram.It isnow fairly well known from independent lattice studies that the transition fromthehadrontothequark gluonplasma phaseatzerobaryon densityisacrossover [1–3].At non-zerodensity,orequivalently nonzeroquark chemicalpotential

μ

,onehasto faceasignprob- lem: quark determinant is complex. This does not allow for an importancesamplingbasedMonteCarlostudy.Severalwayshave beenadvocatedintherecentyearstocircumventthesignproblem in QCD [4–7]. From perturbative studies of model quantum field theorieswith thesame symmetriesasQCD [8]andchiral model investigations at T <<

μ

[9], a critical end-point is expected in the QCD phase diagram. If present, the critical-end point would resultinthedivergenceofthebaryonnumbersusceptibility.Thus itsTaylorexpansion[7]atfinitebaryondensityasaseriesin

μ

B/T can be used to compute the radius of convergence, and there-

*

Correspondingauthor.CurrentlyinBrookhavenNationalLaboratory.

E-mailaddress:sayantans@quark.phy.bnl.gov(S. Sharma).

fore,an estimate ofthe location ofthecritical end-point[10,11].

FirstsuchestimatesoftheradiusofconvergenceoftheTaylorse- ries have predicted the critical end-point to be at TE/Tc=⁰.94 and

μ

B/TE=¹.8(1) [11]. Recently, astudy ona finerlattice has suggested the continuum limit to be around T_E/T_c =⁰.94(1),

μ

B/TE=¹.68(5) [12]. In the heavy-ion experiments, the fluctu- ations of the net proton number could act as a proxy for the netbaryon number.TheSTARexperimentatBrookhavenNational Laboratory hasreportedthemeasurements forthefluctuationsof the net proton number for a wide range of center of mass en- ergy √

s, of the colliding heavy ions between 7.7 and 200 GeV.

At √

s=¹⁹.6 GeV theexperimentaldataareobserved[13] tode- viate from the predictions of the proton fluctuations for models which donot have acritical end-point,andis similarto the lat- tice QCD-based predictions [14] for a critical point, signaling its possiblepresence.Itwouldbethusimportanttohaveathorough understanding of the systematics of the lattice QCD results and makethemasmuchreliableaspossible.

In addition to the usual suspects, such as continuum extrap- olation oreffects duetothe finiteness ofthe lattice spacing, the scale-setting,andthe statisticalprecisionofthe measurements,a key new important factoris that the radius of convergenceesti- mate requires ratios of as many higher orders of quark number susceptibilities(QNSs)aspossible.Currentlythestateoftheartis theeighthorderQNS[10,11].Itisveryimportanttoverifywhether the existing results are stableif ratiosof further higherorder of QNS are takenintoaccount. In orderto calculatethe QNSof or- derm, onehastotake themth-derivativeofthefreeenergywith respecttothequarkchemicalpotential.Sincethepopularmethod http://dx.doi.org/10.1016/j.physletb.2015.07.036

0370-2693/©^{2015TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCCBYlicense}(http://creativecommons.org/licenses/by/4.0/).Fundedby SCOAP³.

ofincorporating the chemical potential on the lattice is through exp(±

μ

a)factorsmultiplyingtheforwardandthebackwardtem- poral gauge links respectively of the fermion operator [15,16], thereis an everincreasing proliferationof terms ofvarying sign asm increases.Their largenumber aswellasthe large cancella- tionsamongstthemataspecificordermakeitdifficulttoincrease m beyondeight atpresent.Introducing thechemical potential by a

μ

N-term,whereN isthecorrespondingconservedcharge,leads tobothmuchfewertermsandlessercancellationsatthesamem [17],therebyreducingthecomputationalcostupto60%ateighth order;more savings oughtto accrue by going toeven higheror- ders.Notonlywillthisimprovetheprecisionofthelocationofthe criticalpoint butmorepreciseTaylorcoefficientsandmoreterms intheTaylorexpansioncanpotentiallyalsoleadtoabettercontrol oftheQCDequationofstateatfinitebaryondensitywhichwillbe neededfortheanalysisoftheheavy-ion datafromthe beamen- ergyscan at RHICas well asthefuture experimentsatFAIR and NICA.

Inthispaper,wediscusswhethersuch alinearin

μ

approach isviable orhasunsurmountableproblemsby comparingwiththe usualexponentialin

μ

method.InSection2,werevisit thenum- berdensityfornon-interactingfermionsinthecontinuumusinga cut-offregulator.Wepointout thatdivergencesappearalreadyin thecontinuum freetheory whenthe cut-offregulator istakento infinity contraryto the commonknowledge. We then discuss an approachtotacklethisdivergenceinthefreetheory.Byperform- ingcontinuum extrapolationof thesecond andfourthorderQNS forquenched QCD for the linear in

μ

method,we validate it in Section3.Thisisthemostimportantresultofourpaper.Wedis- cussitspossibleconsequencesandtheextensionstohigherorder QNS.

2. Thermodynamicsofnon-interactingfermions

QCD thermodynamics can be derived from its partition func- tion,writteninthepathintegralformalism[18]as,

Z

=

DA_μD

ψ ¯

D

ψ

×

^e01/Tdτd3x

−¹/2 Tr(F²μ,ν)− ¯ψ(γμ(∂μ−^igAμ)+^m−μγ4)ψ

,

(1) whereψ,ψ¯ and Aμrepresentthequark,anti-quarkandthegluon fieldsrespectively,whosethecolorindicesarenotwrittenexplic- itlyabove.

μ

isthe chemicalpotential forthe netquark number withthe correspondingconserved chargebeing

d³xψ¯

γ

4ψ.Gen- eralizationstovariousconservedflavor numbersisstraightforward.

Forsimplicity,wewillconsideronlyasingleflavor withthebary- onic chemical potential

μ

B=³

μ

q. Appropriate derivatives of Z leadto various thermodynamicalquantities, e.g., thequark num- ber density, or equivalently (1/3) the baryon number density, is definedas,

n

=

^T V

∂

lnZ

∂ μ ^|

^{T=fixed (2)}

Earlierattempts to discretize the above theory to investigatethe finite baryon density physics on a space–time lattice revealed

μ

-dependentquadraticdivergencesinthenumberdensityandthe energydensitywhen the chemical potential is introduced in the quarkDiracoperatorbymultiplyingitwiththecorrespondingcon- servedchargeonthelattice.Thesedivergences,whichappearasa

μ

/a² terminthe expression forthelattice numberdensitywith a asthelatticespacing, arepresentevenifthegauge interactions areabsent.Throughexplicitcalculationofthenumberdensityfor non-interactingfermionsonthelattice,itwasthenshown[15,16,

19]thatsuitablemodificationofthe

μ

N termintheaction,elimi- natesthesedivergenttermsonthelattice,andyieldsafinitea→⁰ continuum limit.NumericalstudiesoftheQNSfortheinteracting theory subsequently confirmed that once the free theory diver- gences are thus eliminated, no further divergences arise [20,21].

Asuccinct waytodescribe all thevarious actions istointroduce functions f(

μ

a)[^g(

μ

a)]^asthemultiplying factorsfortheforward (backward)timelikegaugefieldsonthelattice.Whileforthenaive discretization, f=¹+^a

μ

andg=¹−^a

μ

leadstoadivergentbary- onicsusceptibilityinthecontinuumlimit,thechoice f=^exp(a

μ

) andg=^exp(−^a

μ

)doesnot.

Clearlysinceallderivativesof f andg arenonzerofortheex- ponentialcase,whereasonlythefirstderivativeisnonzeroforthe linearcase,higherorderQNSare alotsimplerforthe latter.Fur- thermore, for fermions with better chiral properties such as the Overlap fermions or the Domain Wall fermions, the exponential formleads to a loss [22] ofthe exactchiral symmetry onlattice fornonzero

μ

. Indeedthe onlychiral symmetry preserving form these fermions have for finite

μ

and a is the linear form [23].

This motivatesusto revisit theissue ofthe natureandorigin of thesedivergenceswhen thechemicalpotential enters linearlyin- steadoftheexponentialform.Asweshowbelow,thedivergences arepresentforthecontinuumfreefermionsaswell,andthelattice regulatorsimplyfaithfullyreproducesthem.Whileonecanemploy thefreedomoflatticeactiontoeliminatethem,itisnotnecessary.

Indeed, one can perhaps employ simpler subtractionmethods to eliminatethem,aswedemonstrateinthispaper.

2.1. Continuumfreefermions

Results for the continuum free fermions are easily found in textbooks [18]. We review them below solely with the idea of pointing out explicitly the

μ

-dependent divergences present in them. Forsimplicity, we consideronly massless fermions though thisderivationcanbe easilyextendedforfinitemass.Theexpres- sion for the number densityfor free fermions is easily obtained fromthedefinitionsaboveas

n

=

^4iT

^∞

j=−∞

d³p

(

2

π )

³

( ω

j

+

ⁱ

μ )

p²

+ ( ω

j

+

i

μ )

²

≡

^4iT

∞ j=−∞

F

( ω

_j

, μ ) ,

(3)

where p²=^p2₁+^p2₂+^p2₃ ^and

ω

j=(2 j+¹)

π

T . Here we choose the gammamatricesto be all Hermitian asiscommon inlattice studies.Thecontinuum conventionfollowedinthestandard texts hasonly

γ

4 asHermitianandtheothergammamatricesareanti- Hermitian.TheexpressioninEq.(3)canbeevaluatedbytheusual trick of converting the sum over energy states to a contour in- tegral.TheMatsubarafrequencieslieonthereal

ω

-axis.Following [18]again,onecanemployaninfinitesimallysmallcontouraround theeachpoleonthereal

ω

axistorepresentthe

ω

j-sum,andob- tain

2

π

T

j

F

( ω

j

, μ )

=

Lt_→0

⎡

⎢ ⎣

∞+

i

−∞+ⁱ

F

( ω , μ )

d

ω

e^iω/T

+

¹

+

−∞−

i

∞−ⁱ

F

( ω , μ )

d

ω

e^iω/T

+

¹

⎤

⎥ ⎦ .

(4)

ThelineintegralsinEq.(4)caninturnbewrittenintermsofcon- toursinthe upperandlower complex

ω

planes.Using the exact identity,

Fig. 1. Thecontourdiagramforcalculatingthenumberdensityforfreefermionsat zerotemperature.Pdenotesthepole.

F

( ω , μ )

e^iω/T

+

1

=

^F

( ω , μ ) −

^F

( ω , μ )

e^−iω/T

+

1

,

(5)

the line integral in the Im

ω

>0 plane at infinity can be made convergent. The numberdensity expression atfinite temperature anddensityis,

n

=

²ⁱ

π

⎡

⎣

Imω<0

F

( ω , μ )

d

ω

e^iω/T

+

¹

−

Imω>0

F

( ω , μ )

d

ω

e^−iω/T

+

¹

+

∞

−∞

F

( ω , μ )

d

ω

⎤

⎦ .

(6) We note that the last term of the above expression contributes to number density at all temperatures, T and

μ

and has only

μ

-dependence. The first two terms yield the usual Fermi–Dirac distribution functions. These haveno ultraviolet divergences [24]

sincetheultravioletmodesareexponentiallysuppressed.Inorder toexaminethelasttermindetail,letuswriteitexplicitly:

n

=

⁴ⁱ

∞

−∞

d

ω

2

π

d³p

(

2

π )

³

( ω +

ⁱ

μ )

p²

+ ( ω +

ⁱ

μ )

²

.

(7) Underavariabletransformation

ω

+ⁱ

μ

=

ω

,itcanberecastas

n

=

⁴ⁱ

∞+

iμ

−∞+iμ d

ω

2

π

d³p

(

2

π )

³

ω

p²

+ ω

²

^.

⁽⁸⁾

Inordertocomputetheintegral carefully,weimposeacut-off on all thefour momenta anduse thecontour in the complex

ω

planeinFig. 1,leadingto n

=

2i

d³p

(

2

π )

³

×

⎡

⎣ −

ⁱ

( μ ₋

p

) −

2

+

4

+

1

d

ω π

ω

p²

+ ω

²

⎤

⎦ .

(9)

Whilethefirsttermarisesfromtheresidueofthepoleofthe integrand enclosed by the contour, the last three terms in the Eq.(9)arise fromclosingthecontour. Thelineintegral1isiden- tically zerobecausetheintegrandisan oddfunction.The sumof thelineintegrals2and4is

2+⁴

= −

¹ 2

π

^ln

p²

+ ( +

ⁱ

μ )

² p²

+ ( −

ⁱ

μ )

²

.

(10)

Since>>

μ

,beingthecut-off,expandingthelogarithm and retainingtheleadingterm,

2+4

= −

⁴ⁱ

μ

2

π (

²

+

p²

) .

(11)

It is straightforwardto dothe remaining momentum integralsin Eq. (9).The first term leads to the usual

μ

³ term. However, the sum of the two line integrals 2 and4 in Eq. (11)yields a

μ

² divergenceintheexpressionfornumberdensityasbelow,

−

∞ 0

dp 2

π

³

4

μ

p²

²

+

p²

= −

⁴

μ

0

dp 2

π

³

1

−

²

²

+

p²

= −

²

μ

²

π

³

1

− π

4

.

Notethattheleadingdiverging³-contribution,presentfor

μ

=^0, isthesameforthelineintegrals2and4,anddoescancel.Itisthe non-leading

μ

-dependent termwhichleads tothe² divergence above. Thisdivergenceshowsupifone usesPauli–Villarsmethod as well. One then has to introduce additional Pauli–Villars fields to cancelthis

μ

² divergencefromthefree energydistinctfrom thosewhicharerequiredtocanceltheusual⁴ divergence.

Similarly, the lattice as a cut-off regulator forthe free theory also leads to

μ

a⁻² divergence.Using nonzero T as theregulator, followingthemethodoutlined in[18],yields thecontributionsof onlythetwotermsinEq.(6).Thischoiceoftheregulatordoesnot permit the T -independent termof Eq.(6). Such a choice ofreg- ulator is,however,not feasibleforlattice QCDcomputations,and indeed,manyother interactingtheories.Inordertoobtain physi- callymeaningfulresult,thecontributionofthelineintegrals2and 4 has to be subtracted off. This free theory subtraction, though

μ

-dependentisanalogoustothesubtractionfromthepressureat T=^{0,commonlyusedonthelatticeintheequationofstatecom-} putations.Weshowinthenextsectionthatnofurtherdivergences are observed once the free theory divergence is subtracted from thequarknumbersusceptibilities.

3. Quenchedresultsonthelattice

The analytical proof of Refs. [15,16,25] for the lack of diver- gences in the quark number susceptibilityin the exponential

μ

case, outlined briefly abovein theprevious section, was fornon- interactingfermions.Noequivalentproofexistsfortheinteracting fermions even for the exponential case. On the other hand, it is easy to check that for staggered quarks, the chiral symmetry is maintained for

μ

=^{0.Forthelinear}

μ

case,oneeven dealswith conservedbaryoniccurrentsonthelattice,a

μ

beingthecoefficient oftheconservedbaryon numberonthelattice.Onethereforeex- pectsnofurtherdivergencestoarise,andnoextrarenormalization needed, afterswitching on the gauge interactions in eithercase.

Thiswasexplicitlycheckedbynumericalsimulationsintheexpo- nential(indeed,genericallyforany f ·^g=^{1) caseforQCD inthe} quenchedapproximation[21,27].It wasproposed inRefs. [17,26], andagaindemonstratedforthenon-interactingfermions,thatthe spurious anddivergent termsarising inthe linear

μ

casecan be evaluated andexplicitlysubtracted.Cleartestsofsuch aproposal forthe interactingcasearethat the continuumlimit ofa→^{0 of} the so-subtracted quark number susceptibilities should i) exhibit no additionaldivergencesandii)yieldthe sameresultasforthe exponential

μ

form.Inthissection,wereportresultsofourthese numerical testsfor the linear

μ

casein quenchedQCD andver- ify that it passes both the tests, as expected. The choice of the quenched approximation was governed by the fact that the cor- responding publishedavailable results[21,27] forthe exponential casemakeitsimplertocompare.

OnanN³×^NT lattice,thetemperatureisgivenbyT=¹/(NTa), where a the lattice spacing is governed by the gauge coupling β=⁶/g². We employ standard Wilson plaquette action for the gauge fields and use the staggered quarks for our susceptibility determinationswiththecorresponding Diracmatrix D(a

μ

)given by

D

( μ )

_xy

=

3 i=¹

η

_iU_i

(

x

)δ

_x,y−ˆi

− η

_iU_i^†

(

y

)δ

_x,y+ˆi

− (

1

− μ

a

) η

4U₄^†

(

y

)δ

_x,y+ˆ4

+ (

1

+ μ

a

) η

₄U₄

(

x

)δ

_x,y−ˆ4

+

^ma

δ

_x,y

.

(12) Herema isthequarkmassand

η

’saretheusualstaggeredfermion phases.OurquenchedQCDconfigurationsweregeneratedbyusing theCabibbo–Marinaripseudo-heatbathalgorithmwiththreeSU(2) subgroupupdateper sweep usingtheKennedy–Pendleton updat- ing method.We chose to simulateat two differenttemperatures andquarkmassesonavarietyoflatticesizes,aslistedinTable 1, byselecting suitable β values [21,27]such that T was heldcon- stant as NT (or equivalently a) was increased (decreased). This enabledusto make a continuum limit extrapolationat both the temperatures.Wequote thetemperaturesintheunitsofthecrit- icaltemperatureTc correspondingto thefirstordertransitionfor SU(3),definedbyusingtheorderparameter,thePolyakovloop.Al- thoughwedonotneeditexplicitlyanywherebelow,wemention thatTc=²⁷⁶(2)MeV[28]inthecontinuumlimit,usingthestring tensionvaluetobe425MeVtosetthescale.NotingfromEq.(12) that only the first derivative withrespect to a

μ

, D,is nonzero, thequark numbersusceptibilityforthislinearchemical potential actionis:

χ

20

=

^T V

tr

(−

^D−1DD−1D

) +

^tr

(

D⁻¹D

)

²

(13)

Thisexpressionissimilar tothat inRef. [10]butwithout the D termwhichis identicallyzerohere, as arefurther higherderiva- tiveswith respect to

μ

. Our notation is same as in Ref. [10] to facilitatecomparisonofournumericalresultsbelow.Thetracesin theaboveexpressionwerecomputedstochasticallyusingGaussian randomvectors.FromtheMonteCarlotimeevolutionofthediffer- entoperators that enterthe susceptibilitycomputation,including thetwotermsaboveseparately,we estimatethattheautocorrela- tionlengthismuchlessthan1000 sweeps.Inordertoensurethat ourmeasurementsarestatisticallyindependent,theyweredoneon configurations¹ separatedby1000 heatbathsweeps andexcluding thefirst5000–10000 sweepsforthermalization.SuchN_configs con- figurations, whichvaried from24–100, were employed to obtain thethermodynamicaverages.Thedetailsofthenumberofrandom vectorsandnumberofconfigurationsusedateachquarkmassand temperature,aresummarizedinTable 1.

As discussed in the previous section, the choice of D(a

μ

) in Eq. (12) leads to a QNS with a term ∝¹/a² for the free case.

Ateachvalue oflatticecut-off,we computednumericallythe co- efficient for this 1/a²-term for non-interacting fermions on the correspondingN³×∞^{latticeandsubtracteditfromthecomputed} valuesofthesusceptibilityintheinteractingcase. Ifthere areno additionaldivergences inthe interactingtheory, one expects the continuum extrapolation in 1/N²_T ∼^{a2 performedon these sub-} tracted values of the susceptibilities to have a smooth limit. On theother hand,if furtherdivergences do exist inthe interacting theory thenthe 1/a² or equivalently fora fixed temperaturethe

1 TheconfigurationsgeneratedwererotatedtothezeroPolyakovloopsectorto makethemsimilartofullQCD.

Table 1

Theparametersforthelatticesimulations.

T m/Tc β N NT Nconfigs Nrvec

1.25Tc 5.788 16 4 100 500

6.21 24 8 50 500

0.1 6.36 30 10 60 500

6.505 36 12 24 500

1.25Tc 5.788 16 4 100 500

0.01 6.21 24 8 48 500

6.36 30 10 68 500

6.505 36 12 24 500

6.0609 16 4 100 400

6.3331 32 6 50 400

2Tc 0.1 6.45 24 8 80 400

6.75 22 10 80 400

Table 2

Thecontinuumextrapolationresultsforthesecondandfourthorderdiagonalsus- ceptibility.

Observable m/Tc T/Tc c1 c3

1.25 0.838(8) 9.1(2)

0.1

χ20/T² 2 0.94(5) 10(1)

0.01 1.25 0.839(8) 9.4(2)

1.25 0.58(1) 14.9(8)

χ40 0.1

2 0.50(2) 12(2)

N²_T dependenttermwouldsurviveandincreaserapidlytoblowup inthecontinuumlimit.

The results for the dimensionless second order susceptibility

χ

20/T², after the subtraction of the free results at 1.25Tc and 2Tc are displayed in Fig. 2 for the same physical quark mass m/Tc=⁰.1,andintheleftpanelofFig. 3forthesmallerm/Tc= 0.01. For a comparison, we also plotted the available data from [27] ofthesame quantityat 1.25Tc calculated usingtheconven- tional exponential method, in the left panel of Fig. 2. The cut- off effects in the exponential method are larger than the linear method.The difference reduces onfiner lattices, thus converging tothesamecontinuumlimit,inagreementwithexpectationsfrom universality.Ongeneralgrounds,weexpect

χ

20/T² shouldbehave as

χ

20

T²

=

c1

(

T

) +

c2

(

T

)

N²_T

+

^c3

(

T

)

N²_T

+

O

(

1

/

N⁴_T

)

(14) Sincethevaluesofthesecondordersusceptibilitiesreducewith increasing NT in all theplots, thereis clearly no signof anydi- vergence in the interacting theory withc₂(T)=^{0. This is so ir-} respective ofthetemperature, andeven afterlowering thequark massbya factorof10.Ourresultsfortheleastsquaresfitofthe dataforthecoefficientsc₁andc₃aresummarizedinTable 2.Fora propercomparisonwiththeearlierdata [27],wealsoincludedthe NT=^{4 pointsatboth the}temperatures.A measureofhow good theleastsquaresfitrepresentsthedataisgivenbyaquantity R². It is definedas theratio of the sumof squaresin the fitmodel tothe totalsumofsquares.The valueof R²=^{1 would therefore} representthebestfit.Thevaluesforourfitat1.25Tc and2Tc for m/Tc=⁰.1 areR²=⁰.9982 andR²=⁰.9625 respectivelyandthat form/T_c=⁰.01 is R²=⁰.9977.Wealsotriedfitswithalogarith- micterminEq.(12)oftheformln(1/N²_T)bothinplaceof,andin additionto,the 1/N²_T term.While theformerseemsstronglyun- favoured, withan increasein

χ

² bya factorofthree to five,the latteris onlymarginally ruled out.Additional latticeswithlarger N_T areneededtobedefinitiveforthelattercase.Thuseitheronly

Fig. 2. Thesecondordersusceptibilityat1.25Tc(leftpanel)and2Tc(rightpanel)form/Tc=⁰.1.Forcomparingthecut-offeffectswealsoincludetheavailabledatafrom [27]forthesamequantity,calculatedwiththeconventionalexponentialmethod.Itisshownbysquaresintheleftpanelofthefigure.

Fig. 3. Thesecondordersusceptibilityat1.25Tc showsnodivergenceinourmethodevenform/Tc=0.01 (leftpanel).Rightpaneldisplaysthecontinuumextrapolated resultsforthediagonalfourthordersusceptibilityform/Tc=⁰.1 attwodifferenttemperatures.

apower lawdivergence oronlyalogarithmic divergenceisruled out fromour present data.As seen in Fig. 2, our continuum ex- trapolatedvaluesareinagreementwiththecorrespondingresults obtainedusingtheexponentialin

μ

action[27]atthesetempera- tures.Wealsoverifiedthatthereisnomassdependentdivergent termintheexpression forthesecond ordersusceptibilityby per- forming a similarcontinuum extrapolationusing a differentbare quark massofm/Tc=⁰.01,asshownin left panelof Fig. 3.The secondorderoff-diagonalsusceptibilityisidenticalinboththelin- earandtheexponentialmethod.Sinceit iszeroforfree fermion, noadditionalsubtractionisnecessaryforthelinearcase.

Calculationsforfreefermionsshowthatthefourthorderquark number susceptibility has no divergent contributions butan ad- ditionalfinite contribution in the continuum limit forthe linear

μ

-case. Adopting the same procedure for it aswell, we subtract the additional obtained free contribution from the correspond- ing (quenched) lattice QCD determination. For the fourth order susceptibility, there are one diagonal and two off-diagonal com- ponents. Following the convention of [10], these can be written as,

χ

40

=

^T V

O1111

+

6O112

+

4O13

+

3O22

+

O4

−

³ O11

+

O2

²

, χ

22

=

^T

V

O1111

+

²O112

+

O22

−

O11

+

O2

²

−

²O11

²

, χ

31

=

^T

V [O1111

+

3O112

+

O13

−

3 O11

+

O2 O11

]

,

(15) where the operators On satisfy the identities On=On+¹ and Oi j=Oi·Oj and so on. The number density is given by n= T O1/V . O2=^tr(−^D−1DD−1D) was the source of the diver- gence in

χ

20, which was cured by employing O2−O₂^{free div, as}

discussed above. In order to be consistent, a subtraction ofsuch aconstantfromO2 intheexpressionsaboveshouldalsobedone.

It can be easily verified that such a substitution in the expres- sions abovedoesnot changethematall sinceall thefreetheory divergence termsarising out ofitcancel ineach expression.This is, of course, consistent with the fact the direct computation of thefree diagonalsusceptibility

χ

40 hasnodivergenceinthecon- tinuumlimit forthelinear

μ

-case either.Itdoeshavea constant a⁰ term as an artifact though in the term coming from O4. In- deed, it is clear that due to dimensional reasons any difference between the linear andexponential case must be a constant for O4. Moreover, all On for n>4 will only differ by terms which vanishinthe continuumlimit, beingoforderaⁿ⁻⁴. Thusthestill highersusceptibilities mustallagreeinthecontinuum limitwith theexponential

μ

-case.Inspiredbythesuccessforthesecondor- der susceptibility,we computed the O4 for free fermions at the same temperature andsubtracted it fromthe value obtainedfor the quenched theory at that temperature in order to eliminate the a⁰ term artifact. Our results for different N_T are displayed in the right panel of Fig. 3. These results for

χ

40 also show a convergence withincreasing N_T.Thus there are indeedno addi- tional divergences in the continuum limit, as anticipated. More- over, the subtraction of the unphysical artifact appears to have been done correctly on each lattice size. The continuum results for the exponential caseare not available in the literature to fa- cilitate a comparison unlike Fig. 2. The results do show a con- vergingtrendthough.Theresultsforthecontinuumextrapolation are tabulated in Table 2. In all these fits the NT =^{4 data has} not been considered since it clearly stands out of thetrend. The R² forthese fits are 0.996 and 0.978 at 1.25Tc and 2Tc respec- tively.

Fortheoff-diagonalsusceptibilitiesatfourthorder,thefreethe- oryartifactsareagainzerosonosubtractionareexpected.Indeed,

Fig. 4. The fourth order off-diagonal susceptibilityχ22(left panel) andχ31(right panel) for m/Tc=⁰.1.

weobserve thistobe truefor

χ

22 inFig. 4.For

χ

31 intheright panelofthesamefigure,itappearssomewhatdifficulttodrawany definiteconclusionsalthoughafinitecontinuumlimitissuggested.

4. Summary

InvestigationsofQCDatfinitedensity,andinparticularofthe QCD critical point, gain from using the canonical Lagrange mul- tipliertype linearchemical potential terminthe fermionactions onthelattice.Preservationofexactchiralinvarianceonthelattice seemsfeasible[23]onlyforsuchalineartermfortheoverlapand the domain wall fermions. The higher order quark number sus- ceptibilitiesneededforlocatingthecriticalpoint usingthe Taylor expansion approach are easier to compute in the linear case as well.However,itisknown[15,16,19,25]sincelongthatthelinear termleadstoO(1/a²)divergencesinthebaryonnumbersuscepti- bility.Wehaveshownthatsuchadivergenceexistsalreadyinthe continuumforagasoffreefermions,andtherefore,latticemerely faithfully reproduces it. Using simulations of the quenched QCD withthestaggered fermions, we haveverified that oncethe free theorydivergenceisexplicitlysubtractedout,thesusceptibilityhas noadditionaldivergenceinthecontinuumlimit.Thisisonlytobe expectedsincetheconservedcharge,ornumberdensity,doesnot getrenormalizedinaninteractingtheory.Furthermore,itsextrap- olated value inthe continuum agrees very well at two different temperatureswith thesimilar continuum resultsfor thefermion actionusingtermsexponentialin

μ

,whichbyconstructionisfree fromanysuch divergences.Thehigherordersusceptibilities were also shown to be free of divergences in quenched QCD, and it wasargued whythiswas to beexpected. Further workis clearly needed to check that this conclusion hold for the theory with dynamicalquarksaswell,althoughoneexpectsthemaximumpos- sibledifferencetoariseintheimportanceofthelogarithmicterm whosenumericalsignificancemay notbe felt evenwiththe cur- rentbestprecision.

Animportantconsequenceofourstudyisthat thiswoulden- ablecomputationofthehigherorderQNSinasignificantlyshorter timeandwithbettercontroloferrors,therebyenablinguseofra- tiosofstillhigherorderQNSinlocatingtheQCDcriticalpointand foramorepreciseequationofstateatfinitebaryondensity.These couldalso be exciting fortheheavy ionexperiments whichhave alreadyreportedpreliminaryhintsofapossiblecriticalpointand arebeginningtoprobethefinitedensityregionintheongoingand futureprograms.

Acknowledgements

Thisworkwascompleted duringthevisitofoneofus(R.V.G.) to the Universität Bielefeld, Germany. He is very happy to ac- knowledge thekindhospitalityofits PhysicsDepartment, inpar- ticular that of Profs. Frithjof Karsch andEdwin Laermann. R.V.G.

also gratefully acknowledges the financial support from the J.C.

BoseNationalFellowship ofDepartmentofScience& Technology, Government ofIndia(grantNo.SR/S2/JCB-12/2009).S.S.wouldlike to thank Robert Pisarski, Marc Sangel and Sören Schlichting for manyhelpfuldiscussions.

References

[1]C.Bernard,etal.,MILCCollaboration,Phys.Rev.D71(2005)034504.

[2]M.Cheng,N.Christ,S.Datta,J.vanderHeide,C.Jung,etal.,Phys.Rev.D74 (2006)054507.

[3]Y.Aoki,G.Endrodi,Z.Fodor,S.Katz,K.Szabo,Nature443(2006)675.

[4]Z.Fodor,S.Katz,J.HighEnergyPhys.0203(2002)014.

[5]Ph.deForcrand,O.Philipsen,Nucl.Phys.B642(2002)290;

M.D’Elia,M.-P.Lombardo,Phys.Rev.D67(2003)014505.

[6]K.-F.Liu,Int.J.Mod.Phys.B16(2002)2017;

A.Li,A.Alexandru,K.-F.Liu,Phys.Rev.D84(2011)071503.

[7]R.V.Gavai,S.Gupta,Phys.Rev.D68(2003)034506;

C.R.Allton,etal.,Phys.Rev.D68(2003)014507.

[8]R.D.Pisarski,F.Wilczek,Phys.Rev.D29(1984)338.

[9]M.Asakawa,K.Yazaki,Nucl.Phys.A504(1989)668;

A.Barducci,R.Casalbuoni,S.DeCurtis,R.Gatto,G.Pettini,Phys.Lett.B231 (1989)463.

[10]R.V.Gavai,S.Gupta,Phys.Rev.D71(2005)114014.

[11]R.V.Gavai,S.Gupta,Phys.Rev.D78(2008)114503.

[12]S. Datta,R.V.Gavai,S. Gupta,in:ProceedingsofQuarkMatter2012, arXiv:

1210.6784[hep-lat].

[13]L.Adamcyzk,etal.,STARCollaboration,Phys.Rev.Lett.112(2014)032302.

[14]R.V.Gavai,S.Gupta,Phys.Lett.B696(2011)459.

[15]P.Hasenfratz,F.Karsch,Phys.Lett.B125(1983)308.

[16]J.Kogut,etal.,Nucl.Phys.B225(1983)93.

[17]R.V.Gavai,S.Sharma,Phys.Rev.D81(2010)034501.

[18]J.Kapusta,C.Gale,ThermalFieldTheory,CambridgeUniversityPress,2006.

[19]N.Bilic,R.V.Gavai,Z.Phys.C23(1984)77.

[20]S. Gottlieb,W.Liu,D. Toussaint,R.L.Renken,R.L.Sugar, Phys. Rev.Lett. 59 (1987)2247.

[21]R.V.Gavai,S.Gupta,Phys.Rev.D67(2003)034501.

[22]D.Banerjee,R.V.Gavai,S.Sharma,Phys.Rev.D78(2008)014506.

[23]R.Narayanan,S.Sharma,J.HighEnergyPhys.1110(2011)151;

R.V.Gavai,S.Sharma,Phys.Lett.B716(2012)446.

[24]D.Gross,R.Pisarski,L.Yaffe,Rev.Mod.Phys.53(1981)43.

[25]R.V.Gavai,Phys.Rev.D32(1985)519.

[26]R.V.Gavai,S.Sharma,Phys.Rev.D85(2012)054508.

[27]S.Mukherjee,Phys.Rev.D74(2006)054508.

[28]F.Karsch,Nucl.Phys.B,Proc.Suppl.83(2000)14.