(* This notebook computes the data used for the two plots that ended up in the final version of the paper JHEP 01, (2022) 192.

  I am letting r = x/H  *)

Vctilde = ((x-xc)/(xc-xb))*((x+xc+xb)/(xc+2*xb))^(xb/(2*xc)) *(xb/x)^(xb*(2*xc+xb)/2)

Ubtilde = -((x-xb)/(xc-xb))*((x+xc+xb)/(2*xc+xb))^(xc/(2*xb)) *(xc/x)^(xc*(xc+2*xb)/2)

DelUbtilde = (Ubtilde/.{x->x1}) - (Ubtilde/.{x->x2})

DelVctilde = (Vctilde/.{x->x1}) - (Vctilde/.{x->x2})

(* Gbar(x1,x2) = 2 pi G(x1,x2) - T(kappa_b+kappa_c) + Log[omega_0^2/(kappa_b*kappa_c)] + 2*gamma_E  *)

Gbar = - Log[Abs[DelUbtilde*DelVctilde]]

(* Now solve for xc in terms of xb using the identity xb^2 + xb*xc + xc^2 = 1 *)

A = xb^2 + xb*xc+xc^2-1

sol = Solve[A==0,xc]

xcc = Part[sol,2,1,2]

Plot[xcc,{xb,0,1}]

Gbar = Gbar/.{xc->xcc}

(* Now locate the cosmological horizon for various values of xb. *)

xcxbp1 = N[xcc/.{xb->1/10}]

xcxbp2 = N[xcc/.{xb->2/10}]

xcxbp3 = N[xcc/.{xb->3/10}]

xcxbp4 = N[xcc/.{xb->4/10}]

xcxbp5 = N[xcc/.{xb->5/10}]

The following is the plot used in Fig. 1.

p1 = DensityPlot[Gbar/.{xb->2/10,xc->xcxbp2},{x1,1/100,2},{x2,1/100,2}, PlotLegends->Automatic,Frame->True,FrameLabel->{H Subscript[r,1],H Subscript[r,2]}]

(* Now I'll put in the velocity correlation function plot. *)

(* So xb goes from 0 to 1/Sqrt[3]  and xc then goes from 1 to 1/Sqrt[3]  *)

xstar = Log[Abs[x-xb]/(xc-xb)]/(2*kapb) -Log[Abs[x-xc]/(xc-xb)]/(2*kapc) + Log[(x+xb+xc)/(xc+2*xb)]/(2*kapN) - (xc/(4*xb*kapb))*Log[(2*xc+xb)/(xc+2*xb)]-(xb*xc/(2*(xc-xb)))*Log[xb/xc]

h = Log[Abs[x-xb]/(xc-xb)]/(2*kapb) +Log[Abs[x-xc]/(xc-xb)]/(2*kapc) + (1/2)*(xc/(xb*kapb)-1/kapN)*Log[(x+xb+xc)/(xc+2*xb)] - (xb*xc/(2*(xc-xb)))*Log[x^2/(xb*xc)] + (xc/(4*(xb*kapb)))*Log[(xc+2*xb)/(2*xc+xb)]

(* Let's check these.  *)

u = Expand[T-h-xstar]

v = Expand[T-h+xstar]

uHc = Expand[u/.{x->xc}]

Assuming[xc>xb>0,FullSimplify[uHc]]

vHb = Expand[v/.{x->xb}]

Assuming[xc>xb>0,FullSimplify[vHb]]

kapb = (xc-xb)*(xc+2*xb)/(2*xb)

kapc = (xc-xb)*(2*xc+xb)/(2*xc)

kapN = (2*xc+xb)*(xc+2*xb)/(2*(xc+xb))

(* Now solve for xc in terms of xb using the identity xb^2 + xb*xc + xc^2 = 1 *)

A = xb^2 + xb*xc+xc^2-1

sol = Solve[A==0,xc]

xcc = Part[sol,2,1,2]

Plot[xcc,{xb,0,1}]

N[1/Sqrt[3]]

N[xcc/.{xb->1/Sqrt[3]}]

(* So xb goes from 0 to 1/sqrt[3] and xc goes from 1 to 1/sqrt[3]  *)

u1 = u/.{x->x1}

u2 = u/.{x->x2}

v1 = v/.{x->x1}

v2 = v/.{x->x2}

(* Every quantity depends on either u2 - u1 or v2 - v1 and once the derivatives are computed we set T2 = T1 = T so 
the T dependence goes away.  Further they are symmetric in u2-u1 or v2-v1.   *)

Delu = u2-u1

Delv = v2 - v1

(* Now I'll compute the parts of \partial_T \partial_T' G  that depend on Delu and Delv separately.  I'll use lt for less than
in the names and gt for greater than.  I'll start with 0 < x1 < x2 < xb  *)

GTTUbxxxb = Assuming[{x1>0,x2>0,xb>0},Simplify[-kapb^2/(8*Pi*Sinh[kapb*Delu/2]^2)]]

(* This is for one point less than xb and the other greater than xb. *)

GTTUbxxbx = Assuming[{x1>0,x2>0,xb>0} ,Simplify[kapb^2/(8*Pi*Cosh[kapb*Delu/2]^2)]]

(* This is for both points greater than xb. *)

GTTUbxbxx = Assuming[{x1>0,x2>0,xb>0} ,Simplify[-kapb^2/(8*Pi*Sinh[kapb*Delu/2]^2)]]

(* This is for both points less than xc. *)

GTTVcxxxc = Assuming[{x1>0,x2>0,xc>0},Simplify[-kapc^2/(8*Pi*Sinh[kapc*Delv/2]^2)]]

This is the component used for Figure 2.

(* This is for one point less than and one point greater than xc. *)

GTTVcxxcx = Assuming[{x1>0,x2>0,xc>0},Simplify[kapc^2/(8*Pi*Cosh[kapc*Delv/2]^2)]]

(* This is for both points greater than xc. *)

GTTVcxcxx = Assuming[{x1>0,x2>0,xc>0},Simplify[-kapc^2/(8*Pi*Sinh[kapc*Delv/2]^2)]]

(* Now locate the cosmological horizon for various values of xb. *)

xcxbp1 = N[xcc/.{xb->1/10}]

xbdivxcp1 = 1/(10*xcxbp1)

xcdivxbp1 = 1/xbdivxcp1

xcxbp1p1 = N[xcc/.{xb->11/100}]

xcxbp1p2 = N[xcc/.{xb->12/100}]

xbdivxcp1p2 = 12/(100*xcxbp1p2)

xcdivxbp1p2 = 1/xbdivxcp1p2

xcxbp1p5 = N[xcc/.{xb->15/100}]

xcxbp2 = N[xcc/.{xb->2/10}]

xbdivxcp2 = 2/(10*xcxbp2)

xcdivxbp2 = 1/xbdivxcp2

xcxbp3 = N[xcc/.{xb->3/10}]

xbdivxcp3 = 3/(10*xcxbp3)

xcdivxbp3 = 1/xbdivxcp3

xcxbp4 = N[xcc/.{xb->4/10}]

xbdivxcp4 = 4/(10*xcxbp4)

xcdivxbp4 = 1/xbdivxcp4

xcxbp5 = N[xcc/.{xb->5/10}]

xbdivxcp5 = 5/(10*xcxbp5)

xcdivxbp5 = 1/xbdivxcp5

(* If one point is on xb then GTTUb = 0 so long as the other point is elsewhere so I just need to plot GTTVc.  Also in order to plot more than one value of xb on the same plot I will change variables to x = xc*y so that y = 1 is always the location of the cosmological horizon.  The
reason for this is that I get a peak when the other point is outside the cosmological horizon.   *)

GTTVcxxxcy = GTTVcxxxc/.{x1->xc*y1,x2->xc*y2}

GTTVcxxcxy = GTTVcxxcx/.{x1->xc*y1,x2->xc*y2}

(* I'll reproduce p1 here so both the plots are at the end. *)

p1

This is the second plot for the published version of the paper.  However, that actual plot was produced by a different program from the data that is computed below.  The original plot was labeled p7.

p2 = LogPlot[{(GTTUbxbxx/.{xb->1/10,xc->xcxbp1,x1->0.9*xcxbp1})+(GTTVcxxcx/.{xb->1/10,xc->xcxbp1,x1->0.9*xcxbp1}),(GTTUbxbxx/.{xb->1/10,xc->xcxbp1,x1->0.95*xcxbp1})+(GTTVcxxcx/.{xb->1/10,xc->xcxbp1,x1->0.95*xcxbp1}),(GTTUbxbxx/.{xb->1/10,xc->xcxbp1,x1->0.98*xcxbp1})+(GTTVcxxcx/.{xb->1/10,xc->xcxbp1,x1->0.98*xcxbp1}),(GTTUbxbxx/.{xb->1/10,xc->xcxbp1,x1->0.99*xcxbp1})+(GTTVcxxcx/.{xb->1/10,xc->xcxbp1,x1->0.99*xcxbp1})},{x2,xcxbp1,30},PlotRange->{10^(-5),2*10^(-3)},Frame->True, FrameLabel->Subscript[x,2]]

(* This is the data that I used for the plot in Fig. 2 of the paper.  The file is called Fig_p7.png.  *)

T7 = Table[{N[i/200],N[(GTTUbxbxx/.{xb->1/10,xc->xcxbp1,x1->0.9*xcxbp1,x2->i/200})+(GTTVcxxcx/.{xb->1/10,xc->xcxbp1,x1->0.9*xcxbp1,x2->i/200})],N[(GTTUbxbxx/.{xb->1/10,xc->xcxbp1,x1->0.95*xcxbp1,x2->i/200})+(GTTVcxxcx/.{xb->1/10,xc->xcxbp1,x1->0.95*xcxbp1,x2->i/200})],N[(GTTUbxbxx/.{xb->1/10,xc->xcxbp1,x1->0.98*xcxbp1,x2->i/200})+(GTTVcxxcx/.{xb->1/10,xc->xcxbp1,x1->0.98*xcxbp1,x2->i/200})],N[(GTTUbxbxx/.{xb->1/10,xc->xcxbp1,x1->0.99*xcxbp1,x2->i/200})+(GTTVcxxcx/.{xb->1/10,xc->xcxbp1,x1->0.99*xcxbp1,x2->i/200})]},{i,6000}];

(* Export["C:\\Users\\anderson\\Tex\\Research - Current\\BH in Cosmology\\2D Stress Tensor\\Numerical_data_for_plots\\Plot7.dat",T7] *)