(* This notebook calculates the Noise Kernel for hot flat space using the expression derived by Phillips and Hu *)

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

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

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

$Assumptions = {r≥ 0,κ>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,t]];

(* Green function for hot flat space. *)

Gu=1/(8*Pi^2)*(κ/r)*Sinh[κ*r]/(Cosh[κ*r]-Cosh[κ*t])

(* Zero temperature Green function; I use this to cancel divergences when I want a regularized expression. *)

Gr=Series[Gu,{κ,0,0}]//Normal//Simplify

(* Full expression for Green function.  *)

G=Gu(*-Gr*);

(* 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,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]};
coordrule={x->Sqrt[r^2-y^2-z^2]};
coord2rule={Sqrt[r^2-y^2-z^2]->x};

(* Begin calculation of the components of the noise kernel.  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]/.srule/.delrule/.coordrule,{n3,4},{n4,4}]//Expand,{n1,4},{n2,4}]/.coord2rule//Expand

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

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

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

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

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

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

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

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

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

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

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

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

NKcomp=Table[1/72*Sum[N0[[m1,n1,n2,n3,n4]],{m1,13}]+1/288*g[[n1,n2]]*Sum[N2p[[m2,n3,n4]],{m2,11}]+1/288*g[[n3,n4]]*Sum[N2[[m3,n1,n2]],{m3,11}]+1/288*g[[n1,n2]]*g[[n3,n4]]*Sum[N4[[m4]],{m4,10}]//Expand,{n1,4},{n2,4},{n3,4},{n4,4}]

(* 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}]

sis={x,y,z,junk};

ugrule={r->Sqrt[x^2+y^2+z^2]}

(* Testing for conservation.  *)

NKcons=Table[-D[NKcomp[[4,n2,n3,n4]]/.ugrule,t]+D[NKcomp[[1,n2,n3,n4]]/.ugrule,x]+D[NKcomp[[2,n2,n3,n4]]/.ugrule,y]+D[NKcomp[[3,n2,n3,n4]]/.ugrule,z]//Simplify,{n2,1,4},{n3,1,4},{n4,1,4}]

(* Simplifying the ttt't' component to generate a displayable expression.  *)

Nsum1=NKcomp[[4,4,4,4]]//Expand

(* Testing zero temp. limit. *)

Series[Nsum1,{κ,0,1}]//Normal

(* A few substitution rules to help simplify.  *)

sinhrule={Sinh[κ r]^2-> Cosh[κ r]^2-1,Sinh[κ r]^3-> Sinh[κ r]*(Cosh[κ r]^2-1),Sinh[κ r]^4-> (Cosh[κ r]^2-1)^2,
Sinh[κ r]^5-> Sinh[κ r]*(Cosh[κ r]^2-1)^2,
Sinh[κ r]^6-> (Cosh[κ r]^2-1)^3,Sinh[t κ]^2-> Cosh[t κ]^2-1,Sinh[t κ]^4-> (Cosh[t κ]^2-1)^2};
sinh2rule={Cosh[κ r]^2-1->Sinh[κ r]^2, Cosh[t κ]^2-1->Sinh[t κ]^2};

c1rule={Sinh[κ r]->s1,Cosh[κ r]->c1,Sinh[t κ]->s2,Cosh[t κ]->c2};
c2rule={s1->Sinh[κ r],c1->Cosh[κ r],s2->Sinh[t κ],c2->Cosh[t κ]};

(* Splitting the expression up into terms base on dependence on κ.  *)

Nsum1Ks=Table[Coefficient[Nsum1,κ,n]*κ^n/.sinhrule/.c1rule//Simplify,{n,0,8}]/.c2rule/.sinh2rule

Nsum1Ks//Simplify//Expand//Together

8*Nsum1Ks[[4]]/(Sinh[κ r]*κ^3)//Simplify//Expand//Together

8*Nsum1Ks[[5]]/(κ^4)//Simplify//Expand//Together

8*Nsum1Ks[[6]]/(Sinh[κ r]*κ^5)//Simplify//Expand//Together

-8*Nsum1Ks[[7]]/κ^6//Simplify//Expand//Together

Nsum2Ks=Sum[Nsum1Ks[[n]],{n,1,9}]//Simplify

Nsum3Ks=Series[Nsum2Ks,{r,0,0}]//Normal//Simplify