Settings

Plotting options

SetOptions[ListPlot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
SetOptions[Histogram, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->False];
SetOptions[Plot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
SetOptions[ListPlot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
SetOptions[LogLogPlot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
SetOptions[LogPlot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
SetOptions[ListLogPlot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
SetOptions[ListLogLogPlot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
SetOptions[ListLogLinearPlot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
SetOptions[LogLinearPlot, BaseStyle -> {FontFamily -> "Times", FontSize -> 16},Frame->True];
Needs["PlotLegends`"]

Options[legendMaker]=Join[FilterRules[Options[Framed],Except[{ImageSize,FrameStyle,Background,RoundingRadius,ImageMargins}]],{FrameStyle->None,Background->Directive[Opacity[.7],LightGray],RoundingRadius->10,ImageMargins->0,PlotStyle->Automatic,PlotMarkers->None,"LegendLineWidth"->35,"LegendLineAspectRatio"->.3,"LegendMarkerSize"->20,"LegendGridOptions"->{Alignment->Left,Spacings->{.4,.1}}}];


legendMaker::usage="Create a Graphics object with legends given by the list passed as the first argument. The options specify any non-deafult line styles (using PlotStyle -> {...}) or plot markers (using PlotMarkers -> {...}). For more options, inspect Options[legendMaker]";

legendMaker[textLabels_,opts:OptionsPattern[]]:=Module[{f,lineDirectives,markerSymbols,n=Length[textLabels],x},lineDirectives=((PlotStyle/.{opts})/.PlotStyle|Automatic:>Map[ColorData[1],Range[n]])/.None->{None};
markerSymbols=Replace[((PlotMarkers/.{opts})/.Automatic:>(Drop[Normal[ListPlot[Transpose[{Range[3]}],PlotMarkers->Automatic][[1,2]]][[1]],-1]/.Inset[x_,i__]:>x)[[All,-1]])/.{Graphics[gr__],sc_}:>Graphics[gr,ImageSize->("LegendMarkerSize"/.{opts}/.Options[legendMaker,"LegendMarkerSize"]/.{"LegendMarkerSize"->100})],PlotMarkers|None:>Map[Style["",Opacity[0]]&,textLabels]]/.None|{}->Style["",Opacity[0]];
lineDirectives=PadRight[lineDirectives,n,lineDirectives];
markerSymbols=PadRight[markerSymbols,n,markerSymbols];
f=Grid[MapThread[{Graphics[{#1/.None->{},If[#1==={None}||(PlotStyle/.{opts})===None,{},Line[{{-.1,0},{.1,0}}]],Inset[#2,{0,0},Background->None]},AspectRatio->("LegendLineAspectRatio"/.{opts}/.Options[legendMaker,"LegendLineAspectRatio"]/.{"LegendLineAspectRatio"->.2}),ImageSize->("LegendLineWidth"/.{opts}/.Options[legendMaker,"LegendLineWidth"]/.{"LegendLineWidth"->35}),ImagePadding->{{1,1},{0,0}}],Text[#3,FormatType->TraditionalForm]}&,{lineDirectives,markerSymbols,textLabels}],Sequence@Evaluate[("LegendGridOptions"/.{opts}/.Options[legendMaker,"LegendGridOptions"]/.{"LegendGridOptions"->{Alignment->Left,Spacings->{.4,.1}}})]];
Framed[f,FilterRules[{Sequence[opts,Options[legendMaker]]},FilterRules[Options[Framed],Except[ImageSize]]]]];

extractStyles::usage="returns a tuple {\"all line style directives\", \"all plot markers\"} found in the plot, in the order they are drawn. The two sublists need not have the same length if some lines don't use markers ";
extractStyles[plot_]:=Module[{lines,markers,points,extract=First[Normal[plot]]},(*In a plot,the list of lines contains no insets,so I use this to find it:*)lines=Select[Cases[Normal[plot],{___,_Line,___},Infinity],FreeQ[#1,Inset]&];
points=Select[Cases[Normal[plot],{___,_Point,___},Infinity],FreeQ[#1,Inset]&];
(*Most plot markers are inside Inset,except for Point in list plots:*)markers=Select[extract,!FreeQ[#1,Inset]&];
(*The function returns a list of lists:*){(*The first return value is the list of line plot styles:*)Replace[Cases[lines,{c__,Line[__],___}:>Flatten[Directive@@Cases[{c},Except[_Line]]],Infinity],{}->None],(*Second return value:marker symbols*)Replace[Join[Cases[markers,{c__,Inset[s_,pos_,d___],e___}:>If[(*markers "s" can be strings or graphics*)Head[s]===Graphics,(*Append scale factor in case it's needed later;
default 0.01*){s,Last[{.01,d}]/.Scaled[f_]:>First[f]},If[(*For strings,add line color if no color specified via text styles:*)FreeQ[s,CMYKColor|RGBColor|GrayLevel|Hue],Style[s,c],s]],Infinity],(*Filter out Pointsize-legends don't need it:*)Cases[points,{c___,Point[pt__],___}:>{Graphics[{c,Point[{0,0}]}]/.PointSize[_]:>PointSize[1],.01},Infinity]],{}->None]}]

autoLegend::usage="Simplified legending for the plot passed as first argument, with legends given as second argument. Use the option Alignment -> {horizontal, vertical} to place the legend in the PlotRegion in scaled coordinates. For other options, see Options[legendMaker] which are used by autoLegend.";
Options[autoLegend]=Join[{Alignment->{Right,Top},Background->White,AspectRatio->Automatic},FilterRules[Options[legendMaker],Except[Alignment|Background|AspectRatio]]];
autoLegend[plot_Graphics,labels_,opts:OptionsPattern[]]:=Module[{lines,markers,align=OptionValue[Alignment]},{lines,markers}=extractStyles[plot];
Graphics[{Inset[plot,{-1,-1},{Left,Bottom},Scaled[1]],Inset[legendMaker[labels,PlotStyle->lines,PlotMarkers->markers,Sequence@@FilterRules[{opts},FilterRules[Options[legendMaker],Except[Alignment]]]],align,Map[If[NumericQ[#],Center,#]&,align]]},PlotRange->{{-1,1},{-1,1}},AspectRatio->(OptionValue[AspectRatio]/.Automatic:>(AspectRatio/.Options[plot,AspectRatio])/.Automatic:>(AspectRatio/.AbsoluteOptions[plot,AspectRatio]))]]

Parameters used for tabulation (a, k, (Subscript[v, bc]))

(* Scale factors *)
ainit = 1./1011.;
aend    = 1./21.;
Na        = 100;
atab = Table[ainit Exp[(n-1.)Log[aend/ainit]/(Na-1.)], {n, 1, Na}];
logatab = Table[Log[ainit] + (n-1.)Log[aend/ainit]/(Na-1.), {n, 1, Na}];

(* Wavenumbers for the small-scale fluctuations *)
kmin = 10.;
kmax = 2700.;
Nk = 50;
ktab  = Table[kmin Exp[(i-1) Log[kmax/kmin]/(Nk-1.)], {i, 1,Nk}];
Dlogk = Log[kmax/kmin]/(Nk-1.);

(* Relative velocity in km/s *)
Nvbc =  45;
vbcmax = 90.;  
vbctab = Table[(j-1.)*vbcmax/(Nvbc-1.), {j, 1, Nvbc}];

Background cosmology

Cosmological parameters

(*** Independent variables ***)
h = 0.678;
Ωb0 = 0.0456;
Ωc0=0.227;  (* Total dark matter: CDM + massive neutrinos *) 
Tcmb0 = 2.726;
YHe = 0.24;     (* He fraction by mass *) 
Neff = 3.046;
As = 2.196 (10^(-9));
ns = 0.96;
kpivot = 0.05; 
(* As = 2.42 (10^(-9));  (* Hirata *)
ns = 1; *)
(* As = 2.46 (10^(-9));
ns = 0.968;
kpivot = 0.002;  *) (* WMAP *)

(*** Derived constants ***)
ΩR0 = 4.48162687719 (10^(-7))* Tcmb0^4/h^2*(1. + 0.2271073 Neff);
ΩΛ0 = 1.-Ωc0-Ωb0-ΩR0;
fHe = YHe/(1.-YHe)/3.97153; (* nHe/nH, with mHe/mH = 3.97153 *)
H0   = 100. h  ;  (* H0 in km/s/Mpc *)
μ     = 1./(YHe/3.97153  + (1.-YHe)); (* mean molecular weight, for xe << 1 *)   
nh0 = 11.223846333047 (10^(-6))Ωb0 h^2*(1.-YHe); (* H density today in (cm^(-3))*)
kmperMpc = 1./(3.086*^19)  (* km/Mpc *)

(*** Decoupling redshift, using the fitting function of Eisenstein & Hu 1998 ***)
Ωmh2 = (Ωb0 + Ωc0)(h^2);
b1 = 0.313 ((Ωmh2)^(-0.419))(1. + 0.607 ((Ωmh2)^0.674));
b2 = 0.238 ((Ωmh2)^0.223);
zd = 1291. (((Ωmh2)^0.251)/(1. + 0.659 ((Ωmh2)^0.828)))(1. + b1 ((Ωmh2)^b2))

Additional parameters for the small-scale power spectrum

(* Fraction of dark matter in neutrinos today *)
fν = 0.;

Background equations

Tcmb[a_]:=(Tcmb0/a);
nhbar[a_] := nh0 (a^(-3));  (* mean hydrogen density in (cm^(-3)) *)
Hubble[a_] := H0 (Sqrt[(Ωc0+Ωb0) (a^(-3)) +ΩΛ0+ΩR0 (a^(-4))]);  (* in km/s/Mpc *)
(* Importing z, xe, Tm/Tr for Planck cosmological parameters, obtained from HyRec *)

zxeTm = Import["~/Dropbox/21cm Dark Ages/data/xePlanck.dat"];
Nzdata = Dimensions[zxeTm][[1]]
xe = Interpolation[Table[{(((1.+zxeTm[[i, 1]]))^(-1)), zxeTm[[i, 2]]}, {i, 1, Nzdata}]];
Tgas = Interpolation[Table[{(((1.+zxeTm[[i, 1]]))^(-1)), Tcmb[(((1.+zxeTm[[i, 1]]))^(-1))]zxeTm[[i, 3]]}, {i, 1, Nzdata}],InterpolationOrder->3];

cs2 [a_] := (0.008263473/μ) Tgas[a];
 
(* Isothermal sound speed (i.e. without the γ) squared in (km/s)^2 *)

ΓCompton[a_] := ((0.01516667  ((Tcmb[a])^4))/(1.+(1. + fHe)/xe[a]));(* Compton heating coefficient in km/s/Mpc. The numerical coefficient is 8 ar σT/(3 me c)/(kmperMpc) *)

Precompute functions for recombination perturbation

kBoltz = 8.617343 (10^(-5));   (* Boltzmann constant in eV per Kelvin *)
EI = 13.598286071938324;  (* Ionization energy in eV *)
Λ2s = 8.2206; (* 2s-1s 2-photon decay rate in (sec^(-1)) *) 

SahaCoeff = 3.016103031869581 (10^21) ;  (* Used for detailed balance, with Tr in eV *)

(* Exact effective case-B recomination coefficient, obtained from tables in HyRec *)

Alphatab = Import["~/Dropbox/21cm Dark Ages/data/Alpha_inf.dat", "Table"];
TRMIN = 0.004;
TRMAX = 0.4;
NTR = 100;
TMovTRMIN = 0.1;
TMovTRMAX = 1.;
NTM = 40;
logTRtab = Table[Log[TRMIN] + (i-1.)Log[TRMAX/TRMIN]/(NTR-1.), {i, 1, NTR}];
TMovTRtab = Table[TMovTRMIN + (i-1.)*(TMovTRMAX - TMovTRMIN)/(NTM -1.), {i, 1, NTM}];
Alpha2sInterp = Flatten[Table[{{logTRtab[[i]], TMovTRtab[[j]]}, Alphatab[[(i-1)NTM + j, 1]] }, {i, 1, NTR}, {j, 1, NTM}], 1];
Alpha2pInterp = Flatten[Table[{{logTRtab[[i]], TMovTRtab[[j]]}, Alphatab[[(i-1)NTM + j, 2]] }, {i, 1, NTR}, {j, 1, NTM}], 1];
Alpha2slog = Interpolation[Alpha2sInterp,Method->"Spline"];
Alpha2plog = Interpolation[Alpha2pInterp,Method->"Spline"];
αB[Tr_, Tm_] := Alpha2slog[Log[kBoltz Tr], Tm/Tr]+Alpha2plog[Log[kBoltz Tr], Tm/Tr]; 
βB[Tr_] := (SahaCoeff/4.) (((kBoltz Tr))^(3/2))Exp[-((0.25EI)/(kBoltz Tr))]αB[Tr, Tr];

(* Effective 3-level atom model for the recombination perturbations *)
RLyα[a_] :=  1.5109848 (10^(-4)) ((Hubble[a] )/(nhbar[a] (1.-xe[a])));(* Lyman-α escape rate for H in km/s/Mpc and nh in (cm^(-3))*)
CPeebles[a_] :=1.- ((βB[Tcmb[a]])/(βB[Tcmb[a]] + 0.75*RLyα[a] + 0.25 Λ2s)); 
dlnCdlnRanalytic[a_] :=((-(1.-CPeebles[a]))/(1. + ((Λ2s)/(3 RLyα[a])))) ;
dlnCdlnRLog = Interpolation[Table[{logatab[[i]], dlnCdlnRanalytic[Exp[logatab[[i]]]]}, {i, 1, Na}]];
dlnCdlnR[a_] := dlnCdlnRLog[Log[a]];
LogLinearPlot[{-dlnCdlnRanalytic[1./zp1], - dlnCdlnR[1./zp1]}, {zp1, 201, 1000}, PlotRange->Full]

coeffTm1log = Interpolation[Table[{logatab[[i]], SeriesCoefficient[Log[αB[Tcmb[Exp[logatab[[i]]]], Tgas[Exp[logatab[[i]]]](1+ϵ)]], {ϵ, 0, 1}]}, {i, 1, Na}]];

coeffTm1log2 = Interpolation[Table[{Log[logatab[[i]]], SeriesCoefficient[Log[αB[Tcmb[Exp[Log[logatab[[i]]]]], Tgas[Exp[Log[logatab[[i]]]]](1+ϵ)]], {ϵ, 0, 1}]}, {i, 1, Na}]];

dlnAdlnTm[a_] := coeffTm1log[Log[a]]; (* dlog(αB)/dlog(Tgas) *)

LogLinearPlot[{dlnAdlnTm[1./zp1]}, {zp1, 21., 1000}]

Coupled ODEs for baryons, CDM, temperature and ionization fraction

Transfer functions from CAMB at z = 1020 and 1000

data1020 = Import["~/Dropbox/21cm Dark Ages/data/mys_transfer_outz1020.dat","Table"];
data1000 = Import["~/Dropbox/21cm Dark Ages/data/mys_transfer_outz1000.dat","Table"];

dim =Dimensions[data1000]
NkCAMB = dim[[1]];
comk = Table[h *data1000[[i,1]],{i,1,NkCAMB}]; 
Tδc1020=Interpolation[Table[{comk[[i]],data1020[[i,2]]},{i,1,NkCAMB}], Method->"Spline"];
Tδb1020 = Interpolation[Table[{comk[[i]],data1020[[i,3]] },{i,1,NkCAMB}], Method->"Spline"];
Tδc1000=Interpolation[Table[{comk[[i]],data1000[[i,2]]},{i,1,NkCAMB}], Method->"Spline"];
Tδb1000 = Interpolation[Table[{comk[[i]],data1000[[i,3]] },{i,1,NkCAMB}], Method->"Spline"];

Initial conditions at z = 1010

δbinit[k_] := (k^2) ((Tδb1020[k]+Tδb1000[k])/2.) Sqrt[As (((k/kpivot))^(ns-1))];
δcinit[k_] :=k^2 ((Tδc1020[k]+Tδc1000[k])/2.)Sqrt[As (((k/kpivot))^(ns-1))];
θbinit[k_] :=  ((Hubble[ainit])/ainit) k^2 ((Tδb1020[k]-Tδb1000[k])/20.) Sqrt[As (((k/kpivot))^(ns-1))];
θcinit[k_] := ((Hubble[ainit])/ainit) k^2 ((Tδc1020[k]-Tδc1000[k])/20.) Sqrt[As(((k/kpivot))^(ns-1))];

ratio[k_] := 2.42*^-9/(As(((k/kpivot))^(ns-1)));

e=LogLinearPlot[{-ainit θbinit[k]/k, -ainit θcinit[k]/k, ainit( θbinit[k] - θcinit[k])/k, 0}, {k, 0.001, 5.},FrameLabel->{"𝓀 [(Mpc^(-1))]","v [km/s]"},FrameStyle->Directive[Thickness[0.002]],PlotStyle->{{Black},{Red},{Green},{Blue}},ImageSize->700,PlotRange->{{(10^(-3)),1},{-30,40}},AspectRatio->.4, Axes->None];
paramWithLegend5=autoLegend[e,{"baryons","CDM","relative"},Background->Directive[White],Alignment->{.45,-.4}];
Magnify[Show[paramWithLegend5,ImageSize->600],1]

Pvbc[k_] :=ainit^2  ((((θcinit[k] - θbinit[k]))^2)/(k^2)) ;

σ3d = Sqrt[NIntegrate[Pvbc[k]/k,{k,0.0001,10.}]];
σ1d = σ3d/Sqrt[3.];

Solve ODE system and tabulate solutions

Clear[k, vcb, a];

ODEsystem[k_, vbc_, a_] = { (* vbc is vbc.hat{k} here *)
δc'[a]== (( 1.)/(a Hubble[a]))(I((k vbc ainit)/(a^2)) δc[a]-θc[a]), 
θc'[a]==I ((k vbc ainit θc[a])/((a^3) Hubble[a]))-1.5(((H0^2)(Ωc0 δc[a]+Ωb0 δb[a]))/( (a^4)Hubble[a] ))-2.((θc[a])/a),
δb'[a]== ((-θb[a])/(a Hubble[a])),
θb'[a]== -1.5(((H0^2)(Ωc0 δc[a]+Ωb0 δb[a]))/((a^4)Hubble[a] ))-2.((θb[a])/a)+((cs2[a](k^2))/( (a^3) Hubble[a]))(δb[a]+δT[a]),
δT'[a] ==(2./3.) δb'[a]+((ΓCompton[a])/(a Hubble[a]))(( ((Tcmb[a])/(Tgas[a]))-1.)(δe[a]- δT[a])-  δT[a]) , δe'[a] ==((xe'[a])/(xe[a] ))(δe[a] + δb[a] + dlnCdlnR[a](((θb[a])/(3. Hubble[a])) - δb[a]) +dlnAdlnTm[a] δT[a]), 
δb[ainit]==δbinit[k],δc[ainit]==δcinit[k],
θb[ainit]==θbinit[k],θc[ainit]==θcinit[k],
δT[ainit]==0., δe[ainit] == 0.};

δcshorttab = Table[0., {i, 1, Nk}, {j, 1, Nvbc}];
δbshorttab = Table[0., {i, 1, Nk}, {j, 1, Nvbc}];
δTshorttab = Table[0., {i, 1, Nk}, {j, 1, Nvbc}];
θbshorttab = Table[0., {i, 1, Nk}, {j, 1, Nvbc}];
δeshorttab = Table[0., {i, 1, Nk}, {j, 1, Nvbc}];

Clear[k, vbc, a];

Do[k=ktab[[i]]; vbc =vbctab[[j]];(* vbc is vbc.hat{k} here *)
sol=Flatten[NDSolve[ODEsystem[k, vbc, a],{δb,δc,θb,θc, δT, δe},{a,ainit,aend}]];
δcshorttab[[i,j]] = δc/.sol;
δbshorttab[[i,j]] = δb/.sol;
δTshorttab[[i,j]] = δT/.sol; 
θbshorttab[[i,j]] = θb/.sol;
δeshorttab[[i,j]] = δe/.sol;
,{j,1,Nvbc},{i,1,Nk}]

Checking against the result of Tseliakhovich & Hirata for ((Δm)^2)

fb =  ((Ωb0)/(Ωc0+Ωb0))

δmtab = Table[0., {i, 1, Nk}, {j, 1, Nvbc}, {n, 1, Na}];
Do[a = atab[[n]]; δmtab[[i,j,n]] =(1.-fb)δcshorttab[[i,j]][atab[[n]]] + fb δbshorttab[[i,j]][atab[[n]] ];, {i, 1, Nk}, {j, 1, Nvbc}, {n, 1, Na}];

Δc2novbc = Interpolation[Flatten[Table[{{ktab[[i]], atab[[n]]},Abs[δmtab[[i,1, n]]]^2}, {i, 1, Nk}, {n, 1, Na}], 1], Method->"Spline"];

Prob1d[vdotk_] := Sqrt[2./Pi]/σ1d Exp[-0.5 (((vdotk/σ1d))^2)];

Dvbc = (vbcmax/(Nvbc-1.));

Prob1dtab = Table[ Prob1d[(j-1)Dvbc]Dvbc, {j, 1, Nvbc}]; 
Prob1dtab[[1]] = 0.5 Prob1dtab[[1]];
Prob1dtab[[Nvbc]] = 0.5 Prob1dtab[[Nvbc]];

Δc2averaged = Interpolation[Flatten[Table[{{ktab[[i]], atab[[n]]},Sum[Abs[δmtab[[i,j, n]]]^2 Prob1dtab[[j]], {j, 1, Nvbc}]}, {i, 1, Nk}, {n, 1, Na}], 1], Method->"Spline"];

e=LogLinearPlot[{Δc2novbc[k, 1./41.] ratio[k], Δc2averaged[k, 1./41.] ratio[k]}, {k, kmin, kmax},FrameLabel->{"𝓀 [(Mpc^(-1))]","((Subscript[Δ, m])^2)(k)"},FrameStyle->Directive[Thickness[0.002]],PlotStyle->{{Black},{Red},{Green},{Blue}},ImageSize->700,PlotRange->All,AspectRatio->.4, Axes->None];
paramWithLegend5=autoLegend[e,{"No (Subscript[𝓋, bc])","TH w (Subscript[𝓋, bc])"},Background->Directive[White],Alignment->{.45,-.4}];
Magnify[Show[paramWithLegend5,ImageSize->600],1]