Friction and Diffusion Coefficients for Sigmaplus-Sigmaminus Polarized Fields - Energy and Velocity Distributions

This notebook has expressions needed in the appendix of Chapter 18.  This notebook is for (SubPlus[σ])-(SubMinus[σ]) polarization and a G=1 to H=2 transition

Program

<<c:\pberman\mathematica\5.2\clebsh.m

<<c:\pberman\mathematica\5.2\threej.m

<<c:\pberman\mathematica\5.2\ninej.m

<<c:\pberman\mathematica\5.2\sixj.m

<<c:\pberman\mathematica\5.2\prate.m

sq[x_]:=Sqrt[2 x+1]

B[K_]:=Sqrt[16 Pi/3]kron[K,0]+Sqrt[8 Pi/15] kron[K,2]

(Subscript[γ, HG])=2 γ

kron[G_,G1_]:=If[G==G1,1,0]

A[K_,Kp_,Kt_,Kb_,Kbp_,G_,H_]:=Sum[(2 x+1)sixj[{K,Kbp},{Kp,Kb},{Kt,x} ] nj[{K,x,Kb},{G,H,1},{G,H,1}]nj[{x,Kp,Kbp},{H,G,1},{H,G,1}],{x,Abs[K-Kb],K+Kb}]

T[K_,Kp_,Kt_,Kb_,Kbp_,G_,H_]:=3 (-1)^(1+Kb+Kp) B[Kb] sq[H]^2 sq[K] sq[Kp] sq[Kt] sq[Kb] sq[Kbp](Subscript[γ, HG]) χ^2((γ^2)+(δ^2))^(-1) A[K,Kp,Kt,Kb,Kbp,G,H]

S[K_,Kp_,Kb_,H_,G_]:=3 (-1)^(1+Kb+K)  sq[K] sq[Kp]  sq[Kb] (Subscript[γ, HG]) χ^2((γ^2)+(δ^2))^(-1) nj[{K,Kp,Kb},{H,G,1},{H,G,1}]

Sp[K_,Kp_,Kb_,G_,H_]:=-3 (-1)^(H-G+Kb)  sq[K] sq[Kp]  sq[Kb]  χ^2(γ-I δ)^(-1) sixj[{K,G},{Kp,G},{Kb,G}]sixj[{1,G},{1,G},{Kb,H}]

Sm[K_,Kp_,Kb_,G_,H_]:=-3 (-1)^(H-G+Kp+K)  sq[K] sq[Kp]  sq[Kb]  χ^2(γ+I δ)^(-1) sixj[{K,G},{Kp,G},{Kb,G}]sixj[{1,G},{1,G},{Kb,H}]

λ[th_,ph_,Kb_,Kbp_,j_,jp_,Kt_,Qt_]:=(-1)^(Kb+Kbp+Qt)sq[Kt] Sum[T3[{Kb,Qt-Qbp},{Kbp,Qbp},{Kt,-Qt} ] SphericalHarmonicY[Kb,Qt-Qbp,th,ph]ϵ[Kbp,Qbp,j,jp],{Qbp,-Kbp,Kbp}]

interm[K_,Q_,G_,H_,th_,ph_,p_,pp_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] T[K,Kp,Kt,Kb,Kbp,G,H]λ[th,ph,Kb,Kbp,j,jp,Kt,Qp-Q]ρ[G,p+(k-k[j]),pp+k-k[jp],Kp,Qp],{Kt,Abs[K-Kp],K+Kp} ,{Kb,0,2}, {Kbp,0,2},{Kp,0,2 G},{Qp,-Kp,Kp},{j,1,2},{jp,1,2}]

interm[K_,Q_,Kp_,Qp_,G_,H_,th_,ph_,p_,pp_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] T[K,Kp,Kt,Kb,Kbp,G,H]λ[th,ph,Kb,Kbp,j,jp,Kt,Qp-Q]ρ[G,p+(k-k[j]),pp+k-k[jp],Kp,Qp],{Kt,Abs[K-Kp],K+Kp} ,{Kb,0,2}, {Kbp,0,2},{j,1,2},{jp,1,2}]

outterm[K_,Q_,Kp_,Qp_,G_,H_,p_,pp_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] ϵ[Kt,Qp-Q,j,jp](Sp[K,Kp,Kt,G,H]ρ[G,p,pp+k[j]-k[jp],Kp,Qp]+Sm[K,Kp,Kt,G,H]ρ[G,p-k[j]+k[jp],pp,Kp,Qp]),{Kt,Abs[K-Kp],K+Kp} ,{j,1,2},{jp,1,2}]

intermavg[K_,Q_,Kp_,Qp_,G_,H_,p_,pp_]:=(3/(8 Pi))Integrate[Sin[th]interm[K,Q,Kp,Qp,G,H,th,ph,p,pp],{th,0,Pi},{ph,0,2 Pi}]

infp[K_,Q_,Kp_,Qp_,G_,H_,th_,ph_,z_,p_,hb_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] T[K,Kp,Kt,Kb,Kbp,G,H]λ[th,ph,Kb,Kbp,j,jp,Kt,Qp-Q]ρ[G,z,p+hb (k-(k[jp]+k[j])/2),Kp,Qp] Exp[I(k[j]-k[jp])z],{Kt,Abs[K-Kp],K+Kp} ,{Kb,0,2}, {Kbp,0,2},{j,1,2},{jp,1,2}]

infpav[K_,Q_,Kp_,Qp_,G_,H_,z_,p_,hb_]:=(3/(8 Pi))Integrate[Sin[th]infp[K,Q,Kp,Qp,G,H,th,ph,z,p,hb],{th,0,Pi},{ph,0,2 Pi}]

outfp[K_,Q_,Kp_,Qp_,G_,H_,z_,p_,hb_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] ϵ[Kt,Qp-Q,j,jp](Sp[K,Kp,Kt,G,H]ρ[G,z,p+hb(k[j]-k[jp])/2,Kp,Qp]Exp[I(k[j]-k[jp])z]+Sm[K,Kp,Kt,G,H]ρ[G,z,p-hb(k[j]-k[jp])/2,Kp,Qp]Exp[I(k[j]-k[jp])z]),{Kt,Abs[K-Kp],K+Kp} ,{j,1,2},{jp,1,2}]

**  sigmaplus, sigma minus polarization ** **

ϵ[K_,Q_,1,1]:=-(kron[K,0] /(Sqrt[3])-kron[K,1]/(Sqrt[2])+kron[K,2] /(Sqrt[6]))kron[Q,0]

ϵ[K_,Q_,2,2]:=-(kron[K,0] /(Sqrt[3])+kron[K,1]/(Sqrt[2])+kron[K,2] /(Sqrt[6]))kron[Q,0]

ϵ[K_,Q_,1,2]:=-kron[K,2] kron[Q,-2]

ϵ[K_,Q_,2,1]:=-kron[K,2] kron[Q,2]

**population hbar =0 terms**

Simplify[Sum[outfp[0,0,Kp,Qp,1,2,z,p,0]+infpav[0,0,Kp,Qp,1,2,z,p,0] ,{Kp,0,2},{Qp,-Kp,Kp}]]

** population hbar  terms needed in Eq. (18.97) **

infphb[K_,Q_,Kp_,Qp_,G_,H_,th_,ph_,z_,p_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] T[K,Kp,Kt,Kb,Kbp,G,H]λ[th,ph,Kb,Kbp,j,jp,Kt,Qp-Q](k0 Cos[th]-((k[j])/2)-((k[jp])/2))Exp[I(k[j]-k[jp])z],{Kt,Abs[K-Kp],K+Kp} ,{Kb,0,2}, {Kbp,0,2},{j,1,2},{jp,1,2}]

infpavhb[K_,Q_,Kp_,Qp_,G_,H_,z_,p_]:=(3/(8 Pi))Integrate[Sin[th]infphb[K,Q,Kp,Qp,G,H,th,ph,z,p],{th,0,Pi},{ph,0,2 Pi}]

outfphb[K_,Q_,Kp_,Qp_,G_,H_,z_,p_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] ϵ[Kt,Qp-Q,j,jp](Sp[K,Kp,Kt,G,H](1/2) (k[j]-k[jp])Exp[I(k[j]-k[jp])z]-Sm[K,Kp,Kt,G,H](1/2) (k[j]-k[jp])Exp[I(k[j]-k[jp])z]),{Kt,Abs[K-Kp],K+Kp} ,{j,1,2},{jp,1,2}]

Simplify[Sum[(outfphb[0,0,Kp,Qp,1,2,z,p,0]+infpavhb[0,0,Kp,Qp,1,2,z,p,0]) prho[Kp,Qp],{Kp,0,2},{Qp,-Kp,Kp}]]

population hbar  squared  terms needed in Eq. (18.97)

infphbsq[K_,Q_,Kp_,Qp_,G_,H_,th_,ph_,z_,p_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] T[K,Kp,Kt,Kb,Kbp,G,H]λ[th,ph,Kb,Kbp,j,jp,Kt,Qp-Q]((1/8) (((2 k0 Cos[th]-k[j]-k[jp]))^2))Exp[I(k[j]-k[jp])z],{Kt,Abs[K-Kp],K+Kp} ,{Kb,0,2}, {Kbp,0,2},{j,1,2},{jp,1,2}]

infpavhbsq[K_,Q_,Kp_,Qp_,G_,H_,z_,p_]:=(3/(8 Pi))Integrate[Sin[th]infphbsq[K,Q,Kp,Qp,G,H,th,ph,z,p],{th,0,Pi},{ph,0,2 Pi}]

outfphbsq[K_,Q_,Kp_,Qp_,G_,H_,z_,p_]:=Sum[(-1)^Qp T3[{K,Q},{Kp,-Qp},{Kt,Qp-Q} ] ϵ[Kt,Qp-Q,j,jp](Sp[K,Kp,Kt,G,H](1/8) (((k[j]-k[jp]))^2) Exp[I(k[j]-k[jp])z]+Sm[K,Kp,Kt,G,H](1/8) (((k[j]-k[jp]))^2) Exp[I(k[j]-k[jp])z]),{Kt,Abs[K-Kp],K+Kp} ,{j,1,2},{jp,1,2}]

Simplify[outfphbsq[0,0,1,0,1,2,z,p]]

Simplify[Sum[(outfphbsq[0,0,Kp,Qp,1,2,z,p]+infpavhbsq[0,0,Kp,Qp,1,2,z,p])rhop2[Kp,Qp],{Kp,0,2},{Qp,-Kp,Kp}]]

hbar = 0  other terms for equation (18.98)

Simplify[Sum[outfp[1,0,Kp,Qp,1,2,z,p]+infpav[1,0,Kp,Qp,1,2,z,p] ,{Kp,0,2},{Qp,-Kp,Kp}]]

Simplify[Sum[outfp[2,0,Kp,Qp,1,2,z,p]+infpav[2,0,Kp,Qp,1,2,z,p] ,{Kp,0,2},{Qp,-Kp,Kp}]]

Simplify[Sum[outfp[2,2,Kp,Qp,1,2,z,p]+infpav[2,2,Kp,Qp,1,2,z,p] ,{Kp,0,2},{Qp,-Kp,Kp}]]

hbar  terms  only interested in terms related to rho00

Simplify[Sum[(outfphb[1,0,Kp,Qp,1,2,z,p]+infpavhb[1,0,Kp,Qp,1,2,z,p])rhop[Kp,Qp] ,{Kp,0,2},{Qp,-Kp,Kp}]]

Simplify[Sum[(outfphb[2,0,Kp,Qp,1,2,z,p]+infpavhb[2,0,Kp,Qp,1,2,z,p])rho[Kp,Qp] ,{Kp,0,2},{Qp,-Kp,Kp}]]

Simplify[Sum[(outfphb[2,2,Kp,Qp,1,2,z,p]+infpavhb[2,2,Kp,Qp,1,2,z,p]) rho[Kp,Qp],{Kp,0,2},{Qp,-Kp,Kp}]]

steady state equations for lowest order solution

sseqs={(g+I x) a22 +Sqrt[2/3] g a20/5+I g d a10/(g12 5 Sqrt[2])==g/6,(g-I x) a2m2 +Sqrt[2/3] g a20/5-I g d a10/(g12 5 Sqrt[2])==g/6,0==-a20+(15/11) (1/(3 Sqrt[6])-Sqrt[2/3] (a22+a2m2)/5),0==-a10+I d (a2m2-a22)/(g12 Sqrt[2])}

Solve[sseqs,{a22,a2m2,a10,a20}]

m1={{1/5,0,I (δ/γ)/(5 Sqrt[2]),-I (δ/γ)/(5 Sqrt[2])},{0,11/15,Sqrt[6]/15,Sqrt[6]/15},{I (δ/γ)/(5 Sqrt[2]),Sqrt[6]/15,1+I y,0},{-I (δ/γ)/(5 Sqrt[2]),Sqrt[6]/15,0,1-I y}}

MatrixForm[m1]

lam={0, Sqrt[2/3]/6, 1/6, 1/6}Sqrt[3] ro

MatrixForm[lam]

***(Subscript[D, 0])****

Simplify[Inverse[m1].lam]

Simplify[407 ro +28 Sqrt[2] ((5 ro (5 (3+5 (y^2)) (γ^2)+3 (δ^2)))/((Sqrt[2]) (5 (51+55 (y^2)) (γ^2)+51 (δ^2))))-22 Sqrt[3]  ((15 (Sqrt[3]) ro ((5) (γ^2)+(δ^2)))/(10 (51+55 (y^2)) (γ^2)+102 (δ^2)))]/450

((3 ro 766(5 (1+(855/766) (y^2)) (γ^2)+ (δ^2)))/(50  51(5 (1+55/51 (y^2)) (γ^2)+1 (δ^2))))

Simplify[5 855/766]

dif0[y_]:=((383  (5 (1+((855 (y^2))/766)) (γ^2)+(δ^2)))/(425 (5 (1+((55 (y^2))/51)) (γ^2)+(δ^2))))

***(Subscript[D, a])****

lamt={-(Sqrt[(3/2)]),0,-((I δ)/(5 (Sqrt[3]) γ)),((I δ)/(5 (Sqrt[3]) γ))}rhop

MatrixForm[lamt]

Simplify[MatrixForm[Inverse[m1]]]

Simplify[Inverse[m1].lamt]

Simplify[h k Γc rhop ((((Sqrt[12]) (25 (51+55 (y^2)) (γ^2)+34 (δ^2)))/(10 (51+55 (y^2)) (γ^2)+102 (δ^2)))+I (δ/(5γ)) ((13 (51 I+55 y) γ δ)/(2 (Sqrt[3]) (5 (51+55 (y^2)) (γ^2)+51 (δ^2))))-I (δ/(5γ))((13 (-51 I+55 y) γ δ)/(2 (Sqrt[3]) (5 (51+55 (y^2)) (γ^2)+51 (δ^2)))))]/Sqrt[3]

((h k Γc (125 (1+55 (y^2)/51) (γ^2)- (δ^2)))/(5 (5 (1+55 (y^2)/51) (γ^2)+ (δ^2))))

difa[y_]:=(( (125 (1+55 (y^2)/51) (γ^2)- (δ^2)))/(5 (5 (1+55 (y^2)/51) (γ^2)+ (δ^2))))

***(Subscript[D, b])****

m2=-{{0,Sqrt[3]/5,-1/(5 Sqrt[2]),-1/(5 Sqrt[2])},{7/(5 Sqrt[3]),0,-I (δ/γ)Sqrt[2/3]/(10),I (δ/γ)Sqrt[2/3]/(10)},{1/(5 Sqrt[2]),I Sqrt[6](δ/γ)/30,0,0},{1/(5 Sqrt[2]),-I Sqrt[6](δ/γ)/30,0,0}}

MatrixForm[m2]

Simplify[Inverse[m1].m2.Inverse[m1].lam]/ro

Simplify[ (((5 (Sqrt[(3/1)]) (6875 (y^4) (γ^4)+25 (y^2) (255 (γ^4)+2 (γ^2) (δ^2))+51 (5 (γ^2) (δ^2)+(δ^4))))/(((5 (51+55 (y^2)) (γ^2)+51 (δ^2)))^2))+I (δ/(5γ))((5 γ δ (5 (-153 I-390 y+735 I (y^2)+550 (y^3)) (γ^2)-3 (51 I+130 y) (δ^2)))/(2 (Sqrt[3]) (((5 (51+55 (y^2)) (γ^2)+51 (δ^2)))^2)))-I (δ/(5γ))((5 γ δ (5 (153 I-390 y-735 I (y^2)+550 (y^3)) (γ^2)+3 (51 I-130 y) (δ^2)))/(2 (Sqrt[3]) (((5 (51+55 (y^2)) (γ^2)+51 (δ^2)))^2))))]/Sqrt[3]

difb[y_]:=((34375 (y^4) (γ^4)+75 (y^2) (425 (γ^4)-13 (γ^2) (δ^2))+306 (5 (γ^2) (δ^2)+(δ^4)))/(((5 (51+55 (y^2)) (γ^2)+51 (δ^2)))^2))

dif[y_]:=dif0[y]+difa[y]+difb[y]

Simplify[dif[y]]

** *extra term from problem 10 ***

Simplify[Inverse[m1].lam]

ggg[y_]:={-((75 (Sqrt[(3/2)]) ro y γ δ)/(5 (51+55 (y^2)) (γ^2)+51 (δ^2))),((5 ro (5 (3+5 (y^2)) (γ^2)+3 (δ^2)))/((Sqrt[2]) (5 (51+55 (y^2)) (γ^2)+51 (δ^2)))),((15 (Sqrt[3]) ro ((5-5 I y) (γ^2)+(δ^2)))/(10 (51+55 (y^2)) (γ^2)+102 (δ^2))),((15 (Sqrt[3]) ro ((5+5 I y) (γ^2)+(δ^2)))/(10 (51+55 (y^2)) (γ^2)+102 (δ^2)))}

(∂_y) ggg[y]

Simplify[Inverse[m1].lam]

Simplify[Inverse[m1]]

invm1[y_]:={{((25 (51+55 (y^2)) (γ^2))/(5 (51+55 (y^2)) (γ^2)+51 (δ^2))),((50 (Sqrt[3]) y γ δ)/(5 (51+55 (y^2)) (γ^2)+51 (δ^2))),-((5 (51 I+55 y) γ δ)/((Sqrt[2]) (5 (51+55 (y^2)) (γ^2)+51 (δ^2)))),-((5 (-51 I+55 y) γ δ)/((Sqrt[2]) (5 (51+55 (y^2)) (γ^2)+51 (δ^2))))},{((50 (Sqrt[3]) y γ δ)/(5 (51+55 (y^2)) (γ^2)+51 (δ^2))),((75 (5 (1+(y^2)) (γ^2)+(δ^2)))/(5 (51+55 (y^2)) (γ^2)+51 (δ^2))),((5 (Sqrt[6]) (5 I (I+y) (γ^2)-(δ^2)))/(5 (51+55 (y^2)) (γ^2)+51 (δ^2))),((5 (Sqrt[6]) ((-5-5 I y) (γ^2)-(δ^2)))/(5 (51+55 (y^2)) (γ^2)+51 (δ^2)))},{-((5 (51 I+55 y) γ δ)/((Sqrt[2]) (5 (51+55 (y^2)) (γ^2)+51 (δ^2)))),((5 (Sqrt[6]) (5 I (I+y) (γ^2)-(δ^2)))/(5 (51+55 (y^2)) (γ^2)+51 (δ^2))),((10 (53-55 I y) (γ^2)+55 (δ^2))/(10 (51+55 (y^2)) (γ^2)+102 (δ^2))),((20 (γ^2)+55 (δ^2))/(10 (51+55 (y^2)) (γ^2)+102 (δ^2)))},{-((5 (-51 I+55 y) γ δ)/((Sqrt[2]) (5 (51+55 (y^2)) (γ^2)+51 (δ^2)))),((5 (Sqrt[6]) ((-5-5 I y) (γ^2)-(δ^2)))/(5 (51+55 (y^2)) (γ^2)+51 (δ^2))),((20 (γ^2)+55 (δ^2))/(10 (51+55 (y^2)) (γ^2)+102 (δ^2))),((10 (53+55 I y) (γ^2)+55 (δ^2))/(10 (51+55 (y^2)) (γ^2)+102 (δ^2)))}}

Simplify[invm1[y].m2.(∂_y) invm1[y].lam]

Simplify[ (((375 (Sqrt[(3/1)]) ro y (γ^2) (850 (51+55 (y^2)) (γ^4)+165 (51+65 (y^2)) (γ^2) (δ^2)-51 (δ^4)))/(((5 (51+55 (y^2)) (γ^2)+51 (δ^2)))^3))+I (δ/(5γ)) ((375 (Sqrt[3]) ro γ δ (25 (-51+357 I y+385 (y^2)-55 I (y^3)) (γ^4)+5 (-102+357 I y+385 (y^2)) (γ^2) (δ^2)-51 (δ^4)))/(2 (((5 (51+55 (y^2)) (γ^2)+51 (δ^2)))^3)))-I (δ/(5γ))((375 (Sqrt[3]) ro γ δ (25 (-51-357 I y+385 (y^2)+55 I (y^3)) (γ^4)+5 (-102-357 I y+385 (y^2)) (γ^2) (δ^2)-51 (δ^4)))/(2 (((5 (51+55 (y^2)) (γ^2)+51 (δ^2)))^3))))]/Sqrt[3]

c1=275/51

((750 ro y (γ^2) (425 (51+55 (y^2)) (γ^4)+5 (663+1100 (y^2)) (γ^2) (δ^2)-204 (δ^4)))/(51^3(((5 (1+55 (y^2)/51) (γ^2)+ (δ^2)))^3)))

extra[y_]:=((250 ro y (γ^2) (425 (51+55 (y^2)) (γ^4)+5 (663+1100 (y^2)) (γ^2) (δ^2)-204 (δ^4)))/(44217 (((5 (1+((55 (y^2))/51)) (γ^2)+(δ^2)))^3)))

g[d_]:=1/(5+4d^2)

extraT[z_,d_]:=((250 g[d]^2 (21675 (1+z/(5 g[d])) +20 (663+1100 z/(c1 g[d])) (d^2)-3264 (d^4)))/(c1 44217 (((1+z))^3)))

Energy

sigmaplus - sigmaminus  integral needed in Eq. (18.134)

Integrate[a[d] (1+x)/((1+x)^2 (b+e)+(f[d] +c x )(1+x)+(g1[d]+h[d] x+k x^2))/2,{x,0, z^2}]

g[d_]:=1/(5+4d^2)

a[d_]:=-120 g[d]/17

b=383/425

c=((4275/766)/(275/51) -1)383/425

e=5

f[d_]:=-104 g[d] d^2/5

g1[d_]:=1224 g[d] d^2/2601

g1[d]

h[d_]:= g[d] (425-52 d^2)/187

k=13005/(11 2601)

Simplify[g1[d]]

inn[d_,z_]:=-((4982-45464 (d^2)+49600 (d^4)) ArcTanh[((131897+364108 (d^2)+124000 (d^4))/(5 (Sqrt[168646201+1440450712 (d^2)+5440910288 (d^4)+3611951360 (d^6)+615040000 (d^8)])))]+(Sqrt[168646201+1440450712 (d^2)+5440910288 (d^4)+3611951360 (d^6)+615040000 (d^8)]) Log[264 (209+116 (d^2))])/(119442 (1+4 (d^2)) (Sqrt[168646201+1440450712 (d^2)+5440910288 (d^4)+3611951360 (d^6)+615040000 (d^8)]))+((4982-45464 (d^2)+49600 (d^4)) ArcTanh[((131897+124000 (d^4)+119442 (z^2)+4 (d^2) (91027+119442 (z^2)))/(5 (Sqrt[168646201+1440450712 (d^2)+5440910288 (d^4)+3611951360 (d^6)+615040000 (d^8)])))]+(Sqrt[168646201+1440450712 (d^2)+5440910288 (d^4)+3611951360 (d^6)+615040000 (d^8)]) Log[55176+131897 (z^2)+124000 (d^4) (z^2)+59721 (z^4)+4 (d^2) (7656+91027 (z^2)+59721 (z^4))])/(119442 (1+4 (d^2)) (Sqrt[168646201+1440450712 (d^2)+5440910288 (d^4)+3611951360 (d^6)+615040000 (d^8)]))

***energy distribution given in Eq. (18.133)***

temp[a_,d_]:=NIntegrate[z^2Exp[a inn[d,z]/2],{z,-Infinity,Infinity}]/NIntegrate[ Exp[a inn[d,z]/2],{z,-Infinity,Infinity}]

** does not seem to work for arbitrary detuning - calculate for specific detunngs of 1 and 5 and in problem 11 ***

coef[d_]:=4 c1 g[d]/d^2

diff[d_,x_]:=((1+x)^2 (b+e)+(f[d] +c x )(1+x)+(g1[d]+h[d] x+k x^2))

ingr1[d_,x_]:=a[d] (1+x)/((1+x)^2 (b+e)+(f[d] +c x )(1+x)+(g1[d]+h[d] x+k x^2))/2

****fgh is 1/2 integral in Eq. (18.134)***

fgh[d_,z_]:=Integrate[ingr1[1.00001 d,x],{x,0,z^2},Assumptions->z>0&&d>0]

fgh[1,z]

fgh[5,z]

**%127 is fgh[1,z] and %128 is fgh[5,z]***

***Now evaluate Eq. (18.133) for energy distribution for δ=1 and 5***

Plot[ a^2 coef[1]NIntegrate[ z^2 Exp[a %127],{z,0,Infinity}]/NIntegrate[ Exp[a%127],{z,0,Infinity}],{a,25,125},PlotRange->{0,2000},AxesLabel->{"(Subscript[a, c])","<E>/(Subscript[E, R]), δ=1"}]

Plot[ a^2 coef[5]NIntegrate[ z^2 Exp[a %128],{z,0,Infinity}]/NIntegrate[ Exp[a %128],{z,0,Infinity}],{a,300,1750},PlotRange->{0,400},AxesLabel->{"(Subscript[a, c])","<E>/(Subscript[E, R]), δ=5"}]

***velocity distribution****

NIntegrate[ z^2 Exp[800 %128],{z,0,Infinity}]/NIntegrate[ Exp[800 %128],{z,0,Infinity}]

Plot[{Exp[-z^2/(2 .0293)],Exp[800 %128]},{z,-1,1},AxesLabel->{x,"Velocity distributions- Actual distribution narrower than Gaussian "}]

Export["prob18s13.jpg",%]

NIntegrate[ z^2 Exp[400 %128],{z,0,Infinity}]/NIntegrate[ Exp[400 %128],{z,0,Infinity}]

Plot[{Exp[-z^2/(.25)],Exp[400 %128]},{z,-1.5,1.5},AxesLabel->{x,"Velocity distributions - Actual distribution narrower than Gaussian"}]

***contribution from extra term of problem 10  **

((250 g[d]^2 (21675 (1+z/(5 g[d])) +20 (663+1100 z/(c1 g[d])) (d^2)-3264 (d^4)))/(c1 44217 (((1+z))^3)))

extra[d_,z_]:=((250 g[d]^2 (21675 (1+z/(5 g[d])) +20 (663+1100 z/(c1 g[d])) (d^2)-3264 (d^4)))/(c1 44217 (((1+z))^3)))

coef[d_]:=4 c1 g[d]/d^2

ingrand[d_,x_]:=a[d] (1+x)/((1+x)^2 (b+e)+(f[d] +c x )(1+x)+(g1[d]+h[d] x+k x^2))/2

fgh[d_,z_]:=Integrate[ingrand[1.00001 d,x],{x,0,z^2},Assumptions->z>0&&d>0]

ingr1[d_,x_]:=a[d] (1+x)/((1+x)^2 (b+e)+(f[d] +c x )(1+x)+(g1[d]+h[d] x+k x^2))/2

fgh1[d_,z_]:=Integrate[ingr1[1.00001 d,x],{x,0,z^2},Assumptions->z>0&&d>0]

fgh2[d_,z_]:=Integrate[extra[1.0000001 d,x]/((1+x)^2 (b+e)+(f[d] +c x )(1+x)+(g1[d]+h[d] x+k x^2))/2,{x,0,z^2},Assumptions->z>0&&d>0]

fgh1[5,z]

fgh2[5,z]

fgh1[1,z]

fgh2[1,z]

Series[%,{z,0,2}]

Plot[a^2 coef[5]NIntegrate[ z^2 Exp[a %167],{z,0,Infinity}]/NIntegrate[ Exp[a %167],{z,0,Infinity}],{a,400,2000}]

Plot[a^2 coef[5]NIntegrate[ z^2 Exp[a %167],{z,0,Infinity}]/NIntegrate[ Exp[a %167],{z,0,Infinity}],{a,400,2000}]

Plot[a^2 coef[5]NIntegrate[ z^2 Exp[a %167 -%169],{z,0,Infinity}]/NIntegrate[ Exp[a %167-%169],{z,0,Infinity}],{a,400,2000},MaxRecursion->1000000]

Show[%,%173]

****results identical***

fgh1[1,z]

Series[%634,{z,0,2}]

fgh2[1,z]

Series[%635,{z,0,2}]

Plot[a^2 4 5.4/(1 9)NIntegrate[ z^2 Exp[a %634],{z,0,Infinity}]/NIntegrate[ Exp[a %634],{z,0,Infinity}],{a,25,100}]

Plot[a^2 4 5.4/(1 9)NIntegrate[ z^2 Exp[a %634 -%635],{z,0,Infinity}]/NIntegrate[ Exp[a %634-%635],{z,0,Infinity}],{a,25,100}]

Show[%,%641]

***results similar***