Axial QNMs of NS

for Viscous Polytropic Stars using Parametrization A

Quit[]

<<EDCRGTCcode.m

The Star

densities={0.0007429420505200595`,0.0011144130757800893`, 0.001485884101040119`,0.0018573551263001487`, 0.0022288261515601786`, 0.0026002971768202083`,0.002971768202080238`,0.0033432392273402677`, 0.0037147102526002974`};
L=2;
ρ0=SetPrecision[densities[[5]], 70*2];
L0=R;

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, dinv},
ruleP={WorkingPrecision-> 2*prec-2, AccuracyGoal->prec, PrecisionGoal-> prec, Method->"StiffnessSwitching"};
eos[ρ_]:=κ (ρ^(1+(1/n)));
inv[p_]:=(p/κ)^(n/(n+1));
dinv[p_]:=D[inv[p], p];
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]],
"invcs2"->Function[x, N[dinv[p[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],
"(Subscript[cs2, 0])"->N[κ (1+(1/n))((ρ0)^(1/n)), prec]
|>;
Return[tov];
];

Regularity and BCs

TOV=TOVsolverPoly[1, 100,ρ0];

etahat = SetPrecision[0.00012385502154825165/((TOV["(Subscript[p, 0])"]+ρ0) TOV["(Subscript[R, s])"]TOV["(Subscript[cs2, 0])"]) , 50]
viscous$params={ηh->69029706/100000000,τh->10};

Some checks regarding the causality of the  frame

(*Check that this*)
cs20=TOV["(Subscript[cs2, 0])"]
(*is smaller than*)
N[((τh)/(2+τh))/.viscous$params]

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π(ρ[r]-p[r]));


a1=((-1(*(E^(ν[r]/2))*))/r)f[r](D[Log[η[r]], r]+(1/2)D[ν[r],r]);
cn=Sqrt[((η[r])/(τ[r](p[r]+ρ[r])))];
U=Exp[ν[r]](((L(L+1))/(r^2))-((2m[r])/(r^3))+8 π (2p[r]+ρ[r]));
c2=((f[r]^2)/(2r))(Exp[λ[r]](1+8π r^2p[r])-1-2((r η'[r])/(η[r])));
c3=-I ω (p[r]+ρ[r]) Exp[((ν[r])/2)]((1/(η[r]))+16π τ[r])+2f[r]^2((η'[r])/(r η[r]));
c4=r f[r](I ω+(((p[r]+ρ[r]))/(η[r]))(Exp[((ν[r])/2)]-I ω τ[r]));
c5=f[r]((((p[r]+ρ[r])Exp[((ν[r])/2)])/(η[r]))-I (ω/2)(-7+Exp[λ[r]](1+8π r^2 p[r]))+I (ω/(η[r]))(r η'[r]-(p[r]+ρ[r])τ[r]));

eqin1 = (f[r]D[f[r]ψ'[r] ,r]+((ω^2)-V[r])ψ[r]-16 Pi η[r] Exp[ν[r]/2](-I ω ψ[r]+a1 β[r]))//Simplify

eqin2=(E^(λ[r]-ν[r]))(f[r]D[f[r]β'[r], r]+(cn^(-2) ω^2 -U)β[r]-c2  β'[r]-c3 β[r]-c4 ψ'[r]-c5 ψ[r])//Simplify

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])"];

((p0+ρ0)R cs20)

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];
ηnew=Interpolation[Join[{{0,((p0+ρ0) L0 cs20 ηh)/. viscous$params}},Table[{r,((TOV["P"][r]+TOV["ρ"][r]) L0 (1/(TOV["invcs2"][r])) ηh)/. viscous$params},{r,grid}]],InterpolationOrder->3];
τnew=( L0 ηh τh)/. viscous$params;

nOrd=2;
mseries=Chop[Normal@Series[mnew[r],{r,0,2*nOrd-1}], 10^(-6)]
ν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)]
ηseries=Chop[Normal@Series[ηnew[r], {r,0,2*nOrd-2}], 10^(-8)]
τseries=τnew

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

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

p[r_]:=TOV["P"][r];
ρ[r_]:=TOV["ρ"][r];
m[r_]:=TOV["M"][r];
ν[r_]:=TOV["ν"][r];
λ[r_]:=TOV["λ"][r];
cs2[r_]:=(1/(TOV["invcs2"][r]));
η[r_]:=cs2[r](p[r]+ρ[r])L0 ηh /.viscous$params; 
τ[r_]:=L0 ηh τh/.viscous$params;

cs20(p0+ρ0)L0

Plot[cs2[r], {r, 0, R}]

η[r]/.r->R

(( (p[r]+ρ[r]) )/(η[r]))/.r->(R-10^(-4))

cs2[r]/.r->(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/(2M))-1))^(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[ϕ, ψ]
phi2NB=Coefficient[phieq,ϕ''[v]]//Together//Expand
phi1NB=Coefficient[phieq,ϕ'[v]]//Together//Expand
phi0NB=Coefficient[phieq,ϕ[v]]//Together//Expand

α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_:14/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]]
]
Clear[ψ]

Integration and QNMs

M=TOV["(Subscript[M, s])"];
Z=Sqrt[1-2M/R];
surface$condition[{psi_,psip_},{beta_,betap_}, ω_]:=-Z I ω beta+(M/R)ηh f[R]betap+((M(5M-2R)ηh)/(R^3))beta-I ω((M ηh)/Z)psi+R f[R]psip+(Z^2)psi

M=TOV["(Subscript[M, s])"];
CS=(M/R);
SB1[ω_]:=((-4 CS^2ηh + 2CS ηh + R ω (ηh τh R ω + I Sqrt[1-2 CS]))/(ηh R^2 (1-2CS)^2));
SB2=(CS/(2 R (1-2 CS) ));
SB3[ω_]:=((I(I(1-2 CS)^(3/2) + ηh R ω(- 2CS  τh + τh + CS) ))/((1-2CS)^2 ηh R));
SB4[ω_]:=((-1)/(Sqrt[1-2 CS]ηh))+((I τh R ω)/(1-2 CS));
NEWsurface$condition[{psi_,psip_},{beta_,betap_}, ω_]:=(SB1[ω] beta+SB2 betap+SB3[ω] psi+SB4[ω] psip) /. viscous$params

int[ω0_?NumericQ, a0_:14/10, rmin_:(10^(-6)), prec_:14]:=Module[
{params,r0, surf1, surf2,ψin1,ψin2, SubsRules1, SubsRules2, ic1, ic2, βin1, βin2, ψsur, ψpsur, b, out, ψout, ψpout, ruleP, ruleP2, eps},
eps=0;
r0=rmin*R;
params={ω->SetPrecision[ω0, 2*prec], a-> a0*R};
SubsRules1=Flatten@{(Subscript[an, 0])->1,(Subscript[bn, 0])->0};
SubsRules2=Flatten@{(Subscript[an, 0])->0,(Subscript[bn, 0])->1};
ruleP={WorkingPrecision-> prec+8, AccuracyGoal->prec, PrecisionGoal-> prec, Method->"StiffnessSwitching"};
ruleP2={WorkingPrecision-> prec, AccuracyGoal->prec, PrecisionGoal-> prec, Method->"StiffnessSwitching"};
ic1={{ψ[rmin*R], β[rmin*R]}== ini1/.SubsRules1, {ψ'[rmin*R], β'[rmin*R]}==ini2/.SubsRules1}/.params/.r->rmin*R;
ic2={{ψ[rmin*R], β[rmin*R]}== ini1/.SubsRules2, {ψ'[rmin*R], β'[rmin*R]}==ini2/.SubsRules2}/.params/.r->rmin*R; 
{ψin1,βin1}=NDSolveValue[Flatten@{{eqin1==0,eqin2==0}/.params/.SubsRules1,ic1},{ψ,β},{r,rmin*R,(1-eps)R}, ruleP]; 
{ψin2,βin2}=NDSolveValue[Flatten@{{eqin1==0,eqin2==0}/.params/.SubsRules2,ic2},{ψ,β},{r,rmin*R,(1-eps)R}, ruleP];


surf1=(NEWsurface$condition[{ψin1[(1-eps)R],ψin1'[(1-eps)R]},{βin1[(1-eps)R],βin1'[(1-eps)R]},ω]/.params/. viscous$params);
	surf2=(NEWsurface$condition[{ψin2[(1-eps)R],ψin2'[(1-eps)R]},{βin2[(1-eps)R],βin2'[(1-eps)R]},ω]/.params/. viscous$params);
b=-surf1/surf2;
{ψsur, ψpsur} = SetPrecision[{ψin1[(1-eps)R]+b ψin2[(1-eps)R], ψin1'[(1-eps)R]+b ψin2'[(1-eps)R]}, Precision[b]]; 

{ψout,ψpout}=NDSolveValue[Flatten@{eqout==0/. params,ψ[(1-eps) R]==ψsur,ψ'[(1-eps) R]==ψpsur},{ψ[a0*R],ψ'[a0*R]},{r,(1-eps) R,a0*R},ruleP2];


Return[((a0*R (ψpout - ψout * ((χ'[a0*R])/(χ[a0*R]))))/(ψout))]/.params];


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

qnmTimed[ω0_?NumericQ,a0_:14/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->40],om]];
{res,time}]

R/M

{w,tsec}=qnmTimed[0.19315817322988202-0.09845929530589I]

w=qnm[0.155-0.047I]

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

M
R