Preamble

Presets

SetOptions[EvaluationNotebook[],StyleDefinitions->Notebook[{Cell[StyleData[StyleDefinitions->"Default.nb"]],Cell[StyleData["Notebook"],Background->LightGray],Cell[StyleData["Input"],Background->RGBColor[0.87,0.94,1]],Cell[StyleData["Output"],Background->GrayLevel[1]],Cell[StyleData["Graphics"],FontFamily->"CMU Serif",Black,FontSize->14]},StyleDefinitions->"PrivateStylesheetFormatting.nb"]]


colours=ColorData[10,"ColorList"];
colourdashed=Transpose[{Dashed ConstantArray[1,11],ColorData[10,"ColorList"]}];
colourthin=Transpose[{Thin ConstantArray[1,11],ColorData[10,"ColorList"]}];
colourthick=Transpose[{Thick ConstantArray[1,11],ColorData[10,"ColorList"]}];
blue=RGBColor[0,0.5,1]

(*SetOptions[EvaluationNotebook[],Background->White]*)

SetOptions[{Plot,LogPlot,LogLogPlot,ListLogPlot,ParametricPlot,LogLinearPlot,ListPlot,ListLogLinearPlot},PlotStyle->colours,Background->White,Frame->True,FrameStyle->Black,FrameTicksStyle->Directive[FontFamily->"CMU Serif",FontColor->Black],PlotRangePadding->None,FrameTicks->Automatic];

SetOptions[{LogLogPlot,ListLogLogPlot},PlotStyle->colours,Frame->True,FrameStyle->Black,FrameTicksStyle->Directive[FontFamily->"CMU Serif",FontColor->Black],FrameTicks->{{Table[{10^i,Superscript[10,i]},{i,-50,20}],None},{Table[{10^i,Superscript[10,i]},{i,-50,20}],None}}];

SetDirectory[NotebookDirectory[]]
<<MaTeX`

SetOptions[MaTeX,"Preamble"->{"\\usepackage{color,txfonts}"}];

MaTeX

<<MaTeX`
SetOptions[MaTeX,"Preamble"->{"\\usepackage{color,txfonts}"}];
texttoMat[text_]:=MaTeX[text,Magnification->1.3]
numtoMat[num_]:=MaTeX[ToString[num],Magnification->1.3]
mat[x_]:=MaTeX[x,Magnification->1.3]

colours=ColorData[10,"ColorList"];
colourdashed=Transpose[{Dashed ConstantArray[1,11],ColorData[10,"ColorList"]}];
colourlist={ColorData[10][1],ColorData[10][3],ColorData[10][5],ColorData[10][7],ColorData[10][9],ColorData[10][11]};

$Assumptions=cm>0 && g>0&&s>0&&&&Ke>0&&erg>0&&Element[{cm,g,s,Ke,erg},Reals]

Polytrope model

Unit conversions and relevant constants

GeVtogram=1.78266*10^(-24)g;
GeVtocmm1=1/(1.97*10^(-16)m/.m->10^2cm);
GeV4togcm3=(GeVtogram*GeVtocmm1^3);
GeVcm3togcm3=GeVtogram/cm^3;
μe=1.143;κ=0.01;
α=2.32;
Rsun=6.957*10^(10)cm;
Msun= 1.98855*10^30*kg;
LsunGeV=3.827*10^33erg/s/.erg->624.15064799632GeV;
LsunSI=3.827*10^26W/.W->1kg*m^2/s^3/.kg->10^3g/.m->10^2cm;
Lsunerg=3.827*10^33erg/s;
σB=5.670374419  (10^(-5))(g (cm^2)/(s^2)) (cm^(-2))(s^(-1))(K^(-4));
kB=1.38064852*10^−23*kg*m^2/(s^2)/K;

μ1=1.24;
μ1modGen[xh_]=(1/μ1+3/2*(xh(1-xh))/(2-xh))^(-1);
μ1mod=μ1modGen[0.51/2];(* Model D in 1607.04338 table 1 *)
Rgas=8.314462618 J/Ke/.J->kg* m^2/s^2/.kg->10^3g/.m->10^2cm;
G=6.6743*10^(-11)  m^3/s^2/kg/.kg->10^3g;
msun=1.989*10^33g;
cPoly=10^13cm^4g^(-2/3)/s^2;
me=9.1093837*10^(-28)g;
hbar=1.05457182*10^(-34)m^2 kg/s/.m->10^2cm/.kg->10^3g;
Na=6.02214076*10^(23)/g;
E0=6.73857*10^(49)erg;
mp=1.67262192×(10^(-24)) g;

Central quantities and luminosities

γ(ψ)

gp[ψ_]=-(5/16) ψ Log[1+(E^(-1/ψ))]+(15/8) (ψ^2) (((π^2)/3)+PolyLog[2,-(E^(-1/ψ))])

Central Temperature in units of 10^6 K. 

M here and in the following in units of Msun.

Tc[M_,ψ_]=7.68097*10^2*M^(4/3)*μe^(8/3)*ψ/(1+gp[ψ]+α*ψ)^2

Central Temperature taking Coulomb pressure into account. Default is the Debye-Huckel approximation, for an ionized gas. Wigner-Seitz approximation for a zero temperature gas.

X=0.75;

electronCharge=4.8032*10^(-10)cm^(3/2) g^(1/2)/s;
Kc=-3/10(4*Pi*z2bar/3)^(1/3)*electronCharge^2*(Na/μe)^(4/3)/.z2bar->(X*(1^2)+(1-X)*(2^2))/.X->0.75//Simplify;

Pcoulomb[ρ_,T_]=-(1/3) (electronCharge^3) (Sqrt[(Pi/(kB T))])*((((((X(1+1)+(((1-X)(2+4))/4))*ρ)/mp)))^(3/2))/.ρ->g/cm^3*ρ/.T->T*K/.kg->(10^3)g/.m->(10^2)cm/.g->1/.s->1/.cm->1//Simplify
TcCoulomb[M_,ψ_]=(((((centralPressure[M,ψ])/(centralPressure[M,ψ]+Pcoulomb[centralDensity[M,ψ],(10^6)Tc[M,ψ]]))))^2)Tc[M,ψ]

PcoulombWS[ρ_,T_]=(Kc*ρ^(4/3)/.ρ->ρ*g/cm^3//Simplify)/.g->1/.cm->1/.s->1

TcCoulombWS[M_,ψ_]=(((((centralPressure[M,ψ])/(centralPressure[M,ψ]+PcoulombWS[centralDensity[M,ψ],(10^6)Tc[M,ψ]]))))^2)Tc[M,ψ]

Radius in units of Rsun.

Radius[M_,ψ_]=(2.80858*10^9cm/Rsun)*M^(-1/3)*μe^(-5/3)*(1+gp[ψ]+α*ψ)

Central density in g/cm3.

centralDensity[M_,ψ_]=1.28412*10^5M^2*μe^(5)/(1+gp[ψ]+α*ψ)^(3)

Central pressure in units of g/cm/s^2.

centralPressure[M_,ψ_]=3.26763*10^(21)M^(10/3)*μe^(20/3)/(1+gp[ψ]+α*ψ)^4

Luminosity in hydrogen burning luminosity in units of Lsun

LHB[M_,ψ_]=7.33*10^16*M^(11.977)*ψ^(6.0316)/(1+gp[ψ]+α*ψ)^(16.466);

Luminosity in dark matter burning in units of Lsun including gravitational focusing. Note that v is in cm/s.

LDM[ρ_,v_,M_,ψ_]=(((Sqrt[((8*Pi)/3)]*(R0^2)*(((μe)^(-10/3))/(M^(2/3)))*(((1+gp[ψ]+α*ψ))^2)*ρ*v)/Lsun)(1+(3/2)((2 G M msun/(R0 *((μe)^(-5/3))(1+gp[ψ]+α*ψ)(M^(-1/3))))/(v^2))))/.R0->2.80858*10^9cm/.ρ->ρ*(GeVtogram/(cm^3))*c^2/.c->3*10^10cm/s/.v->v*cm/s/.Lsun->3.827*10^26W/.W->1kg*m^2/s^3/.kg->10^3g/.m->10^2cm

Quantities depending on the photosphere phase transition

Quantities depend on the parameters b1 and ν which depend on the model for the phase transition, see 

G. Chabrier, D. Saumon, W. B. Hubbard, and J. I. Lunine, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 391, no. 2, June 1, 1992, p. 817-826. 391, 817 (1992).

Teffmodel[b1_,ν_,ρe_,ψ_]=b1 (10^6)((ρe)^(2/5))(ψ^ν)K;

ρemodel[b1_,ν_]=ρe/.Solve[((ρe NA KB Teffmodel[b1,ν,ρe,ψ])/(μ2))==((11.2193 ((10^6) g (cm^(-1))(s^(-2)))(*to convert to cgs*))/(κR))(((MinMsun))^(5/3))(((μE)^(10/3))/(((1+γ+αa ψ))^2)),ρe][[1]]

(* μ2 is the mean molecular weight assuming a certain masss fraction hydrogen/helium *)
TeffModel[b1_,ν_]=Simplify[Teffmodel[b1,ν,ρemodel[b1,ν]((((g/(cm^3))))^(-5/7)),ψ]/.NA->Na/.KB->kB/.μE->μe/.μ2->((((0.75/2)+(0.25/4)))^(-1))/.kg->(10^3)g/.cm->(10^(-2))m,Assumptions->{g>0,m>0}]

Rad[MinMsun_,ψ_]=(2.80858*(10^9)cm)*(MinMsun^(-1/3))*((μe)^(-5/3))*(1+γ+αa*ψ);

LsurfModel[b1_,ν_]=Simplify[4π σB ((Rad[MinMsun,ψ])^2)(((TeffModel[b1,ν])^4)/LsunGeV)/.cm->(10^(-2))m/.g->(10^(-3))kg/.kg->((299792458)^2/(1.602176565*10^−19))*eV/.
s->((1.602176565*10^−19)/(1.054571800*10^−34))*eV^-1/.
m->((1.602176565*10^−19)/(1.054571800*10^−34)/(299792458))*eV^-1/.eV->(10^(-9))GeV,Assumptions->{MinMsun>0,b1>0}]

General model

Surface luminosity in units of Lsun

LsurfMod[M_,ψ_]=LsurfModel[b1,ν]/.κR->κ/.MinMsun->M/.γ->gp[ψ]/.αa->α
LsurfMod1[M_,ψ_,b1_,ν_]=LsurfMod[M,ψ];

Effective Temperature in units of 10^4 K

TeffMod[M_,ψ_]=((TeffModel[b1,ν])/((10^4)K))/.κR->κ/.MinMsun->M/.γ->gp[ψ]/.αa->α

Cooling differential equation - L is any source of burning.

ψODEMod[L_,ψ_,M_]=ψ'[t]==μ1mod/E0/1.5/M^(7/3)/μe^(8/3)*(1+gp[ψ]+α*ψ)^(2)*(L-LsurfMod[M,ψ]*Lsunerg/.s->1)/.erg->1

model D, used in most of the plots in the paper

Surface luminosity in units of Lsun.

Lsurf[M_,ψ_]=0.4147/κ^(1.1424)*M^(1.239)*ψ^(4.5797)/(1+gp[ψ]+α*ψ)^(0.2848);

Central Temperature in units of 10^6 K.

Tc[M_,ψ_]=7.68097*10^2*M^(4/3)*μe^(8/3)*ψ/(1+gp[ψ]+α*ψ)^2;

Effective Temperature in units of 10^4 K.

Teff[M_,ψ_]=2.57881*M^(0.4764)/κ^(0.2856)*ψ^(1.1456)/(1+gp[ψ]+α*ψ)^(0.5712);

Radius in units of Rsun.

Radius[M_,ψ_]=(2.80858*10^9cm/Rsun)*M^(-1/3)*μe^(-5/3)*(1+gp[ψ]+α*ψ)

Cooling differential equation - L is any source of burning.

ψODE[L_,ψ_,M_]=ψ'[t]==μ1mod/E0/1.5/M^(7/3)/μe^(8/3)*(1+gp[ψ]+α*ψ)^(2)*(L-Lsurf[M,ψ]*Lsunerg/.s->1)/.erg->1

Cooling curves

Time evolution in model D

Hydrogen burning only (SM)

The right hand side of the ODE is

ψODESMRHS[M_,ψ_]:=(ψODE[LHB[M,ψ]*Lsunerg/.s->1,ψ,M]/.erg->1)[[2]]

With derivative

ψODEderiv[M_,ψ_]:=D[ψODESMRHS[M,ψ],ψ]

Solve the ODE

ψODESMStarRHSYEARS[M_,ψ_]:=ψODESMRHS[M,ψ]/s/.s->3.17098*10^(-8)yr/.yr->1

ψsolSM[M_?NumericQ,IC_?NumericQ]:=NDSolve[{ψ'[t]==ψODESMStarRHSYEARS[M,ψ[t]],ψ[0]==IC},ψ,{t,0,10^12},AccuracyGoal->11,PrecisionGoal->11,WorkingPrecision->20]

Hydrogen burning + DM annihilation

The right hand side of the ODE is

ψODEhybridStarRHS[M_,ψ_,ρ_,v_]:=(ψODE[(LHB[M,ψ]+LDM[ρ,v,M,ψ])*Lsunerg/.s->1,ψ,M]/.erg->1)[[2]]

With derivative

ψODEderivhybridStar[M_,ψ_,ρ_,v_]:=D[ψODEhybridStarRHS[M,ψ,ρ,v],ψ]

ψODEhybridStarRHSYEARS[M_,ψ_,ρ_,v_]:=ψODEhybridStarRHS[M,ψ,ρ,v]/s/.s->3.17098*10^(-8)yr/.yr->1

Solve the ODE

ψsolhybridStar[M_?NumericQ,ρ_?NumericQ,v_?NumericQ,IC_]:=NDSolve[{ψ'[t]==ψODEhybridStarRHSYEARS[M,ψ[t],ρ,v],ψ[0]==IC},ψ,{t,0,3*10^10},AccuracyGoal->10,PrecisionGoal->10]

Test the late-time behavior converges to the theoretical value

ψEquilibriumhybridStar[M_?NumericQ,IC_?NumericQ,ρ_?NumericQ,v_?NumericQ]:=FindRoot[ψODEhybridStarRHS[M,ψ,ρ,v]==0,{ψ,IC},MaxIterations->Infinity]

(*theoretical value*)
ψ/.ψEquilibriumhybridStar[0.1,1,10^5,50*10^5]

(*late time behaviour*)
Evaluate[ψ[3 10^10]/.ψsolhybridStar[0.1,10^5,50*10^5,10][[1]]]

Stability in model D

HB only (SM)

There are two equilibria when the mass is larger than the minimum

ψEquilibriumSM[M_?NumericQ,IC_?NumericQ]:=FindRoot[ψODESMRHS[M,ψ]==0,{ψ,IC}]

{ψEquilibriumSM[0.075,0.02],ψEquilibriumSM[0.075,1]}

To evaluate stability we need to evaluate the derivative of the RHS as a function of mass. We start with the branch where ψ<ψmin with the latter being the value corresponding to the minimum mass.

Plot[Evaluate[ψ/.ψEquilibriumSM[M,0.02]],{M,0.075,0.1}]

We can see this branch is unstable

Plot[Evaluate[ψODEderiv[M,ψ]/.ψEquilibriumSM[M,0.02]],{M,0.075,0.1}]

Now plot the large ψ branch

Plot[Evaluate[ψ/.ψEquilibriumSM[M,1]],{M,0.075,0.1}]

This branch is stable:

Plot[Evaluate[ψODEderiv[M,ψ]/.ψEquilibriumSM[M,1]],{M,0.075,0.1}]

HB + DM

We can see that hybrid stars are stable

Plot[Evaluate[ψODEderivhybridStar[M,ψ,10^5,50*10^5]/.ψEquilibriumhybridStar[M,6,10^5,50*10^5]],{M,0.03,0.2},PlotRange->All]

Time evolution plots

Benchmark initial condition for ψ

icbench=.3;

Clear[plotLuminosity]

plotLuminosity[r_,mass_,colour_]:=Module[{ρgrid,LumPlot,linecolors,legend,SMLumPlot,label},
ρgrid=10^r;
LumPlot=ListLogPlot[
Table[{Log10[t],Lsurf[mass,ψ[t]/.ψsolhybridStar[mass,ρgrid,50*10^5,icbench]]}/.t->(10^tt)//Flatten,{tt,6.5,10.5,.2}],Joined->True,
PlotStyle->colour];
label={mat[StringJoin[ToString[mass],"{\\rm M}_\\odot"]],mat[HoldForm[ρ=10^r"GeV/(cm^3)"]]};
SMLumPlot=ListLogPlot[
Table[{Log10[t],Lsurf[mass,ψ[t]/.ψsolSM[mass,icbench]]}/.t->(10^tt)//Flatten,{tt,6.5,10.5,.2}],Joined->True,PlotStyle->{colour,Dashed},
Frame->True,PlotRange->{(10^(-6.8)),5 (10^(-2))},PlotRangePadding->None,FrameTicks->{{Automatic,Table[{10^i,""},{i,-20,20}]},{Table[{i,Superscript[10,i]},{i,0,12}],Table[{i,""},{i,0,10,1}]}},
FrameStyle->Directive[Black,FontSize->15],FrameLabel->{mat["t/{\\rm years}"],mat["L/{\\rm L}_\\odot"],label}];
Show[SMLumPlot,LumPlot]
]//Quiet

plotHBratioLog[r_,mass_,colour_]:=Module[{ρgrid,LumPlot,linecolors,legend,SMLumPlot,label},
ρgrid=10^r;
LumPlot=ListLogPlot[
Table[{Log10[t],((LHB[mass,ψ[t]/.ψsolhybridStar[mass,ρgrid,50*10^5,icbench]])/(Lsurf[mass,ψ[t]/.ψsolhybridStar[mass,ρgrid,50*10^5,icbench]]))}/.t->(10^tt)//Flatten,{tt,6.5,10.5,.2}],Joined->True,
PlotStyle->colour];
label={mat[StringJoin[ToString[mass],"{\\rm M}_\\odot"]],mat[HoldForm[ρ=10^r"GeV/(cm^3)"]]};
SMLumPlot=ListLogPlot[
Table[{Log10[t],((LHB[mass,ψ[t]/.ψsolSM[mass,icbench]])/(Lsurf[mass,ψ[t]/.ψsolSM[mass,icbench]]))}/.t->(10^tt)//Flatten,{tt,6.5,10.5,.2}],Joined->True,PlotStyle->{colour,Dashed},
FrameTicks->{{Automatic,Table[{10^i,""},{i,-20,20}]},{Table[{i,Superscript[10,i]},{i,0,12}],Table[{i,""},{i,0,10,1}]}},
Frame->True,PlotRange->{(10^(-6.8)),2 (10^0)},PlotRangePadding->None,
FrameStyle->Directive[Black,FontSize->15],FrameLabel->{mat["t/{\\rm years}"],mat["L_{\\rm HB}/L_{\\rm total}"],label}];
Show[SMLumPlot,LumPlot]
]//Quiet

plotRad[r_,mass_,colour_]:=Module[{ρgrid,LumPlot,linecolors,legend,SMLumPlot,label},
ρgrid=10^r;
LumPlot=ListPlot[
Table[{Log10[t],Radius[mass,ψ[t]/.ψsolhybridStar[mass,ρgrid,50*10^5,icbench]]}/.t->(10^tt)//Flatten,{tt,5.9,10.5,.1}],Joined->True,
PlotStyle->colour];
label={mat[StringJoin[ToString[mass],"{\\rm M}_\\odot"]],mat[HoldForm[ρ=10^r"GeV/(cm^3)"]]};
SMLumPlot=ListPlot[
Table[{Log10[t],Radius[mass,ψ[t]/.ψsolSM[mass,icbench]]}/.t->(10^tt)//Flatten,{tt,6.5,10.5,.2}],Joined->True,PlotStyle->{colour,Dashed},
Frame->True,PlotRange->{.07,.23},PlotRangePadding->None,
FrameStyle->Directive[Black,FontSize->15],FrameTicks->{{Automatic,Table[{i,""},{i,0,10,.05}]},{Table[{i,Superscript[10,i]},{i,0,12}],Table[{i,""},{i,0,10,1}]}},
FrameTicksStyle->Black,FrameLabel->{mat["t/{\\rm years}"],mat["R/{\\rm R}_\\odot"],label}];
Show[SMLumPlot,LumPlot]
]//Quiet

plotTeff[r_,mass_,colour_]:=Module[{ρgrid,LumPlot,linecolors,legend,SMLumPlot,label},
ρgrid=10^r;
LumPlot=ListPlot[
Table[{Log10[t],Teff[mass,ψ[t]/.ψsolhybridStar[mass,ρgrid,50*10^5,icbench]]}/.t->(10^tt)//Flatten,{tt,5.9,10.5,.1}],Joined->True,
PlotStyle->colour];
label={mat[StringJoin[ToString[mass],"{\\rm M}_\\odot"]],mat[HoldForm[ρ=10^r"GeV/(cm^3)"]]};
SMLumPlot=ListPlot[
Table[{Log10[t],Teff[mass,ψ[t]/.ψsolSM[mass,icbench]]}/.t->(10^tt)//Flatten,{tt,6.5,10.5,.2}],Joined->True,PlotStyle->{colour,Dashed},GridLines->None,PlotRange->{0.01,.42},
Frame->True,PlotRangePadding->None,FrameTicks->{{Automatic,Table[{i,""},{i,0,10,.1}]},{Table[{i,Superscript[10,i]},{i,0,12}],Table[{i,""},{i,0,10,1}]}},
FrameStyle->Directive[Black,FontSize->15],FrameLabel->{mat["t/{\\rm years}"],mat["T/10^4{\\rm K}"],label}];
Show[SMLumPlot,LumPlot]
]//Quiet

plotTc[r_,mass_,colour_]:=Module[{ρgrid,LumPlot,linecolors,legend,SMLumPlot,label},
ρgrid=10^r;
LumPlot=ListPlot[
Table[{Log10[t],Tc[mass,ψ[t]/.ψsolhybridStar[mass,ρgrid,50*10^5,icbench]]}/.t->(10^tt)//Flatten,{tt,5.9,10.5,.1}],Joined->True,
PlotStyle->colour];
label={mat[StringJoin[ToString[mass],"{\\rm M}_\\odot"]],mat[HoldForm[ρ=10^r"GeV/(cm^3)"]]};
SMLumPlot=ListPlot[
Table[{Log10[t],Tc[mass,ψ[t]/.ψsolSM[mass,icbench]]}/.t->(10^tt)//Flatten,{tt,6.5,10.5,.2}],Joined->True,PlotStyle->{colour,Dashed},GridLines->None,PlotRange->{0.01,3.9},
Frame->True,PlotRangePadding->None,FrameTicks->{{Automatic,Table[{i,""},{i,0,10,.1}]},{Table[{i,Superscript[10,i]},{i,0,12}],Table[{i,""},{i,0,10,1}]}},
FrameStyle->Directive[Black,FontSize->15],FrameLabel->{mat["t/{\\rm years}"],mat["T_c/10^6{\\rm K}"],label}];
Show[SMLumPlot,LumPlot]
]//Quiet

plotTcCoulomb[r_,mass_,colour_]:=Module[{ρgrid,LumPlot,linecolors,legend,SMLumPlot,label},
ρgrid=10^r;
LumPlot=ListLogLinearPlot[
Select[Table[{t,TcCoulomb[mass,ψ[t]/.ψsolhybridStar[mass,ρgrid,50*10^5,icbench]]}/.t->(10^tt)//Flatten,{tt,5,10.4,.2}],#[[2]]>0&],Joined->True,
PlotStyle->colour];
label={mat[StringJoin[ToString[mass],"{\\rm M}_\\odot"]],mat[HoldForm[ρ=10^r"GeV/(cm^3)"]]};
SMLumPlot=ListLogLinearPlot[
Table[{t,TcCoulomb[mass,ψ[t]/.ψsolSM[mass,icbench]]}/.t->(10^tt)//Flatten,{tt,6,10.4,.2}],Joined->True,PlotStyle->{colour,Dashed},
Frame->True,PlotRange->{0,4.7},PlotRangePadding->None,
FrameStyle->Directive[Black,FontSize->15],FrameTicks->Automatic,FrameTicksStyle->Black,FrameLabel->{mat["t/{\\rm years}"],mat["T_c/10^6\\rm K"],label}];
Show[SMLumPlot,LumPlot]
]//Quiet

List of masses to be plotted (in solar masses) -- users can define their own list

mlist={.03,.05,.08,.1};

plotLum4=Show[Table[plotLuminosity[4,mlist[[i]],colourlist[[i]]],{i,1,Length[mlist]}],
LogLogPlot[(10^(-9)),{x,(10^6),(10^11)},PlotLegends->Placed[LineLegend[colourlist,Table[mat[ToString[mlist[[i]]]],{i,1,Length[mlist]}],Spacings->0,LegendLabel->mat["M/\\rm M_\\odot"]],{Left,Bottom}]],Epilog->Inset[mat["\\rho = 10^4 \\rm \\, GeV \\, cm^{-3}"],{9.4,Log[(10^(-1.8))]}]
]

plotHBrat4=Show[Table[plotHBratioLog[4,mlist[[i]],colourlist[[i]]],{i,1,Length[mlist]}]
]

plotRad4=Show[Table[plotRad[4,mlist[[i]],colourlist[[i]]],{i,1,Length[mlist]}]]

mlist={.05,.08,.1};

plotTcCoulomb5=Show[Table[plotTcCoulomb[5,mlist[[i]],colourlist[[i+1]]],{i,1,Length[mlist]}],
LogLinearPlot[(10^(-9)),{x,(10^6),(10^11)},PlotLegends->Placed[LineLegend[colourlist[[2;;]],Table[mat[ToString[mlist[[i]]]],{i,1,Length[mlist]}],Spacings->0,Background->White,LegendLabel->mat["M/\\rm M_\\odot \\quad\\quad\\quad\\quad\\quad\\quad"],LegendLayout->{"Column",4}],{Left,Bottom}]],Epilog->{Background->White,Inset[mat["\\rho = 10^5 \\rm \\, GeV \\, cm^{-3}"],{Log[(10^9.5)],0.4}]}
]

Steady state solutions

Hydrogen burning only

Balance the surface luminosity and hydrogen burning luminosity:

Clear[eqnNoDMmod]
Lsurfmod1[M_,ψ_]=11.019690367991773` (b1^3) (ψ^(3 ν)) ((((( (ψ^(-ν)))/(b1 (((1.` +2.32` ψ-(5/16) ψ Log[1+(E^(-1/ψ))]+(15/8) (ψ^2) (((π^2)/3)+PolyLog[2,-(E^(-1/ψ))])))^2)))))^(1/7));

eqnNoDMmod=Lsurfmod1[M,ψ]==(LHB[M,ψ]/(M^1.238))

Mathematica can solve this:

Clear[MsolNoDMmod]
MsolNoDMmod[b1_,ν_]=Solve[eqnNoDMmod,M][[1]]

Which shows that there is a minimum value of M for a solution to exists, here an example for model D.

Plot[{M/.MsolNoDMmod[2,1.6],M/.MsolNoDM},{ψ,0,0.17},AxesStyle->Directive[Black,FontSize->15],AxesLabel->{mat["\\psi"],mat["M/{\\rm M}_\\odot"]},PlotRange->{Automatic,{0.07,0.11}}]//Quiet

DM annihilation only

Analytic solution not including gravitational focusing -- equation (8) in the paper

Clear[LsurfModel1]
(* Stripped of M dependence *)
LsurfModel1[b1_,ν_]=(( L0 (b1^3) (ψ^(3 ν)) ((((( (ψ^(-ν)))/(b1 κR (((1 +γ+αa ψ))^2)))))^(1/7)))/(κR))/.κR->κ/.αa->α;
(* L0 = 0.05707620723398145` *)

MeqnDarkStarMod=(M^1.239)*(M^(2/3))==((Sqrt[((8 π)/3)])*ρ*v*f*(R0^2)*((μe)^(-10/3))(((1+γ+α*ψ))^2))((LsurfModel1[b1,ν])^(-1))

Rewriting this

MeqnDarkStarMod=(M^(1.9056666666666668`))==((0.009601908920332507` f (R0^2) v ρ (ψ^(-3 ν)) (((1+γ+α ψ))^(2+2/7.)))/((b1^(3-1/7.)) L0 ((((ψ^(-ν))))^(1/7))))

This enables us to solve for M in a simple form

(* Using v = 50 km/s *)
MsolDarkStarMod[b1_,ν_]=Solve[MeqnDarkStarMod,M][[1]];
MsolDarkStarModbench[b1_,ν_]=MsolDarkStarMod[b1,ν]/.R0->2.80858*10^9cm (* R0 *)/.f->1/.ρ->ρ*(GeVtogram/(cm^3))*c^2/.c->3*10^10cm/s/.v->v*((50 km)/s)/.L0->
0.0578Lsun/.Lsun->3.827*10^26W/.W->1kg*m^2/s^3/.kg->10^3g/.km->(10^3)m/.m->10^2cm

eqnDarkStarMod[ρ_,v_,M_,ψ_,b1_,ν_]=(LsurfMod1[M,ψ,b1,ν])-(LDM[ρ,v,M,ψ])==0

Including gravitational focusing (no closed form solution for (Subscript[M, 0]))

Clear[LsurfModel1]
(* Stripped of M dependence *)
LsurfModel1[b1_,ν_]=(( L0 (b1^3) (ψ^(3 ν)) ((((( (ψ^(-ν)))/(b1 κR (((1 +γ+αa ψ))^2)))))^(1/7)))/(κR))/.κR->κ/.αa->α;
(* L0 = 0.05707620723398145` *)

MeqnDarkStarModGF=(M^1.239)*(M^(2/3))==((Sqrt[((8 π)/3)])(1+(3/2)(((2 Gn M M0/(R0 *((μe)^(-5/3))(1+γ+α*ψ)(M^(-1/3))))/(v^2))))*ρ*v*f*(R0^2)*((μe)^(-10/3))(((1+γ+α*ψ))^2))((LsurfModel1[b1,ν])^(-1))

MeqnDarkStarModGF=(M^(1.9056666666666668`))==((0.009601908920332507` f (R0^2) v ρ (ψ^(-3 ν)) (((1+γ+2.32` ψ))^(2+2/7)) (1+((3.7485643625284077` Gn (M^(4/3))M0)/(R0 (v^2) (1+γ+2.32` ψ)))))/((b1^(3-1/7)) L0 ((((ψ^(-ν))))^(1/7))))

When does gravitational focussing become important?

MeqnDarkStarModbench[ρ_]=MeqnDarkStarModGF/.M0->Msun/.γ->gp[ψ]/.Gn->G/.R0->2.80858*10^9cm/.f->1/.ρ->ρ*(GeVtogram/(cm^3))*c^2/.c->3*10^10cm/s/.v->(50 (km/s))/.L0->
0.0578Lsun/.Lsun->3.827*10^26W/.W->1kg*m^2/s^3/.kg->10^3g/.km->(10^3)m/.m->10^2cm/.ψ->.3/.ν->1.6/.b1->2

Table[{(10^logrho),M/.FindRoot[MeqnDarkStarModbench[(10^logrho)],{M,(10^0.)}]},{logrho,0,5,.2}];
Show[ListLogLogPlot[%,Joined->True],LogLogPlot[M/.MsolDarkStarModbench[2,1.6]/.γ->gp[ψ]/.ψ->.3/.v->1,{ρ,1,(10^6)},PlotStyle->Red]]

eqnDarkStarMod[ρ_,v_,M_,ψ_,b1_,ν_]=(LsurfMod1[M,ψ,b1,ν])-(LDM[ρ,v,M,ψ])==0

eqnDarkStarMod[ρ,v,M,ψ,2,1.6]

Hydrogen burning + DM annihilation

Steady state plot

Equations for stars that burn hydrogen only, DM annihilation only (dark stars), and both, in model D.

eqnDM[ρ_,v_,M_,ψ_]=Lsurf[M,ψ]-LDM[ρ,v,M,ψ]-LHB[M,ψ]==0

eqnDarkStar[ρ_,v_,M_,ψ_]=eqnDarkStarMod[ρ,v,M,ψ,2,1.6]
eqnNoDM=eqnNoDMmod/.b1->2/.ν->1.6

tableψ=Table[i,{i,(10^(-5)),1,((1-(10^(-5)))/150)}];
tableM=Table[i,{i,(*(10^(-5))*)3.8 (10^(-3)),1.2,((1.1-3.8 (10^(-3)))/500)}];

Plot of the steady state solutions for different DM densities

rhoDM=.;
makelogcontourplot[rhoDM_]:=Module[{noDM,onlyDM,hybrid,tableM,plot,DMequalH},
noDM=Table[{M/.FindRoot[Evaluate[eqnNoDM],{M,.1}],ψ},{ψ,tableψ}]//N//Quiet;
onlyDM=Table[{M/.FindRoot[Evaluate[eqnDarkStar[rhoDM,50*(10^5),M,ψ]],{M,1}],ψ},{ψ,tableψ}]//N//Quiet;
tableM=Table[(10^i),{i,Log10[onlyDM[[-1,1]]],Log10[1.1],((Log10[1.1]-Log10[onlyDM[[-1,1]]])/50)}]; 
hybrid=Select[Table[{M,ψ/.FindRoot[Evaluate[eqnDM[rhoDM,50*(10^5),M,ψ]],{ψ,1}]},{M,tableM}]//N//Quiet,#[[1]]>0&]//Quiet;

plotstyle={{Thick,Orange},{Thick,blue},{Dashing[.03],Black,Thickness[.007]}};
legend=LineLegend[{Orange,blue,{Black,Thickness[.008]}},{mat["\\text{H burning only}"],mat["\\text{DM only}"],mat["\\text{hybrid star}"]},Spacings->0,Background->Directive[White,Opacity[.8]],LegendLabel->mat[StringJoin["\\rho = 10^",ToString[Log10[rhoDM]],"\\, \\rm GeV \\, cm^{-3}"]]];
DMequalH=(Graphics`Mesh`FindIntersections@ListPlot[{noDM,onlyDM},Joined->True])[[1,1]];
plot=ListLogLinearPlot[{noDM,onlyDM,hybrid},Joined->True,FrameStyle->Black,
GridLines->{{DMequalH},{}},GridLinesStyle->Directive[Black, Dashed],
PlotRange->{{(10^(-4)),1},{0,1}(*{(10^(-2)),1}*)},PlotRangePadding->Scaled[.02],ImageSize->350,AspectRatio->2/3,
PlotStyle->plotstyle,PlotLegends->Placed[legend,{Left,Top}],
FrameLabel->{mat["M ({\\rm M}_\odot)"],mat["\\psi"]},
Epilog->{
Inset[Rotate[mat["\\color[rgb]{1.,.5,0.}{\\text{Red dwarf}}\\downarrow"],90Degree],{Log[1.5 DMequalH],.7}],
Inset[Rotate[mat["\\color[rgb]{0.,.5,1.}{\\text{Dark dwarf}}\\uparrow"],90Degree],{Log[.75 DMequalH],.7}]}];
(*SelectFirst[onlyDM[[;;,2]]-noDM[[;;2/]]];*)
plot
]

{makelogcontourplot[(10^0)],makelogcontourplot[(10^5)]}

Min mass for H burning -- parameter dependence

eqnDMmod[ρ_,v_,M_,ψ_,b1_,ν_]=LsurfMod1[M,ψ,b1,ν]-LDM[ρ,v,M,ψ]-LHB[M,ψ]==0

eqnDarkStarMod[ρ_,v_,M_,ψ_,b1_,ν_]=(LsurfMod1[M,ψ,b1,ν])-(LDM[ρ,v,M,ψ])==0

tableψ=Table[i,{i,(10^(-5)),1,((1-(10^(-5)))/250)}];

findminM2mod[rhoDM_,b1_,ν_]:=Module[{tableM,hybrid,onlyDM},
onlyDM=Table[{M/.FindRoot[Evaluate[eqnDarkStarMod[rhoDM,50*(10^5),M,ψ,b1,ν]],{M,1}],ψ},{ψ,tableψ}]//N//Quiet;
tableM=Table[(10^i),{i,Log10[onlyDM[[-1,1]]],Log10[1.1],((Log10[1.1]-Log10[onlyDM[[-1,1]]])/350)}]; 
hybrid=Select[Table[{M,ψ/.FindRoot[Evaluate[eqnDMmod[rhoDM,50*(10^5),M,ψ,b1,ν]],{ψ,1}]},{M,tableM}]//N//Quiet,#[[1]]>0&]//Quiet;

{rhoDM,SelectFirst[Table[{hybrid[[i,1]],LDM[rhoDM,50*(10^5),hybrid[[i,1]],hybrid[[i,2]]]-LHB[hybrid[[i,1]],hybrid[[i,2]]]},{i,1,Length[hybrid]}],#[[2]]<0&][[1]]}
]

Note that evaluating the next cell takes some time.

Note also that “not found” indicates no solution is found within the range of masses searched for.

rholist=Join[Table[i,{i,0,3}],Table[i,{i,3.25,4.,.25}],Table[i,{i,4.1,7.5,.1}]];

minMlistmodA=Monitor[Table[findminM2mod[(10^i),2.87,1.58],{i,rholist}],i]
minMlistmodB=Monitor[Table[findminM2mod[(10^i),2.70,1.59],{i,rholist}],i]
minMlistmodC=Monitor[Table[findminM2mod[(10^i),2.26,1.59],{i,rholist}],i]
minMlistmodD=Monitor[Table[findminM2mod[(10^i),2,1.6],{i,rholist}],i]
minMlistmodE=Monitor[Table[findminM2mod[(10^i),1.68,1.61],{i,rholist}],i]
minMlistmodF=Monitor[Table[findminM2mod[(10^i),1.29,1.59],{i,rholist}],i]
minMlistmodG=Monitor[Table[findminM2mod[(10^i),0.6,1.44],{i,rholist}],i]
minMlistmodH=Monitor[Table[findminM2mod[(10^i),0.4,1.3],{i,rholist}],i]

legendminH=LineLegend[colours,{mat["\\text{model A}"],mat["\\text{model B}"],mat["\\text{model C}"],mat["\\text{model D}"],mat["\\text{model E}"],mat["\\text{model F}"],mat["\\text{model G}"],mat["\\text{model H}"]},Spacings->0]

ListLogLinearPlot[{minMlistmodA,minMlistmodB,minMlistmodC,minMlistmodD,minMlistmodE,minMlistmodF,Select[minMlistmodG,#[[2]]>0&],Select[minMlistmodH,#[[2]]>0&]},Joined->True,PlotStyle->colours,Frame->True,PlotRange->{{(10^1),(10^8)},{0,1.}},FrameLabel->{mat["\\rho \\left( \\textrm{GeV cm}^{-3} \\right)"],mat["\\textrm{Minimum mass for H burning (M}_\\odot)"]},AspectRatio->.8,ImageSize->350,PlotLegends->Placed[legendminH,{Left,Top}]]

Lithium depletion

LaneEmden[n_]=(1/(ξ^2))(∂_ξ)((ξ^2)(∂_ξ)θ[ξ])+((θ[ξ])^n)==0;

θLE32[ξ_]=θ[ξ]/.NDSolve[{LaneEmden[3/2],θ[0.001]==1,θ'[0.001]==0},θ,{ξ,0.001,10}, MaxSteps->1000000];
Plot[θLE32[ξ],{ξ,0.01,10},PlotRange->{{0.001,10},{-1,1}}]
θLE32root=ξ/.FindRoot[θLE32[ξ]==0,{ξ,3}][[1]]//Chop

Lithium burning rate from G. R. Caughlan and W. A. Fowler, Atomic Data and Nuclear Data Tables 40, 283 (1988)

T9=(T/(10^9));
Q34=17.346;
T9A=T9/(1.0+0.759*T9);
T9A13=T9A^(1.0/3.0);
T9A56=T9A^(5.0/6.0);

(* this is the rate in ((((cm^3)/moles))^(N-1))/s where N is the number of particles involved in the reaction *)
Clear[f34]
f34[T_]=1.096*10^9/T9^(2/3)*Exp[-8.472/T9^(1/3)]-4.830*10^8*T9A56/T9^(3/2)*Exp[-8.472/T9A13]+1.06*10^10/T9^(3/2)*Exp[-30.442/T9];

(* this is the reverse rate *)
Clear[r34]
r34[T_]=f34*4.69*Exp[-201.291/T9];

ODE following https://arxiv.org/pdf/astro-ph/9612155 (with rate above)

Finding the RHS numerically in the polytrope model

icbench=.3;
Mbench=.065;
Rsuncm=7 (10^10);

RHSlitODEnum[M_,Tcen_ ?NumericQ,rhocen_ ?NumericQ,R_ ?NumericQ]:=-((4π X )/(M mp ))NIntegrate[(r^2)(((rhocen ((θLE32[((θLE32root r)/R)])^(3/2))))^2)((f34[Tcen θLE32[((θLE32root r)/R)]])/((Na g))),{r,((0.001 R)/(θLE32root)),R}](g^2)((((g/(cm^3))))^2)(((cm^3)/(s g)))((1/g))((1/Msun))((cm^3))(s)/.kg->(10^3)g

initLNF=Log[(10^(-2.))]; (* initial ln of the lithium fraction, since this is a diff eq for log f this should not matter *)
Mbench=.065;

lnfODEnum[M_,rhoDM_]:=Module[{RHSlitODEnumlist,RHSlitODEnumint},
RHSlitODEnumlist=Table[Flatten[{tt,Block[{t=(10^tt)},Evaluate[RHSlitODEnum[M,
(10^6)TcCoulomb[M,ψ[t]/.ψsolhybridStar[M,rhoDM,50*(10^5),icbench]][[1]],
centralDensity[M,ψ[t]/.ψsolhybridStar[M,rhoDM,50*(10^5),icbench]][[1]],
Radius[M,ψ[t]/.ψsolhybridStar[M,rhoDM,50*(10^5),icbench]][[1]]Rsun/cm] ]]}],{tt,3,10,.5}];
RHSlitODEnumint[t_]=Interpolation[RHSlitODEnumlist,InterpolationOrder->1][Log10[t]];
NDSolve[{lnf'[t]==(60^2)24 365RHSlitODEnumint[t] ,lnf[(10^3)]==initLNF},lnf,{t,(10^3),3 (10^10)}]]//Quiet

lnfODEnumlist={lnfODEnum[Mbench,0],lnfODEnum[Mbench,(10^5)],lnfODEnum[Mbench,(10^6)]}

LogLinearPlot[Evaluate[E^lnf[t]/.lnfODEnumlist],{t,(10^0),3 (10^10)},FrameTicks->Automatic,FrameLabel->{mat["t/\\rm years"],mat["f_{\\rm lithium}"]},PlotStyle->colours,PlotRange->{0,1.1 (E^initLNF)},GridLines->{None,{(E^initLNF)}},PlotLegends->Placed[LineLegend[colours,{mat["\\rm SM"],mat["\\rho_{\\rm DM} = 10^5 \\rm \\, GeV \\, cm^{-3}"],mat["\\rho_{\\rm DM} = 10^6 \\rm \\, GeV \\, cm^{-3}"]},Spacings->0,LegendLabel->mat[StringJoin["M_*/\\rm M_\\odot =",ToString[Mbench]]]],{Left,Bottom}]]

Note that the following takes some time to evaluate.

initLNF=Log[(10^(-2.))]; (* initial ln of the lithium fraction, since this is a diff eq for log f this should not matter *)

lnfODEnumGyr[M_,rhoDM_]:=Module[{RHSlitODEnumlist,RHSlitODEnumint,lnfsol,lnf,survGyr},
RHSlitODEnumlist=Table[Flatten[{tt,Block[{t=(10^tt)},Evaluate[RHSlitODEnum[M,
(10^6)TcCoulomb[M,ψ[t]/.ψsolhybridStar[M,rhoDM,50*(10^5),icbench]][[1]],
centralDensity[M,ψ[t]/.ψsolhybridStar[M,rhoDM,50*(10^5),icbench]][[1]],
Radius[M,ψ[t]/.ψsolhybridStar[M,rhoDM,50*(10^5),icbench]][[1]]Rsun/cm] ]]}],{tt,3,10,.5}];
RHSlitODEnumint[t_]=Interpolation[RHSlitODEnumlist,InterpolationOrder->1][Log10[t]];
lnfsol=NDSolve[{lnf'[t]==(60^2)24 365RHSlitODEnumint[t] ,lnf[(10^3)]==initLNF},lnf,{t,(10^3),3 (10^10)}][[1]];
survGyr=Evaluate[((E^lnf[(10^9)])/(1(E^initLNF)))/.lnfsol];
PrintTemporary[{M,survGyr}];
survGyr
]//Quiet

smlistGyr=Monitor[Table[{Mi,lnfODEnumGyr[Mi,0]},{Mi,0.045,0.08,0.0025}],Mi]

(*lnfODEnumlist={lnfODEnum[Mbench,0],lnfODEnum[Mbench,(10^5)],lnfODEnum[Mbench,(10^6)]}*)

rho1d3Gyr=Monitor[Table[{Mi,lnfODEnumGyr[Mi,(10^3)]},{Mi,0.045,0.07,0.0025}],Mi]

rho1d4Gyr=Monitor[Table[{Mi,lnfODEnumGyr[Mi,(10^4)]},{Mi,0.045,0.07,0.0025}],Mi]

rho1d5Gyr=Monitor[Table[{Mi,lnfODEnumGyr[Mi,(10^5)]},{Mi,0.045,0.075,0.0025}],Mi]

rho3d5Gyr=Monitor[Table[{Mi,lnfODEnumGyr[Mi,3 (10^5)]},{Mi,0.05,0.09,0.0025}],Mi]

rho1d6Gyr=Monitor[Table[{Mi,lnfODEnumGyr[Mi,(10^6)]},{Mi,0.065,0.11,0.0025}],Mi]

rho1d7Gyr=Monitor[Table[{Mi,lnfODEnumGyr[Mi,(10^7)]},{Mi,0.07,0.125,0.0025}],Mi]

rho1d8Gyr=Monitor[Table[{Mi,lnfODEnumGyr[Mi,(10^8)]},{Mi,0.09,0.1375,0.0025}],Mi]

rho3d5Gyr=DeleteDuplicates[Join[{{0.045,1}},rho3d5Gyr]];
rho1d6Gyr=DeleteDuplicates[Join[{{0.045,1},{0.05,1},{0.06,1}},rho1d6Gyr]];
rho1d7Gyr=DeleteDuplicates[Join[{{0.045,1},{0.05,1},{0.06,1}},rho1d7Gyr]];
rho1d8Gyr=DeleteDuplicates[Join[{{0.045,1},{0.05,1},{0.06,1},{0.07,1},{0.08,1}},rho1d8Gyr]];

coloursplus=Join[{{Black,Dashed}},colours[[;;10;;2]],{Purple,Pink}];

smintGyr[m_]=10^Interpolation[smlistGyr/.{x_,y_}:>{x,Log10[y]}][m];
rho1d3intGyr[m_]=10^Interpolation[rho1d3Gyr/.{x_,y_}:>{x,Log10[y]}][m];
rho1d4intGyr[m_]=10^Interpolation[rho1d4Gyr/.{x_,y_}:>{x,Log10[y]}][m];
rho1d5intGyr[m_]=10^Interpolation[rho1d5Gyr/.{x_,y_}:>{x,Log10[y]}][m];
rho3d5intGyr[m_]=10^Interpolation[rho3d5Gyr/.{x_,y_}:>{x,Log10[y]}][m];
rho1d6intGyr[m_]=10^Interpolation[rho1d6Gyr/.{x_,y_}:>{x,Log10[y]}][m];
rho1d7intGyr[m_]=10^Interpolation[rho1d7Gyr/.{x_,y_}:>{x,Log10[y]}][m];
rho1d8intGyr[m_]=10^Interpolation[rho1d8Gyr/.{x_,y_}:>{x,Log10[y]}][m];

LogPlot[{smintGyr[m],rho1d3intGyr[m],rho1d4intGyr[m],rho1d5intGyr[m],rho3d5intGyr[m],rho1d6intGyr[m],rho1d7intGyr[m],rho1d8intGyr[m]},{m,0.045,0.14},FrameLabel->{mat["M/\\rm M_\\odot"],mat["\\text{$^7$Li survival after 1 Gyr}"]},PlotStyle->coloursplus,PlotRange->{{0.045,0.14},{(10^(-4)),2 (10^0)}},ImageSize->400,
PlotLegends->Placed[LineLegend[coloursplus,{mat[" 0"],
mat["10^3 "],mat["10^4"],mat["10^5"],mat["3 \\times 10^5"],
mat["10^6"],
mat["10^7 "],mat["10^8 "]},Spacings->0,Background->White,LegendLabel->mat["\\rho_{\\rm DM}\\left(\\rm GeV \\, cm^{-3}\\right)"],LegendFunction->(Framed[#1, FrameMargins -> 0] & ),LegendLayout->{"Column",4},LegendMarkerSize->20],{Right,Bottom}]]

Validity of the polytrope model

Solving the steady state equation

eqnDM[ρ_,v_,M_,ψ_]=Lsurf[M,ψ]-LDM[ρ,v,M,ψ]-LHB[M,ψ]==0;
ψsolutionDM[M_?NumericQ,r_?NumericQ,IC_]:=ψ/.FindRoot[eqnDM[r,50*10^5,M,ψ],{ψ,IC},WorkingPrecision->30,AccuracyGoal->10,PrecisionGoal->10,MaxIterations->Infinity]

PcouloverPc[M_,ρ_]=-(Pcoulomb[centralDensity[M,ψsolutionDM[M,ρ,0.5]],(10^6)Tc[M,ψsolutionDM[M,ρ,0.5]]]/centralPressure[M,ψsolutionDM[M,ρ,0.5]])/.g->1/.s->1/.cm->1//Quiet;

Pcoultable=Table[Table[{M,PcouloverPc[M,(10^(logρ))]},{M,0.0005,0.1,0.001}],{logρ,{1,3,5}}]//Chop//Quiet;

P_Coulomb/P_polytrope. The n=1.5 polytropic approximation breaks down then the lines cross unity. Interestingly, the n=1.5 approximation is better for higher DM densities.

linecolors={Black,Red,Blue};
legend=
LineLegend[linecolors,{mat["10{\\rm GeV/cm}^3"],mat["10^3{\\rm GeV/cm}^3"],mat["10^5{\\rm GeV/cm}^3"]},Spacings->0,LegendLabel->mat["\\rho_{\\rm DM}"]];

Show[ListLogLogPlot[Pcoultable,Joined->True,PlotRange->{{0.005,.065},{6 (10^(-2)), 5}},Frame->True,FrameTicks->Automatic,FrameLabel->{mat["M/{\\rm M}_\odot"],mat["|P_{\\rm C}|/P_c"]},PlotStyle->linecolors,FrameStyle->Directive[Black,FontSize->15],PlotLegends->Placed[legend,{Left,Bottom}]],LogLogPlot[1,{M,0.001,1},PlotStyle->{Gray,Dashed}]]//Quiet

Breakdown of the approximation in the SM

Pcoulint[m_]=10^Interpolation[Transpose[{Pcoultable[[1,;;,1]],Log10[Pcoultable[[1,;;,2]]]}]][m];
m/.NSolve[Pcoulint[m]==1,m]