This is a modification of a program originally written by Zachary Scofield.  The modification is by Paul Anderson.  In this version I am using fits to get p^b_w for large values of T.   The data read in at the beginning for the scattering parameters and the mode functions was computed by Zachary Scofield.

SetDirectory["C:\\Users\\anderson\\Research\\BH_cosmology\\Massive Field in 2D\\Zachs_Final_Notebooks_2022"]

<<"ELFLvalues.m";

<<"ERFRvalues.m";

<<"chifunctions5.m";

Notes

This program is computing PBW and PCW for 3 values of omega (0.1, 1, and 10).  The small omega approximation for the massive case is not the one we should be using, I have not had time to switch everything over. We should instead be  extrapolating to omega = 0 from the value at omega = 10^-4. However, the contribution of this piece to the integral is very small so there is not a large impact on the result. 

Most of the computation in this program is to find the values for PBW and PCW between omega = 10^-4 and 100. The plots at the end are created in real and imaginary parts. The massive plots are created by first subtracting the massless small omega contribution from the massive small omega contribution, and then the integral for PBW/PCW is added on along with the analytic value of PBW/PCW.

The massless plots are analytic solutions for the massless Unruh modes at different points.

Definitions

xb = 1/10;

xc = ((-xb)/2)+Sqrt[1-(3*xb^2)/4];

x0=(101*xb)/100;

x1=(((10^9)-1)*xc)/(10^9);

(* r1 and r2 are the values of x that we want to compute the mode functions *)

r1 = 3/10;

r2 = 6/10;

m=1;

f = -(1/x)*(x-xc)*(x-xb)*(x+xc+xb);

fp = D[f,x];

κb=(1/(2*xb))*(xc-xb)*(xc+(2*xb));

κc=(1/(2*xc))*(xc-xb)*((2*xc)+xb);

κN=(1/(2(xc+xb)))*((2*xc)+xb)*(xc+(2*xb));

xstar = (1/(2*κb))*Log[(x-xb)/(xc-xb)]-(1/(2*κc))*Log[-(x-xc)/(xc-xb)]+(1/(2*κN))*Log[(x+xc+xb)/(xc+(2*xb))]-(xc/(4*xb*κb))*Log[((2*xc)+xb)/(xc+(2*xb))]-(xb*xc)/(2*(xc-xb))*Log[(xb/xc)];

xstar1 = N[xstar/.{x-> x1}];

xstarp = D[xstar,x];

xstarp1 = N[D[xstar,x]/.{x->x1}];

(* He is first solving the equations with omega = 0. *)

omega0 = 0;

(* ct is the value of chi at x = x1 which is near the cosmological horizon. It is equal to 1 since we are taking the limit of chi^c -> cos (w xstar).  The derivative is zero in the limit w -> 0.  *)

ct=(α+β*xstar)/.{α-> 1, β-> 0};

cdt=β/.{β-> 0};

sol1 =NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((omega0^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1] == ct, chi'[x1]==cdt}, chi, {x, x1,x0}];

(* Now he is doing the limit of sin(w xstar).  So chi/w -> xstar and
chi'/w -> 0.  Prime is d/d xstar *)

ct2=(α+β*xstar1)/.{α-> 0, β-> 1};

cdt2=β/.{β-> xstarp1};

sol2 =NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((omega0^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1] == ct2, chi'[x1]==cdt2}, chi, {x, x1,x0}];

These are the values for the scattering coefficients B and D.  They will be fitted so they can be extended to w’ = 0.

Dtable = Table[{ReERfin[[i,1]], ReERfin[[i,2]]-ReFRfin[[i,2]]}, {i,1,991}];

Btable = Table[{ImERfin[[i,1]], -(ImERfin[[i,1]]*(ImERfin[[i,2]]-ImFRfin[[i,2]]))},{i,1,991}];

ScriptD = Fit[Dtable,{1,x,x^2,x^3,x^4,x^5,x^6,x^7, x^8, x^9, x^10},x];

ScriptB = Fit[Btable,{1,x,x^2,x^3,x^4,x^5,x^6,x^7, x^8, x^9, x^10},x];

DValue = ScriptD/.{x-> 0}

BValue = ScriptB/.{x-> 0}

(* Let's test the fits for accuracy *)

ReERfin[[991,1]]

ListPlot[Btable]

TestB = Table[{Btable[[i,1]],Abs[(Btable[[i,2]]-(ScriptB/.{x->Btable[[i,1]]}))/Btable[[i,2]]]},{i,991}];

ListPlot[TestB,PlotRange->All]

ListPlot[Dtable,Joined->True]

TestD = Table[{Dtable[[i,1]],Abs[(Dtable[[i,2]]-(ScriptD/.{x->Dtable[[i,1]]}))/Dtable[[i,2]]]},{i,991}];

ListPlot[TestD,PlotRange->All]

These are the Bogolubov coefficients.

αcw = (1/(2*Pi))*Sqrt[omegap/omega]*(κc^(-1-(I*(omegap/κc))))*((Gamma[δ-(I*(omegap/κc))])/(((-I*omega)+ϵ)^(-I*(omegap/κc))));

βcw = (1/(2*Pi))*Sqrt[omegap/omega]*(κc^(-1+(I*(omegap/κc))))*((Gamma[δ+(I*(omegap/κc))])/(((-I*omega)+ϵ)^(I*(omegap/κc))));

αbw = (1/(2*Pi))*Sqrt[omegap/omega]*(κb^(-1-(I*(omegap/κb))))*((Gamma[δ-(I*(omegap/κb))])/(((-I*omega)+ϵ)^(-I*(omegap/κb))));

βbw =(1/(2*Pi))*Sqrt[omegap/omega]*(κb^(-1+(I*(omegap/κb))))*((Gamma[δ+(I*(omegap/κb))])/(((-I*omega)+ϵ)^(I*(omegap/κb))));

hx = (1/(2*κb))*Log[Abs[(x-xb)]/(xc-xb)]+(1/(2*κc))*Log[Abs[(x-xc)]/(xc-xb)]+(1/2)*((xc)/(xb*κb)-(1/κN))*Log[(x+xc+xb)/(xc + 2*xb)]-((xb*xc)/(2*(xc-xb)))*Log[(x^2)/(xc*xb)]+(xc/(4*xb*κb))*Log[(xc+2*xb)/(2*xc+xb)];

u=(T-hx)-xstar;

v=(T-hx)+xstar;

N[hx/.{x-> r1}]

omegatable1 = Table[i,{i,1/10000,1/100,1/100000}];

omegatable2 = Table[i,{i,1/100,1,1/1000}];

omegatable3 = Table[i,{i,1,6,5/1000}];

omegatable4 = Table[i,{i,6,5000/100,4/100}];

omegatable5 = Table[i,{i,5000/100,10000/100,5/100}];

chi10 = (chi[x]/.sol1)

chi20 = (chi[x]/.sol2)

omegavalues = Table[10^(i), {i,-1,1,1}]

(* Tvalues = Table[i,{i,-1,100,1}]; *)

Tvalues = Table[5*i,{i,0,20}]

Length[Tvalues]

Omega Prime Delta Values:

Delta1 = omegatable1[[2]]-omegatable1[[1]]

Delta2 = omegatable2[[2]]-omegatable2[[1]]

Delta3 = omegatable3[[2]]-omegatable3[[1]]

Delta4 = omegatable4[[2]]-omegatable4[[1]]

Delta5 = omegatable5[[2]]-omegatable5[[1]]

Pbw Calculations:

Small omega prime calculations:

• Massive small omega prime contribution:

m1swr1 = Table[-((chi10/.{x-> r1})*(10^(-4)))/(Pi^(3/2)*(omegavalues[[i]])^(1/2)*(BValue)),{i,1,3,1}]

Exact integrand for p^b_w  at r = r1 for the massless case.

Int1r1 =  Table[{(I/(2*Pi*Sqrt[4*Pi*omega]))*(1/omegap)*(Exp[-I*omegap*u]*Gamma[1-I*omegap/κb]*(-I*omega/κb)^(I*omegap/κb) - 1)/.{omega-> omegavalues[[i]],T-> Tvalues[[k]],x-> r1}},{k,1,1,1},{i,1,3,1}];

TableI1r1 = Table[{NIntegrate[Int1r1[[k,i]],{omegap,-(10^-4),10^-4}]},{k,1,1,1},{i,1,3,1}];

m0swr1 = Table[{omegavalues[[i]],(1/(2*Sqrt[4*Pi*omegavalues[[i]]]) + TableI1r1[[k,i]])}, {k,1,1,1},{i,1,3,1}];

Massless pbw exact

Ub = (-(1/κb)*Exp[-κb*u]);

Pbwm0r1 = Table[{((1/Sqrt[4*Pi*omega])*Exp[-I*omega*Ub])/.{x-> r1,T-> Tvalues[[k]],omega-> omegavalues[[i]]}},{k,1,1,1},{i,1,3,1}];

Pbwm0r2 = Table[{((1/Sqrt[4*Pi*omega])*Exp[-I*omega*Ub])/.{x-> r2,T-> Tvalues[[k]],omega-> omegavalues[[i]]}},{k,1,1,1},{i,1,3,1}];

Integrand calculations:

Here I am going to deviate from Zach and make separate tables for each w and T so that I can plot the integrands.  I’ll start with w = 1/10 and T = 0
bBand1int is the integrand for the band 1 values of w’.

bBand1int = Table[{N[omegatable1[[j]]],((((1/(2*Pi))*Sqrt[N[omegatable1[[j]]]/omega]*(κb^(-1-(I*(N[omegatable1[[j]]]/κb))))*(Gamma[δ-(I*(N[omegatable1[[j]]]/κb))]/((-I*omega)+ϵ)^(-I*(N[omegatable1[[j]]]/κb))))*(((Exp[-I*N[omegatable1[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable1[[j]]]])*(1/(ReERfin[[j,2]]+I*(ImERfin[[j,2]])))*(Part[chi[x]/.Part[chiR1[[j]],1],1]))-((Exp[-I*N[omegatable1[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable1[[j]]]])))+(((1/(2*Pi))*Sqrt[N[omegatable1[[j]]]/omega]*(κb^(-1+(I*(N[omegatable1[[j]]]/κb))))*(Gamma[δ+(I*(N[omegatable1[[j]]]/κb))]/((-I*omega)+ϵ)^(I*(N[omegatable1[[j]]]/κb))))*(((Exp[I*N[omegatable1[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable1[[j]]]])*(1/((ReERfin[[j,2]])-I*(ImERfin[[j,2]])))*(Conjugate[(Part[chi[x]/.Part[chiR1[[j]],1],1])]))-((Exp[I*N[omegatable1[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable1[[j]]]]))))},{j,1,991}];

OK we can plot this as a function of T if we fix the value of w’ as the lowest value.

bBand1Twp1m4 = bBand1int[[1,2]]/.{omega->1/10, ϵ->0,δ->0,x->r1};

The explanation is that I am subtracting the massless modes and there is a pole at w’ = 0.  That is why it has such a large amplitude.  But I have to be careful because this is going to cause a lot of oscillation in  w’ for the integral over w’ and so I may not end up computing the integral that accurately.  I’ll proceed cautiously!

bBand1intw0p1T0 = bBand1int/.{omega->omegavalues[[1]],T->Tvalues[[1]],x->r1,ϵ->0,δ->0};

ImbBand1intw0p1T0 = Table[{bBand1intw0p1T0[[i,1]], Im[bBand1intw0p1T0[[i,2]]]},{i,1,991}];

Now we’ll do time T = 10

bBand1intw0p1Tm10 = bBand1int/.{omega->omegavalues[[1]],T->Tvalues[[2]],x->r1,ϵ->0,δ->0};

ImbBand1intw0p1Tm10 = Table[{bBand1intw0p1Tm10[[i,1]], Im[bBand1intw0p1Tm10[[i,2]]]},{i,1,991}];

bBand1intw0p1Tm100 = bBand1int/.{omega->omegavalues[[1]],T->Tvalues[[3]],x->r1,ϵ->0,δ->0};

ImbBand1intw0p1Tm100 = Table[{bBand1intw0p1Tm100[[i,1]], Im[bBand1intw0p1Tm100[[i,2]]]},{i,1,991}];

bBand2int = Table[{N[omegatable2[[j]]],((((1/(2*Pi))*Sqrt[N[omegatable2[[j]]]/omega]*(κb^(-1-(I*(N[omegatable2[[j]]]/κb))))*(Gamma[δ-(I*(N[omegatable2[[j]]]/κb))]/((-I*omega)+ϵ)^(-I*(N[omegatable2[[j]]]/κb))))*(((Exp[-I*N[omegatable2[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable2[[j]]]])*(1/(ReERfinb[[j,2]]+I*(ImERfinb[[j,2]])))*(Part[chi[x]/.Part[chiR2[[j]],1],1]))-((Exp[-I*N[omegatable2[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable2[[j]]]])))+((1/(2*Pi))*Sqrt[N[omegatable2[[j]]]/omega]*(κb^(-1+(I*(N[omegatable2[[j]]]/κb))))*(Gamma[δ+(I*(N[omegatable2[[j]]]/κb))]/((-I*omega)+ϵ)^(I*(N[omegatable2[[j]]]/κb))))*(((Exp[I*N[omegatable2[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable2[[j]]]])*(1/((ReERfinb[[j,2]])-I*(ImERfinb[[j,2]])))*(Conjugate[(Part[chi[x]/.Part[chiR2[[j]],1],1])]))-((Exp[I*N[omegatable2[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable2[[j]]]])))},{j,1,991}];

bBand2intw0p1Tm1 = bBand2int/.{omega->omegavalues[[1]],T->Tvalues[[1]],x->r1,ϵ->0,δ->0};

ImbBand2intw0p1Tm1 = Table[{bBand2intw0p1Tm1[[i,1]], Im[bBand2intw0p1Tm1[[i,2]]]},{i,1,991}];

bBand2intw0p1Tm10 = bBand2int/.{omega->omegavalues[[1]],T->Tvalues[[2]],x->r1,ϵ->0,δ->0};

ImbBand2intw0p1Tm10 = Table[{bBand2intw0p1Tm10[[i,1]], Im[bBand2intw0p1Tm10[[i,2]]]},{i,1,991}];

bBand2intw0p1Tm100 = bBand2int/.{omega->omegavalues[[1]],T->Tvalues[[3]],x->r1,ϵ->0,δ->0};

ImbBand2intw0p1Tm100 = Table[{bBand2intw0p1Tm100[[i,1]], Im[bBand2intw0p1Tm100[[i,2]]]},{i,1,991}];

OK, I have enough values of w’ for T = -1 and -10 but not for -100.  I’ll go ahead for now anyway.

bBand3int = Table[{N[omegatable3[[j]]],((((1/(2*Pi))*Sqrt[N[omegatable3[[j]]]/omega]*(κb^(-1-(I*(N[omegatable3[[j]]]/κb))))*(Gamma[δ-(I*(N[omegatable3[[j]]]/κb))]/((-I*omega)+ϵ)^(-I*(N[omegatable3[[j]]]/κb))))*(((Exp[-I*N[omegatable3[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable3[[j]]]])*(1/(ReERfinc[[j,2]]+I*(ImERfinc[[j,2]])))*(Part[chi[x]/.Part[chiR3[[j]],1],1]))-((Exp[-I*N[omegatable3[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable3[[j]]]])))+((1/(2*Pi))*Sqrt[N[omegatable3[[j]]]/omega]*(κb^(-1+(I*(N[omegatable3[[j]]]/κb))))*(Gamma[δ+(I*(N[omegatable3[[j]]]/κb))]/((-I*omega)+ϵ)^(I*(N[omegatable3[[j]]]/κb))))*(((Exp[I*N[omegatable3[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable3[[j]]]])*(1/((ReERfinc[[j,2]])-I*(ImERfinc[[j,2]])))*(Conjugate[(Part[chi[x]/.Part[chiR3[[j]],1],1])]))-((Exp[I*N[omegatable3[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable3[[j]]]])))},{j,1,1001}];

bBand3intw0p1Tm1 = bBand3int/.{omega->omegavalues[[1]],T->Tvalues[[1]],x->r1,ϵ->0,δ->0};

ImbBand3intw0p1Tm1 = Table[{bBand3intw0p1Tm1[[i,1]], Im[bBand3intw0p1Tm1[[i,2]]]},{i,1,991}];

bBand3intw0p1Tm10 = bBand3int/.{omega->omegavalues[[1]],T->Tvalues[[2]],x->r1,ϵ->0,δ->0};

ImbBand3intw0p1Tm10 = Table[{bBand3intw0p1Tm10[[i,1]], Im[bBand3intw0p1Tm10[[i,2]]]},{i,1,991}];

bBand3intw0p1Tm100 = bBand3int/.{omega->omegavalues[[1]],T->Tvalues[[3]],x->r1,ϵ->0,δ->0};

ImbBand3intw0p1Tm100 = Table[{bBand3intw0p1Tm100[[i,1]], Im[bBand3intw0p1Tm100[[i,2]]]},{i,1,991}];

OK, I might be able to proceed with T = -10 but definitely no point in doing T = -100.  Not nearly enough values of w’.  I am going to have to use fits to go to large negative T.

bBand4int = Table[{N[omegatable4[[j]]],((((1/(2*Pi))*Sqrt[N[omegatable4[[j]]]/omega]*(κb^(-1-(I*(N[omegatable4[[j]]]/κb))))*(Gamma[δ-(I*(N[omegatable4[[j]]]/κb))]/((-I*omega)+ϵ)^(-I*(N[omegatable4[[j]]]/κb))))*(((Exp[-I*N[omegatable4[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable4[[j]]]])*(1/(ReERfind[[j,2]]+I*(ImERfind[[j,2]])))*(Part[chi[x]/.Part[chiR4[[j]],1],1]))-((Exp[-I*N[omegatable4[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable4[[j]]]])))+((1/(2*Pi))*Sqrt[N[omegatable4[[j]]]/omega]*(κb^(-1+(I*(N[omegatable4[[j]]]/κb))))*(Gamma[δ+(I*(N[omegatable4[[j]]]/κb))]/((-I*omega)+ϵ)^(I*(N[omegatable4[[j]]]/κb))))*(((Exp[I*N[omegatable4[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable4[[j]]]])*(1/((ReERfind[[j,2]])-I*(ImERfind[[j,2]])))*(Conjugate[(Part[chi[x]/.Part[chiR4[[j]],1],1])]))-((Exp[I*N[omegatable4[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable4[[j]]]])))},{j,1,1101}];

bBand4intw0p1Tm1 = bBand4int/.{omega->omegavalues[[1]],T->Tvalues[[1]],x->r1,ϵ->0,δ->0};

ImbBand4intw0p1Tm1 = Table[{bBand4intw0p1Tm1[[i,1]], Im[bBand4intw0p1Tm1[[i,2]]]},{i,1,1101}];

bBand4intw0p1Tm10 = bBand4int/.{omega->omegavalues[[1]],T->Tvalues[[2]],x->r1,ϵ->0,δ->0};

ImbBand4intw0p1Tm10 = Table[{bBand4intw0p1Tm10[[i,1]], Im[bBand4intw0p1Tm10[[i,2]]]},{i,1,1101}];

bBand4intw0p1Tm100 = bBand4int/.{omega->omegavalues[[1]],T->Tvalues[[3]],x->r1,ϵ->0,δ->0};

ImbBand4intw0p1Tm100 = Table[{bBand4intw0p1Tm100[[i,1]], Im[bBand4intw0p1Tm100[[i,2]]]},{i,1,1101}];

This really doesn’t even work for T = -10 and barely works well for T = -1.

bBand5int = Table[{N[omegatable5[[j]]],((((1/(2*Pi))*Sqrt[N[omegatable5[[j]]]/omega]*(κb^(-1-(I*(N[omegatable5[[j]]]/κb))))*(Gamma[δ-(I*(N[omegatable5[[j]]]/κb))]/((-I*omega)+ϵ)^(-I*(N[omegatable5[[j]]]/κb))))*(((Exp[-I*N[omegatable5[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable5[[j]]]])*(1/(ReERfine[[j,2]]+I*(ImERfine[[j,2]])))*(Part[chi[x]/.Part[chiR5[[j]],1],1]))-((Exp[-I*N[omegatable5[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable5[[j]]]])))+((1/(2*Pi))*Sqrt[N[omegatable5[[j]]]/omega]*(κb^(-1+(I*(N[omegatable5[[j]]]/κb))))*(Gamma[δ+(I*(N[omegatable5[[j]]]/κb))]/((-I*omega)+ϵ)^(I*(N[omegatable5[[j]]]/κb))))*(((Exp[I*N[omegatable5[[j]]]*(T - hx)]/Sqrt[4*Pi*N[omegatable5[[j]]]])*(1/((ReERfine[[j,2]])-I*(ImERfine[[j,2]])))*(Conjugate[(Part[chi[x]/.Part[chiR5[[j]],1],1])]))-((Exp[I*N[omegatable4[[j]]]*((T-hx)-xstar)])/Sqrt[4*Pi*N[omegatable5[[j]]]])))},{j,1,1001}];

bBand5intw0p1Tm1 = bBand5int/.{omega->omegavalues[[1]],T->Tvalues[[1]],x->r1,ϵ->0,δ->0};

ImbBand5intw0p1Tm1 = Table[{bBand5intw0p1Tm1[[i,1]], Im[bBand5intw0p1Tm1[[i,2]]]},{i,1,991}];

bBand5intw0p1Tm10 = bBand5int/.{omega->omegavalues[[1]],T->Tvalues[[2]],x->r1,ϵ->0,δ->0};

ImbBand5intw0p1Tm10 = Table[{bBand5intw0p1Tm10[[i,1]], Im[bBand5intw0p1Tm10[[i,2]]]},{i,1,991}];

bBand5intw0p1Tm100 = bBand5int/.{omega->omegavalues[[1]],T->Tvalues[[3]],x->r1,ϵ->0,δ->0};

ImbBand5intw0p1Tm100 = Table[{bBand5intw0p1Tm100[[i,1]], Im[bBand5intw0p1Tm100[[i,2]]]},{i,1,991}];

From here I will deviate from Zach.  What I want to do is to integrate over the bands using using fits for the radial mode functions.  That is the only way I can get to large positive and negative values of T.

I’ll just fit the relevant part which is the radial mode function plus scattering coefficients.

chib1T = Table[{N[omegatable1[[j]]],((1/(ReERfin[[j,2]]+I*(ImERfin[[j,2]])))*(Part[chi[x]/.Part[chiR1[[j]],1],1]))/.{x->r1}},{j,1,991}];

fit1 = Fit[chib1T,{wp,wp^2,wp^3,wp^4,wp^5,wp^6,wp^7,wp^8,wp^9,wp^10},wp]

dif1 = Table[{N[omegatable1[[j]]],(Re[chib1T[[j,2]]]-Re[(fit1/.{wp->omegatable1[[j]]})])/Re[chib1T[[j,2]]]},{j,1,991}];

dif1Im = Table[{N[omegatable1[[j]]],(Im[chib1T[[j,2]]]-Im[(fit1/.{wp->omegatable1[[j]]})])/Im[chib1T[[j,2]]]},{j,1,991}];

dif1[[3]]

chib2T = Table[{N[omegatable2[[j]]],((1/(ReERfinb[[j,2]]+I*(ImERfinb[[j,2]])))*(Part[chi[x]/.Part[chiR2[[j]],1],1]))/.{x->r1}},{j,1,991}];

N[omegatable2[[991]]]

fit2 = Fit[chib2T,{wp,wp^2,wp^3,wp^4,wp^5,wp^6,wp^7,wp^8,wp^9,wp^10,wp^11,wp^12,wp^13,wp^14,wp^15},wp]

dif2 = Table[{N[omegatable2[[j]]],(Re[chib2T[[j,2]]]-Re[(fit2/.{wp->omegatable2[[j]]})])/Re[chib2T[[j,2]]]},{j,1,991}];

dif2Im = Table[{N[omegatable2[[j]]],(Im[chib2T[[j,2]]]-Im[(fit2/.{wp->omegatable2[[j]]})])/Im[chib2T[[j,2]]]},{j,1,991}];

ListPlot[{dif2},PlotRange->All]

ListPlot[dif2Im]

chib3T = Table[{N[omegatable3[[j]]],((1/(ReERfinc[[j,2]]+I*(ImERfinc[[j,2]])))*(Part[chi[x]/.Part[chiR3[[j]],1],1]))/.{x->r1}},{j,1,1001}];

N[omegatable3[[1001]]]

fit3 = Fit[chib3T,{1,1/wp,1/wp^2,1/wp^3,1/wp^4,1/wp^5,1/wp^6,1/wp^7,1/wp^8,1/wp^9,1/wp^10,1/wp^11,1/wp^12,1/wp^13,1/wp^14,1/wp^15},wp]

dif3 = Table[{N[omegatable3[[j]]],(Re[chib3T[[j,2]]]-Re[(fit3/.{wp->omegatable3[[j]]})])/Re[chib3T[[j,2]]]},{j,1,1001}];

dif3Im = Table[{N[omegatable3[[j]]],(Im[chib3T[[j,2]]]-Im[(fit3/.{wp->omegatable3[[j]]})])/Im[chib3T[[j,2]]]},{j,1,1001}];

dif3[[1]]

ListPlot[dif3,PlotRange->All]

ListPlot[dif3Im]

chib4T = Table[{N[omegatable4[[j]]],((1/(ReERfind[[j,2]]+I*(ImERfind[[j,2]])))*(Part[chi[x]/.Part[chiR4[[j]],1],1]))/.{x->r1}},{j,1,1101}];

N[omegatable4[[1]]]

N[omegatable4[[1101]]]

fit4 = Fit[chib4T,{1,1/wp,1/wp^2,1/wp^3,1/wp^4,1/wp^5,1/wp^6,1/wp^7,1/wp^8,1/wp^9,1/wp^10,1/wp^11,1/wp^12,1/wp^13,1/wp^14,1/wp^15,1/wp^16,1/wp^17,1/wp^18,1/wp^19,1/wp^20},wp]

dif4 = Table[{N[omegatable4[[j]]],(Re[chib4T[[j,2]]]-Re[(fit4/.{wp->omegatable4[[j]]})])/Re[chib4T[[j,2]]]},{j,1,1101}];

dif4Im = Table[{N[omegatable4[[j]]],(Im[chib4T[[j,2]]]-Im[(fit4/.{wp->omegatable4[[j]]})])/Im[chib4T[[j,2]]]},{j,1,1101}];

dif4[[1]]

ListPlot[dif4]

ListPlot[dif4Im]

chib5T = Table[{N[omegatable5[[j]]],((1/(ReERfine[[j,2]]+I*(ImERfine[[j,2]])))*(Part[chi[x]/.Part[chiR5[[j]],1],1]))/.{x->r1}},{j,1,1001}];

N[omegatable5[[1]]]

N[omegatable5[[1001]]]

fit5 = Fit[chib5T,{1,1/wp,1/wp^2,1/wp^3,1/wp^4,1/wp^5,1/wp^6,1/wp^7,1/wp^8,1/wp^9,1/wp^10,1/wp^11,1/wp^12,1/wp^13,1/wp^14,1/wp^15,1/wp^16,1/wp^17,1/wp^18,1/wp^19,1/wp^20},wp]

dif5 = Table[{N[omegatable5[[j]]],(Re[chib5T[[j,2]]]-Re[(fit5/.{wp->omegatable5[[j]]})])/Re[chib5T[[j,2]]]},{j,1,1001}];

dif5Im = Table[{N[omegatable5[[j]]],(Im[chib5T[[j,2]]]-Im[(fit5/.{wp->omegatable5[[j]]})])/Im[chib5T[[j,2]]]},{j,1,1001}];

dif5[[1]]

ListPlot[dif5]

ListPlot[dif5Im]

Save["Fits_for_chib.m",fit1,fit2,fit3,fit4,fit5]