(* In this version I will compute N1, N2, and N3 and see if I get a more accurate answer than doing the whole thing by brute force.  
 This is for 2D.  
Here I compute N_{jn} for the case with no effective potential.  This version has the same parameters as Ray's 4D version.  It has j = 1 and ep = 1/10.
I'll set v_0 = 0 and kappa = 1/4.  I'll set v0 = 0 and kappa = 1/4, i.e. M = 1 which is what Ray used.
Therefore vH = -4M = -4
I'm effectively scaling M out of the equations by defining \bar(w) = w/M and bar(v0) = M v0  so that v0 w  = bar(v0) bar(w).  
  I'll use in this version the variable z = log[w'].   *)

(* This is the exact beta except without the exp(-i wp vH) because that cancels out when the particle number is computed.  *)

bwwp = -Exp[-I*w*vH]*(Sqrt[w*wp]/(2*Pi*kap*(w+wp)))*(-I*(w+wp)/kap)^(-I*w/kap) * Gamma[I*w/kap]

(*  What happened is that for some reason I don't get the full coverage in z unless I use the form here which factors out a wp from the w+wp factors. Note z = (1/kap)*Log(wp/kap)  so wp = kap*exp(kap*z) *)

bwwpz = -Exp[-I*w*vH]*(Sqrt[w]*Exp[-kap*z/2]/(2*Pi*kap^(3/2)*(1+(w/kap)*Exp[-kap*z])))*Exp[-I*w*z]*Exp[-Pi*w/(2*kap)]*(1+(w/kap)*Exp[-kap*z])^(-I*w/kap) * Gamma[I*w/kap]

(* We need B(w,wp) and its complex conjugate.  *)

B = -Exp[-I*w*vH]*(Sqrt[w]/(2*Pi*kap*(1+(w/kap)*Exp[-kap*z])))*Exp[-Pi*w/(2*kap)]*(1+(w/kap)*Exp[-kap*z])^(-I*w/kap) * Gamma[I*w/kap]

Bs = -Exp[I*w*vH]*(Sqrt[w]/(2*Pi*kap*(1+(w/kap)*Exp[-kap*z])))*Exp[-Pi*w/(2*kap)]*(1+(w/kap)*Exp[-kap*z])^(I*w/kap) * Gamma[-I*w/kap]

BH = - Exp[-I*w*vH]*Sqrt[w]*Exp[-Pi*w/(2*kap)]*Gamma[I*w/kap]/(2*Pi*kap)

BHs =  - Exp[I*w*vH]*Sqrt[w]*Exp[-Pi*w/(2*kap)]*Gamma[-I*w/kap]/(2*Pi*kap)

(* It works better if I put the a factor of kap Exp[kap*z/2] here rather than inserting
a factor of kap*Exp[kap*z] directly into the integral over z.  This comes from dwp = kap^2*exp[kap*z] dz.  *)

betajnzint = (1/Sqrt[ep])*Exp[2*Pi*I*w*n/ep]*bwwpz*kap*Exp[kap*z/2]

(* This is what we need for N3b.  I think I must have changed B0 to BH after doing this?  There are no more B0's here. *)

bwwpzBH = Exp[-kap*z/2]*Exp[-I*w*z]*B0/Sqrt[kap]

betajnzintBH = (1/Sqrt[ep])*Exp[2*Pi*I*w*n/ep]*bwwpzB0*kap*Exp[kap*z/2]

(* Next we'll try integrating over various ranges of z up to Lam = 10 which corresponds to Z = kap*Exp[Lam/kap]*)

Z = kap*Log[Lam/kap]

ZLam10 =N[Z/.{kap->1/4,Lam->10},15]

(* The stationary phase point is at  ()

wps = kap*Exp[2*Pi*n*kap/ep]

zwps =kap*Log[wps/kap]

N[wps/.{kap->1/4,ep->1/10,n->1}]

N[zwps/.{kap->1/4,ep->1/10,n->1}]

N[zwps/.{kap->1/4,ep->1/10,n->10}]

Ta = Table[{i,NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-Infinity,-100},WorkingPrecision ->15]},{i,1,10}]

Tb = Table[{i,NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-100,-10},WorkingPrecision ->15]},{i,1,10}]

Tc = Table[{i,NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10,0},WorkingPrecision ->15]},{i,1,10}]

Td = Table[{i,NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,0,ZLam10},WorkingPrecision ->15]},{i,1,10}]

TN1 = Table[{i,Ta[[i,2]]+Tb[[i,2]]+Tc[[i,2]]+Td[[i,2]]},{i,10}]

Tn2N1 = Table[{i,i^2*TN1[[i,2]]},{i,10}]

ListPlot[Tn2N1]

(* Now let's do N2.  We need the integrand for just using BH.  *)

N2int = (kap/ep)*Exp[I*(2*Pi*n/ep-z)*(w1-w2)]*((B/.{w->w1})*(Bs/.{w->w2})-(BH/.{w->w1})*(BHs/.{w->w2}))

ZLam10

(* The problem is that in integrating over z the imaginary part which is really zero is slowing things down.  Taking the real part of the integrand doesn't help.  So I'll just do two NIntegrates as I do above.  *)

T2a = Table[{i,NIntegrate[Re[NIntegrate[N2int/.{kap->1/4,ep->1/10,n->i,vH->-4,z->zz},{w1,1/10,2/10},{w2,1/10,2/10}]],{zz,ZLam10,10}]},{i,10}]

T2an2 = Table[{i,i^2*T2a[[i,2]]},{i,10}]

ListPlot[T2an2]

(* T2b refuses to evaluate in a reasonable time for 10 to 100.  So I probably need to look at the integrand in this range for n = 10.  *)

T2bintegrand = Table[{N[k], Re[NIntegrate[N2int/.{kap->1/4,ep->1/10,n->10,vH->-4,z->k},{w1,1/10,2/10},{w2,1/10,2/10},WorkingPrecision->15]]},{k,10,50}]

ListPlot[T2bintegrand]

(* I simply cannot get the triple integral below to work!  It takes too long. *)

T2b = Table[{i,NIntegrate[Re[NIntegrate[N2int/.{kap->1/4,ep->1/10,n->i,vH->-4,z->zz},{w1,1/10,2/10},{w2,1/10,2/10},WorkingPrecision->15]],{zz,10,20}]},{i,10}]

(* What I will try to do is to use Simpson's rule on a somewhat more dense table than the one generated above.  I'll start with the one generated above. *)

Simpsonn10 =  T2bintegrand[[1,2]]/3 +T2bintegrand[[41,2]]/3 + (4/3)*Sum[T2bintegrand[[j,2]],{j,2,40,2}]+(2/3)*Sum[T2bintegrand[[j,2]],{j,3,39,2}]

(* I'll go from z = 10 to 40.
    Even this with k from 0 to 2 took too long! *)

T2bint = Table[{i,N[10+k/10], Re[NIntegrate[N2int/.{kap->1/4,ep->1/10,n->i,vH->-4,z->10+k/10},{w1,1/10,2/10},{w2,1/10,2/10},WorkingPrecision->15]]},{k,0,2},{i,10}]

(* OK, what I'll do is use the fact that B is proportional to BH and compare the value of the resulting factor to see how different they are compared to the leading order approximation for this factor for z = 10.  If I can use the leading order value then it should shorten things dramatically. *)

facdif = (1+(w/kap)*Exp[- kap*z])^(-1-I*w/kap) - 1

facdifleading = -(w/kap)*(1+ I*w/kap)*Exp[-kap*z]

N[facdif/.{kap->1/4,w->2/10,z->10}]

N[facdifleading/.{kap->1/4,w->2/10,z->10}]

(* Ok they are not that close.   *)

(* I think the only reasonable solution is to compute the two integrals separately and then subtract them! So I'll try that next for the previous range and see how well it works. *)

T2aexactmodes = Table[{i,NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,ZLam10,10},WorkingPrecision ->15]},{i,1,10}]

T2aHawking = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,ZLam10,10},WorkingPrecision ->15]},{i,1,10}]

T2anew = Table[{i,T2aexactmodes[[i,2]]-T2aHawking[[i,2]]},{i,10}]

T2a

(* OK it works!  Now let's try the next band *)

T2bexactmodes = Table[{i,NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,10,10^2},WorkingPrecision ->15]},{i,1,10}]

T2bHawking = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,10,10^2},WorkingPrecision ->15]},{i,1,10}]

T2b = Table[{i,T2bexactmodes[[i,2]]-T2bHawking[[i,2]]},{i,10}]

T2cexactmodes = Table[{i,NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,10^2,10^3},WorkingPrecision ->15]},{i,1,10}]

T2cHawking = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,10^2,10^3},WorkingPrecision ->15]},{i,1,10}]

T2c = Table[{i,T2cexactmodes[[i,2]]-T2cHawking[[i,2]]},{i,10}]

TN2 = Table[{i,T2a[[i,2]]+T2b[[i,2]]+T2c[[i,2]]},{i,10}]

Tn2N2 = Table[{i,i^2*TN2[[i,2]]},{i,10}]

ListPlot[Tn2N2]

T3a = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-Infinity,-10^4},WorkingPrecision ->15]},{i,1,10}]

T3b = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10^4,-10^3},WorkingPrecision ->15]},{i,1,10}]

T3c = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10^3,-10^2},WorkingPrecision ->15]},{i,1,10}]

T3d = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10^2,-10},WorkingPrecision ->15]},{i,1,10}]

T3e = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10,0},WorkingPrecision ->15]},{i,1,10}]

T3f = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,0,ZLam10},WorkingPrecision ->15]},{i,1,10}]

(* I forgot a minus sign.  I'll put it in here. *)

TN3 = Table[{i,-(T3a[[i,2]]+T3b[[i,2]]+T3c[[i,2]]+T3d[[i,2]]+T3e[[i,2]]+T3f[[i,2]])},{i,10}]

Tn1N3 = Table[{i,i*TN3[[i,2]]},{i,10}]

p1 = ListPlot[Tn1N3]

(*  Now let's compare to the asympotic analytic value for n N_{3b} *)

Njn1 = -(kap/(2*Pi* n))*((B0/.{w->(j+1)*ep})*(B0s/.{w->(j+1)*ep})+(B0/.{w->j*ep})*(B0s/.{w->j*ep}))

nNjn1 = n*Njn1

nNjn1plot = Re[N[nNjn1/.{j->1,ep->1/10,kap->1/4},15]]

p2 = Plot[nNjn1plot,{n,0,10}]

Show[p1,p2]

(* Let's see if we can get just a bit more accurate. *)

T3a1 = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-Infinity,-10^5},WorkingPrecision ->15]},{i,1,10}]

T3a2 = Table[{i,NIntegrate[Abs[NIntegrate[betajnzintBH/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10^5,-10^4},WorkingPrecision ->15]},{i,1,10}]

T3atot = Table[{i,T3a1[[i,2]]+T3a2[[i,2]]},{i,10}]

T3a

TN3better = Table[{i,-(T3atot[[i,2]]+T3b[[i,2]]+T3c[[i,2]]+T3d[[i,2]]+T3e[[i,2]]+T3f[[i,2]])},{i,10}]

Tn1N3better = Table[{i,i*TN3better[[i,2]]},{i,10}]

p3 = ListPlot[{Tn1N3,Tn1N3better}]

Show[p3,p2]

Save["C:\\Users\\anderson\\Research\\BH_evaporation\\Particle_Production\\2019\\N1N2N32Dresults.m",{TN1,TN3better,nNjn1,nNjn1plot,TN2}]

Export["C:\\Users\\anderson\\Research\\BH_evaporation\\Particle_Production\\2019\\N12D.dat",TN1]

Export["C:\\Users\\anderson\\Research\\BH_evaporation\\Particle_Production\\2019\\N22D.dat",TN2]

Export["C:\\Users\\anderson\\Research\\BH_evaporation\\Particle_Production\\2019\\N32D.dat",TN3better]