Faraday conversion & rotation in uniformly magnetized relativistic plasmas (Huang & Shcherbakov, 2011, MNRAS submitted)

Definitions

Note, that the numerical factor ((e^2)/((Subscript[m, e])c))is consistently omitted and density of particles is set to unity: (Subscript[n, e])=1

Rfin=(ω^2) ((ξ^2)-(((-2+(ξ^2) ((Ω0)^2)+2 Cos[ξ Ω0]) ((Sin[θ])^2))/((Ω0)^2)));
(tdot={{-1,0,0,0},{0,Cos[Ω0 τ],-Cos[θ]Sin[Ω0 τ],-Sin[θ]Sin[Ω0 τ]},{0,Cos[θ]Sin[Ω0 τ],((Sin[θ])^2)+((Cos[θ])^2)Cos[Ω0 τ],-Sin[θ]Cos[θ]+Sin[θ]Cos[θ]Cos[Ω0 τ]},{0,Sin[θ]Sin[Ω0 τ],-Sin[θ]Cos[θ]+Sin[θ]Cos[θ]Cos[Ω0 τ],((Cos[θ])^2)+((Sin[θ])^2)Cos[Ω0 τ]}});
Π1=((Sin[R p0])/(R^3))-p0 ((Cos[R p0])/(R^2));Π2=3((Sin[R p0])/(R^5))-p0 ((3Cos[R p0])/(R^4))-(p0^2)((Sin[R p0])/(R^3));Π3=γ ((Sin[R p0])/R);Π4=γ ((Sin[R p0])/(R^3))-γ p0 ((Cos[R p0])/(R^2));
XX=(Sin[Ω0 ξ]-Ω0 ξ);YY=(1-Cos[Ω0 ξ]);
resLei(*(α^ij) in equation 12*)=4π Exp[I γ ξ ω]((((tdot/.τ->ξ)(I ω ξ Π1+Π3))[[2;;3,2;;3]])-(I ω ξ Π2+Π4)(((ω^2)((Sin[θ])^2))/((Ω0)^2)){{-(YY^2),-XX YY Cos[θ]},{XX YY Cos[θ],(XX^2)((Cos[θ])^2)}});
AALei(*(α^ij) in equation 13*)=4π Exp[I γ ξ ω]((((tdot/.τ->ξ)Π1)[[2;;3,2;;3]])-Π2(((ω^2)((Sin[θ])^2))/((Ω0)^2)){{-(YY^2),-XX YY Cos[θ]},{XX YY Cos[θ],(XX^2)((Cos[θ])^2)}});
Rsub={p0->(Sqrt[(γ^2)-1]),R->(Sqrt[Rfin])};

Sample integration for a monoenergetic distribution

sub={γ->500.,ω->1,Ω0->(10^(-4.)),θ->N[45°]};num=25;
xsub=Join[sub,Rsub/.sub];ttx=TimeUsed[];
rec=8;accF=0;accI=8;acc=0;ang=(π/4.);
SetOptions[NIntegrate,Method->{"GlobalAdaptive",Method->{"LobattoKronrodRule","Points"->5}},MinRecursion->rec,MaxRecursion->rec+5,MaxPoints->+∞,WorkingPrecision->MachinePrecision,AccuracyGoal->6];
ttx=TimeUsed[];
expr=Re[resLei[[1,2]](Cos[α]+I Sin[α])/.xsub/.ξ->ks (Cos[α]+I Sin[α])/.α->ang];
i=1;res={};ma=0;While[i<(10^6),ne=Abs[expr/.ks->i];If[ne>ma,ma=ne];If[(ne/ma)<(10^(-100.)),Break[]];i*=1.1];
Print["Outer = ",Round[ksinf=i]];kin=Log[ksinf+(Sqrt[1+(ksinf^2)])];
ans1=NIntegrate[Re[resLei[[1,2]]Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];(*Faraday rotation*)
ans2=NIntegrate[Im[(resLei[[1,1]]-resLei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];(*Faraday conversion*)
ans3=NIntegrate[Im[resLei[[1,2]]Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->I num+Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];(*V emission*)
ans4=NIntegrate[Re[(resLei[[1,1]]-resLei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->I num+Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];
ans5=NIntegrate[Re[(resLei[[1,1]]+resLei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->I num+Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];
ang=0.5;
ans1x=NIntegrate[Re[resLei[[1,2]]Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];(*Faraday rotation*)
ans2x=NIntegrate[Im[(resLei[[1,1]]-resLei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];(*Faraday conversion*)
ans3x=NIntegrate[Im[resLei[[1,2]]Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->I num+Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];(*V emission*)
ans4x=NIntegrate[Re[(resLei[[1,1]]-resLei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->I num+Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];
ans5x=NIntegrate[Re[(resLei[[1,1]]+resLei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->I num+Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];
If[Abs[(cur1=(ans1/ans1x)-1)]<(10^(-2.)),Print[ans1,";  ρV relative error is ",cur1],Print[ans1,"; ρV diverged ", cur1]]
If[Abs[(cur2=(ans2/ans2x)-1)]<(10^(-2.)),Print[ans2,";  ρQ relative error is ",cur2],Print[ans2,"; ρQ diverged ", cur2]]
If[Abs[(cur3=(ans3/ans3x)-1)]<(10^(-2.)),Print[ans3,";  ηV relative error is ",cur3],Print[ans3,"; ηV diverged ", cur3]]
If[Abs[(cur4=(ans4/ans4x)-1)]<(10^(-2.)),Print[ans4,";  ηQ relative error is ",cur4],Print[ans4,"; ηQ diverged ", cur4]]
If[Abs[(cur5=(ans5/ans5x)-1)]<(10^(-2.)),Print[ans5,";  ηI relative error is ",cur5],Print[ans5,"; ηI diverged ", cur5]]
Print[TimeUsed[]-ttx," seconds"]

Accuracy of fitting formulae

sub={γ->40.,ω->1,Ω0->(10^(-4.)),θ->N[45°]};num=25;ruru={me->9.1 (10^(-28)),ee->4.8 (10^(-10)),c->3 (10^10)};XA=(Sqrt[(((Sqrt[2])Sin[θ]Ω0)/((10^(-4))ω))]);ttx=TimeUsed[];
rec=8;accF=0;accI=8;acc=0;ang=(π/4.);
SetOptions[NIntegrate,Method->{"GlobalAdaptive",Method->{"LobattoKronrodRule","Points"->5}},MinRecursion->rec,MaxRecursion->rec+5,MaxPoints->+∞,WorkingPrecision->MachinePrecision,AccuracyGoal->6];
ttx=TimeUsed[];

ρQacc(*accurate Faraday conversion*):=(xsub=Join[sub,Rsub/.sub];expr=Re[resLei[[1,2]](Cos[α]+I Sin[α])/.xsub/.ξ->ks (Cos[α]+I Sin[α])/.α->ang];
i=1;res={};ma=0;While[i<(10^6),ne=Abs[expr/.ks->i];If[ne>ma,ma=ne];If[(ne/ma)<(10^(-100.)),Break[]];i*=1.1];ksinf=i;kin=Log[ksinf+(Sqrt[1+(ksinf^2)])];
ans2=-2π NIntegrate[Re[(resLei[[1,1]]-resLei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];ang=0.5;
ans2x=-2π NIntegrate[Re[(resLei[[1,1]]-resLei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];
If[Abs[(cur2=(ans2/ans2x)-1)]<(10^(-2.)),conv=True ;ans2,conv=False;0]);
ρVacc(*accurate Faraday rotation*):=(xsub=Join[sub,Rsub/.sub];expr=Re[resLei[[1,2]](Cos[α]+I Sin[α])/.xsub/.ξ->ks (Cos[α]+I Sin[α])/.α->ang];
i=1;res={};ma=0;While[i<(10^6),ne=Abs[expr/.ks->i];If[ne>ma,ma=ne];If[(ne/ma)<(10^(-100.)),Break[]];i*=1.1];ksinf=i;kin=Log[ksinf+(Sqrt[1+(ksinf^2)])];
ans2=4π NIntegrate[Im[resLei[[1,2]]Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];ang=0.5;
ans2x=4π NIntegrate[Im[resLei[[1,2]]Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];
If[Abs[(cur2=(ans2/ans2x)-1)]<(10^(-2.)),conv=True ;ans2,conv=False;0]);
ρQBacc(*boundary term for accurate Faraday conversion*):=(xsub=Join[sub,Rsub/.sub];expr=Re[resLei[[1,2]](Cos[α]+I Sin[α])/.xsub/.ξ->ks (Cos[α]+I Sin[α])/.α->ang];
i=1;res={};ma=0;While[i<(10^6),ne=Abs[expr/.ks->i];If[ne>ma,ma=ne];If[(ne/ma)<(10^(-100.)),Break[]];i*=1.1];ksinf=i;kin=Log[ksinf+(Sqrt[1+(ksinf^2)])];
ans2=-2π NIntegrate[Re[(AALei[[1,1]]-AALei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];ang=0.5;
ans2x=-2π NIntegrate[Re[(AALei[[1,1]]-AALei[[2,2]])Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];
If[Abs[(cur2=(ans2/ans2x)-1)]<(10^(-2.)),conv=True ;ans2,conv=False;0]);
ρVBacc(*boundary term for accurate Faraday rotation*):=(xsub=Join[sub,Rsub/.sub];expr=Re[resLei[[1,2]](Cos[α]+I Sin[α])/.xsub/.ξ->ks (Cos[α]+I Sin[α])/.α->ang];
i=1;res={};ma=0;While[i<(10^6),ne=Abs[expr/.ks->i];If[ne>ma,ma=ne];If[(ne/ma)<(10^(-100.)),Break[]];i*=1.1];ksinf=i;kin=Log[ksinf+(Sqrt[1+(ksinf^2)])];
ans2=4π NIntegrate[Im[AALei[[1,2]]Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];ang=0.5;
ans2x=4π NIntegrate[Im[AALei[[1,2]]Cosh[kex](Cos[α]+I Sin[α])/.xsub/.ξ->Sinh[kex] (Cos[α]+I Sin[α])/.α->ang],{kex,0,kin}];
If[Abs[(cur2=(ans2/ans2x)-1)]<(10^(-2.)),conv=True ;ans2,conv=False;0]);

γmin=1.001;γmax=(10^4);nn=120;γtab=Table[(*z==Log[γ-1]*)γ=Exp[z]+1.,{z,Log[γmin-1],Log[γmax-1],Log[((γmax-1)/(γmin-1))]/(nn-1)}];
Clear[γ,γx];ga=Log[10,XA γ];ga=XA γ;sub={γ->γx,ω->1,Ω0->(10^(-4.)),θ->N[45°]};
HX[x_](*equation 52*):=Piecewise[{{9.29 (10^(-9))(Sqrt[1-(γ^(-1))])(x^3.036),x<40},{-0.000203 (x^0.4343)-0.0013 Cos[0.5646Log[x]-4.03]+0.002 Exp[-((((Log[x]-4.2137))^2)/0.5429)]+0.00083Exp[-((((Log[x]-4.2137))^2)/0.2121)],x≥40}}];fX[x_]:=1-0.4 Exp[-((((Log[x]-9.21))^2)/11.93)]-0.05Exp[-((((Log[x]-5.76))^2)/1.33)]+0.075Exp[-((((Log[x]-4.03))^2)/0.65)];
hB[x_]:=Piecewise[{{4.67 (10^(-9))(((1-(γ^(-1))))^(3/2))(x^3.84),x<40},{0.864-0.2082((Log[x])^2)+0.0175 ((Log[x])^4)-0.000626 ((Log[x])^6)+1.0175 (10^(-5))((Log[x])^8)-7.686 (10^(-8))((Log[x])^10)
-0.01Exp[-((((Log[x]-4.0755))^2)/0.0763)],x≥40}}];fB[x_]:=1-0.0045(x^0.52);
ρQapp:=8(π^2)(*((ee^2)/(me cc))*)XA HX[ga/.sub]/.sub;ρVapp:=8(π^2) Ω0 Cos[θ] Log[((1+(Sqrt[1-(γ^(-2))]))/(1-(Sqrt[1-(γ^(-2))])))]fX[ga/.sub]/.sub;
ρQBapp:=8(π^2) hB[ga/.sub]/.sub;
ρVBapp:=8(π^2) Ω0 Cos[θ] (γ Log[((1+(Sqrt[1-(γ^(-2))]))/(1-(Sqrt[1-(γ^(-2))])))]-2(Sqrt[(γ^2)-1]))fB[ga/.sub]/.sub;
Dynamic[k]
Quiet[ρQapptab=Table[γx=γtab[[k]];{γx-1,ρQapp},{k,nn}];
ρVapptab=Table[γx=γtab[[k]];{γx-1,ρVapp},{k,nn}];
ρQBapptab=Table[γx=γtab[[k]];{γx-1,ρQBapp},{k,nn}];
ρVBapptab=Table[γx=γtab[[k]];{γx-1,ρVBapp},{k,nn}];

ρQacctab=Table[γx=γtab[[k]];{γx-1,ρQacc},{k,nn}];
ρVacctab=Table[γx=γtab[[k]];{γx-1,ρVacc},{k,nn}];
ρQBacctab=Table[γx=γtab[[k]];{γx-1,ρQBacc},{k,nn}];
ρVBacctab=Table[γx=γtab[[k]];{γx-1,ρVBacc},{k,nn}];]

{lpl=ListLogLogPlot[{Abs[ρQapptab],Abs[ρQacctab]},Joined->True,Frame->True,PlotStyle->{Red,Blue},ImageSize->300],
ListLogLogPlot[Transpose[{γtab-1,Abs[((ρQapptab[[All,2]])/(ρQacctab[[All,2]]))]}],Joined->True,Frame->True,PlotStyle->{Red,Blue},ImageSize->300]}
{lpl=ListLogLogPlot[{Abs[ρVapptab],Abs[ρVacctab]},Joined->True,Frame->True,PlotStyle->{Red,Blue},ImageSize->300],
ListLogLogPlot[Transpose[{γtab-1,Abs[((ρVapptab[[All,2]])/(ρVacctab[[All,2]]))]}],Joined->True,Frame->True,PlotStyle->{Red,Blue},ImageSize->300]}
{lpl=ListLogLogPlot[{Abs[ρQBapptab],Abs[ρQBacctab]},Joined->True,Frame->True,PlotStyle->{Red,Blue},ImageSize->300],
ListLogLinearPlot[Transpose[{γtab-1,Abs[((ρQBapptab[[All,2]])/(ρQBacctab[[All,2]]))]}],Joined->True,Frame->True,PlotStyle->{Red,Blue},ImageSize->300]}
{lpl=ListLogLogPlot[{Abs[ρVBapptab],Abs[ρVBacctab]},Joined->True,Frame->True,PlotStyle->{Red,Blue},ImageSize->300],
ListLogLogPlot[Transpose[{γtab-1,Abs[((ρVBapptab[[All,2]])/(ρVBacctab[[All,2]]))]}],Joined->True,Frame->True,PlotStyle->{Red,Blue},ImageSize->300]}
(*left column - quantity, right column - ratio of approximate to exact - should be 1,
first row - Faraday conversion;
second row - Faraday rotation; 
third row - Faraday conversion boundary term;
fourth row - Faraday rotation boundary term;*)

Test integration over thermal distribution

Clear[γ];θe=6;
ρQinter=Interpolation[ρQapptab,InterpolationOrder->2][γ-1];
ρQexinter=Interpolation[ρQacctab,InterpolationOrder->2][γ-1];

NIntegrate[ρQinter(1/(4π θe BesselK[2,((θe)^(-1))]))Exp[-(γ/(θe))],{γ,γmin,γmax}]
NIntegrate[ρQexinter(1/(4π θe BesselK[2,((θe)^(-1))]))Exp[-(γ/(θe))],{γ,γmin,γmax}]

NIntegrate[(1/(4π θe BesselK[2,((θe)^(-1))]))Exp[-(γ/(θe))]4π γ(Sqrt[(γ^2)-1]),{γ,1,+∞}](*normalization of thermal distribution*)

Clear[γ];ρVinter=Interpolation[ρVapptab,InterpolationOrder->2][γ-1];
θe=6;
NIntegrate[ρVinter(1/(4π θe BesselK[2,((θe)^(-1))]))Exp[-(γ/(θe))],{γ,γmin,γmax}]

XR=θe(Sqrt[(Sqrt[2])Sin[θ]((10^3)((Ω0)/ω))])/.sub;fX=2.011 Exp[-((XR^1.035)/4.7)]-Cos[(XR/2)]Exp[-((XR^1.2)/2.73)]-0.011Exp[-(XR/47.2)];gX=1-0.11Log[1+0.035XR];
4π(*((e^2)/((Subscript[m, e])c))*)((Ω0)/(ω^2))gX(((BesselK[0,((θe)^(-1))])/(BesselK[2,((θe)^(-1))])))Cos[θ]/.sub
4π(*((e^2)/((Subscript[m, e])c))*)(1/2) (((Ω0)^2)/(ω^3))fX(((BesselK[1,((θe)^(-1))])/(BesselK[2,((θe)^(-1))]))+6θe)((Sin[θ])^2)/.sub