(*  This program computes the time dependend |beta(u)|^2 for the paper for the Bunch Davies state in pure Sitter space with flat spatial sections.   *)

gamma2 = m^2 - 1/4

gamma = Sqrt[gamma2]

(* This is the Bunch-Davies solution with H = 1. Note that it is H1 because
   we have -eta = Exp[-u] in flat spatial setions.  Also we are assuming
that gamma is real, i.e. m^2 > 1/4.  *)
f = (1/2)*Sqrt[Pi]*Exp[-Pi*gamma/2]*Exp[I*Pi/4]*HankelH1[I*gamma,k*Exp[-u]]

(* Note that because of the way Hankel functions are defined so that H1 = J+i N and H2 = J-iN, to take the complex conjugate of H1 when it has an imaginary order you must compute H2 with the complex conjugate of the imaginary order.  I tried this and it works.  *)

fs = (1/2)*Sqrt[Pi]*Exp[-Pi*gamma/2]*Exp[-I*Pi/4]*HankelH2[-I*gamma,k*Exp[-u]]

fdot = D[f,u]

fdots = D[fs,u]

(* Now check the Wronskian. *)
Wron = Expand[f*fdots-fs*fdot]

N[Wron/.{k->1,m->1,u->10}]

N[Wron/.{k->1,m->1,u->0}]

(* Next compute the WKB terms.  *)

W0 = om

W2 = om+3*omp^2/(8*om^3) - om2p/(4*om^2)

V1 = -omp/om

(* First I'll compute |beta|^2 using W0 *)

beta20 = (1/(2*W0))*(fdot+(I*W0-V1/2)*f)*(fdots+(-I*W0-V1/2)*fs);

(* Next I'll compute it using W2 *)

beta2 = (1/(2*W2))*(fdot+(I*W2-V1/2)*f)*(fdots+(-I*W2-V1/2)*fs)

om = Sqrt[k^2*Exp[-2*u]+gamma2]

omp = D[om,u]

om2p = D[omp,u]

beta20 = Expand[beta20];

beta2 = Expand[beta2];

(* Let's do gamma = 1 so m^2 = 5/4  *)

beta20g1k1 = Expand[beta20/.{m->Sqrt[5/4],k->1}];

beta2g1k1 = Expand[beta2/.{m->Sqrt[5/4],k->1}];

p1 = Plot[{beta2g1k1,beta20g1k1},{u,-10,10}]

(* This is where the particle creation event should occur according to (5.7) when k = 1 and gamma = 1 *)

upc = N[-Log[1.00454]]

(* Now we'll compute the superadiabatic particle number using results from Emil's notebook  *)

w=Sqrt[k^2Exp[-2u] + gam^2]

(* x = -sigma_k.  Fk0 = F_k(0)  *)

x= Sqrt[2gam/Pi] * ( (w/gam) - Log[w + gam] + Log[k] - u)

Fk0 = gam*Pi

(* Nk is the superadiabatic particle number. *)

Nk = Erfc[x]^2*Exp[-2*Fk0]/4

p2 = Plot[Nk/.{gam->1,k->1},{u,-10,10}]

(* In the plot below blue is zeroth order, orange is second order and green is the superadiabatic *)

p3 = Plot[{Re[beta20g1k1],Re[beta2g1k1],Nk/.{gam->1,k->1}},{u,-2,4}]

(* The time the Stokes curve crosses the real axis is supposed to be at *)

ukgam = Log[k/(gam*0.662743)]

N[ukgam/.{k->1,gam->1}]

(* So it looks like this works.  *)

(* To make the particle creation event look like it takes less time we can plot farther out so that *)

p4 = Plot[{Re[beta20g1k1],Re[beta2g1k1],Nk/.{gam->1,k->1}},{u,-4,6}]

p5= Plot[{Re[beta20g1k1],Re[beta2g1k1],Nk/.{gam->1,k->1}},{u,-5,10}]

(* We need to scale Nk by the constant 1/(1-exp(-2pi))  *)

C1 = 1/(1-Exp[-2*Pi*gam])

N[C1/.{gam->1}]

N[Exp[-2*Pi]]

p6= Plot[{Re[beta20g1k1],Re[beta2g1k1],(C1/.{gam->1})*(Nk/.{gam->1,k->1})},{u,-5,10}]

T1 = Table[{N[-10+i/100],N[Re[beta20g1k1/.{u->-10+i/100}]]},{i,0,2001}];

T1[[2001,2]]

ListPlot[T1,Joined->True]

T2 = Table[{N[-10+i/100],N[Re[beta2g1k1/.{u->-10+i/100}]]},{i,0,2001}];

ListPlot[T2,Joined->True]

T3 = Table[{N[-10+i/100],N[(C1/.{gam->1})*(Nk/.{gam->1,k->1,u->-10+i/100})]},{i,0,2001}];

T3[[2001,2]]

ListPlot[T3,Joined->True]

(* Exact |B|^2  with gamma = 1 *)
N[1/(Exp[2*Pi]-1)]

N[Exp[-2*Pi]]

(* I must have saved this data but have no record of doing so on mathematica.  However, the saved data files exist and are the same as the above tables.  *)

T3[[2001,1]]

T3[[1000,1]]

T3[[1000,2]]

T3[[990,1]]

T3[[990,2]]

T2[[1000,1]]

T2[[1000,2]]

T2[[990,1]]

T2[[990,2]]

T1[[1000,1]]

T1[[1000,2]]

T1[[990,2]]

T1[[990,1]]

(* These numbers above all agree with the data files. *)