This file generates the equation of state and MR plots from tabulated particle numbers (use the ones from tabluated.nb or create your own using EoS.nb).The parameter muMatch gives the quark number chemical potential (in GeV) where the quark and hadronic equations of states are matched.

Set path for hadronic equations of state (these files are also provided on the webpage):

PathAPR="~/Downloads/av18_uix_dvb_conv.dat";
PathKaons="~/Downloads/schaf1_conv.dat";
PathHyperons="~/Downloads/V18UIXNSC89_conv.dat";

-Pick a value for the matching point (quark number chemical mu_B/3 potential in GeV):

muMatch=0.987/3

Hint : See Fig. 8 of our article (pdf version provided on the webpage) to get an idea of muMatch for caseI and case II, respectively

-Choose one of the hadronic Equations of state :

APR :

HadronicEoS=0;

Hyperonic :

HadronicEoS=1;

Kaons :

HadronicEoS=2;

-Choose quark equation of state (from tabutaned.nb or create your own Ntable using EoS.nb) :

Central values

No CSC

Ntable=LambdaMSbar0378ms0010Δ0X1;

Ntable=LambdaMSbar0378ms0010Δ0X2;

Ntable=LambdaMSbar0378ms0010Δ0X4;

With CSC

Ntable=LambdaMSbar0378ms0010Δ01X4;

Ntable=LambdaMSbar0378ms0010Δ01X2;

Ntable=LambdaMSbar0378ms0010Δ01X1;

Favoured values

No CSC

Ntable=LambdaMSbar0346ms0007Δ0X1;

Ntable=LambdaMSbar0346ms0007Δ0X2;

Ntable=LambdaMSbar0346ms0007Δ0X4;

With CSC

Ntable=LambdaMSbar0346ms0007Δ01X4;

Ntable=LambdaMSbar0346ms0007Δ01X2;

Ntable=LambdaMSbar0346ms0007Δ01X1

Other Selected values

No CSC

Ntable=LambdaMSbar0378ms0007Δ0X4;

Run This

Hadronic EoS:

gevtoinvfm=5.06842372022301;

APR

nucl=Import[PathAPR];
nucl=Delete[nucl,1];

nuclp=Transpose[{nucl[[All,1]]/3,nucl[[All,2]]/gevtoinvfm^3}];
nucleps=Transpose[{nucl[[All,1]]/3,nucl[[All,3]]/gevtoinvfm^3}];
nucln=Transpose[{nucl[[All,1]]/3,nucl[[All,4]]/gevtoinvfm^3}];

Show[ListPlot[nuclp,PlotStyle-> RGBColor[1,0,0]],ListPlot[nucleps]]

Hyperons

hyperon=Import[PathHyperons];
hyperon=Delete[hyperon,1];

hyperonp=Transpose[{hyperon[[All,1]]/3,hyperon[[All,2]]/gevtoinvfm^3}];
hyperoneps=Transpose[{hyperon[[All,1]]/3,hyperon[[All,3]]/gevtoinvfm^3}];
hyperonn=Transpose[{hyperon[[All,1]]/3,hyperon[[All,4]]/gevtoinvfm^3}];

Show[ListPlot[hyperonp,PlotStyle-> RGBColor[1,0,0]],ListPlot[hyperoneps]]

Kaonic

kaon=Import[PathKaons];
kaon=Delete[kaon,1];

kaonp=Transpose[{kaon[[All,1]]/3,kaon[[All,2]]/gevtoinvfm^3}];
kaoneps=Transpose[{kaon[[All,1]]/3,kaon[[All,3]]/gevtoinvfm^3}];
kaonn=Transpose[{kaon[[All,1]]/3,kaon[[All,4]]/gevtoinvfm^3}];

Show[ListPlot[kaonp,PlotStyle-> RGBColor[1,0,0]],ListPlot[kaoneps]]

hadronicp=Which[HadronicEoS==0,nuclp,HadronicEoS==1,hyperonp,HadronicEoS==2,kaonp];

hadroniceps=Which[HadronicEoS==0,nucleps,HadronicEoS==1,hyperoneps,HadronicEoS==2,kaoneps];

hadronicn=Which[HadronicEoS==0,nucln,HadronicEoS==1,hyperonn,HadronicEoS==2,kaonn];

Quark EoS :

mu0=0.31

NOTE: "mu0" is just a new name for the constant mu0 (Lambdabar) appearing in Eq.(70) of our paper 
http : // arxiv.org/abs/0912.1856
In principle, the value of mumax is fixed by requiring the total number density to vanish at mu=mu0. However, for automation purposes it is much easier to set it to one (low) value since all the difference can be absorbed into the constant "B". Thus, the values for B generated below are dependent on the value chosen for mu0, and hence only of limited usefulness.

n0f3[mu_]:=(((mu^3) Nf)/(π^2))/.Nf-> 3
p0[mu_]:=Nf /4/Pi^2 mu^4/.Nf-> 3
eps0[mu_]:=p0[mu]+mu n0f3[mu]

Neutraln2total=Transpose[{Ntable[[All,1]],(Ntable[[All,2]]+Ntable[[All,3]]+Ntable[[All,4]]+Ntable[[All,5]])/n0f3[Ntable[[All,1]]]}];
Neutraln2quark=Transpose[{Ntable[[All,1]],(Ntable[[All,2]]+Ntable[[All,3]]+Ntable[[All,4]])/n0f3[Ntable[[All,1]]]}];

neutraln2totalI:=Interpolation[Neutraln2total]
neutraln2quark:=Interpolation[Neutraln2quark]
Ns:=Interpolation[Transpose[{Ntable[[All,1]],(Ntable[[All,2]])/n0f3[Ntable[[All,1]]]}]]
Nu:=Interpolation[Transpose[{Ntable[[All,1]],(Ntable[[All,3]])/n0f3[Ntable[[All,1]]]}]]
Nd:=Interpolation[Transpose[{Ntable[[All,1]],(Ntable[[All,4]])/n0f3[Ntable[[All,1]]]}]]
Ne:=Interpolation[Transpose[{Ntable[[All,1]],(Ntable[[All,5]])/n0f3[Ntable[[All,1]]]}]]
Mue[mu_]:=((3π^2)Ne[mu]*n0f3[mu])^(1/3)
dMuedmu[mu_]:=Evaluate[D[Mue[mu],mu]]
MuQ[mu_]:=mu-1/3 Mue[mu]
MuS[muq_]:=mu/.FindRoot[mu-1/3Mue[mu]==muq,{mu,muq}]

neutralp2total[mu_,B_]:=If[mu<MuS[muMatch],-B,NIntegrate[Ns[mup]n0f3[mup]+Nd[mup]n0f3[mup]+(1-dMuedmu[mup])Nu[mup]n0f3[mup]+dMuedmu[mup]Ne[mup]n0f3[mup],{mup,mu0,mu}]-B]

neutrale2total[mu_,B_]:=-neutralp2total[mu,B]+(mu Ns[mu]+mu Nd[mu]+(mu-Mue[mu])Nu[mu]+Mue[mu]Ne[mu])*n0f3[mu]

Matching

FindRoot[Interpolation[hadronicp][muMatch] == neutralp2total[MuS[muMatch], B], {B, 0.1}]
Bmatch = B /. %

matchedeos = {}
Quiet[Do[mmu = hadronicp[[i, 1]]; If[mmu > muMatch, matchedeos = Append[matchedeos, {mmu, neutralp2total[MuS[mmu], Bmatch], neutrale2total[MuS[mmu], Bmatch]}], matchedeos = Append[matchedeos, {mmu, hadronicp[[i, 2]], hadroniceps[[i, 2]]}]], {i, 1, Length[hadroniceps]}]]

Mass - Radius curves from Tolmann - Oppenheimer - Volkov

G = 1/(1.2211* 10^(19))^2;

GeVtokg = 1.7827 *10^(-27);

ooGeVtokm = 1.9733 10^(-16)/1000;

newNUQSPress[EOST_, LAST_] := Module[{rho, rhou, M, dp, p, r, dr, dr2, cnt},
  (*LAST=Length[EOST];*)
  (*LAST=30;*)
  p = EOST[[LAST, 1]];
  rho = EOST[[LAST, 2]];
  dp = (EOST[[LAST - 1, 1]] - EOST[[LAST, 1]]);
  dr = Sqrt[-dp/4/Pi/G/(p + rho)/(p + rho/3)];
  r = dr;
  M = 4*Pi*dr^3/3*rho;
  cnt = 0;
  While[LAST - cnt > 2, cnt = cnt + 1;
   p = EOST[[LAST - cnt, 1]];
   rho = EOST[[LAST - cnt, 2]];
   rhou = EOST[[LAST - cnt - 1, 2]];
   dp = (EOST[[LAST - 1 - cnt, 1]] - EOST[[LAST - cnt, 1]]);
   dr = -dp*1/N[G*(p + rho)*(M + 4*Pi*p*r^3)/(r - 2*G*M)/r]; p = EOST[[LAST - cnt - 1, 1]];
   dr2 = -dp*1/N[G*(p + rhou)*(M + 4 Pi rhou (r + dr)^2 dr + 4*Pi*p*(r + dr)^3)/(r + dr - 2*G*(M + 4 rhou Pi (r + dr)^2 dr))/(r + dr)];
   dr = (dr + dr2)/2;
   M = N[M + 4*Pi*(rho + rhou)/2*(r + dr/2)^2*dr];
   r = N[r + dr];]; {r*ooGeVtokm, M*GeVtokg/1.989/10^(30)}]

hybrid = {};
Do[hybrid = Append[hybrid, {matchedeos[[Length[matchedeos] - i, 2]] , matchedeos[[Length[matchedeos] - i, 3]]}], {i, 0, Length[matchedeos] - 1}]

HybridStar=Table[newNUQSPress[hybrid,n],{n,2,Length[hybrid]-2,1}];

hadroneos=Reverse[Transpose[{hadronicp[[All,2]],hadroniceps[[All,2]]}]];

HadronStar=Table[newNUQSPress[hadroneos,n],{n,2,Length[hadroneos]-2,1}];

Results

Some information:

Print["Matching at rho=",gevtoinvfm^3Interpolation[hadronicn][muMatch],"/fm^3 and mu_B=",3*muMatch" GeV and B=(",Sign[Bmatch],")*(",Abs[Bmatch]^(1/4),")^4 GeV"]

Print["You have chosen to match at mu_B=",3*muMatch,"For case I, matching is possible up to (approximately) mu_B=",3*MuQ[mu/.Quiet[FindRoot[Interpolation[hadroniceps][MuQ[mu]]==neutrale2total[mu,0],{mu,0.31,0.35}]]]]; Print["You may consider adjusting muMatch"]

Pressure (normalized by free quarks)

Blue dots : hardonic EoS
Thin line : Quark EoS
Red dots : Hybrid EoS

Show[ListPlot[Transpose[{hadronicp[[All,1]],hadronicp[[All,2]]/p0[matchedeos[[All,1]]]}],AxesLabel-> {"(Subscript[μ, B])/3","P / (Subscript[P, 0])"}],ListPlot[Transpose[{matchedeos[[All,1]],matchedeos[[All,2]]/p0[matchedeos[[All,1]]]}],PlotStyle-> RGBColor[1,0,0]],Plot[neutralp2total[MuS[mu],Bmatch]/p0[mu],{mu,muMatch,1.7}],PlotRange-> {{0.31,1.7},{0,1}}]

Energy density

Blue dots : hardonic EoS
Thin line : Quark EoS
Red dots : Hybrid EoS

Quiet[Show[ListPlot[Transpose[{hadroniceps[[All,1]],hadroniceps[[All,2]]/eps0[matchedeos[[All,1]]]}],AxesLabel-> {"(Subscript[μ, B])/3","ϵ / (Subscript[ϵ, 0])"}],ListPlot[Transpose[{matchedeos[[All,1]],matchedeos[[All,3]]/eps0[matchedeos[[All,1]]]}],PlotStyle-> RGBColor[1,0,0]],Plot[neutrale2total[MuS[mu],Bmatch]/eps0[mu],{mu,muMatch,1.7}],PlotRange-> {{0.31,1.7},{0,1}}]]

Mass Radius relationship

ListPlot[HybridStar,PlotRange-> {{0,40},{0,3}},AxesLabel-> {"R[km]","M / (Subscript[M, sol])"}]

ListPlot[HadronStar,PlotRange-> {{0,40},{0,3}},AxesLabel-> {"R[km]","M / (Subscript[M, sol])"},PlotStyle->RGBColor[1,0,0]]

Show[%,%%]

Tabulate EoS in nuclear physics units

Gives the EoS in tabulated format EoSlist[[n]]={Baryon mu[GeV], Pressure[GeV/fm^3], Energy Dens.[GeV/fm^3], Baryon Dens.[1/fm^3]}

EoSlist={};
Quiet[Do[mmu=hadronicp[[i,1]];
If[mmu>muMatch,EoSlist=Append[EoSlist,{3*mmu,neutralp2total[MuS[mmu],Bmatch]*gevtoinvfm^3,neutrale2total[MuS[mmu],Bmatch]*gevtoinvfm^3,neutraln2quark[MuS[mmu]]*gevtoinvfm^3*n0f3[MuS[mmu]]/3}],EoSlist=Append[EoSlist,{3*mmu,hadronicp[[i,2]]*gevtoinvfm^3,hadroniceps[[i,2]]*gevtoinvfm^3,hadronicn[[i,2]]*gevtoinvfm^3}]],{i,1,Length[hadroniceps]}]]

EoSlist