Axial QNMs of NS

for Polytropic Non - Viscous Stars

<<EDCRGTCcode.m

Quit[]

The star

densities={0.0007429420505200595`,0.0011144130757800893`, 0.001485884101040119`,0.0018573551263001487`, 0.0022288261515601786`, 0.0026002971768202083`,0.002971768202080238`,0.0033432392273402677`, 0.0037147102526002974`};

L=2; (*mode*)
ρ0=SetPrecision[densities[[5]], 70*2]; (*sets compacticity*)

TOVsolverPoly[n_,κ_,ρ0_,ri_:(10^(-6)),prec_:50,rmax_:100,pmin_:(10^(-12))]:=Module[
{p0,initialconds,eqs,ruleP,R,M,nuR,PolyDifSol,difsol,m,p,ρ,ν,tov, eos, inv},
ruleP={WorkingPrecision-> 2*prec-2, AccuracyGoal->prec, PrecisionGoal-> prec, Method->"StiffnessSwitching"};
eos[ρ_]:=κ (ρ^(1+(1/n)));
inv[p_]:=(p/κ)^(n/(n+1));
p0=eos[ρ0];
eqs={ν'[r]==Re[((2 m[r]+8 π (r^3) p[r])/((r^2)-2 r m[r]))],
m'[r]==Re[4 π (r^2) ρ[r]],
p'[r]==Re[-(((m[r]+4 π (r^3) p[r]) (p[r]+ρ[r]))/(r (r-2 m[r])))],
WhenEvent[p[r]<p0*pmin,R=r;"StopIntegration"]}/.ρ[r]->(Abs[y]/. First@Solve[eos[y]==z,y]/.z->p[r]);
initialconds={m[ri]==Re@N[4/3*Pi*(ri^3)*ρ0,prec*2],ν[ri]==Re@N[1,prec*2],p[ri]==Re@N[p0,prec*2](*m[ri]==Re@N[0,prec*2],ν[ri]==Re@N[1,prec*2],p[ri]==Re@N[p0,prec*2]*)
};

difsol=First@NDSolve[Flatten[{eqs,initialconds}],
{ν[r],m[r],p[r]},{r,ri,rmax},ruleP];

M=m[r]/.difsol/.{r->R};
nuR=ν[r]/.difsol/.{r->R};

tov = <|
	"M"->Function[x,N[ m[r]/.difsol/.r->x,prec]],
	"P"-> Function[x, N[p[r]/.difsol/.r->x,prec]],
(*"ρ"-> Function[x,N[(y/. First@Solve[eos[y]==z,y]/.z->p[r])/.difsol/.r->x,prec]],*)
"ρ"->Function[x,N[inv[p[r]]/. difsol/. r->x,prec]],
	"ν"-> Function[x,N[ν[r]-nuR+Log[1-2M/R]/.difsol/.r->x,prec]],
"λ"->Function[x, N[-Log[1-2((m[r])/r)]/.difsol/.r->x,prec]],
	"(Subscript[R, s])"-> N[R,prec],
	"(Subscript[M, s])"->N[M,prec],
"(Subscript[ν, 0])"->N[1-nuR+Log[1-2M/R],prec],
"(Subscript[ρ, 0])"->N[ρ0,prec],
"(Subscript[p, 0])"->N[p0,prec],
"(Subscript[λ, 0])"->N[-Log[1-8/3*Pi*(ri^2)*ρ0], prec],
"(Subscript[m, 0])"->N[4/3*Pi*(ri^3)*ρ0,prec]
|>;
Return[tov];
];

Regularity and BCs

TOV=TOVsolverPoly[1, 100,  ρ0]; (*sets the family of polytropic stars and solves TOV*)

Inside

Clear[eqin1, eqin2, ψ, f]
f[r_]:= Exp[ν[r]/2-λ[r]/2]
V[r_] := Exp[ν[r]](((L(L+1))/(r^2))-((6m[r])/(r^3))-4 π (p[r]-ρ[r]));


eqin = f[r]D[f[r]ψ'[r] ,r]+((ω^2)-V[r])ψ[r];

p0=TOV["(Subscript[p, 0])"];
M=TOV["(Subscript[M, s])"]
R=TOV["(Subscript[R, s])"]
ρ0=TOV["(Subscript[ρ, 0])"];
ν0=TOV["(Subscript[ν, 0])"];
λ0=TOV["(Subscript[λ, 0])"];
m0=TOV["(Subscript[m, 0])"];

grid=Join[Range[10^-6,R-10^-3,10^-3],{R}];
mnew=Interpolation[Join[{{0,0}},Table[{r,TOV["M"][r]},{r,grid}]],InterpolationOrder->3];
νnew=Interpolation[Join[{{0,ν0}},Table[{r,TOV["ν"][r]},{r,grid}]],InterpolationOrder->3];
λnew=Interpolation[Join[{{0,λ0}},Table[{r,TOV["λ"][r]},{r,grid}]],InterpolationOrder->3];
pnew=Interpolation[Join[{{0,p0}},Table[{r,TOV["P"][r]},{r,grid}]],InterpolationOrder->3];
ρnew=Interpolation[Join[{{0,ρ0}},Table[{r,TOV["ρ"][r]},{r,grid}]],InterpolationOrder->3];

nOrd=2;
mseries=Chop[Normal@Series[mnew[r],{r,0,2*nOrd-1}], 10^(-4)]
νseries= Chop[Normal@Series[νnew[r],{r,0,2*nOrd-2}], 10^(-5)]
λseries=Chop[Normal@Series[λnew[r],{r,0,2*nOrd-2}], 10^(-6)]
pseries=Chop[Normal@Series[pnew[r],{r,0,2*nOrd-2}], 10^(-7)]
ρseries=Chop[Normal@Series[ρnew[r], {r,0,2*nOrd-2}], 10^(-6)]

seriesRule=Flatten@{FuncRepRules[m[r],mseries], FuncRepRules[ν[r],νseries], FuncRepRules[λ[r],λseries], FuncRepRules[p[r],pseries], FuncRepRules[ρ[r],ρseries]};

ψ[r_]:=(r^(L+1)) φ[r];
φ[r_]:=Sum[(Subscript[an, n]) (r^(2n)),{n,0,nOrd}]/. (Subscript[an, 0])->1 ; 
eqreduced={(eqin/(r^(L+1)))}/.seriesRule;
Chop[Series[eqreduced,{r,0,2nOrd-2}],(10^(-15))];
coefs=First@Solve[%==0,Flatten@Table[(Subscript[an, n]),{n,1,nOrd}]]//Simplify
{ini1, ini2} = {ψ[r]/.coefs, ψ'[r]/.coefs}//Simplify;
Precision[%]
Clear[ψ, φ, coefs];

p[r_]:=TOV["P"][r];
ρ[r_]:=TOV["ρ"][r];
m[r_]:=TOV["M"][r];
ν[r_]:=TOV["ν"][r];
λ[r_]:=TOV["λ"][r];

Outside

Clear[M];

eqout = (1-2(M/r))D[(1-2(M/r))ψ'[r],r]+((ω^2)-Vout[r])ψ[r];
Vout[r_] := (1-2(M/r))(((L(L+1))/(r^2))-6(M/(r^3)));

χ[r_] := ((((r/1)-2M))^(2M I ω)) Exp[I ω r]; ψ[r_] := χ[r] ϕ[1-(a/r)];
phieq=Collect[(((a^3) (((a+2 M (-1+v)))^2))/( χ[r](-1+v)^2 (((a+2 M (-1+v)))^3) ))eqout  /. r-> ((-a)/(v-1)), {ϕ[v], ϕ'[v], ϕ''[v]},FullSimplify]
 {f0,f1,f2,f3}=CoefficientList[Coefficient[phieq, ϕ''[v]],v]
 {d0,d1,d2}=CoefficientList[Coefficient[phieq, ϕ'[v]],v]
{e0,e1}=CoefficientList[Coefficient[phieq, ϕ[v]],v]
Clear[ϕ, ψ]

αn = f0 n (n+1);
βn = d0 n + f1 n (n-1);
γn = e0 + d1 (n-1) + f2 (n-1)(n-2);
δn = e1 + d2(n-2) + f3 (n-2)(n-3);

elim[ω0_?NumericQ, a0_:18/10,N0_:30]:=Module[
{params,alpha,beta,gamma,delta,alp,bet,gam,Rp},
params = {ω->ω0, a->a0*R};
alpha[n0_] := αn/.params/.n->n0;
beta[n0_]:=βn/.params/.n->n0;
gamma[n0_]:=γn/.params/.n->n0;
delta[n0_] := δn/.params/.n->n0;

Module[{k},{
For[k=0, k<N0+2,{
alp[k]= alpha[k];
If[k<2,
bet[k]=beta[k]; gam[k]=gamma[k];,
bet[k] = beta[k]-(( delta[k] alp[k-1])/(gam[k-1])); 
gam[k] = gamma[k]-(( delta[k]bet[k-1])/(gam[k-1] ));
];
k++;
}];
}];

Rp[N0+1]=0;
Do[
Rp[n] =- (( gam[n+1])/(bet[n+1]+alp[n+1]Rp[n+1]));,
{n,N0,0,-1}
];

Return[Simplify@Rp[0]]
]

Integration and QNMs

M=TOV["(Subscript[M, s])"];

int[ω0_?NumericQ, a0_:18/10, rmin_:(10^(-6)), prec_:16]:=Module[
{params,ic,sol1,ψsur,ψ2sur,sol2, ψ1,ψ2, ruleP, ruleP2},
params={ω->SetPrecision[ω0, 2*prec], a-> a0*R};
ruleP={WorkingPrecision-> prec+8, AccuracyGoal->prec, PrecisionGoal-> prec, Method->"StiffnessSwitching"};
ruleP2={WorkingPrecision-> prec, AccuracyGoal->prec, PrecisionGoal-> prec, Method->"StiffnessSwitching"};
ic={ψ[rmin*R]== ini1, ψ'[rmin*R]==ini2}//.params //. r->rmin*R;
sol1=First@NDSolve[Flatten@{eqin==0/.params,ic}, ψ,{r,rmin*R,R}, ruleP];
{ψsur, ψ2sur} = {ψ[R], ψ'[R]}/.sol1;
sol2 = First@NDSolve[Flatten@{eqout==0/.params,ψ[R]==ψsur, ψ'[R]==ψ2sur}, ψ,{r,R,a0*R}, ruleP2];
{ψ1, ψ2} = {ψ[a0*R], ψ'[a0*R]}/.sol2;
Return[((a0*R (ψ2 - ψ1 * ((χ'[a0*R])/(χ[a0*R]))))/(ψ1))/.params]
];


qnm[ω0_?NumericQ, a0_:18/10, rmin_:(10^(-6)), N0_:30]:=om/.Monitor[FindRoot[int[om,a0,rmin]-elim[om,a0,N0],{om,ω0}],om]

qnmTimed[ω0_?NumericQ,a0_:18/10,rmin_:10^-6,N0_:30]:=Module[{res,time},{time,res}=AbsoluteTiming[om/. Monitor[FindRoot[int[om,a0,rmin]-elim[om,a0,N0],{om,ω0,ω0+0.01+0.005 I},Method->"Secant",MaxIterations->25],om]];
{res,time}]

w=qnm[0.198-0.08I]

M

R

M/R

int[w,18/10,(10^(-6))]-elim[w,18/10,30] (*check*)