Axion cosmology, lattice QCD and the dilute instanton gas

Contents lists available atScienceDirect

Physics Letters B

www.elsevier.com/locate/physletb

Axion cosmology, lattice QCD and the dilute instanton gas

Sz. Borsanyi

^a

, M. Dierigl

^b

, Z. Fodor

^a,c,d

, S.D. Katz

^d,e,∗

, S.W. Mages

^f,c

, D. Nogradi

^d,e,g

, J. Redondo

^h,i

, A. Ringwald

^b

, K.K. Szabo

^a,c

aDepartmentofPhysics,WuppertalUniversity,Gaussstrasse20,D-42119Wuppertal,Germany bDeutschesElektronen-SynchrotronDESY,Notkestrasse85,D-22607Hamburg,Germany cIAS/JSC,ForschungszentrumJülich,D-52425Jülich,Germany

dInstituteforTheoreticalPhysics,EötvösUniversity,PázmányPetersétany1/A,H-1117Budapest,Hungary eMTA-ELTELendületLatticeGaugeTheoryResearchGroup,Budapest,Hungary

fUniversityofRegensburg,D-93053Regensburg,Germany

gKavliInstituteforTheoreticalPhysics,UniversityofCalifornia,SantaBarbara,CA93106-4030,USA hDepartamentodeFísicaTeórica,UniversidaddeZaragoza,PedroCerbuna12,E-50009Zaragoza,Spain iMax-Planck-InstitutfürPhysik(Werner-Heisenberg-Institut),FöhringerRing6,D-80805München,Germany

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

Articlehistory:

Received3September2015

Receivedinrevisedform6November2015 Accepted9November2015

Availableonline11November2015 Editor:B.Grinstein

Keywords:

Axiondarkmatter QCDonthelattice Instantons

Axionsare one ofthe mostattractive dark mattercandidates.The evolutionoftheir numberdensity intheearlyuniversecanbedeterminedbycalculatingthetopologicalsusceptibilityχ(T) ofQCDasa functionofthetemperature.LatticeQCDprovidesanabinitiotechniquetocarryoutsuchacalculation.

A full result needs two ingredients:physical quark massesand acontrolledcontinuum extrapolation fromnon-vanishingtozerolatticespacings.Wedetermineχ(T)inthequenchedframework(infinitely large quarkmasses)and extrapolateitsvaluestothe continuumlimit. Theresultsare comparedwith the prediction of the dilute instanton gas approximation (DIGA). A nice agreement is found for the temperaturedependence,whereastheoverallnormalizationoftheDIGAresultstilldiffersfromthenon- perturbativecontinuumextrapolatedlatticeresultsbyafactoroforderten.Wediscusstheconsequences ofourfindingsforthepredictionoftheamountofaxiondarkmatter.

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

1. Introduction

Oneofthegreatest puzzlesinparticlephysicsisthenatureof darkmatter. A prominentparticle candidate forthe latter is the axion A [1,2]: apseudoNambu–Goldstonebosonarising fromthe breakingofa hypotheticalglobalchiral U(1)extension [3]ofthe StandardModelatan energyscale f_A muchlargerthantheelec- troweakscale.

A key input for the prediction of the amount of axion dark matter [4–6] is its potential as a function of the temperature, V(A,T).ItisrelatedtothefreeenergydensityinQCD, F(θ,T)≡

−^lnZ(θ,T)/V,via V

(

A

,

T

) ≡ −

_V^1ln

Z

(θ,

T

)

Z

(

0

,

T

)

|

θ=A/fA

,

(1) where V ^{is the Euclidean space–time volume. Here A}(x) is the axion field and fA is the axion decay constant. The axion field

*

Correspondingauthor.

E-mailaddress:katz@bodri.elte.hu(S.D. Katz).

hasmass dimensionone, therefore A(x)/f_A isdimensionlessand canbe interpretedasan xdependentθ value.Theangleθ enters theEuclideanQCDLagrangianviatheadditionalterminvolvingthe topologicalchargedensityq(x),

−

ⁱ

θ

q

(

x

) ≡ −

ⁱ

θ α

s

16

π

_μνρσF^a_μν

(

x

)

F^a_ρσ

(

x

),

(2) with Fμν^a beingthe gluonicfield strength and

α

s≡^g2s/(4

π

) the finestructureconstantofstronginteractions.

Ongeneralgrounds,thefreeenergydensityandthustheaxion potential hasanabsoluteminimumatθ=^A/fA=^{0.In fact,this} is thereason whyin thisextension of theStandard Model there isnostrong CPproblem[3].Thecurvature aroundthisminimum determinestheaxionmassm_Aatfinitetemperature,

m²_A

(

T

) ≡ ∂

²V

(

A

,

T

)

∂

A²

|

A=⁰

= χ (

T

)

f_A²

,

(3)

in terms of the topologicalsusceptibility, i.e. the variance of the θ=0 topologicalchargedistribution,

http://dx.doi.org/10.1016/j.physletb.2015.11.020

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

χ (

T

) ≡

d⁴x

^q

(

x

)

q

(

0

)

T

|

θ=0

=

^lim

V→∞

^Q2

T

|

θ=0

V

,

(4)

whereQ≡

d⁴xq(x).Similarly,self-interactiontermsinthepoten- tial,e.g.theA⁴ termoccurringintheexpansionofthefreeenergy densityaroundθ=^A/fA=^0,

V

(

A

,

T

) =

¹ 2

χ (

T

)θ

²

1

+

^b2

(

T

)θ

²

+ . . .

|

_θ=A/f_A

,

(5) aredeterminedbyhighermoments,

b2

(

T

) = −

^Q4

T

−

³

^Q2

²_T

12

Q²

T

|

θ=⁰

.

(6)

These non-perturbative quantities enter in a prediction of the lower bound on the fractional contribution of axions to the ob- servedcolddarkmatterasfollows¹

R_A

¹⁰⁵

θ

²

^f

(θ

²

)

^GeV3 g_∗S

(

T_osc

)

T_osc³

χ (

T_osc

)

χ (

0

) ,

(7) whereθ^{2isthevariance ofthespatial} distributionoftheinitial values ofthe axion field A/fA before the formation of the dark mattercondensatebythemisalignmentmechanism,whichoccurs, when the Hubble expansion rate gets of the order of the axion mass,mA(Tosc)=^3H(Tosc),i.e.atatemperature

Tosc

^{50 GeV} g¹_∗^/_R⁴

(

T_osc

)

mA

μeV

1/2

χ (

Tosc

) χ (

0

)

1/4

.

(8)

Thefunction f(θ²)istakinginto accountanharmonicityeffects arisingfromtheself-interactiontermsinV(A,T)anddependson thespecificformofthepotential.Thefunctionsg_∗R(T)andg_∗S(T) denotetheeffectivenumberofrelativisticenergyandentropyde- greesoffreedom,respectively.

What is urgently needed for axion cosmology is thus a pre- cisedeterminationofthetopologicalsusceptibilityandhighermo- mentsofthetopologicalchargedistribution.Inthiscontext,most predictions have been entirely based on the semi-classical ex- pansionofthe Euclideanpath integral offinite temperatureQCD aroundadilutegasofinstantons–finiteactionminimaoftheEu- clideanactionwithunittopologicalcharge–seee.g.Ref.[8]foran earlyexhaustingstudyandRef.[9]forarecentupdateconcerning thequarkmasses.Acomparativestudyofthesepredictionsbased onthediluteinstantongasapproximation(DIGA)hasbeencarried outinRef.[10],wherealsoananalysisintermsofaphenomeno- logical instanton liquid model(IILM) [11] is presented. However, up to now, in all DIGA investigations of the topological suscep- tibilityonlytheone-loop expressionintheexpansion aroundthe instantonbackgroundfieldwasused.Thisresultsinastrongrenor- malization scale dependence and thus large uncertainties which wereneglected intheprevious DIGAbasedpredictions.Infact,in thetemperaturerangeTosc∼^{GeV ofinterest,oneexpectsalarge} uncertaintyin theoverall normalization dueto theneglection of higherorderloopeffects,sinceinthisregion

α

s(T_osc)isnotsmall.

We will exploit in this letter both the one-loop DIGA result as wellasitstwo-looprenormalizationgroupimproved(RGI)version inordertostudythetheoreticaluncertainties arisingfromhigher loopcorrections. Mostimportantly,wecomparethesepredictions

1 Thisisarewritingofequation(2.10)ofRef.[7]wherewehaveusedχ(0)= (mAfA)²3.6×10⁻⁵GeV⁴ from the chiralLagrangian toexpress fA interms ofmA,thezerotemperaturemass,thatisitselfafunctionofTosc andtheratio χ(T)/χ(0)=(mA(T)/mA)²throughtheconditionfortheonsetoftheoscillations, mA(Tosc)=^3H(Tosc)withH²=⁸π³GNg_∗R(T)T⁴/90 theHubbleexpansionrate.

Fig. 1.Demonstrationthatvolumeissufficientlylargetohavenegligiblefinitevol- umecorrectionsonQ².ThedataforT/Tc=2 ismultipliedby3forbettervisibility ofthecomparison.

with the outcome of our lattice based fullynon-perturbative re- sults.

Actually, there have been a number of lattice calculations of

χ

(T)andb2(T)attemperaturesbeloworslightlyabovetheQCD phase transition, mostly inquenched QCD,see e.g. [12–16].Here we go beyond thoselattice calculations. We extend theavailable temperaturerangeandcarry outacontrolled continuumextrapo- lationforthisextendedrange.Inaddition,notethattherehasbeen noquantitativeinvestigationwhetherandwherethelatticeresults turn intotheDIGA results.Wewill presentournewhigh-quality latticedataforthetopologicalsusceptibilityinquenchedQCD(i.e.

neglecting the effectsof lightquarks) andcompare themquanti- tatively totheDIGA resultspecializedto thecaseofn_f =^{0 light} quarks.

2. Axionpotentialcoefficientsfromthelattice

Onthelattice,thetopologicalsusceptibilityismeasuredonthe torusasthesecondmomentofthedistributionoftheglobaltopo- logicalcharge

χ

t

=

^Q2

/

V

,

where Q isanyrenormalizeddiscretizationoftheglobaltopolog- ical charge, andV ^{is the} four-volume of the lattice. There are a lot of different fermionic and gluonic definitions of Q available.

We chooseagluonicdefinitionbasedonthe Wilsonflow[17,18], which hasthe correctcontinuum limit similarly tothe fermionic definitionsbutisnumericallyalotcheaper.Inparticularweevolve ourgaugefieldconfigurationswiththeWilsonplaquetteactionto a flow time t anddefine the globaltopologicalcharge as thein- tegraloverthe cloverdefinitionofthetopologicalcharge density.

Thisdefinitiongivesaproperlyrenormalizedobservablewhenthe flowtimetisfixedinphysicalunits.

Weuseatree-levelSymanzikimprovedgaugeaction.Ourtem- peratures range frombelow Tc up to 4Tc.Here Tc is thecritical temperature,whichisthequantityusedforscalesetting.Thecrit- icaltemperaturesfordifferentlatticespacingsweredeterminedin earlierwork[19,20].Forthewholetemperaturerangewekeepthe spatial lattice size approximately at L=²/T_c. We checked with dedicatedhighvolume runsat1.5Tc and2Tc that thisvolumeis sufficiently large to have negligible finite volume corrections on Q².ThisisshowninFig. 1.OurspatialgeometryisL×^L×^2Lto enable testsof subvolumemethods whichwill be reportedsepa- rately–hereweonlyusethefullvolume.Foralltemperatureswe

Table 1

Listofsimulationpoints:temperatures,latticesizesNt=¹/(aT),Ns=^L/a,number ofsweepsNsweeps,andestimatedintegratedautocorrelationtimetint,Q aregiven.

T/Tc Nt Ns Nsweeps tint,Q

0.9 12 5 32 K 3.3(2)

12 6 48 K 7.1(4)

16 8 170 K 36.4(32)

1.0 12 5 48 K 4.6(2)

12 6 64 K 12.6(9)

16 8 180 K 53.9(52)

1.1 12 5 48 K 5.3(3)

16 6 160 K 12.7(8)

20 8 330 K 65.3(71)

1.3 16 5 64 K 5.7(3)

16 6 220 K 15.4(9)

24 8 550 K 70.2(87)

1.5 16 5 96 K 5.8(2)

20 6 210 K 14.5(10)

24 8 660 K 63.8(75)

1.7 20 5 420 K 6.9(2)

20 6 1300 K 18.9(6)

28 8 8200 K 88.1(42)

2.0 20 5 440 K 7.5(2)

24 6 1900 K 18.0(6)

32 8 8400 K 90.5(51)

2.3 24 5 740 K 8.0(3)

28 6 2400 K 18.7(5)

36 8 8500 K 60.4(31)

2.6 28 5 960 K 6.9(2)

32 6 3500 K 17.8(5)

44 8 8000 K 51.7(35)

3.0 32 5 1500 K 6.2(2)

36 6 5700 K 16.0(5)

48 8 11 000 K 54.3(46)

3.5 36 5 2200 K 5.3(1)

44 6 5200 K 15.6(6)

56 8 12 000 K <45.0

4.0 40 5 2600 K <5.0

48 6 5900 K 15.6(8)

64 8 12 000 K <45.0

havethreelattice spacings(aT)⁻¹=⁵,6,8 to beable toperform an independent continuum extrapolation for every temperature.

The local heatbath/overrelaxation algorithm is used for the up- date,onesweepconsistsof1heatbathand4overrelaxationsteps.

Wefound that theautocorrelationtime of thetopological charge depends weakly on (aT), i.e. if the temperature is increased by decreasingthelattice spacing. The numberofupdate sweeps be- tweenmeasurementswaschosen inaccordancewiththeautocor- relationtime.Table 1liststhesimulationpointswiththenumber ofsweeps.

We integrated the Wilson flow numerically to a maximum flow-timeofabout8t≈¹/(2T_c²) forall temperatures.Fig. 2gives thedependenceofthesusceptibilityontheflowtimeforT=^2Tc. While in the continuum limit the result is independent of the choice of the flow time t, different choices have very different lattice artefacts. For small flow times the different lattice spac- ingsgiveverydifferentresults.Forlargerflowtimestheexpected plateaubehaviorcan beobservedforeach latticespacingandthe latticeartefacts alsodecreasesignificantly. Thechoice oftheflow time brings in some arbitrariness into the analysis, however the continuumresultshouldnotdependonthischoiceoncet isfixed inphysicalunits.Butthisiscertainlyasubleadingsourceoferror comparedtothestatisticalerrorduetotheraretopologytunneling

Fig. 2.DemonstrationthatWilsonflow/latticerenormalizationareundercontroland thatthedependenceoftheresultsonthechoiceoftheflowtimedecreaseswhen thecontinuumlimitisapproached.

Fig. 3.LatticedataonthetopologicalsusceptibilityatNt=5,6,8 andlatticecon- tinuumextrapolationtogetherwithfitofsimplepowerlaw.

eventsathightemperatures.Inthisanalysiswechooseatemper- aturedependentflowtimefortheevaluationof Q²as

8t

=

1

/(

1

.

5Tc

)

²

,

T

<

1

.

5

1

/

T²

,

T

≥

¹

.

5

.

(9)

For low/high temperatures this means a temperature indepen- dent/dependent flow time. This choice is safely in the expected plateauregionforalltemperatures.

TheresultingvaluesforthesusceptibilityareplottedinFig. 3.

Thisplotalsogivestheresultofaglobalcontinuum extrapolation usingasetoftemperaturesandthe6-parameterpowerlawansatz

χ

_t

₌ ( χ

₀

₊ χ

₀a²

)

^T T₀

+

^T₀^a2

b+^ba2

,

(10)

where

χ

0,T₀,andbarefitparametersgivingthecontinuumlimit.

χ

₀,T₀,andbarefitparametersdescribingthedeviationfromthe continuumlimit.Thepowerlawformofthefitismotivatedbythe expectedhightemperaturebehavior ofthesusceptibility.Fortem- peratures closeto T_c thisisonlyan empiricalfitwhichseemsto describe thelattice resultsquitewell. Thefit parameter T0 isin- cludedasaconsistencycheckandshouldgive1inunitsofTc.This issatisfiedbythefitresult.Thevariationbetweendifferentchoices forthestartingtemperatureofthefitrangeT_min/T_c=¹.3,1.5,1.7

Fig. 4.Latticedataontheanharmonicitycoefficientb2oftheaxionpotentialcom- paredtoitsDIGAprediction.Thedatapointsareshiftedabithorizontallyforbetter visibility.

givesanestimateofthesystematicerroroftheresult.Thebestfit parametersare

χ

₀

₌

0

.

11

(

2

)(

1

),

b

= −

⁷

.

1

(

4

)(

2

),

T₀

=

¹

.

02

(

4

)(

2

),

(11) where thefirst error is the statistical, the second is the system- atic. The point-wise continuum extrapolation is consistent with theglobalfit,evidentlythelatterhassmallererrorsforlargetem- peratures.Notethat though acontrolled continuumextrapolation ispossible usingthree lattice spacings,estimatingthe systematic uncertaintyofthisextrapolationwouldrequireatleastonemore latticeresolution.

Inarecentanalysis[15] thetopologicalsusceptibilitywas cal- culated using thetechniques [21,22].The calculation was carried outattwotemporalextensions,correspondingtotwolatticespac- ingsateach temperature. The exponent b= −⁵.64(4) was found whichdiffersfromourvalue.Notehowever,thatourtemperature rangeislarger,thusweareclosertotheapplicabilityrangeofthe DIGA.Furthermore,thetwolatticespacingswerenotsufficientfor acontrolledcontinuumextrapolationthusuncertaintiesrelatedto thisfinalsteparenotincludedintheresultof[15].

Wehavealsodeterminedthesecondimportantcoefficientb2of theaxionpotential,characterizingitsanharmonicity,bymeasuring theobservable

b2,t

= −

^Q4

−

³

^Q2

² 12

^Q2

.

TheresultisplottedinFig. 4.

3. ComparisonbetweenlatticeandDIGAresults

Inthissection,weconfrontthelatticeresultswiththeonesob- tainedfromtheDIGAframework.Forthesakeofcompletenesslet uscollectfirsttheavailableformulasforthelatter[23–30].Atvery hightemperatures, far above the QCD phase transition, it makes senseto infertheθ dependenceofQCDfromtheDIGA, inwhich thepartitionfunction,foranyn_f,canbewrittenas[24],

Z

(θ,

T

)

nI,n_¯I

1 n_I

!

ⁿ¯I

!

^Z

n_I+ⁿ¯I

I

(

T

)

exp

i

θ (

n_I

−

ⁿ¯I

)

,

(12)

where Z_I= ^Z¯I isthe contribution arising from the expansion of thepathintegral aroundasingle instanton I (anti-instanton ¯I).It followsdirectlythatthepotentialhastheform

V

(

A

,

T

) χ (

T

) (

1

−

^cos

θ ) |

θ=A/fA

,

(13)

fromwhichoneinfers b2

(

T

) −

¹

12

.

(14)

Thiscanbeconfrontedright-awaywithourlatticeresults,cf.Fig. 4.

Similar toRef.[13]we findthat thepredictionfromtheDIGAfor b2 isreachedalreadyatsurprisinglylowvaluesofT/Tc^1.

The whole temperature dependence of the axion potential arisesintheDIGAthroughthetopologicalsusceptibility,whichin thiscaseisexplicitlygivenby

χ (

T

)

^ZI

(

T

) +

^Z¯I

(

T

)

V

=

²

∞

0

d

ρ

D

( ρ )

G

( πρ

T

),

(15)

in terms of the instanton size distribution at zero temperature, D(

ρ

),andafactorG(

πρ

T)takingintoaccount finitetemperature effects. Theformeris knowninthe frameworkofthe semiclassi- cal expansion around theinstanton forsmall

α

s(

μ

r)ln(

ρ μ

r)and

ρ

mi(

μ

r),where

α

s isthestrongcoupling,

μ

r istherenormaliza- tion scaleandm_i(

μ

r)are therunningquark masses.Toone-loop accuracy,itisgivenby²

D

( ρ ) ≡

^dMS

ρ

⁵

2

π α

_MS

( μ

_r

)

6

exp

−

²

π α

_MS

( μ

_r

)

(16)

× ( ρ μ

r

)

^β0

1

+

O

( α

_MS

( μ

r

)) ,

with d_MS

=

^e5/6

π

^{2 e}

−⁴.534122

; β

0

=

11

.

(17)

Atfinitetemperature,electricDebyescreeningprohibitstheex- istenceoflarge-scale coherentfieldsintheplasma,leading tothe factor[26,27],

G

(

x

) ≡

^exp

−

^2x2

−

^18A

(

x

)

,

(18)

with A

(

x

) −

¹

12ln

1

+ ( πρ

T

)

²

/

3

+ α

1

+ γ ( πρ

T

)

^−3/2

₋8

,

(19)

and

α

₌0.01289764 and

γ

₌0.15858, in Eq. (15). This factor cuts off the integration over the size distribution in Eq. (15) at x=

πρ

T∼^{1 andensuresthevalidityoftheDIGAatlargetemper-} atures,atwhich

α

s(

π

T)^1.

Collecting all the factors, the topological susceptibility, in the one-loopDIGA,reads

χ (

T

)

^2d_MS

( π

T

)

⁴

μ

r

π

T

11

I 2

π α

_MS

( μ

r

)

6

×

^exp

−

²

π

α

_MS

( μ

r

)

¹

+

O

( α

_MS

( μ

r

))

,

(20)

with I

=

∞ 0

dx x⁶G

(

x

) =

⁰

.

267271

.

(21)

This result, however, still suffers from a sizeable dependence on therenormalizationscale

μ

r,reflecting theimportance ofthe neglected two-loop and higher order contributions. In fact, it is

2 ForquenchedQCDthenumberoflightquarksisnf=0;thegeneralformula canbefoundexplicitlyine.g.Ref.[30],whichcontainsapioneeringconfrontation ofcooledlatticedataonD(ρ)withthetwo-loopRGimprovedDIGAresult.

Fig. 5.PredictionofthetopologicalsusceptibilityintheDIGA:comparisonbetween one-loopandtwo-loopRGIresults.Weusedthefour-loopexpressionfortherun- ningcouplinginthemodifiedminimalsubtractionschemeasgivenintheappendix ofRef.[31]andthecentralvalueofTc/⁽_MS^nf=0)=1.26(7)asdeterminedfromthe latticeinRef.[20].

Table 2

Temperatureslopesofthetopologicalsusceptibilitypredictedinthetwo-loopRGI DIGA,forarangeofrenormalizationscalesaccordingtoEq.(24).

T/Tc 1.5 2 3 4 5

b(κ=⁰.6) −⁶.04 −⁶.26 −⁶.43 −⁶.50 −⁶.55 b(κ=¹⁾ −⁶.37 −⁶.46 −⁶.55 −⁶.59 −⁶.62 b(κ=2) −6.55 −6.59 −6.64 −6.67 −6.69

tamedbytakingintoaccount theultravioletpartofthe two-loop correction.ThelatterhasbeencalculatedinRef.[29]andshownto haveexactlytheformthatthegauge couplingbecomesa param- eterrunningaccordingto therenormalizationgroup (RG). There- fore,theultimate, allorderresultforthetopologicalsusceptibility becomesindependentof

μ

r,for

μ

r→ ∞^{.Attwoloop,thecorrec-} tionsamounttoafactor

( ρμ

r

)

^{(β1−12β0)αMS(μr)/(4π)}

; β

1

=

102

,

(22) inD(

ρ

).Therefore,thisRGimprovementcanbetakenintoaccount byreplacingthefactorI inEq.(20)by

˜

I

= μ

r

π

T

₋30α_MS(μr)/(4π)

∞ 0

dx x^{6−30αMS(μr)/(4π)}G

(

x

).

(23)

Infact, exploiting the two-loop RG improvement,the

μ

r depen- denceisheavilyreduced,asisobviousfromFig. 5,wherewehave usedasanaturalrenormalizationscale

μ

r

= κ / ρ

m

= κπ

T

/

1

.

2

,

(24) with

ρ

mbeingapproximatelythemaximumoftheintegrandof˜I, andvariedtheremainingfreeparameter

κ

between0.6and2.The renormalizationscaledependenceappearstobehighlyreducedin theregimes>3T_c and<T_c.However,inan intermediateregion,

∼^Tc–2Tc,itiscomparableinsizetotheoneatone-loop.

WepresentinTable 2thepower-lawbehaviorpredictedbythe two-loop RGI DIGA at various temperatures, which can be com- paredtothefit(11)tothecontinuum latticeresult.Asfarasthe overall normalizationofthe DIGA result for

χ

is concerned, one stillexpectsalargeuncertaintyinthetemperaturerangeavailable from the lattice. In fact, at these temperatures,

α

s is not small, seeTable 3.Apartfromtheultravioletpart,there willbe afinite partofthetwo-loopcorrectionwhichwillaffectmainlytheover- all normalizationof

χ

and willdepend on the temperatureonly

Table 3

Strongcouplingconstantatμr=¹/ρm forthetemperaturerangecoveredbythe lattice.

T/Tc 1 2 3 4 5

α_MS(1/ρm) 0.36 0.23 0.19 0.17 0.16

Fig. 6.Rescaledone-loopandtwo-loopRGIDIGAresultscomparedtolatticecontin- uumextrapolation.TheDIGAresultsshowninFig. 5werescaledbyafactor Kof ordertensuchthattheycoincideatT/Tc=^{2 withthecentralvalueofthelattice} continuumdata.

logarithmically. Unfortunately, this finite part is not known, yet.

Therefore,when comparingtothe continuumextrapolatedlattice results,weallowa multiplicativefactor K to accountforthisun- certainty,i.e.weabsorbtheremaininghigherloopuncertaintiesby replacing

1

+

O

( α

_MS

( μ

r

= κπ

T

/

1

.

2

))

→

^K

(

T

/

T_c

)

(25) inEq.(20)andthecorrespondingtwo-loopRGIexpression.Clearly, the K-factorshouldapproachunityatverylarge T/Tc.

Fig. 6 nicely illustrates the agreement between the DIGA and thelatticeresult,ifa K-factorofordertenisincluded.³Morepre- cisely, fitting the lattice continuum data withthe rescaled DIGA expressioninthetemperaturerangeT/Tc≥^1,onefinds

K

=

⁸

.

9

±

⁰

.

7

,

at 95% CL

,

(26)

whileafitinthetemperaturerangeT/Tc≥^{2 yields}

K

=

⁷

.

9

±

³

.

3

,

at 95% CL

.

(27)

Apparently,inthetemperaturerangeaccessibletothelattice, the higherordercorrectionstothepre-factoroftheDIGAarestillap- preciable, but there are indications of a trend that the K-factor getssmaller,asexpected,towardslargervaluesofT/Tc.

The K-factorstrongly dependsonthe valueof T_c//^(nf=0)

MS : it reduces toone for Tc//^(nf=0)

MS ¹.03.However, the lattervalue isabout3sigmabelowthecentralvalue determinedinRef.[20], Tc//⁽_MS^nf=0)1.26(7).

4. Conclusions

ThispaperpresentslatticeandDIGAcalculationsoftheθ and temperaturedependenceofthefreeenergydensityofQCDinthe

3 Kfactorsofordertentofiftyarenotuncommonevenatnext-to-leadingorder inordinaryperturbativeQCD,seee.g.Ref.[32].

Fig. 7.Consequenceofourfindingsforaxiondarkmatterfromquenched(top)and fullQCD(bottom).ThedarkyellowregionisexcludedbecauseRA>1 justfrom themisalignmentmechanism,cf.Eq.(7),assumingaflatdistributionofinitialmis- alignmentanglesintheobservableuniverseθ∈ [−π,π]^{.Thelightyellowregion} indicates RA>1 whenstring contributions accordingtoRef.[33] areincluded.

Inthedarkgreenregionaxionsevenincludingstringsgiveasmallcontribution (RA<0.1)todarkmatter.ThelightgreenregionindicatesRA<0.1 justfromthe misalignmentmechanism.Inordertocomparethequenchedresultswiththeaxion darkmatterscenariowehadtotransformdimensionlessquantitiesintodimension- fulones.Forillustrativepurposesweuse, inthe quenchedcase, Tc=^{294 MeV.}

Ourquenchedlatticeresultsareshownbythebluepointsandthetwo-loopRGI DIGApredictionbythegraybandinthetoppanel.SincethetransitioninfullQCD isacrossover [34],Tc isnotunambiguouseveninthe unquenchedcase. Using thetransitiontemperaturedefinedbythechiralsusceptibility,Tc=147(2)(3)MeV [35–37], weneedK=9.22±0.6 inorderthatthetwo-loop RGIDIGAmatches χ(Tc)/χ(0)=^{1.Theblueandgraybands inthe lowerpanelcorrespondtothe} two-loopRGIDIGApredictionswithK=^{1 andK}=⁹.22±⁰.6,respectively.The dashedredlineshowstheIILMpredictionofRef.[11].Thedashedblacklinescor- respondtofixedaxionmasses(inunitsofμeV).UsingK=⁹.22 andagainassuming atleasttenpercentaxioncontributiontodarkmatterwecanreadofffromthe lowerpaneltherange40 μeVmA930 μeV,whileusingK=1 wewouldget 50 μeVmA1100 μeV,i.e.onlya∼20% correctionforanO(10)Kfactoruncer- tainty.(Forinterpretationofthereferencestocolorinthisfigurelegend,thereader isreferredtothewebversionofthisarticle.)

quenched limit, i.e. for infinite quark masses. In the lattice ap- proach the temperatureranges from 0.9Tc to 4Tc and thus sig- nificantlyextendsformerresults.Thecomparisonofthequenched continuum extrapolatedlattice dataandtheDIGAframework has ledtotwohighlynon-trivialfindings.Thefirstoneistheobserva- tionthat theexponentofthetemperaturedependenceof

χ

(T)is correctlydescribed bytheDIGA, evenallthewaydownto T_c (at leastwithinourerrorbars),cf.Fig. 7(top).Thesecondnon-trivial findingisabouttheprefactorofthistemperaturedependence:the one-loop based DIGA prefactor is off by a large factor of about ten(wecallthisa K-factor,afamiliarnotioninperturbativeQCD predictions). This is expected, since in the considered tempera- ture range the strong coupling is still sizeable. However, such a K-factoruncertaintyhasnotbeenaccountedforinanyofthepre- vious studiesoftheaxion darkmatterabundance,e.g. Refs. [8,9].

Nevertheless,thereisapieceofphysicsinformation,whichessen- tiallyfixesthisK-factor.Asexpected,thetopologicalsusceptibility changesonlyalittlefromT=^{0 toT}=^Tcanddropsfurtheronfor T>Tc.Actually, thequencheddata isvery well describedif one takes

χ

(T)/

χ

(0)=^{1 atT}/T_c∼^{1 inordertofixtheK-factor4and} thenusestheexponentoftheDIGAframeworkfortherestofthe temperaturedependenceinthehightemperatureplasmaphase,cf.

Fig. 7(top).

Unquenchedlatticesimulationswithphysicalquarkmassesare abouttwo-to-three ordersofmagnitudemoreCPUintensivethan quenchedones. Inaddition oneexpectsmuch smallertopological susceptibilities andlarger cutoffeffects.Ifonecombines allthese costfactorsoneendsupwithaveryCPUdemandingproject.Thus, it is of extremeimportance to give an estimate ofthe ranges in temperatureandtopologicalsusceptibilityofinterestforaxioncos- mology.

This canbe done by exploitingthe two highly non-trivialob- servations we discussed above: we set

χ

(T_c)/

χ

(0)=^{1 and use} the exponentsuggestedby theDIGA framework,which isreadily available alsoforfullQCD, forthe temperaturedependence.⁵ Our finding isshownin thelower panel ofFig. 7,wherealsothe re- sultwiththe DIGAprefactorandthat oftheinteractinginstanton liquidmodelaredepicted. Fortunately, theestimatedvalue ofthe K-factordoesnotaffectverystronglytheextractionofthevalues oftheaxionmassofinterestforaxiondarkmatter,becauseofthe steepness ofthefalloffof

χ

(T),cf.Fig. 7(bottom).

Therearetwoimportantmessages,whichdominantlyinfluence anyfuturelatticestudy.Thefirstoneisthatoneneedssomewhat larger T/Tc values than inthe quenchedcase making thecalcu- lationquiteCPUdemanding.Theother importantmessageisthat using dynamical QCD one has to determine about one order of magnitudesmallervaluesof

χ

(T)/

χ

(0)thaninthequenchedcase.

Since already

χ

(0)isfarlessintheunquenchedcasethaninthe quenched one,thesmallnessoftheneeded

χ

(T)/

χ

(0)makesthe task even more CPUdemanding. The planning ofany future un- quenchedprojectshouldconsidertheestimateslistedabove.

Acknowledgements

The authors wish to thank G. Moore for useful discussions.

Computations were performed on JUQUEEN at FZ-Jülich and on GPUclustersatWuppertalandBudapest.Thisprojectwas funded bytheDFGgrantSFB/TR55.byOTKAundergrantOTKA-NF-104034

4 WeobtainK=¹¹.7⁺₋⁵₃^._.⁵₁bytakingχ(Tc)/χ(0)=^{1,wheretheerrorscomefrom} varyingκ inthe(0.5,2)interval,and K=⁷.9±³.3,bytakingχ(T)/χ(0)=^{1 at} T/Tc=0.92,forκ=1.

5 Qualitatively,boththeapproximateconstancyofχ(T)/χ(0)for T/Tc1 as wellastherapidfalloffforT/Tc^{1 havebeenobservedalreadyinpioneering} latticestudiesinfullQCD[12,14].