Tidal deformability in massive scalar tensor gravity

The tidal deformability in ST theory is useful for: i) inspiral phase GW phase study and compare with GR; ii) study the f_peak universal relation and tidal deformabilty (possibly broken)

ClearAll[X,ω,r,g,x,ds,K,m]
$Assumptions={r≥0,θ≥0,θ≤π};
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[θ]}};
h=DiagonalMatrix[{-H0[r]*Exp[ν[r]],H2[r]*(((1-2((M[r])/r)))^(-1)),(r^2)*K[r],(r^2)*X[θ]*K[r]}]*P2[θ]*ϵ;
gdd0=g0+h;
rule1[Z_]:=Z/.Union[{((X'[θ])^2)->4 X[θ](1-X[θ]),X''[θ]->2 (1-2X[θ])}];

gdd0

1. Spacetime

gdd=Simplify[rule1[Normal[Series[gdd0,{ϵ,0,order}]]]];
guu=Simplify[rule1[Normal[Series[Inverse[gdd],{ϵ,0,order}]]]];

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[ds[i]=x[[i]];,{i,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[ig[i,j]*Gdd[j,k],{j,1,dim}],{ϵ,0,order}]]]];,{i,1,4},{k,1,4}]

Normal[Series[Gud[1,1],{ϵ,0,1}]]

2. matter (calculated in the Jordan frame, metric and four velocity are in the Einstein frame, p and rho are in the Jordan frame)

Clear[ϕ];
ϕ[r_]:=φ0[r]+ϵ*ϕ2[r]P2[θ];
u={(1/(Sqrt[-g[1,1]])),0,0,0};
Δp=ϵ*δp[r]*P2[θ];
Δe=dedp[r]*Δp;

Do[Tuu[i,j]=Simplify[(1/((A[ϕ[r]])^2))*rule1[Normal[Series[(p[r]+e[r]+Δp+Δe)u[[i]]u[[j]]+ig[i,j] (p[r]+Δp),{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[Tdd[i,j]=Simplify[rule1[((A[ϕ[r]])^4)*Normal[Series[(Sum[Tuu[a,b]g[a,i]g[b,j],{a,1,dim},{b,1,dim}]),{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
(*Do[Tud[i,k]=Simplify[rule1[Normal[Series[Sum[ig[i,j]*Tdd[j,k],{j,1,dim}],{ϵ,0,order}]]]];,{i,1,dim},{k,1,dim}];*)
T=Simplify[rule1[Normal[Series[Sum[Tdd[i,j]((A[ϕ[r]])^(-2)) ig[i,j],{i,1,dim},{j,1,dim}],{ϵ,0,order}]]]];

Do[Tmunu[i,j]=Simplify[rule1[Normal[Series[2*D[ϕ[r],ds[i]]*D[ϕ[r],ds[j]]-(1/2)g[i,j]V[ϕ[r]]-g[i,j]*Sum[ig[a,b]*D[ϕ[r],ds[a]]*D[ϕ[r],ds[b]],{a,1,dim},{b,1,dim}]+8π*Tdd[i,j]*((A[ϕ[r]])^2),{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}]
Do[Einstein[i,j]=Simplify[rule1[Normal[Series[Gdd[i,j]- Tmunu[i,j],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[DivT[i]=Normal[Series[Sum[ig[a,j]ig[b,i](D[Einstein[a,b],ds[j]]-Sum[Chr[l,a,j]Einstein[l,b]+Chr[l,b,j]Einstein[l,a],{l,1,dim}]),{a,1,dim},{b,1,dim},{j,1,dim}],{ϵ,0,order}]];,{i,1,dim}]

Einstein equation and energy conservation

(*Do[Einstein[i,j]=Simplify[rule1[Normal[Series[Ricci[i,j]-2*D[ϕ[r],ds[i]]*D[ϕ[r],ds[j]]-(1/2)g[i,j]*V[ϕ[r]]-8π*((A[ϕ[r]])^2) Tdd[i,j]+4π g[i,j]*((A[ϕ[r]])^4)T,{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[DivT[i]=Simplify[rule1[Normal[Series[Sum[D[((A[ϕ[r]])^6)Tuu[i,j],ds[j]]+((A[ϕ[r]])^6)*Sum[Chr[i,a,j]*Tuu[a,j]+Chr[j,a,j]Tuu[a,i],{a,1,dim}]-T ig[i,j]*dAdϕ[r]*D[ϕ[r],ds[j]]/A[ϕ[r]],{j,1,dim}],{ϵ,0,order}]]]];,{i,1,dim}]*)

(*Do[Einstein[i,j]=Simplify[rule1[Normal[Series[Gdd[i,j]- Tmunu[i,j],{ϵ,0,order}]]]];,{i,1,dim},{j,1,dim}];
Do[DivT[i]=Normal[Series[Sum[ig[a,j]ig[b,i](D[Einstein[a,b],ds[j]]-Sum[Chr[l,a,j]Einstein[l,b]+Chr[l,b,j]Einstein[l,a],{l,1,dim}]),{a,1,dim},{b,1,dim},{j,1,dim}],{ϵ,0,order}]];,{i,1,dim}]*)

Scalar function

DET=Simplify[rule1[Normal[Series[Det[gdd],{ϵ,0,order}]]]];
LHS=Simplify[rule1[Normal[Series[(1/(Sqrt[-DET]))Sum[D[(Sqrt[-DET])D[ϕ[r],ds[i]]*ig[i,j],ds[j]],{i,1,dim},{j,1,dim}],{ϵ,0,order}]]]];

RHS=Simplify[rule1[Normal[Series[(1/4)V'[ϕ[r]]-4*π*A'[ϕ[r]]*((A[ϕ[r]])^3)*T,{ϵ,0,order}]]]];
scalar=Simplify[(LHS-RHS)];

scalar

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];
eqs=scalar;

SetDirectory[NotebookDirectory[]];
DumpSave["matter.mx",{eq11,eq12,eq13,eq14,eq22,eq23,eq24,eq33,eq34,eq44,div1,div2,div3,div4,eqs}];

Quit[]

Clear[r,X];
SetDirectory[NotebookDirectory[]];
<<"matter.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,eqs};

0th order equations

SS[eq11,0]
SS[eq22,0]
SS[div2,0]
SS[eqs,0]

syst0=Solve[{SS[eq11,0]==0,SS[eq22,0]==0,SS[div2,0]==0,SS[eqs,0]==0},{M'[r],ν'[r],p'[r],φ0''[r]}][[1]]//FullSimplify
syst0=Union[syst0,D[syst0,r],D[syst0,{r,2}]];

1st order equations

systl2=Solve[{(SS[div3,1]//.syst0)==0,(SS[(eq44/((Sin[θ])^2))-eq33,1]//.syst0)==0,(SS[eq23,1]//.syst0)==0},{δp[r],H2[r],K'[r]}][[1]]//FullSimplify
systl2=Union[systl2,D[systl2,r],D[systl2,{r,2}]];

kdr=K'[r]/.systl2

eqsol=Solve[(SS[eqs,1]//.Union[syst0,systl2])==0,(ϕ2)''[r]][[1]]//FullSimplify
systS=Union[eqsol,D[eqsol,r],D[eqsol,{r,2}]];

masterϕ2=(ϕ2)''[r]/. eqsol

d1=-Simplify[Coefficient[masterϕ2,(ϕ2)'[r]]]
d0=-Simplify[Coefficient[masterϕ2,ϕ2[r]]]
ds=Simplify[Coefficient[masterϕ2+d1*(ϕ2)'[r]+d0*ϕ2[r],H0[r]]]

ksol=Solve[(SS[eq22,1]//.Union[syst0,systl2,systS])==0,K[r]][[1]]//FullSimplify

kd=Simplify[D[K[r]/.ksol,r]//.Union[{e'[r]->p'[r]*dedp[r]},syst0,systl2,systS]]

h0sol=Solve[kdr==kd,H0''[r]][[1]]//FullSimplify

masterh = H0''[r]/.h0sol

c1=-FullSimplify[Coefficient[masterh,H0'[r]]]
c0=-FullSimplify[Coefficient[masterh,H0[r]]]
cs=FullSimplify[Coefficient[masterh+c1*H0'[r]+c0*H0[r],ϕ2[r]]]

coe0=Simplify[Coefficient[c0,dedp[r]]]

TeXForm[Simplify[coe0/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)'}]]

Expand[c0-Coefficient[c0,dedp[r]]*dedp[r]]

coe1 =Simplify[Coefficient[Expand[c0-Coefficient[c0,dedp[r]]*dedp[r]],(((φ0)'[r])^2)]]

TeXForm[Simplify[coe1/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)'}]]

coe2=Simplify[Coefficient[Expand[c0-Coefficient[c0,dedp[r]]*dedp[r]],(((φ0)'[r])^4)]]

coe3=FullSimplify[c0 - coe0*dedp[r]-coe1*(((φ0)'[r])^2)-coe2*(((φ0)'[r])^4)]

coe4=Coefficient[coe3,V[φ0[r]]]

TeXForm[Simplify[coe4/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)'}]]

coe5=Coefficient[coe3-coe4*V[φ0[r]],((V[φ0[r]])^2)]

TeXForm[Simplify[coe5/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)'}]]

coe6 = FullSimplify[coe3-coe4*V[φ0[r]]-coe5*((V[φ0[r]])^2)]

TeXForm[FullSimplify[coe6/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)'}]]

Simplify[c0- coe0*dedp[r]-coe1*(((φ0)'[r])^2)-coe2*(((φ0)'[r])^4)-coe4*V[φ0[r]]-coe5*((V[φ0[r]])^2)-coe6]

The final form of the equation

c1

c1/d1

cs

TeXForm[Simplify[cs/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)'}]]

Simplify[ds/cs]

Simplify[c0]

TeXForm[Simplify[c0/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)'}]]

d0

doe1=Simplify[Coefficient[d0,dedp[r]]]

TeXForm[Simplify[doe1/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)'}]]

doe2=Simplify[d0-Coefficient[d0,dedp[r]]dedp[r]]

TeXForm[Simplify[doe2/.{M[r]->m,A[φ0[r]]->A,V[φ0[r]]->V,e[r]->e,p[r]->p,dedp[r]->dedp,V'[φ0[r]]->V',A'[φ0[r]]->A*α,(φ0)'[r]->φ',(φ0)'[r]->(φ0)',V''[φ0[r]]->V'',A''[φ0[r]]->A''}]]

cs

Simplify[Coefficient[cs,A'[φ0[r]]]]

cs1=Simplify[cs-Coefficient[cs,A'[φ0[r]]]A'[φ0[r]]]

Simplify[Coefficient[cs1,((A[φ0[r]])^4)]]

cs2=Simplify[cs1-Coefficient[cs1,((A[φ0[r]])^4)]((A[φ0[r]])^4)]

Simplify[c1/d1]
Simplify[cs/ds]

d0

Simplify[Coefficient[d0,dedp[r]]]

Simplify[d0-Coefficient[d0,dedp[r]]*dedp[r]]

Simplify[(24-96 π (r^2) ((A[φ0[r]])^2) (P[r]-ρ[r]) ((A'[φ0[r]])^2)+16 r (r-2 M[r]) (((φ0)'[r])^2)+16 π (r^2) ((A[φ0[r]])^3) (-3 P[r]+ρ[r]) A''[φ0[r]]+(r^2) V''[φ0[r]])/(-4)]

Simplify[c1/.{p[r]->0,e[r]->0}]
Simplify[c0/.{p[r]->0,e[r]->0}]
Simplify[cs/.{p[r]->0,e[r]->0}]

FullSimplify[(24 (r^2)+(r^4) (4 V[φ0[r]]+(r^2) (((V[φ0[r]]-2 (((φ0)'[r])^2)))^2))+8 r M[r] (-6+(r^2) (-2 V[φ0[r]]+(2+(r^2) V[φ0[r]]) (((φ0)'[r])^2)-2 (r^2) (((φ0)'[r])^4))))]

Simplify[d0/.{p[r]->0,e[r]->0}]

doe1=Simplify[Coefficient[d0,dedp[r]]]

Simplify[d0-doe1*dedp[r]]

Final form of equations

RULE=Union[syst0,systl2,systS,h0sol];

E01 = (M'[r])-(M'[r]//.RULE)
E02 = (ν'[r])-(ν'[r]//.RULE)
E03=(P'[r])-(P'[r]//.RULE)
E04=(φ0''[r])-(φ0''[r]//.RULE)
E05 = H0''[r]-(H0''[r]//.RULE)
E06=(ϕ2)''[r]-((ϕ2)''[r]//.RULE)

DumpSave["final_form.mx",{E01,E02,E03,E04,E04,E05,RULE}];

Quit[];

SetDirectory[NotebookDirectory[]];
<<"final_form.mx";

EQs01={E01,E02,E03,E04}

rhoP[x_]:=P2ρ[x]

ORD=6;
RULEc={};
M[r_]:=Sum[Mc[i](r^i),{i,3,ORD}];
ν[r_]:=Sum[νc[i](r^i),{i,0,ORD}];
P[r_]:=Sum[Pc[i](r^i),{i,0,ORD}];
ρ[r_]:=rhoP[P[r]];
φ0[r_]:=Sum[φ0c[i](r^i),{i,0,ORD}];
ss0=Series[EQs01//.RULEc,{r,0,ORD}];
EQc0=Union[Table[SeriesCoefficient[ss0[[1]],i]==0,{i,2,ORD-1}],Table[SeriesCoefficient[ss0[[2]],i]==0,{i,-1,ORD-1}],Table[SeriesCoefficient[ss0[[3]],i]==0,{i,-1,ORD-1}],Table[SeriesCoefficient[ss0[[4]],i]==0,{i,-2,ORD-2}]];
yc0=Union[Table[Mc[i],{i,3,ORD}],Table[νc[i],{i,1,ORD}],Table[Pc[i],{i,1,ORD}],Table[φ0c[i],{i,1,ORD}]];
systc0=Solve[EQc0,yc0][[1]]//Simplify

Mc[3]//.systc0
νc[2]//.systc0
Pc[2]//.systc0
φ0c[2]//.systc0

Table[SeriesCoefficient[ss0[[1]],i]==0,{i,2,ORD}]

Table[SeriesCoefficient[ss0[[2]],i]==0,{i,-1,ORD}]

Simplify[d0-Coefficient[d0,dρdp[r]]dρdp[r]]

Simplify[Coefficient[d0,dρdp[r]]]

Simplify[c1/.{P[r]->0,ρ[r]->0}]

Simplify[c0/.{P[r]->0,ρ[r]->0}]

Simplify[cs/.{P[r]->0,ρ[r]->0}]

Simplify[d1/.{P[r]->0,ρ[r]->0}]

Simplify[d0/.{P[r]->0,ρ[r]->0}]

sol0=Solve[{SS[eq11,0]==0,SS[eq22,0]==0,SS[div2,0]==0,SS[eqs,0]==0},{M'[r],ν'[r],P'[r],φ0''[r]}][[1]]//FullSimplify

ϕdd = (φ0)''[r]/.sol0

Coefficient[ϕdd,(-P[r]+ρ[r])]
Coefficient[ϕdd,(-3P[r]+ρ[r])]
Coefficient[ϕdd,V[φ0[r]]]
Coefficient[ϕdd,V'[φ0[r]]]

dpdr=P'[r]/.sol0

FortranForm[-(1/(4 r (r-2 M[r])))(24+16 r (-π r ((A[φ0[r]])^2) ((6+dρdp[r]) P[r]+(-6+dρdp[r]) ρ[r]) ((A'[φ0[r]])^2)+(r-2 M[r]) (((φ0)'[r])^2)+π r ((A[φ0[r]])^3) (-3 P[r]+ρ[r]) A''[φ0[r]])+(r^2) V''[φ0[r]])]/.{Pi->pi,φ0(r)->φ0}

TeXForm