Interfacial modes of neutron stars

ClearAll["`*"]

1. Spacetime

$Assumptions={r≥0,θ≥0,θ≤π,r>2 M[r]&&r>0,ν[r]∈Reals};
x={t,r,θ,ϕ};
δ={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
dim=4;
order =1;
P2[θ_]:=((3(1-X[θ])-1)/2);
(*g0={{-Exp[ν[r]],0,0,0},{0,(((1-2((M[r])/r)))^(-1)),0,0},{0,0,(r^2),0},{0,0,0,(r^2) X[θ]}};*)
g0={{-Exp[ν[r]],0,0,0},{0,Exp[λ[r]],0,0},{0,0,(r^2),0},{0,0,0,(r^2) X[θ]}};
rule1[Z_]:=Z/.Union[{((X'[θ])^2)->4 X[θ](1-X[θ]),X''[θ]->2 (1-2X[θ])}]
(*ξ0= (r^2){0,W[r](r^(-1))*(((1-2((M[r])/r)))^(1/2))P2[θ]*ϵ,-V[r](r^(-2))D[P2[θ],θ]*ϵ,0}Exp[I*ω*t];*)
ξ0= (r^2){0,W[r](r^(-1))*Exp[-λ[r]/2]*P2[θ]*ϵ,-V[r](r^(-2))D[P2[θ],θ]*ϵ,0}Exp[I*ω*t];
(*ξ0={0,W[r](r^(-2))*Exp[-λ[r]/2]*P2[θ]*ϵ,-V[r](r^(-2))D[P2[θ],θ]*ϵ,0}Exp[I*ω*t];*)
h=-(r^2){{H0[r]*Exp[ν[r]],I*ω*r*H1[r],0,0},{I*ω*r*H1[r],H2[r]*Exp[λ[r]],0,0},{0,0,(r^2)*K[r],0},{0,0,0,(r^2)*X[θ]*K[r]}}*P2[θ]Exp[I*ω*t]*ϵ;
(*h={{H0[r]*Exp[ν[r]],I*ω*r*H1[r],0,0},{I*ω*r*H1[r],H2[r]*Exp[λ[r]],0,0},{0,0,(r^2)*K[r],0},{0,0,0,(r^2)*X[θ]*K[r]}}*P2[θ]Exp[I*ω*t]*ϵ;*)
gdd = g0 +h;
guu = Simplify[rule1[Normal[Series[Inverse[gdd],{ϵ,0,order}]]]];
Do[ξ[i]=ξ0[[i]],{i,1,dim}];
u0={Exp[-ν[r]/2],0,0,0};
ud0={-Exp[ν[r]/2],0,0,0};
Do[u[i]=u0[[i]],{i,1,dim}];
Do[ud[i]=ud0[[i]],{i,1,dim}];
Do[g[i,j]=Simplify[gdd[[i]][[j]]];,{i,1,dim},{j,1,dim}];
Do[ig[i,j]=Simplify[guu[[i]][[j]]];,{i,1,dim},{j,1,dim}];
Do[δg[i,j]=Simplify[h[[i]][[j]]];,{i,1,dim},{j,1,dim}];
Do[ds[i]=x[[i]];,{i,1,dim}];
Do[gud[i,j]=Simplify[Normal[Series[Sum[g[i,k]*ig[k,j],{k,1,dim}],{ϵ,0,order}]]];,{i,1,dim},{j,1,dim}];
Do[Chr[i,j,k]=Simplify[rule1[Normal[Series[(1/2)*Sum[ig[i,s]*(-D[g[j,k],ds[s]]+D[g[j,s],ds[k]]+D[g[s,k],ds[j]]),{s,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim},{k,1,dim}];
Do[Ruddd[i,k,l,m]= Simplify[rule1[Normal[Series[D[Chr[i, k, m], ds[l]]
       - D[Chr[i, k, l], ds[m]]
       + Sum[Chr[i, s, l]*Chr[s, k, m]-Chr[i, s, m]*Chr[s, k, l], {s, 1, dim}], {ϵ, 0, order}]]]]; , {i, 1, dim}, {k, 1, dim}, {l, 1, dim}, {m, 1, dim}];

Do[Ricci[i,k]=Simplify[rule1[Normal[Series[Sum[Ruddd[l,i,l,k],{l,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{k,1,dim}];
Do[Ruu[i,k]=Simplify[rule1[Normal[Series[Sum[Ricci[l,m]*ig[i,l]*ig[k,m],{l,1,dim},{m,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{k,1,dim}];
Rs=Simplify[rule1[Normal[Series[Sum[ig[i,j]*Ricci[i,j],{i,1,dim},{j,1,dim}],{ϵ,0,order}]]]];
Do[Gdd[i,j]=Simplify[rule1[Normal[Series[Ricci[i,j]-(1/2) g[i,j]Rs,{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[Gud[i,k]=Simplify[rule1[Normal[Series[Sum[Gdd[i,j]*ig[j,k],{j,1,dim}],{ϵ,0,order}]]]];,{i,1,4},{k,1,4}]
Do[Lie[i]=rule1[Normal[Series[Sum[u[j]*D[ξ[i],ds[j]]-ξ[j]*D[u[i],ds[j]],{j,1,dim}],{ϵ,0,order}]]];,{i,1,dim}]

Do[q[i,j]=ig[i,j]+u[i]*u[j];,{i,1,dim},{j,1,dim}];
Do[ξd[i]=Simplify[rule1[Normal[Series[Sum[g[i,j]*ξ[j],{j,1,dim}],{ϵ,0,order}]]]];,{i,1,dim}];
(*Lagrangian perturbation of the metric*)
Do[Δg[i,j]=Simplify[rule1[Normal[Series[δg[i,j]+ D[ξd[i],ds[j]]+ D[ξd[j],ds[i]] -2Sum[Chr[k,i,j]*ξd[k],{k,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];

Matter: fluid component

(*restmass fraction Lagrangian perturbation*)
Δnb=FullSimplify[rule1[Normal[Series[((-1)/2)*Sum[q[i,j]*Δg[i,j],{i,1,dim},{j,1,dim}],{ϵ,0,order}]]]];
(*perutrbation of four-velocity*)
Do[qud[i,j]=Simplify[rule1[Normal[Series[Sum[g[i,k]*q[k,j],{k,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[δu[i]=Simplify[rule1[Normal[Series[Sum[qud[i,j]*Lie[j],{j,1,dim}],{ϵ,0,order}]]+Normal[Series[Sum[(1/2)*u[i]*u[j]*u[k]*δg[j,k],{j,1,dim},{k,1,dim}],{ϵ,0,order}]]]];,{i,1,dim}];
Do[δud[i]=Simplify[Normal[Series[Sum[g[i,j]*δu[j],{j,1,dim}],{ϵ,0,order}]]];,{i,1,dim}];
δρ=Simplify[(ρ[r]+p[r])*Δnb-D[ρ[r],r]*ξ[2]];
δp=Simplify[γ*p[r]*Δnb-D[p[r],r]*ξ[2]];
Do[δ1[i,j]=δ[[i]][[j]];,{i,1,dim},{j,1,dim}];
Do[δTud[i,j]=Simplify[rule1[Normal[Series[(δρ+δp)*ud[i]*u[j]+δp*δ1[i,j]+(ρ[r]+p[r])*(u[j]*δud[i]+ud[i]*δu[j]+Sum[u[j]*u[k]*δg[k,i],{k,1,dim}]),{ϵ,0,1}]]]];,{i,1,dim},{j,1,dim}];
Do[Tud[i,j]=Simplify[rule1[Normal[Series[(ρ[r]+p[r])*(ud[i]*u[j])+δ1[i,j]*p[r]+δTud[i,j],{ϵ,0,1}]]]];,{i,1,dim},{j,1,dim}];
Do[Einstein[i,j]=Simplify[rule1[Normal[Series[Gud[i,j]- 8*π*Tud[i,j],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[DivT[i]=Simplify[rule1[Normal[Series[Sum[D[Tud[i,j],ds[j]],{j,1,dim}]+Sum[Chr[j,j,k]*Tud[i,k]-Chr[k,i,j]*Tud[k,j],{k,1,dim},{j,1,dim}],{ϵ,0,order}]]]];,{i,1,dim}];
eq11=Einstein[1,1];
eq12=Einstein[1,2];
eq13=Einstein[1,3];
eq14=Einstein[1,4];
eq22=Einstein[2,2];
eq23=Einstein[2,3];
eq24=Einstein[2,4];
eq33=Einstein[3,3];
eq34=Einstein[3,4];
eq44=Einstein[4,4];
div1=DivT[1];
div2=DivT[2];
div3=DivT[3];
div4=DivT[4];
SetDirectory[NotebookDirectory[]];
DumpSave["fluid.mx",{eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4}];

Matter: solid component

(*restmass fraction Lagrangian perturbation*)
Δnb=FullSimplify[rule1[Normal[Series[((-1)/2)*Sum[q[i,j]*Δg[i,j],{i,1,dim},{j,1,dim}],{ϵ,0,order}]]]];
(*perutrbation of four-velocity*)
Do[qud[i,j]=Simplify[rule1[Normal[Series[Sum[g[i,k]*q[k,j],{k,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[δu[i]=Simplify[rule1[Normal[Series[Sum[qud[i,j]*Lie[j],{j,1,dim}],{ϵ,0,order}]]+Normal[Series[Sum[(1/2)*u[i]*u[j]*u[k]*δg[j,k],{j,1,dim},{k,1,dim}],{ϵ,0,order}]]]];,{i,1,dim}];
Do[δud[i]=Simplify[Normal[Series[Sum[g[i,j]*δu[j],{j,1,dim}],{ϵ,0,order}]]];,{i,1,dim}];
x1=FullSimplify[ -γ*p[r]*Δnb/(r^2) * Exp[ν[r]/2]/P2[θ]/Exp[I*ω*t]]/ϵ;
δρ=Simplify[(ρ[r]+p[r])*Δnb-D[ρ[r],r]*ξ[2]];
(*δp=Simplify[γ*p[r]*Δnb-D[p[r],r]*ξ[2]];*)
δp=Simplify[-X1[r]*(r^2)*Exp[-ν[r]/2]*P2[θ]*Exp[I*ω*t]*ϵ-D[p[r],r]*ξ[2]];
Do[δ1[i,j]=δ[[i]][[j]];,{i,1,dim},{j,1,dim}];
Do[δsud[i,j]=Simplify[rule1[Normal[Series[Sum[(1/2)(qud[i,k]*q[l,j]-(1/3)qud[i,j]*q[k,l])*Δg[k,l],{k,1,dim},{l,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}]
Do[δπud[i,j]=-2*μ[r]*δsud[i,j];,{i,1,dim},{j,1,dim}]
t1=-Simplify[δπud[2,2]/P2[θ]/Exp[I*ω*t]]/ϵ;
t2=Simplify[δπud[2,3]/D[P2[θ],θ]/Exp[I*ω*t]]/ϵ*r;
δπud[2,2]=-T1[r]*P2[θ]*Exp[I*ω*t]*ϵ;
δπud[2,3]=T2[r]*D[P2[θ],θ]*Exp[I*ω*t]*ϵ/r;
δπud[3,2]=T2[r]*D[P2[θ],θ]*Exp[I*ω*t]*ϵ/r * (r^2)*Exp[-λ[r]];
δπud[3,3]=(1/2)*T1[r]*P2[θ]*Exp[I*ω*t]*ϵ+3Exp[I*ω*t]ϵ V[r]* X[θ]* μ[r];
δπud[4,4]=(1/2)*T1[r]*P2[θ]*Exp[I*ω*t]*ϵ-3Exp[I*ω*t]ϵ V[r]* X[θ] *μ[r];
Do[δTud[i,j]=Simplify[rule1[Normal[Series[(δρ+δp)*ud[i]*u[j]+δp*δ1[i,j]+(ρ[r]+p[r])*(u[j]*δud[i]+ud[i]*δu[j]+Sum[u[j]*u[k]*δg[k,i],{k,1,dim}]),{ϵ,0,1}]]+δπud[i,j]]];,{i,1,dim},{j,1,dim}];
Do[Tud[i,j]=Simplify[rule1[Normal[Series[(ρ[r]+p[r])*(ud[i]*u[j])+δ1[i,j]*p[r]+δTud[i,j],{ϵ,0,1}]]]];,{i,1,dim},{j,1,dim}];
Do[Einstein[i,j]=Simplify[rule1[Normal[Series[Gud[i,j]- 8*π*Tud[i,j],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[DivT[i]=Simplify[rule1[Normal[Series[Sum[D[Tud[i,j],ds[j]],{j,1,dim}]+Sum[Chr[j,j,k]*Tud[i,k]-Chr[k,i,j]*Tud[k,j],{k,1,dim},{j,1,dim}],{ϵ,0,order}]]]];,{i,1,dim}];
eq11=Einstein[1,1];
eq12=Einstein[1,2];
eq13=Einstein[1,3];
eq14=Einstein[1,4];
eq22=Einstein[2,2];
eq23=Einstein[2,3];
eq24=Einstein[2,4];
eq33=Einstein[3,3];
eq34=Einstein[3,4];
eq44=Einstein[4,4];
div1=DivT[1];
div2=DivT[2];
div3=DivT[3];
div4=DivT[4];
SetDirectory[NotebookDirectory[]];
DumpSave["solid.mx",{eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4,x1,t1,t2}];

(*restmass fraction Lagrangian perturbation*)
Δnb=FullSimplify[rule1[Normal[Series[((-1)/2)*Sum[q[i,j]*Δg[i,j],{i,1,dim},{j,1,dim}],{ϵ,0,order}]]]];
(*perutrbation of four-velocity*)
Do[qud[i,j]=Simplify[rule1[Normal[Series[Sum[g[i,k]*q[k,j],{k,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[δu[i]=Simplify[rule1[Normal[Series[Sum[qud[i,j]*Lie[j],{j,1,dim}],{ϵ,0,order}]]+Normal[Series[Sum[(1/2)*u[i]*u[j]*u[k]*δg[j,k],{j,1,dim},{k,1,dim}],{ϵ,0,order}]]]];,{i,1,dim}];
Do[δud[i]=Simplify[Normal[Series[Sum[g[i,j]*δu[j],{j,1,dim}],{ϵ,0,order}]]];,{i,1,dim}];
x1=FullSimplify[ -γ*p[r]*Δnb/(r^2) * Exp[ν[r]/2]/P2[θ]/Exp[I*ω*t]]/ϵ;
δρ=Simplify[(ρ[r]+p[r])*Δnb-D[ρ[r],r]*ξ[2]];
(*δp=Simplify[γ*p[r]*Δnb-D[p[r],r]*ξ[2]];*)
δp=Simplify[-X1[r]*(r^2)*Exp[-ν[r]/2]*P2[θ]*Exp[I*ω*t]*ϵ-D[p[r],r]*ξ[2]];
Do[δ1[i,j]=δ[[i]][[j]];,{i,1,dim},{j,1,dim}];
Do[δsud[i,j]=Simplify[rule1[Normal[Series[Sum[(1/2)(qud[i,k]*q[l,j]-(1/3)qud[i,j]*q[k,l])*Δg[k,l],{k,1,dim},{l,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}]
Do[δπud[i,j]=-2*μ[r]*δsud[i,j];,{i,1,dim},{j,1,dim}]

δπud[3,2]
δπud[2,3]

Quit[]

Fluid master equations

SetDirectory[NotebookDirectory[]];
<<"fluid.mx";
X[θ_]:=((Sin[θ])^2);
SS[eq_,ord_]:=Normal[Series[eq,{ϵ,0,ord}]];
EQs={eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4};

syst0=Solve[{SS[eq11,0]==0,SS[eq22,0]==0,SS[div2,0]==0},{λ'[r],ν'[r],p'[r]}][[1]]//FullSimplify
syst0=Union[syst0,D[syst0,r],D[syst0,{r,2}]];
(*syst0=Union[D[syst0,r],D[syst0,{r,2}]];*)

syst1=Solve[{SS[eq12,1]==0,SS[eq13,1]==0,SS[div3,1]==0,SS[eq33-eq44,1]==0},{H1'[r],K'[r],W'[r],H2[r]}][[1]]//FullSimplify
syst11=Union[syst1,D[syst1,r],D[syst1,{r,2}]];

Vdot=Collect[FullSimplify[Solve[(Coefficient[div2,ϵ]//.syst11)==0,V'[r]][[1]]],{H1[r],V[r],H0[r],W[r]}]

verify consistent with Kruger, Ho and Andersson, PRD (2015)

coe1=Simplify[Coefficient[FullSimplify[Coefficient[V'[r]/.Vdot,V[r]]],ν'[r]]]

coe2=Simplify[Coefficient[FullSimplify[Coefficient[V'[r]/.Vdot,V[r]]],p'[r]]]

FullSimplify[Coefficient[V'[r]/.Vdot,V[r]]-coe1*ν'[r]-coe2*p'[r]]

eqH01dot = FullSimplify[Solve[{(Coefficient[eq22,ϵ]//.Union[syst11])==0},{H0'[r]}]];
eqH02dot=Simplify[Solve[{(eq23//.Union[syst11])==0},{H0'[r]}]];
eqH0=Collect[Solve[FullSimplify[(H0'[r]/. First[eqH01dot])-(H0'[r]/. First[eqH02dot])]==0,H0[r]],{H0[r],H1[r],K[r],W[r],V[r]}][[1]]

Union[syst1,Vdot,eqH0]

Solid master equations

Quit[]

Clear[r,X];
SetDirectory[NotebookDirectory[]];
<<"solid.mx";
X[θ_]:=((Sin[θ])^2);
SS[eq_,ord_]:=Normal[Series[eq,{ϵ,0,ord}]];
EQs={eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4,x1,t1,t2};

syst0=Solve[{SS[eq11,0]==0,SS[eq22,0]==0,SS[div2,0]==0},{λ'[r],ν'[r],p'[r]}][[1]]//FullSimplify
(*syst0=Union[];*)
syst0=Union[syst0,D[syst0,r],D[syst0,{r,2}]];

syst1=Solve[{SS[eq12,1]==0,SS[eq13,1]==0,SS[eq23,1]==0,SS[eq33-eq44,1]==0,X1[r]-x1==0,T2[r]-t2==0},{H1'[r],K'[r],H0'[r],H2[r],W'[r],V'[r]}][[1]]//FullSimplify
syst11=Union[syst1,D[syst1,r],D[syst1,{r,2}]];

syst3=List[Solve[{X1[r]-x1==0,T1[r]-t1==0},{W'[r],T1[r]}][[1]][[2]]//FullSimplify]
syst2=Collect[FullSimplify[Solve[{(SS[div3,1]//.syst11)==0},{T2'[r]}][[1]]],{H0[r],V[r],W[r],X[r],T1[r],T2[r]}]

equal1=FullSimplify[(Coefficient[-eq22,ϵ]//.Union[syst11])/(-2+3 ((Sin[θ])^2))*2r /Exp[I*ω*t]==0];
eqH0=Collect[Solve[equal1,H0[r]],{H0[r],H1[r],K[r],X1[r],T1[r],T2[r],W[r]}][[1]]

Union[syst1,syst2,syst3,eqH0]