(* This notebook calculates the Noise Kernel for static de Sitter (without delta function) using the expression derived by Phillips and Hu *)

SetAttributes[κ,Constant];
SetAttributes[α,Constant];

(* This section defines a few shorthand assignments for various functions and derivatives to make the expression more readable. 

The variables i,j,k,l are simply placeholders so that I can keep track of indices. 

s is the three-dimensional biscalar sigma, and r is the geodesic distance given by Sqrt[s/2] *)

r:=R[i,j,k,l];
s:=S[i,j,k,l];
(Subscript[s, i]):=Derivative[1,0,0,0][S][i_,j_,k_,l_];
(Subscript[s, j]):=Derivative[0,1,0,0][S][i_,j_,k_,l_]
(Subscript[s, k]):=Derivative[0,0,1,0][S][i_,j_,k_,l_]
(Subscript[s, l]):=Derivative[0,0,0,1][S][i_,j_,k_,l_]

Format[R[i,j,k,l]]:="r";
Format[(Subscript[δ, ij])[i,j,k,l]]:="(Subscript[δ, ij])";
Format[(Subscript[δ, ik])[i,j,k,l]]:="(Subscript[δ, ik])";
Format[(Subscript[δ, il])[i,j,k,l]]:="(Subscript[δ, il])";
Format[(Subscript[δ, jk])[i,j,k,l]]:="(Subscript[δ, jk])";
Format[(Subscript[δ, jl])[i,j,k,l]]:="(Subscript[δ, jl])";
Format[(Subscript[δ, kl])[i,j,k,l]]:="(Subscript[δ, kl])";
Format[Derivative[1,0,0,0][S][i_,j_,k_,l_]]:="(Subscript[s, i])";
Format[Derivative[0,1,0,0][S][i_,j_,k_,l_]]:="(Subscript[s, j])";
Format[Derivative[0,0,1,0][S][i_,j_,k_,l_]]:="(Subscript[s, k])";
Format[Derivative[0,0,0,1][S][i_,j_,k_,l_]]:="(Subscript[s, l])";
Format[(Subscript[s, x])]:="(Subscript[s, x])";
Format[(Subscript[s, y])]:="(Subscript[s, y])";
Format[(Subscript[s, z])]:="(Subscript[s, z])";

(* Replacement rule for s in terms of r. *)

rsrule={s-> r^2/2};

$Assumptions = {r≥ 0,κ>0, A1≥ 0, A2≥ 0,α≥0,ρ1≥ 0,ρ2≥0, B≥ 0,B1≥ 0, B2≥ 0};

(* Defining all of the second and third derivatives of s. *)

Derivative[m1_,m2_,m3_,m4_][R][i_,j_,k_,l_]:=D[Sqrt[2*s],{i,m1},{j,m2},{k,m3},{l,m4}]/.rsrule//Simplify;
Derivative[1,1,0,0][S][i_,j_,k_,l_]:=(Subscript[δ, ij]);
Derivative[1,0,1,0][S][i_,j_,k_,l_]:=(Subscript[δ, ik]);
Derivative[1,0,0,1][S][i_,j_,k_,l_]:=(Subscript[δ, il]);
Derivative[0,1,1,0][S][i_,j_,k_,l_]:=(Subscript[δ, jk]);
Derivative[0,1,0,1][S][i_,j_,k_,l_]:=(Subscript[δ, jl]);
Derivative[0,0,1,1][S][i_,j_,k_,l_]:=(Subscript[δ, kl]);
Derivative[m1_,m2_,m3_,m4_][(Subscript[δ, ij])][i,j,k,l]:=0;
Derivative[m1_,m2_,m3_,m4_][(Subscript[δ, ik])][i,j,k,l]:=0;
Derivative[m1_,m2_,m3_,m4_][(Subscript[δ, il])][i,j,k,l]:=0;
Derivative[m1_,m2_,m3_,m4_][(Subscript[δ, jk])][i,j,k,l]:=0;
Derivative[m1_,m2_,m3_,m4_][(Subscript[δ, jl])][i,j,k,l]:=0;
Derivative[m1_,m2_,m3_,m4_][(Subscript[δ, kl])][i,j,k,l]:=0;

(* Minkowski metric. *)

g={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,-1}};

(* The del function is a hack which lets me take different derivates depending on which component I'm working with.  *)

del[n_,v_,A_]:=Switch[n,1,D[A,v],2,D[A,v],3,D[A,v],4,D[A,η]];

(* Green function for flat space, but using η rather than t for ease of notation. *)

G=Series[1/(8*Pi^2)*(κ/r)*Sinh[κ*r]/(Cosh[κ*r]-Cosh[κ*η]),{κ,0,0}]//Expand//Normal

(* Derivatives of the Green function.

The position of the 1's tell me which index I'm derivating with respect to.  So, 1000 means derivate once with respect to the first index, and 0110 is derivate with respect to the second and third indexes.  *)

G1000[n1_]:=del[n1,i,G];
G0100[n2_]:=del[n2,j,G];
G0010[n3_]:=del[n3,k,G];
G0001[n4_]:=del[n4,l,G];
G1100[n1_,n2_]:=del[n2,j,G1000[n1]];
G1010[n1_,n3_]:=del[n3,k,G1000[n1]];
G1001[n1_,n4_]:=del[n4,l,G1000[n1]];
G0110[n2_,n3_]:=del[n3,k,G0100[n2]];
G0101[n2_,n4_]:=del[n4,l,G0100[n2]];
G0011[n3_,n4_]:=del[n4,l,G0010[n3]];
G1110[n1_,n2_,n3_]:=del[n3,k,G1100[n1,n2]];
G1101[n1_,n2_,n4_]:=del[n4,l,G1100[n1,n2]];
G1011[n1_,n3_,n4_]:=del[n4,l,G1010[n1,n3]];
G0111[n2_,n3_,n4_]:=del[n4,l,G0110[n2,n3]];
G1111[n1_,n2_,n3_,n4_]:=del[n4,l,G1110[n1,n2,n3]];

(*  Phillips and Hu split the Noise Kernel up into four objects based on their dependence on the metric.  I follow a similar procedure here, then further split each object up by their terms and components, which I assign to a table.

Here, I initialize the table to zero. *)

N4=Table[0,{n,10}];
N2=Table[0,{n,11},{n1,4},{n2,4}];
N2p=Table[0,{n,11},{n3,4},{n4,4}];
N0=Table[0,{n,13},{n1,4},{n2,4},{n3,4},{n4,4}];
NKcomp=Table[0,{n1,4},{n2,4},{n3,4},{n4,4}];

(* Replacement rules for dealing with components and derivatives.  srule and delrule replace the derivatives of s with their values in terms of the coordinates.  coordrule and coord2rule collect terms in order to simplify the expression. *)

delrule={(Subscript[δ, ij])->If[n1==n2,1,0],(Subscript[δ, ik])->If[n1==n3,1,0],(Subscript[δ, il])->If[n1==n4,1,0],(Subscript[δ, jk])->If[n2==n3,1,0],(Subscript[δ, jl])->If[n2==n4,1,0],(Subscript[δ, kl])->If[n3==n4,1,0]};
srule={(Subscript[s, i])->Switch[n1,1,(Subscript[s, x]),2,(Subscript[s, y]),3,(Subscript[s, z]),4,0],(Subscript[s, j])->Switch[n2,1,(Subscript[s, x]),2,(Subscript[s, y]),3,(Subscript[s, z]),4,0],(Subscript[s, k])->Switch[n3,1,(Subscript[s, x]),2,(Subscript[s, y]),3,(Subscript[s, z]),4,0],(Subscript[s, l])->Switch[n4,1,(Subscript[s, x]),2,(Subscript[s, y]),3,(Subscript[s, z]),4,0]};
s1rule={(Subscript[s, i])->Switch[n1,1,(Subscript[s, x]),2,(Subscript[s, y]),3,(Subscript[s, z]),4,0],(Subscript[s, j])->Switch[n2,1,(Subscript[s, x]),2,(Subscript[s, y]),3,(Subscript[s, z]),4,0]};
s2rule={(Subscript[s, k])->Switch[n3,1,(Subscript[s, x]),2,(Subscript[s, y]),3,(Subscript[s, z]),4,0],(Subscript[s, l])->Switch[n4,1,(Subscript[s, x]),2,(Subscript[s, y]),3,(Subscript[s, z]),4,0]};
scrule={(Subscript[s, i])->Switch[n1,1,x,2,y,3,z,4,0],(Subscript[s, j])->Switch[n2,1,x,2,y,3,z,4,0],(Subscript[s, k])->Switch[n3,1,x,2,y,3,z,4,0],(Subscript[s, l])->Switch[n4,1,x,2,y,3,z,4,0]};
sc1rule={(Subscript[s, i])->Switch[n1,1,x,2,y,3,z,4,0],(Subscript[s, j])->Switch[n2,1,x,2,y,3,z,4,0]};
sc2rule={(Subscript[s, k])->Switch[n3,1,x,2,y,3,z,4,0],(Subscript[s, l])->Switch[n4,1,x,2,y,3,z,4,0]};
sfrule={(Subscript[s, x])->x,(Subscript[s, y])->y, (Subscript[s, z])->z};
coordrule={(Subscript[s, x])->Sqrt[r^2-(Subscript[s, y])^2-(Subscript[s, z])^2]};
coord2rule={Sqrt[r^2-(Subscript[s, y])^2-(Subscript[s, z])^2]->(Subscript[s, x])};

(* Begin calculation of the components of the noise kernel in flat space.  Note that for flat space most of the terms in the full expression are identically zero and therefore are not computed. *)

N4[[1]]=2*Sum[g[[n1,n2]]*g[[n3,n4]]*G1010[n1,n3]*G0101[n2,n4]/.srule/.delrule,{n1,4},{n2,4},{n3,4},{n4,4}]/.coordrule//Expand

N2[[1]]=Table[-8*Sum[g[[n3,n4]]*G1010[n1,n3]*G0101[n2,n4]/.sc1rule/.s2rule/.delrule/.coordrule,{n3,4},{n4,4}]//Simplify//Expand,{n1,4},{n2,4}]/.coord2rule/.sfrule//Expand

N2p[[1]]=Table[-8*Sum[g[[n1,n2]]*G1010[n1,n3]*G0101[n2,n4]/.s1rule/.sc2rule/.delrule/.coordrule,{n1,4},{n2,4}]//Simplify//Expand,{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N2[[4]]=Table[4*Sum[g[[n3,n4]]*G0010[n3]*G1101[n1,n2,n4]/.sc1rule/.s2rule/.delrule/.coordrule,{n3,4},{n4,4}]//Simplify//Expand,{n1,4},{n2,4}]/.coord2rule/.sfrule//Expand

N2p[[4]]=Table[4*Sum[g[[n1,n2]]*G1000[n1]*G0111[n2,n3,n4]/.s1rule/.sc2rule/.delrule/.coordrule,{n1,4},{n2,4}]//Simplify//Expand,{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N0[[1]]=Table[4*G0110[n2,n3]*G1001[n1,n4]/.srule/.delrule/.coordrule//Simplify//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N0[[2]]=Table[4*G1010[n1,n3]*G0101[n2,n4]/.srule/.delrule/.coordrule//Simplify//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N0[[3]]=Table[G1100[n1,n2]*G0011[n3,n4]/.srule/.delrule/.coordrule//Simplify//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N0[[4]]=Table[G*G1111[n1,n2,n3,n4]/.srule/.delrule/.coordrule//Simplify//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N0[[5]]=Table[-2*G0100[n2]*G1011[n1,n3,n4]/.srule/.delrule/.coordrule//Simplify//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N0[[6]]=Table[-2*G1000[n1]*G0111[n2,n3,n4]/.srule/.delrule/.coordrule//Simplify//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N0[[7]]=Table[-2*G0001[n4]*G1110[n1,n2,n3]/.srule/.delrule/.coordrule//Simplify//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

N0[[8]]=Table[-2*G0010[n3]*G1101[n1,n2,n4]/.srule/.delrule/.coordrule//Simplify//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]/.coord2rule/.sfrule//Expand

(* Defining the conformal transformation. *)

Ω[η_]:=-α/η;

Ω[η1]^2*Ω[η2]^2

(* Collecting terms and conformalling transforming in order to generate the full expression for the noise kernel.  *)

NKcomp=Table[(1/9*Sum[N0[[m1,n1,n2,n3,n4]],{m1,13}]+1/36*g[[n1,n2]]*Sum[N2p[[m2,n3,n4]],{m2,11}]+1/36*g[[n3,n4]]*Sum[N2[[m3,n1,n2]],{m3,11}]+1/36*g[[n1,n2]]*g[[n3,n4]]*Sum[N4[[m4]],{m4,10}])*Ω[η1]^-2*Ω[η2]^-2//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]

NKcomp[[4,4,4,4]]//Simplify

(* Code to test that the partial traces vanish. *)

NKTrace=Table[Sum[g[[n1,n1]]*NKcomp[[n1,n1,n3,n4]],{n1,1,4}]//Simplify,{n3,1,4},{n4,1,4}]

(* Defining rules and functions for the coord. transform to static de Sitter. *)

deltarule={η->η1-η2,x->x1-x2, y->y1-y2,z->z1-z2};

Arule={α^2-ρ1^2-> A1^2,α^2-ρ2^2->A2^2}
A2rule={A1-> Sqrt[α^2-ρ1^2],A2->Sqrt[α^2-ρ2^2]}
A3rule={ρ1-> Sqrt[α^2-A1^2],ρ2->Sqrt[α^2-A2^2]}

z1c[x_,y_,z_,η_]:=-α x/η;
z2c[x_,y_,z_,η_]:=-α y/η;
z3c[x_,y_,z_,η_]:=-α z/η;
z4c[x_,y_,z_,η_]:= -α^2/η;

z1s[ρ_,θ_,ϕ_,t_]:=ρ Sin[θ]Cos[ϕ];
z2s[ρ_,θ_,ϕ_,t_]:=ρ Sin[θ]Sin[ϕ];
z3s[ρ_,θ_,ϕ_,t_]:=ρ Cos[θ];
z4s[ρ_,θ_,ϕ_,t_]:=Sqrt[α^2-ρ^2]Exp[t/α];

T1[A_,B_]:=Solve[{z1c[A[[1]],A[[2]],A[[3]],A[[4]]]==z1s[B[[1]],B[[2]],B[[3]],B[[4]]],z2c[A[[1]],A[[2]],A[[3]],A[[4]]]==z2s[B[[1]],B[[2]],B[[3]],B[[4]]],z3c[A[[1]],A[[2]],A[[3]],A[[4]]]==z3s[B[[1]],B[[2]],B[[3]],B[[4]]],z4c[A[[1]],A[[2]],A[[3]],A[[4]]]==z4s[B[[1]],B[[2]],B[[3]],B[[4]]]},A][[1]]

c1={x1,y1,z1,η1};
c2={x2,y2,z2,η2};
s1={ρ1,θ1,ϕ1,t1};
s2={ρ2,θ2,ϕ2,t2};

Vs[A_,B_]:=A/.T1[A,B]

V1=Vs[c1,s1]/.Arule//Simplify
V2=Vs[c2,s2]/.Arule//Simplify

Vtemp=Vs[c1,s1]//Simplify

Bsrule={α^2-ρ^2->B*α^2}
csrule={t1->T,ρ1->ρ,θ1->θ,ϕ1->ϕ}
tsrule={E^(-t1/α)->τ}

Vtemp/.tsrule/.csrule/.Bsrule//Simplify

staticrule={c1[[1]]-> V1[[1]],c1[[2]]-> V1[[2]],c1[[3]]-> V1[[3]],c1[[4]]-> V1[[4]],c2[[1]]-> V2[[1]],c2[[2]]-> V2[[2]],c2[[3]]-> V2[[3]],c2[[4]]-> V2[[4]]};

rrule={r->Sqrt[x^2+y^2+z^2/.deltarule/.staticrule//Expand//Simplify]};

(* Jacobian for the transform from comoving cartesian to static spherical coords. *)

Js[A_,B_]:=Table[D[Vs[A,B][[i]],B[[j]]],{i,4},{j,4}]

J1=Js[c1,s1]//Simplify
J2=Js[c2,s2]//Simplify

J1//Simplify//MatrixForm

(* Testing the transform against the metric. *)

gc1={{α^2/η1^2,0,0,0},{0,α^2/η1^2,0,0},{0,0,α^2/η1^2,0},{0,0,0,-α^2/η1^2}};

gc2={{α^2/η2^2,0,0,0},{0,α^2/η2^2,0,0},{0,0,α^2/η2^2,0},{0,0,0,-α^2/η2^2}};

gs1=Table[Sum[gc1[[a1,a2]]*J1[[a1,b1]]*J1[[a2,b2]]/.staticrule,{a1,4},{a2,4}],{b1,4},{b2,4}]/.Arule/.A2rule//Simplify;
gs2=Table[Sum[gc2[[a1,a2]]*J2[[a1,b1]]*J2[[a2,b2]]/.staticrule,{a1,4},{a2,4}],{b1,4},{b2,4}]/.Arule/.A2rule//Simplify;

gs2//Simplify//MatrixForm

(* Changing variables. *)

NKcomp2=Table[NKcomp[[n1,n2,n3,n4]]/.rrule/.deltarule/.staticrule/.Arule//Expand,{n1,4},{n2,4},{n3,4},{n4,4}];

(* Transforming a given component to static coords.  In this case, the ttt't' component. *)

cmp={4,4,4,4};NKstatic=Sum[NKcomp2[[n1,n2,n3,n4]]*J1[[n1,cmp[[1]]]]*J1[[n2,cmp[[2]]]]*J2[[n3,cmp[[3]]]]*J2[[n4,cmp[[4]]]],{n1,4},{n2,4},{n3,4},{n4,4}]//Expand

(* Raising the first and third indicies and starting to collect terms and simplify. *)

NKst1=NKstatic/((-1+ρ1^2/α^2)*(-1+ρ2^2/α^2))/.Arule//Together//Expand

(* More rules for simplifying the exression.  *)

ind1=1/gs1[[cmp[[1]],cmp[[1]]]];
ind2=1/gs1[[cmp[[2]],cmp[[2]]]];
ind3=1/gs2[[cmp[[3]],cmp[[3]]]];
ind4=1/gs2[[cmp[[4]],cmp[[4]]]];

orule={θ2->0,ϕ1->0,ϕ2->0};

dtrule={t2->t1-t};
drrule={ρ2->ρ1-ρ,A2->Sqrt[A1^2-(ρ^2)+2 ρ ρ1]};

trule={t1->α Log[τ1],t2->α Log[τ2],t->α Log[τ]};
t2rule={τ->Exp[t/α],τ1->Exp[t1/α],τ2->Exp[t2/α]};

r2rule={Cos[θ1]->(-r^2+(((α^2) ((ρ1)^2))/((A1^2) ((τ1)^2)))+(((α^2) ((ρ2)^2))/((A2^2) ((τ2)^2))))/(((2 (α^2) ρ1 ρ2 )/(A1 A2 τ1 τ2)))}

Brule = {ρ1 -> α*Sqrt[1-B1], ρ2 -> α*Sqrt[1-B2]}
B2rule={1-B1-> ρ1^2/α^2,1-B2 -> ρ2^2/α^2,-1+B1 -> -ρ1^2/α^2,-1+B2 -> -ρ2^2/α^2,B1-B1^2 ->B1* ρ1^2/α^2,B2-B2^2 -> B2*ρ2^2/α^2}
B3rule={B1->1-ρ1^2/α^2,B2->1- ρ2^2/α^2}

$Assumptions={A1∈Reals,A1>0,A1^2>0,A2∈Reals,A2>0,A2^2>0,B1∈Reals,B1>0,B1^2>0,B2∈Reals,B2>0,B2^2>0,B1*B2>0,ρ1∈Reals,ρ1≥0,ρ1^2≥0,ρ2∈Reals,ρ2≥0,ρ2^2≥0,α∈Reals,α>0,α^2>0}

(* Simplifying the expression.  I set θ2, ϕ1, and ϕ2 to zero so that θ1 becomes the angle between the two points. Also, I replace the time coordinate with its log. *)

NKst2=NKst1/.orule/.dtrule/.trule//Simplify

(* Writing in terms of radial separation rather than radial coordinates. *)

NKst3=NKst2/.drrule//Simplify

(* Time separation. *)

NKtime1 = NKst3/.{ρ->0,θ1->0}//Simplify

NKtime1a=NKtime1/.Brule//Simplify

NKtime1b=NKtime1a/.{((τ^4)/(((-1+τ))^8))->1/(t-2)^4}//Simplify

NKtime1b/.{t->τ,B1->B}//TeXForm

NKtime2=NKtime1/.t2rule/.Brule//Simplify

NKtime3=Series[NKtime2,{t,0,-4}]//Normal

TeXForm[NKtime3]

(* Plot of the relative error of the truncated series versus the exact expression. *)

t/a

B0=0.5;
a0=1;
h=0.01;
d=4*a0;
LogPlot[{(NKtime3-NKtime2)/NKtime2/.{B1->B0,α->a0}},{t,0,d},RegionFunction-> Function[t,t>0.015],PlotStyle->Black,Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},FrameLabel->{Style["|ΔT| / α",Bold,12]}]

(* Radial separation. *)

NKrad1=NKst3/.{τ->1,θ1->0}/.A2rule//Simplify//Expand//Simplify

NKrad2=NKrad1/.Brule//Simplify

$Assumptions={ρ1∈Reals,ρ1≥0,ρ1^2≥0,ρ2∈Reals,ρ2≥0,ρ2^2≥0,α∈Reals,α≥0,α^2≥0,B1∈Reals,B1≥0,B2∈Reals,B2≥0};

NKrad2a = NKrad2/.{ρ->ρ1-ρ2}/.B2rule//Simplify

NKrad2b=NKrad2a/.{B1->1-ρ1^2/α^2}//Simplify

NKrad2c=NKrad2b/.Brule//Simplify

NKrad2d=NKrad2c/.B2rule//Simplify//Expand

NKrad2e=Sum[NKrad2d[[n]]/.{B1->B12/B2}/.{ρ1->α^2*(1-p12)/ρ2}//Simplify,{n,1,7}]

NKrad2f=NKrad2e/.{B12->B1*B2}//Simplify//Expand//Together//Simplify

NKrad2f/.{B1->B,B2->B'}//TeXForm

NKrad3=Series[NKrad2,{ρ,0,-4}]//Normal

NKrad4=NKrad3//Simplify//Expand

NKrad5=NKrad4/.{(Sqrt[1-B1])->(1-B1)*α/ρ1}//Simplify//Expand

NKrad6=Sum[Coefficient[NKrad5*(π^4) ,ρ,n]*ρ^n/Pi^4//Simplify,{n,-8,-4}]

(* Plot of the relative error of the truncated series versus the exact expression. *)

NKerr=(NKrad6-NKrad2)/NKrad2/.{B1->1-ρ1^2/α^2,ρ->ρ1-ρ2};

$Assumptions={};

ρ/α

p0=0.25;
a0=1;
h=0.01;
d=a0;
LogPlot[NKerr/.{α->a0,ρ2->p0},{ρ1,0,d},RegionFunction->Function[ρ1,Not[p0-0.1<ρ1<p0+0.095]],PlotStyle->Black,Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},FrameLabel->{Style["ρ / α",Bold,12]}]

p0=0.75;
a0=1;
h=0.01;
c=4;
q=1;
d=(c*(1-p0))^q*a0;
LogPlot[NKerr/.{α->a0,ρ2->p0},{ρ1,0,d},RegionFunction->Function[ρ1,Not[p0-0.06<ρ1<p0+0.05]],PlotStyle->Black,Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},FrameLabel->{Style["ρ / α",Bold,12]}]

p0=0.95;
a0=1;
h=0.01;
c=20;
q=1;
d=(c*(1-p0))^q*a0;
LogPlot[NKerr/.{α->a0,ρ2->p0},{ρ1,a0-d,a0},RegionFunction->Function[ρ1,Not[p0-0.017<ρ1<p0+0.013]],PlotStyle->Black,Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},FrameLabel->{Style["ρ / α",Bold,12]}]

p0=0.99;
a0=1;
h=0.01;
c=20;
q=1;
d=(c*(1-p0))^q*a0;
LogPlot[NKerr/.{α->a0,ρ2->p0},{ρ1,a0-d,a0},RegionFunction->Function[ρ1,Not[p0-0.002<ρ1<p0+0.002]],PlotStyle->Black,Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},FrameLabel->{Style["ρ / α",Bold,12]}]

Clear[NKerr1]

p0={0.25,0.75,0.95,0.99};
a0=1;
NKerr1:=If[p0[[1]]-0.1<ρ1<p0[[1]]+0.095,0,NKerr/.{α->a0,ρ2->p0[[1]]}];
NKerr2:=If[p0[[2]]-0.06<ρ1<p0[[2]]+0.05,0,NKerr/.{α->a0,ρ2->p0[[2]]}];
NKerr3:=If[p0[[3]]-0.017<ρ1<p0[[3]]+0.013,0,NKerr/.{α->a0,ρ2->p0[[3]]}];
NKerr4:=If[p0[[4]]-0.002<ρ1<p0[[4]]+0.002,0,NKerr/.{α->a0,ρ2->p0[[4]]}];
LogPlot[{NKerr1,NKerr2,NKerr3,NKerr4},{ρ1,0,1},PlotStyle->{{Black,Dotted},{Black,Dashed},{Black,Dashing[0.025]},Black},Frame->True,PlotRange->{0.00000001,1},FrameTicks->{{Automatic,None},{Automatic,None}},FrameLabel->{Style["ρ / α",Bold,12]}]

p0=0.1;
a0=1;
h=1;
d=a0/10;
(NKrad6-NKrad2)/NKrad2/.{ρ->p0,α->a0,B1->1-ρ1^2/a0^2}/.{ρ1->0.999}
LogPlot[{NKrad6/NKrad2-1/.{ρ->p0,α->a0,B1->1-ρ1^2/a0^2}},{ρ1,1-d,1},PlotStyle->Black,Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},FrameLabel->{Style["ρ / α",Bold,12]}]

TeXForm[NKrad6]

NKrerr1=NKrad6/NKrad2-1//Simplify

NKrerr2=NKrerr1/.{ρ1->α*Sqrt[1-B1]}//Simplify

NKrerr3=Series[NKrerr2,{B1,0,2}]//Normal

NKrerr4=Series[NKrerr3,{ρ,0,0}]//Normal

NKrerr4/.{B1->0.1,ρ->1,α->1}

(*Angular separation.*)

NKang1=NKst3/.{τ->1,ρ->0,Cos[2θ1]->Cos[θ1]^2-1}/.{Cos[θ1]->η+1}/.A2rule//Simplify//Expand//Simplify

NKang2 = Series[NKang1,{η,0,2}]//Normal

TeXForm[NKang2/.{ρ1->"ρ"}]

(* Simplifying the exact expression into a displayable form. *)

$Assumptions={A1∈Reals,A1>0,A1^2>0,A2∈Reals,A2>0,A2^2>0,B1∈Reals,B1>0,B1^2>0,B2∈Reals,B2>0,B2^2>0,B1*B2>0,ρ1∈Reals,ρ1≥0,ρ1^2≥0,ρ2∈Reals,ρ2≥0,ρ2^2≥0,α∈Reals,α>0,α^2>0,τ∈Reals,τ>0};

NKst4=NKst1/.orule/.trule/.A2rule/.Brule//Simplify

NKst5=NKst4/.B2rule//Simplify//Expand

NKst6=NKst5/.{B1*B2->B12}/.{B1^n_*B2^n_->B12^n}//Simplify

NKst7=NKst6/.{B1*B2->B12}/.{B1^n_*B2^n_->B12^n}//Simplify

NKst8=NKst7/.t2rule/.dtrule/.trule//Simplify//Expand//Together//Simplify

NKst9=NKst8/.{B12->B1*B2}

NKst9//Simplify

NKst9/.{B1->1,B2->1,ρ1->0,ρ2->0,θ1->0}/.taurule//Simplify

taurule=Solve[τ+1/τ==T,τ][[1]]

NKst10=NKst9/τ^6*((((α^2) (-2 τ+(Sqrt[B1 B2]) (1+(τ^2)))+2 ρ1 ρ2 τ Cos[θ1]))^6)//Simplify

Table[Coefficient[NKst10,τ,n],{n,-10,10}]

NKst11a=(Coefficient[NKst10,τ,1]*τ+Coefficient[NKst10,τ,-1]/τ)/.taurule//Simplify

NKst11b=(Coefficient[NKst10,τ,2]*τ^2+Coefficient[NKst10,τ,0]+Coefficient[NKst10,τ,-2]/τ^2)/.taurule//Simplify//Expand//Together

NKst11=(NKst11a+NKst11b)*12*Pi^4//Simplify//Expand

Sum[Coefficient[NKst11,α,n]*α^n//Simplify,{n,0,10}]

(* Exact expression on the horizon. *)

NKst9h1=NKst9/.{B1->0,ρ1->α}//Simplify

NKst9h2=NKst9h1/.{B2->"B'",ρ2->"ρ'",θ1->γ}

NKst9h2//TeXForm

NKst9h=NKst9/.{B1->0,B2->0,ρ1->α,ρ2->α}//Simplify

NKst9ha=NKst9/.{Cos[2θ1]->Cos[θ1]^2-1}/.{Cos[θ1]->η+1}/.{B1->0,B2->0,ρ1->α,ρ2->α}//Simplify

(* Finding a second component. *)

cmp={4,1,4,1};NKstatic2=Sum[NKcomp2[[n1,n2,n3,n4]]*J1[[n1,cmp[[1]]]]*J1[[n2,cmp[[2]]]]*J2[[n3,cmp[[3]]]]*J2[[n4,cmp[[4]]]],{n1,4},{n2,4},{n3,4},{n4,4}]//Expand

NKtr1=NKstatic2/((-1+ρ1^2/α^2)*(-1+ρ2^2/α^2))/.Arule//Together//Expand

NKtr2=NKtr1/.orule/.dtrule/.trule//Simplify

NKtr4=NKtr2/.A2rule/.Brule//Simplify

NKtr5=NKtr4/.B2rule//Simplify//Expand

NKtr6=NKtr5/.{B1*B2->B12}/.{B1^n_*B2^n_->B12^n}//Simplify

NKtr7=NKtr6/.{B1*B2->B12}/.{B1^n_*B2^n_->B12^n}//Simplify

NKtr9=NKtr7/.{B12->B1*B2}//Simplify

(* Expression on the horizon. *)

NKtr9h1=NKtr9*B1*B2/.{B1->0,ρ1->α}//Simplify

NKtr9h2=NKtr9h1/.{B2->"B'",ρ2->"ρ'",θ1->γ}

NKtr9h2//TeXForm

NKtr9h=NKtr9*B1*B2/.{B1->0,B2->0,ρ1->α,ρ2->α}//Simplify

TeXForm[NKst9]

TeXForm[NKtr9]

TeXForm[NKst9h]

TeXForm[NKtr9h]