Love numbers for General Relativity [ApJ 677 : 1216-1220]

Exit[]

(* This set up the working directory of the notebook. *)
path=SetDirectory[NotebookDirectory[]];
SetDirectory[path];

<<GR`

Metric tensor

(* I define the metric here. This just corresponds to time independent polar perturbations *)

g0={{-A[r],0,0,0},{0,1/B[r],0,0},{0,0,(r^2),0},{0,0,0,(r^2)((Sin[θ])^2)}};

h={{-A[r]H0[r],H1[r],0,0},{H1[r],H2[r]/B[r],0,0},{0,0,(r^2)K[r],0},{0,0,0,(r^2)((Sin[θ])^2)K[r]}}Y[θ,ϕ];

(* ϵ is the usual bookkeping parameter for the perturbation *)
ℊ=g0+ϵ h;

ℊℐ=Normal[Series[Inverse[ℊ],{ϵ,0,1}]];

dg=Series[-Det[g],{ϵ,0,1}];

y={t,r,θ,ϕ};

ℊ//MatrixForm

Field equations and stress-energy tensor

G=Normal[Series[Einstein[ℊ,y],{ϵ,0,1}]];

(* Mixed compoenent of the Einstein tensor *)
𝒢=Normal[Series[Table[Sum[ℊℐ[[l1,l3]]G[[l3,l2]],{l3,1,4}],{l1,1,4},{l2,1,4}],{ϵ,0,1}]];

(* Define the 4-velocity of the fluid and compute the stress-energy tensor *)
u={u0,0,0,0};

u0=Normal[Series[Solve[Sum[ℊℐ[[l1,l2]]u[[l1]]u[[l2]],{l1,1,4},{l2,1,4}]==-1,u0][[2,1,2]],{ϵ,0,1}]];

T=Collect[8π Series[(ρ[r]+p[r]+ϵ ((dpdρ)^(-1))δp[r]Y[θ,ϕ]+ϵ δp[r]Y[θ,ϕ])u⊗u+ℊ (p[r]+ϵ δp[r]Y[θ,ϕ]),{ϵ,0,1}],ϵ,FullSimplify];

(* Again this is the mixed compoent form *)
𝒯=Normal[Series[Table[Sum[ℊℐ[[l1,l3]]T[[l3,l2]],{l3,1,4}],{l1,1,4},{l2,1,4}],{ϵ,0,1}]];

(* Christoffel symbols *)
Γ=Collect[Series[Chr[ℊ,y],{ϵ,0,1}],ϵ];

(* Covariant derivative of the stress enegy tensor *)
DT=Collect[Table[Sum[ℊℐ[[i,j]]D[T[[i,s]],y[[j]]],{j,1,4},{i,1,4}]-Series[Sum[ℊℐ[[i,j]]Γ[[β]][[j,i]]T[[β,s]],{β,1,4},{i,1,4},{j,1,4}]+Sum[ℊℐ[[i,j]]Γ[[β]][[s,j]]T[[i,β]],{β,1,4},{i,1,4},{j,1,4}],{ϵ,0,1}],{s,1,4}],ϵ,FullSimplify];

(* Export everything *)

Export["stressenergy.m",𝒯];

Export["EoMcontinuity.m",DT];

Export["eq_metric.m",𝒢];

Solve field’s equation order by order

𝒯=Import["stressenergy.m"];
𝒢=Import["eq_metric.m"];
𝒟T=Import["EoMcontinuity.m"];

(* Einstein's equations *)
ℰ=𝒢-𝒯;

ϵ = 0 - Background solution. This gives you the usual TOV equations

Coefficient[ℰ[[1,1]],ϵ,0];
dB=Solve[%==0,B'[r]][[1,1,2]]

Coefficient[ℰ[[2,2]],ϵ,0];
dA=Solve[%==0,A'[r]][[1,1,2]]

Coefficient[𝒟T[[2]],ϵ,0];
dp=Solve[%==0,p'[r]][[1,1,2]]

dM=Solve[D[1-((2M[r])/r),r]==(dB//.{B[r]->1-((2M[r])/r)}),M'[r]][[1,1,2]]

dp2=dp//.{A'[r]->dA,B[r]->1-((2M[r])/r)}//Simplify

(* Check for the vacuum case: it must reproduce Schwarzschild *)

dB//.ρ[r]->0;
BS=DSolve[%==B'[r],B[r],r][[1,1,2]]//.C[1]->-2M

dA//.{p[r]->0,B[r]->BS};
As=DSolve[%==A'[r],A[r],r][[1,1,2]]//.C[1]->1//Simplify

ϵ = 1 - Fist order equations in the quadrupolar perturbation

(* From this equation I find H0+H2 = 0 *)
Simplify[Coefficient[ℰ[[3,3]]-ℰ[[4,4]],ϵ,1]]

(* I impose the previous condition *)
H0[r_]:=H[r];
H2[r_]:=-H[r];

(* Find the solution for K *)
Coefficient[ℰ[[2,3]],ϵ,1]//Simplify;
dK=Solve[%==0,K'[r]][[1,1,2]]

(* Find the solution for δp *)
Coefficient[ℰ[[3,3]]+ℰ[[4,4]],ϵ,1]//.{K'[r]->dK,K''[r]->D[dK,r]}//Simplify;
solδp=Solve[%==0,δp[r]][[1,1,2]]

(* Simplify the master equation for H[r] *)
Assuming[((Csc[θ])^2) Derivative[0,2][Y][θ,ϕ]+Cot[θ] Derivative[1,0][Y][θ,ϕ]+Derivative[2,0][Y][θ,ϕ]==-6 Y[θ,ϕ],Simplify[Coefficient[ℰ[[1,1]]-ℰ[[2,2]],ϵ,1]//.{K'[r]->dK,K''[r]->D[dK,r],δp[r]->solδp}]];

Eq=Collect[-((%)/(Y[θ,ϕ]B[r]))//.{B'[r]->dB,A'[r]->dA,A''[r]->D[dA,r],p'[r]->dp},{H''[r],H'[r],H[r],r,dpdρ},FullSimplify]

Matching the solution with the esxterior

(* We need now to find the form of the equation in vacuum to make the match between the interior and exterior solution *) 
EqR=Collect[Eq//.{p[r]->0,ρ[r]->0,A[r]->1-((2M)/r),B[r]->1-((2M)/r)},{H[r],H'[r],H''[r]},FullSimplify]

(* We change the radial variable to r = M(1+x) *) 
Coefficient[EqR,H[r]]H[x]+Coefficient[EqR,H'[r]]((H'[x])/M)+Coefficient[EqR,H''[r]]((H''[x])/(M^2))//.r->M(1+x);
HR2[x_]=Collect[%(-1+(x^2))(M^2),{H[x],H'[x],H''[x]},FullSimplify]

(* Back again to r *) 
solHR=Collect[DSolve[HR2[z]==0,H[z],z][[1,1,2]]//.{Log[1-z]->Log[z-1],z->(r/M)-1},{C[1],C[2]},Simplify]

Series[solHR,{r,∞,3}]

gttACMC=-1+((2M)/r)+((3𝒬)/(r^3))-ℰT (r^2)

(* I match the solution at infinity with ACMC coordinates to extract the quadrupole *) 
c1=Solve[Coefficient[gttACMC,r,2]==Coefficient[Series[(-1+((2M)/r))(1+solHR),{r,∞,3}],r,2],C[1]][[1,1,2]]

c2=Solve[Coefficient[gttACMC,r,-3]==Coefficient[Series[(-1+((2M)/r))(1+solHR),{r,∞,3}],r,-3],C[2]][[1,1,2]]

(* Finally I introduce the euqation for (𝒬 = λ ℰ) and I solve for λ *) 

Eq1=Collect[solHR//.{r->R,C[1]->c1,C[2]->c2,𝒬->-λ ℰT},λ,FullSimplify];

Eq2=Collect[D[solHR,r]//.{r->R,C[1]->c1,C[2]->c2,𝒬->-λ ℰT},λ,FullSimplify];

Eq1

Eq2

(* k2 is defined in terms of the compactness C=M/R and of the varibale 𝓎 = R H'[R]/H[r], i.e. a combination of H and its derivative at the stellar surface *)
k2=3/2/R^5Solve[{Eq1==HR,Eq2==DHR},{ℰT,λ}][[1,2,2]]//.{M->R 𝒞,DHR->((HR 𝓎)/R)}//FullSimplify