This file generates tables of the particle number densities of 2+1 flavor 
Nc=3 quark + electron, electrically neutral matter to ((Subscript[α, s])^2)needed for the thermodynamics. The file takes as its input the value of (Subscript[Λ, (OverBar[MS])]) (LambdaMSbar)  and (Subscript[m, s])(ms0)in MSbar-scheme with renormalization scale at 2 GeV. Also,the user must specify the range of strange quark chemical potentials for which the particle number densities are evaluated. The output, Ntable, is of the form
Ntable[[n]]={(Subscript[μ, s]),(Subscript[n, s]),(Subscript[n, u]),(Subscript[n, d]),(Subscript[n, e])}.  All quantities in units  of  GeV.

Set parameters :

LambdaMSbar=0.378;
ms0=.1;
Δ=0.0;

(OverBar[Λ])=X μ:

X=4;

Range of strange quark chemical potentials to tabulate :

muMin=0.27491;
muMax=2.5;
dmu=0.005;

ln[Z]

These must not be changed unless you know what you are doing. Nc = 3 hardwired in plasmon sum.

CA=3;
Nc=3;
dA=Nc^2-1;
CF=(Nc^2-1)/(2 Nc);
Nf=3;

ln Z^m = 0

lnZm0[vecmu_,g2_,barLhat_]:=1/4/Pi^2Sum[vecmu[[i]]^4(Nc/3-dA g2/(4 Pi)^2+dA (g2/(4 Pi)^2)^2(4/3(Nf-11 CA/2)Log[barLhat/2/vecmu[[i]]]-142/9 CA+17/2 CF+22/9 Nf)),{i,1,Length[vecmu]}]

ln Z^xVM

firstintapr2[t_]:=3 t^3/(1-t)^2(Log[t]-1)(0.83 t+0.06 t(1-t)-0.056 t (1-t)^2+Log[1-t](1.005-0.272 t (1-t)+0.154 t (1-t)^2))

lnZxVM[vecmu_,g2_,barLhat_,mhat_]:=-dA g2^2/(4 Pi)^2 mhat^4/12/Pi^4 Sum[firstintapr2[Simplify[1/(1+mhat/vecmu[[i]])]],{i,1,Length[vecmu]}]

ln Z^m

ubar[x_]:=(Sqrt[1-x^2])
z[x_]:=ubar[x]-x^2Log[(1+ubar[x])/x]

δa5[x_]:=+  0.0123069 x^2+0.00788621 x^2Log[x]+0.00165308 x^2 (Log[x])^2 - 0.0125499 x^4+0.0175865 x^4Log[x]-0.0156622 x^4(Log[x])^2   +  0.00497165 x^4 (Log[x])^3-0.00251205 x^4 (Log[x])^4  ;
δa7[x_]:=0.00157785 x^2+0.000685771 x^2Log[x]-0.0000536257 x^2 (Log[x])^2   - 0.00157794 x^4+0.00244162 x^4Log[x]+0.00291526 x^4 (Log[x])^2  +  0.00128126 x^4 (Log[x])^3+0.00153878 x^4 (Log[x])^4 ;

G1[x_]:=-58.06886100277347` (x^2)+63.5119057094417` (x^4)-121.57641998626995` (x^4) Log[x]+80.44991984500753` (x^4) ((Log[x])^2)-98.27970397499895` (x^4) ((Log[x])^3)+35.8676186564003` (x^4) ((Log[x])^4)
G2[x_]:=-622.9027525956599` (x^2)+619.5515464296311` (x^4)-149.52930985837634` (x^2) Log[x]-1112.4343595202872` (x^4) Log[x]-1.1006583345905372`*^-17 (x^2) ((Log[x])^2)+948.4166421450439` (x^4) ((Log[x])^2)-502.204722548742` (x^4) ((Log[x])^3)+305.1762641712499` (x^4) ((Log[x])^4)
G3[x_]:=(1/24) ((71+8 (π^2)) (x^2) (Sqrt[1-(x^2)])+768 (π^4) δa5[x]-(71+8 (π^2)) (x^4) Log[((1+(Sqrt[1-(x^2)]))/x)])
G4[x_]:=-(-3+(π^2)) (x^2) (Sqrt[1-(x^2)])+32 (π^4) δa7[x]+(-3+(π^2)) (x^4) Log[((1+(Sqrt[1-(x^2)]))/x)]

F1[x_]:=CA/(24π^2)(2 ubar[x]^3-3x^2 z[x])
F2a[x_]:=dA/(4π^2)(2  ubar[x]^4-4x^2z[x]-3 z[x]^2)
F2b[x_]:=-(3 dA x^2 z[x])/(2π^2)
F3a[x_]:=dA/(2π^2)( CA(-11/3Log[x/2]-71/9+G1[x])+CF(17/4+G2[x])+Nf(2/3Log[x/2]+11/9 +G3[x])+G4[x])
F3b[x_]:=dA/(6π^2)(CA(22 ubar[x]^4-185/2 x^2 z[x]-33 z[x]^2)+CF 9/2(16 x^2 ubar[x](1-ubar[x])-3(7 x^2-8ubar[x])z[x]-24 z[x]^2)-Nf(4 ubar[x]^4-13 x^2 z[x]-6 z[x]^2))
F3c[x_]:=-dA x^2/2/Pi^2((11CA-2 Nf)z[x]+18 CF(2z[x]-ubar[x]));

lnZm[g2_,barLhat_,mhat_]:=F1[mhat]+g2/(4 Pi)^2 (F2a[mhat]+F2b[mhat] Log[barLhat/mhat])+(g2/(4 Pi)^2)^2(F3a[mhat]+F3b[mhat] Log[barLhat/mhat]+F3c[mhat]Log[barLhat/mhat]^2)

Plot[lnZm[0.5,3,mhat],{mhat,0.001,0.99}]

lnZ^ring

m0Zring[vecmu_,g2_,barLhat_]:=-dA/4/Pi^2 (vecmu.vecmu)^2 (g2/(4Pi)^2)^2 (2 Log[g2/(4 Pi)^2]-22/3+16/3Log[2](1-Log[2])+AVdelta+2Pi^2/3+I15[vecmu]/(vecmu.vecmu)^2)

AVdelta=-0.856383209;

I13[mhat_]:=0.818274 ubar[mhat]^3+0.31845 ubar[mhat]^6-0.136725 ubar[mhat]^7

I14[x_]:=1.991188165230242` (1-x)-0.9912208629031928` (((1-x))^2)-0.26942499045800095` (1-x) Log[1-x]+0.26026098267857506` (((1-x))^2) Log[1-x]-0.11855128965430835` (((1-x))^2) ((Log[1-x])^2)

I15[vecmu_]:=16/3 Log[2] Sum[vecmu[[i]]^4,{i,1,Length[vecmu]}]-2 vecmu.vecmu Sum[vecmu[[i]]^2 Log[vecmu[[i]]^2/(vecmu.vecmu)],{i,1,Length[vecmu]}]+2/3 Sum[(vecmu[[i]]-vecmu[[j]])^4 1/2Log[((vecmu[[i]]^2-vecmu[[j]]^2)/vecmu[[i]]/vecmu[[j]])^2]+4 vecmu[[i]] vecmu[[j]](vecmu[[i]]^2+vecmu[[j]]^2) Log[(vecmu[[i]]+vecmu[[j]])^2/vecmu[[i]]/vecmu[[j]]]-(vecmu[[i]]^4-vecmu[[j]]^4)Log[vecmu[[i]]/vecmu[[j]]],{i,1,Length[vecmu]},{j,i+1,Length[vecmu]}]

I16help1[x_]:=4.6291936663850946` (1-x)-3.6292695238575305` (((1-x))^2)+0.6064386349232564` (1-x) Log[1-x]+2.0213406785306756` (((1-x))^2) Log[1-x]-0.5658832917823072` (((1-x))^2) ((Log[1-x])^2)

I16[mhat_,vec2_]:=AVdelta+2Log[(vec2^2+2 vec2 I14[mhat]+I14[mhat] I16help1[mhat])/(vec2^2+vec2 I14[mhat])]

I17input1[x_]:=4.8037832033625705` (1-x)-2.2242383349496593` (((1-x))^2)-0.41145665291433065` (1-x) Log[1-x]+0.05618460508642291` (((1-x))^2) Log[1-x]-0.07180597699877932` (((1-x))^2) ((Log[1-x])^2)

I17input2[x_]:=-8.574463854093054` (1-x)+3.0906075080265083` (Sqrt[1-(x^2)])+9.654355344233919` (1-(x^2))+10.861457739326454` (((1-(x^2)))^(3/2))-12.445046357366184` (((1-(x^2)))^2)+0.37400274001840395` (Sqrt[1-(x^2)]) Log[1-x]+8.731314784103287` (((1-(x^2)))^(3/2)) Log[1-x]

I17ai[x_]:=-7.543846397415098` (1-x)+1.1421668724164455` (Sqrt[1-(x^2)])+6.638741765598661` (((1-(x^2)))^(5/2))-2.0053887442078895` (1-x) Log[1-x]+5.865202661442799` (((1-x))^2) ((Log[1-x])^2)-10.62161589244306` (((1-x))^3) ((Log[1-x])^3)-8.880917939218795` (((1-(x^2)))^(3/2)) Log[1+x]

I17bi[x_]:=-8.085131197844229` (((1-x))^3)+6.646189419214665` (((1-x))^5)-0.5845813413276937` (((1-x))^12)+0.3465094002382129` x (Sqrt[1-(x^2)])+1.3593707432490079` (1-(x^2))-3.770523846688646` (((1-(x^2)))^(3/2))+5.395415821785071` (((1-(x^2)))^2)-0.9525244276080395` (((1-(x^2)))^(5/2))

I17cll[mx_]:=I17input1[mx]-7/3I14[mx]-I17ai[mx]

I17chh[mx_]:=I17input2[mx]-7/3I14[mx]-I17ai[mx]

I17[x_,y_]:=((I14[x] ((((1-y))^4) Log[(((1-(1/y)))^2)]+(((1+y))^4) Log[(((1+(1/y)))^2)]-2 Log[(1/(y^2))]))/(6 (y^2)))+I17ai[x]+I17chh[x] y/(y+1)+I17cll[x]/(y^2+1)+I17bi[x] y/(3.38 Sqrt[1-x]+7.96(1-x)+y^3);

I18[x_]:=2.3003340724232086` (1-x)-1.3003455037065073` (((1-x))^2)-0.3476679926538811` (1-x) Log[1-x]+0.6484072905393345` (((1-x))^2) Log[1-x]-0.2409261667958931` (((1-x))^2) ((Log[1-x])^2)

I19help1[x_]:=-3.6643515105370845` (1-x)+2.808009035802549` (((1-x))^2)-0.7305488237960135` (1-x) Log[1-x]-1.2202673500975705` (((1-x))^2) Log[1-x]-0.04367660731823744` (1-x) ((Log[1-x])^2)+0.31153783190979534` (((1-x))^2) ((Log[1-x])^2)-0.02364776435743023` (((1-x))^2) ((Log[1-x])^3)

I19help2[x_]:=0.10416428010990529` (1-x)+0.4643661415965035` (((1-x))^2)+0.43156507764371566` (((1-(x^2)))^(5/2))-0.039399797727379865` (1-x) Log[1-x]-0.9889881831437727` (((1-x))^2) Log[1-x]+0.45508471067027356` (((1-x))^2) ((Log[1-x])^2)

I19[x_,y2_]:=I19help1[x]+I14[x](2 Log[y2]+2 Log[(y2^2+2 y2 I13[x]/I14[x]+I19help2[x])/(y2^2+y2 I13[x]/I14[x])])

I20[x_]:=3.843182933271279` (1-x)-7.9361743472477` (((1-x))^2)+8.062088562951423` (1-(x^2))+3.810951566866307` (1-x) Log[1-x]+9.703546789038324` (((1-x))^2) Log[1-x]-7.414762950214298` (((1-x))^2) ((Log[1-x])^2)

I21[x_]:=-0.11110744187523043` (1-x)+1.1103546415446395` (((1-x))^2)-0.059260849073043395` (1-x) Log[1-x]-2.0711122427813073` (((1-x))^2) Log[1-x]+0.39906137132795194` (((1-x))^2) ((Log[1-x])^2)

I22help1[x_]:=-0.8254720711551927` (((1-(x^2)))^2)-0.03084009925821547` (((1-(x^2)))^(7/2))+0.0016123815590526566` (1-x) Log[1-x]

I22help2[x_]:=1.3126254088011997` (Sqrt[1-(x^2)])+0.4343304030259942` (1-(x^2))+0.2530276354635296` (((1-(x^2)))^(5/2))

I22help3[x_]:=1.8442204296027336` (1-x)-0.8455831178924472` (((1-x))^2)+0.19443516859364646` (1-x) Log[1-x]-0.3920154346890649` (((1-x))^2) Log[1-x]-0.7036871885445087` (((1-x))^2) ((Log[1-x])^2)

I22[x_,y2_]:=I22help1[x]+I13[x](2 Log[y2]+2 Log[(y2^2+y2 I22help2[x]+I22help3[x])/(y2^2+y2 I22help2[x]/2)])

lnZring[vecmu_,g2_,barLhat_,mhat_]:=-dA/512/Pi^6 g2^2((vecmu.vecmu)^2(Log[g2/(4 Pi)^2/Sqrt[E]]+1/2(-19/3+2Pi^2/3 +I15[vecmu]/(vecmu.vecmu)^2+16/3 Log[2](1-Log[2])+I16[mhat,vecmu.vecmu]))+2 Sum[vecmu[[i]]^2 (I14[mhat]Log[g2/(4Pi)^2/Sqrt[E]]+1/2(I17[mhat,vecmu[[i]]]+16/3 Log[2](1-Log[2]) I18[mhat]+I19[mhat,vecmu.vecmu])),{i,1,Length[vecmu]}]+I13[mhat] Log[g2/(4Pi)^2/Sqrt[E]]+1/2(I20[mhat]+16/3Log[2](1-Log[2])I21[mhat]+I22[mhat,vecmu.vecmu]))

lnZring[{1.00001,2},1.5,1,0.5]

Plot[lnZring[{1.00001,2},1,1.5,mhat],{mhat,0,1}]

Sum of all

lnZtotal[vecmu_,g2_,barLhat_,mhat_]:=lnZring[vecmu,g2,barLhat,mhat]+lnZm[g2,barLhat,mhat]+lnZxVM[vecmu,g2,barLhat,mhat]+lnZm0[vecmu,g2,barLhat]

lnZtotal[{1.0000001,1.001},1,0.1,0.01]

Running alpha

beta0:=(11/3 CA-2/3 Nf)
beta1:=17/3 CA^2-CF Nf-5/3 CA Nf

Lfullbar[L_]:=2*Log[L/LambdaMSbar]

alphas0[barL_]:=4Pi/beta0/Lfullbar[barL]
alphas1[barL_]:=4Pi/beta0/Lfullbar[barL] (1-2beta1/beta0^2 Log[Lfullbar[barL]]/Lfullbar[barL])

asprime0[as_]:=-beta0 Pi/2 (as/Pi)^2
asprime1[as_]:=-beta0 Pi/2 (as/Pi)^2-beta1 Pi/4 (as/Pi)^3

Running ms

gamma0=3 CF;
gamma1=CF (97/6 CA+3/2 CF-5/3 Nf);
gamma2=CF (129/2 CF^2-129/4 CF CA+11413/108 CA^2+CF Nf (-23+24 Zeta[3])+CA Nf (-278/27-24 Zeta[3])-35/27 Nf^2);
beta0ver=1/4*beta0;
beta1ver=1/8*beta1;
beta2ver:=1/64*(2857/2-5033/18*Nf+325/54*Nf^2);
gamma0ver=1/4*gamma0;
gamma1ver=1/16*gamma1;
gamma2ver=1/64*gamma2;
A1=-beta1ver gamma0ver/beta0ver^2+gamma1ver/beta0ver;
A2=gamma0ver/beta0ver^2*(beta1ver^2/beta0ver-beta2ver)-beta1ver gamma1ver/beta0ver^2+gamma2ver/beta0ver;
Ialtmsrun1[barL_]:=ms0*(alphas0[barL]/alphas0[2])^(gamma0ver/beta0ver)*(1+A1 (alphas0[barL]/Pi))/(1+A1 (alphas0[2]/Pi))
Ialtmsrun2[barL_]:=ms0*(alphas1[barL]/alphas1[2])^(gamma0ver/beta0ver)*(1+A1 (alphas1[barL]/Pi)+(A1^2+A2) (alphas1[barL]/Pi)^2/2)/(1+A1 (alphas1[2]/Pi)+(A1^2+A2) (alphas1[2]/Pi)^2/2)

massive eos

Thermodynamics :

barLfunc[muS_,mu1_,mu2_,x_]:=x(muS+mu1+mu2)/3

momega0[mu_,vecmu_,as_,barL_,m_]:=-mu^4(1/4/Pi^2Sum[vecmu[[i]]^4(Nc/3),{i,1,Length[vecmu]}]+F1[m/mu])
momega1[mu_,vecmu_,as_,barL_,m_]:=-mu^4(1/4/Pi^2Sum[vecmu[[i]]^4(Nc/3-dA as/(4 Pi)),{i,1,Length[vecmu]}]+F1[m/mu])-mu^4as/(4 Pi) (F2a[m/mu]+F2b[m/mu] Log[barL/m])
momega2[mu_,vecmu_,as_,barL_,m_]:=-mu^4lnZtotal[vecmu,4 Pi as,barL/mu,m/mu]
omegaCFL[muS_,mu1_,mu2_]:=-3Δ^2 ((muS+mu1+mu2)/3)^2/π^2;

sn1[mu_,mu1_,mu2_,as_,barL_,m_,x_]:=Evaluate[-D[omegaCFL[mu1,mu2,mu]+momega1[mu,{mu1/mu,mu2/mu},as,barL,m],mu]]
sn2[mu_,mu1_,mu2_,as_,barL_,m_,x_]:=Evaluate[Evaluate[-D[omegaCFL[mu1,mu2,mu]+momega2[mu,{mu1/mu,mu2/mu},as,barL,m],mu]]]
ln1[mu_,mu1_,mu2_,as_,barL_,m_,x_]:=Evaluate[Evaluate[-D[omegaCFL[mu1,mu2,mu]+momega1[mu,{mu1/mu,mu2/mu},as,barL,m],mu1]]]
ln2[mu_,mu1_,mu2_,as_,barL_,m_,x_]:=Evaluate[Evaluate[-D[omegaCFL[mu1,mu2,mu]+momega2[mu,{mu1/mu,mu2/mu},as,barL,m],mu1]]]

1 - loop

sn1[mu_,mue_,x_]:=Max[Re[sn1[mu,mu,mu-mue,alphas0[barL],barL,Ialtmsrun1[barL],x]/.barL->barLfunc[mu,mu,mu-mue,x]],0]
dn1[mu_,mue_,x_]:=ln1[mu,mu,mu-mue,alphas0[barL],barL,Ialtmsrun1[barL],x]/.barL->barLfunc[mu,mu,mu-mue,x]
un1[mu_,mue_,x_]:=ln1[mu,mu-mue,mu,alphas0[barL],barL,Ialtmsrun1[barL],x]/.barL->barLfunc[mu,mu-mue,mu,x]

sn2[mu_,mue_,x_]:=Max[Re[sn2[1.000001mu,mu,mu-mue,alphas1[barL],barL,Ialtmsrun2[barL],x]/.barL->barLfunc[mu,mu,mu-mue,x]],0]
sn2b[mu_,mue_,x_]:=Re[sn2[1.000001mu,mu,mu-mue,alphas1[barL],barL,Ialtmsrun2[barL],x]/.barL->barLfunc[mu,mu,mu-mue,x]]
dn2[mu_,mue_,x_]:=Max[ln2[1.000001mu,mu,mu-mue,alphas1[barL],barL,Ialtmsrun2[barL],x]/.barL->barLfunc[mu,mu,mu-mue,x],0]
un2[mu_,mue_,x_]:=Max[ln2[1.000001mu,mu-mue,mu,alphas1[barL],barL,Ialtmsrun2[barL],x]/.barL->barLfunc[mu,mu-mue,mu,x],0]

Electric neutrality :

getmue1[mu_,x_]:=mue/.FindRoot[2/3un1[mu,mue,x]-1/3dn1[mu,mue,x]-1/3sn1[mu,mue,x]-mue^3/3/Pi^2,{mue,0.001,0.01}]
getmue2[mu_,x_]:=mue/.FindRoot[2/3un2[mu,mue,x]-1/3dn2[mu,mue,x]-1/3sn2[mu,mue,x]-mue^3/3/Pi^2,{mue,0.001,0.01}]

neutraln1total[mu_,x_]:=(sn1[mu,mue,x]+un1[mu,mue,x]+dn1[mu,mue,x])/.mue->getmue1[mu,x]
neutraln2total[mu_,x_]:=(sn2[mu,mue,x]+un2[mu,mue,x]+dn2[mu,mue,x]+mue^3/(3π^2))/.mue->getmue2[mu,x]
neutralsn1total[mu_,x_]:=(sn1[mu,mue,x])/.mue->getmue1[mu,x]
neutralsn2total[mu_,x_]:=(sn2[mu,mue,x])/.mue->getmue2[mu,x]

n2total[mu_,x_]:={mu,sn2[mu,mue,x],un2[mu,mue,x],dn2[mu,mue,x],mue^3/(3π^2)}/.mue->getmue2[mu,x]

dgetmue2[mu_,x_,eps_]:=(getmue2[mu+eps,x]-getmue2[mu-eps,x])/2/eps;

EoS

If[n2total[muMin,X][[2]]≠0,Print["WARNING: strange quark number density not vanishing at mu=muMin. It is strongly suggested you decrease muMin!"]]

Ntable=Table[ n2total[mu,X], {mu, muMin, muMax, dmu}];