Quit[]

Cosmological Parameters, Time-Redshift Relation

(*Here we define the values of various cosmological parameters, and the relationship between the age of the universe and redshift, neglecting the cosmological constant.*)
melec = 510.998903 10^3;
OmegaM = 0.266;
OmegaRad = 8 10^(-5);
H0 = 1/(4.34602852 10^17);
Time[redshift_,OmegaM_,OmegaRad_,H0_] := 2 (OmegaM Sqrt[OmegaRad + OmegaM/redshift] + 2 OmegaRad^(3/2) redshift - 2 OmegaRad Sqrt[redshift (OmegaM + OmegaRad redshift)])/(3 OmegaM^2 redshift H0)
(*Approximate relationship neglecting dark energy, will fail at low redshifts - difference is <0.1% at 1+z=10. Note 'redshift' here means 1+z.*)
StandardTime[redshift_] := Time[redshift,OmegaM, OmegaRad, H0]
rhoDM = 0.1116 1.05375 10^(-5) ; (*present-day DM density in (GeV/c^2)/cm^3, taken from the PDG*)

Import Data, Set Up Plotting Defaults

(*Here we import the data files we will use in the remainder of the calculation*)
DataFolder = "~/deposition/data/";
ElecDataFile =DataFolder <>"resultsgrid_elec.hdf";
{outputredshift,energy,inputredshift,fractions,fcheck,conversion} = Import[ElecDataFile,{"Datasets",{"output_redshift","log10(energy/eV)","input_redshift","deposition_fractions","f_check","conversion_factor"}}];
PhotDataFile =DataFolder <>"resultsgrid_phot.hdf";
{outputredshiftphot,energyphot,inputredshiftphot,fractionsphot,fcheckphot,conversionphot} = Import[PhotDataFile,{"Datasets",{"output_redshift","log10(energy/eV)","input_redshift","deposition_fractions","f_check","conversion_factor"}}];
(*Annihilation PCs*)
AnnPCAWMAP =DataFolder <>"epsilon_ann_WMAP7_PCs.hdf";
AnnPCAPlanck =DataFolder <>"epsilon_ann_Planck_PCs.hdf";
AnnPCACVL =DataFolder <>"epsilon_ann_CVL_PCs.hdf";
{LVals,AnnRedshiftZ,AnnWMAPEigenvalues,AnnEpsilon,AnnWMAPPCsZ,AnnWMAPPCs,AnnWMAPPCsPerp} = Import[AnnPCAWMAP,{"Datasets",{"l","z","eigenvalues","normfac","zPCs", "ClPCs","ClPCsperp"}}];
{LVals,AnnRedshiftZ,AnnPlanckEigenvalues,AnnEpsilon,AnnPlanckPCsZ,AnnPlanckPCs,AnnPlanckPCsPerp} = Import[AnnPCAPlanck,{"Datasets",{"l","z","eigenvalues","normfac","zPCs", "ClPCs","ClPCsperp"}}];
{LVals,AnnRedshiftZ,AnnCVLEigenvalues,AnnEpsilon,AnnCVLPCsZ,AnnCVLPCs,AnnCVLPCsPerp} = Import[AnnPCACVL,{"Datasets",{"l","z","eigenvalues","normfac","zPCs", "ClPCs","ClPCsperp"}}];
(*Decay PCs*)
DecPCAWMAP =DataFolder <>"epsilon_dec_WMAP7_PCs.hdf";
DecPCAPlanck =DataFolder <>"epsilon_dec_Planck_PCs.hdf";
DecPCACVL =DataFolder <>"epsilon_dec_CVL_PCs.hdf";
{LVals,DecRedshiftZ,DecWMAPEigenvalues,DecEpsilon,DecWMAPPCsZ,DecWMAPPCs,DecWMAPPCsPerp} = Import[DecPCAWMAP,{"Datasets",{"l","z","eigenvalues","normfac","zPCs", "ClPCs","ClPCsperp"}}];
{LVals,DecRedshiftZ,DecPlanckEigenvalues,DecEpsilon,DecPlanckPCsZ,DecPlanckPCs,DecPlanckPCsPerp} = Import[DecPCAPlanck,{"Datasets",{"l","z","eigenvalues","normfac","zPCs", "ClPCs","ClPCsperp"}}];
{LVals,DecRedshiftZ,DecCVLEigenvalues,DecEpsilon,DecCVLPCsZ,DecCVLPCs,DecCVLPCsPerp} = Import[DecPCACVL,{"Datasets",{"l","z","eigenvalues","normfac","zPCs", "ClPCs","ClPCsperp"}}];
(*Sort into more compact arrays*)
AnnPCsZ = {AnnWMAPPCsZ, AnnPlanckPCsZ, AnnCVLPCsZ};
AnnEigenvalues = {AnnWMAPEigenvalues, AnnPlanckEigenvalues, AnnCVLEigenvalues};
AnnPCs = {AnnWMAPPCs, AnnPlanckPCs, AnnCVLPCs};
AnnPCsPerp = {AnnWMAPPCsPerp, AnnPlanckPCsPerp, AnnCVLPCsPerp};
DecPCsZ = {DecWMAPPCsZ, DecPlanckPCsZ, DecCVLPCsZ};
DecEigenvalues = {DecWMAPEigenvalues, DecPlanckEigenvalues, DecCVLEigenvalues};
DecPCs = {DecWMAPPCs, DecPlanckPCs, DecCVLPCs};
DecPCsPerp = {DecWMAPPCsPerp, DecPlanckPCsPerp, DecCVLPCsPerp};

PlotFolder = "~/deposition/plots/";
PrettyExp[x_,MAX_:0]:=If[Abs[x]>MAX,Superscript[10,x],If[x≥0,(10^x),(10.^x)]]
fticks[min_,max_,AMAX_:0]:=Flatten[Table[Prepend[Flatten[Table[{x+Log[10,i]," ",{.01,0}},{i,2,9}],{1}],{x,PrettyExp[x,AMAX],{.03,0}}],{x,Ceiling[min],Floor[max],1}],1]
ticks={fticks[-5,1],Automatic};
ticksb={fticks[0,12],fticks[1,4],None,None}; (*ticks on left and bottom only*)
tickst={fticks[0,12],fticks[10,26],None,None};
ticksz = {fticks[1,4],fticks[1,4],None,None};
log3[{x_,y_,z_}]:={x,y,Log[10,z]};
clabel=Function[{x,y,z},Text[Framed[z],{x,y},Background->White]];

Make Basic Plots

(*Plots showing various integrals over the energy absorption function. Not necessary for later sections*)
(*Here we plot the aggregate fraction of the particle's energy that is deposited by redshift 10*)
Aggregate = Table[{energy[[j]], Log[10,inputredshift[[k]]],Sum[fractions[[k,j,i]],{i,1,Length[outputredshift]}]},{j,1,Length[energy]},{k,1,Length[inputredshift]}];
thisPlot = Show[{ListContourPlot[Flatten[Aggregate,1],  Contours -> {0.1,0.3,0.5,0.7,0.9},PlotRange ->{All,Log[10,{10,3000}],All}]},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, Epilog->{clabel[10,1.9,0.1],clabel[10,2.4,0.3],clabel[10,2.7,0.5],clabel[9.8,3,0.7], clabel[8,3,0.9]}]
Export[PlotFolder<>"totaldep_electrons.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]
AggregatePhot = Table[{energyphot[[j]], Log[10,inputredshiftphot[[k]]],Sum[fractionsphot[[k,j,i]],{i,1,Length[outputredshiftphot]}]},{j,1,Length[energyphot]},{k,1,Length[inputredshiftphot]}];
thisPlotPhot = Show[{ListContourPlot[Flatten[AggregatePhot,1],  Contours -> {0.01,0.1,0.3,0.5,0.7,0.9},PlotRange ->{All,Log[10,{10,3000}],All}]},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, Epilog->{clabel[10,2.6,0.1],clabel[9,2.75,0.3],clabel[9,3.0,0.5],clabel[7.5,3,0.7], clabel[4,3,0.9],clabel[10,2.2,0.01]}]
Export[PlotFolder<>"totaldep_photons.pdf",thisPlotPhot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

(*Account for lack of positron annihilation in pure electron case*)
PositroniumContinuum[eng_] := Module[{x,lg},x = eng/melec; lg = Log[1 - x]; If[0 < eng < melec,x(1-x)/(2 - x)^2 - 2 (1 - x)^2 lg/(2 - x)^3 + (2 - x)/x + 2lg(1 -x)/x^2,0]]; (*dN/dE*)
NormFac = NIntegrate[PositroniumContinuum[x]*x,{x,0,melec}];
PositroniumContinuumNormed[eng_] := PositroniumContinuum[eng]*(3*melec/2)/NormFac;
PositroniumSpectrum[eng_] := PositroniumContinuumNormed[eng] + (1/2)DiracDelta[eng - melec];
(*Want to multiply positronium spectrum by fractionsphot, appropriately, to find the energy deposition history for such photons - it's a good approximation that the annihilation to positronium happens at the same time as injection. Not true for positrons produced in cascades, but we want to keep those included - only the injected particles are changed.*)
PositroniumContinuumFractionTable = Quiet[Table[{Log[10,outputredshiftphot[[i]]],Log[10,inputredshiftphot[[j]]],If[j ≥ i,NIntegrate[(Interpolation[Table[{energyphot[[k]],fractionsphot[[j,k,i]]},{k,1,Length[energyphot]}]])[Log[10,x]]* PositroniumContinuumNormed[x],{x,0,melec}]/(9/4),0]},{i,1,Length[outputredshiftphot]},{j,1,Length[inputredshiftphot]}]];
PositronDeltaFractionTable = Quiet[Table[{Log[10,outputredshiftphot[[i]]],Log[10,inputredshiftphot[[j]]],If[j ≥  i,(Interpolation[Table[{energyphot[[k]],fractionsphot[[j,k,i]]},{k,1,Length[energyphot]}]])[Log[10,melec]],0]},{i,1,Length[outputredshiftphot]},{j,1,Length[inputredshiftphot]}]];
fractionselec =  Table[((melec + 10^energy[[j]])*fractions[[i,j,k]] - melec*(0.75*PositroniumContinuumFractionTable[[k,i,3]] + 0.25*PositronDeltaFractionTable[[k,i,3]]))/10^energy[[j]],{i,1,Length[inputredshift]},{j,1,Length[energy]},{k,1,Length[outputredshift]}];
(*'Aggregate' plot, but showing fraction of KINETIC energy deposited in electron case only (assuming no way to liberate mass energy, and extra electron adds exactly one free electron, so assuming E_kin >> 13.6 eV, should be very subdominant)*)
AggregateElecOnly = Table[{energy[[j]], Log[10,inputredshift[[k]]],Min[Sum[fractionselec[[k,j,i]],{i,1,Length[outputredshift]}],1]},{j,1,Length[energy]},{k,12,Length[inputredshift]}];
thisPlot = Show[{ListContourPlot[Flatten[AggregateElecOnly,1],  Contours -> {0.1,0.3,0.5,0.7,0.9, 0.95, 0.99}]},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, PlotRange ->{ {5.5,12},Log[10,{Min[inputredshift]*3.0,inputredshift[[Length[inputredshift]]]}]},Epilog->{clabel[10,1.9,0.1],clabel[10,2.4,0.3],clabel[10,2.7,0.5],clabel[9.8,3,0.7], clabel[8,3,0.9], clabel[7.7,3.2,0.95],clabel[7.4,2.2,0.99]  }]
Export[PlotFolder<>"totaldep_electrons_only.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

Redshift-dependence functions (for use in examples)

(*Define sample rescaling functions. We are making the approximation that all decaying/oscillating particles are non-relativistic.*)
DecayFunction[redshift_,decaytime_] := Exp[-(StandardTime[redshift]/10^7)/(decaytime/10^7)];
(*Mathematica seems to have some trouble dealing with exponentials of large numbers; dividing by 10^7 in both numerator and denominator improves stability*)
OscillationFunction[redshift_,oscillationtime_] := 1  - Exp[-(StandardTime[redshift]/10^7)/(oscillationtime/10^7)];
DecayAnnConversion[redshift_] := redshift^(-3);
DecaySwitchFunc[redshift_] := If[DecaySwitch, DecayAnnConversion[redshift],1]

Spectrum for electrons from DM annihilation/decay to muons (for use in examples)

MuonElectronSpectrum[eng_,muoneng_] := (12/muoneng)*((1/4)(1 - (eng/muoneng)^2) - (1/9)(1 - (eng/muoneng)^3)) 
(*Spectrum of electrons (dN/dE) produced by decay of muons from DM annihilation/decay. Approximate form okay for mDM >> muon mass*)
MuonElectronSpectrumNormedApprox[eng_,muoneng_,muonmass_] := MuonElectronSpectrum[eng,muoneng]

MuonFrameElectrons[eng_,muonmass_] := eng^2((3/2) (muonmass/2) - eng)
BoostedEng[costheta_, gamma_, eng_] := gamma eng (1 + costheta Sqrt[1 - 1/gamma^2])
BaseEng[costheta_,gamma_,boostedeng_] := boostedeng/(gamma(1 + costheta Sqrt[1 - 1/gamma^2]))
MuonIntegral[gamma_,boostedeng_, muonmass_]:= Module[{eng,int, val},eng[costheta_] := BaseEng[costheta,gamma, boostedeng];int := Assuming[{gamma >  1, boostedeng < gamma muonmass (1 + Sqrt[1 - 1/gamma^2])/2},Integrate[MuonFrameElectrons[eng[costheta],muonmass](eng[costheta]/boostedeng), costheta]]; val = Simplify[(int /. costheta -> 1)] - Simplify[(int /. costheta -> Max[-1,((2 boostedeng / (gamma muonmass) - 1)/Sqrt[1 - 1/gamma^2])])]; val] 
(*Spectrum of electrons (dN/dE) produced by decay of muons from DM annihilation/decay. More precise form accurate for any value of gamma*)
UpperBoundMuonElectronSpectrum[gamma_,muonmass_] := gamma muonmass (1 + Sqrt[1 - 1/gamma^2])/2;
MuonElectronSpectrumNormed[eng_,muoneng_,muonmass_] := 32*MuonIntegral[muoneng/muonmass,eng,muonmass]/muonmass^4

(*Example comparison between approximate and exact spectra*)
Plot[{eng MuonElectronSpectrumNormed[eng,1,105.7/1000], eng * MuonElectronSpectrumNormedApprox[eng,1,105.7/1000]},{eng,0,1},LabelStyle->16, FrameLabel ->{"E (GeV)","dN/dlnE"},  Frame -> {True,True,False,False}]

Choose  and implement general inputs for Examples

DecaySwitch = True; (*switch to False for annihilation, True for decay*)
(*Choose relevant spectrum and rescaling function*)
Spectrum[eparam_,mwimp_] := MuonElectronSpectrumNormedApprox[eparam,mwimp,1057 10^5]
TimeScale = (10.0)^13; (*Input decay lifetime in seconds*)
RescaleFunction[redshift_] := DecayFunction[redshift,TimeScale];
(*Choose inputs for computing CMB constraints*)
ExptLabel =2; (* 1 = WMAP7, 2 = Planck, 3 = CVL *)
nPCs = 10; (*Number of PCs to keep*)

(*Read in principal components in z-space*)
AnnPCTable = Table[Interpolation[Table[{AnnRedshiftZ[[i]],AnnPCsZ[[ExptLabel,j,i]]},{i,1,Length[AnnRedshiftZ]}]],{j,1,nPCs}];
DecPCTable = Table[Interpolation[Table[{DecRedshiftZ[[i]],DecPCsZ[[ExptLabel,j,i]]},{i,1,Length[DecRedshiftZ]}]],{j,1,nPCs}];
(*Determine normalization of principal component functions*)
RescaleAnnPCTable = Table[NIntegrate[(AnnPCTable[[i]])[z]^2,{z,Min[AnnRedshiftZ],Max[AnnRedshiftZ]}],{i,1,nPCs}];
RescaleDecPCTable = Table[NIntegrate[(DecPCTable[[i]])[z]^2/z,{z,Min[DecRedshiftZ],Max[DecRedshiftZ]}],{i,1,nPCs}];
(*Specify whether PCs optimized for annihilation-like or decay-like histories are being used*)
(*Determine information for PC analysis depending on whether scenario is annihilation-like or decay-like*)
PCTable = If[DecaySwitch,DecPCTable,AnnPCTable];
EigenvalueTable = If[DecaySwitch, DecEigenvalues,AnnEigenvalues];
Epsilon = If[DecaySwitch,DecEpsilon,AnnEpsilon];
RescaleTable = If[DecaySwitch,RescaleDecPCTable,RescaleAnnPCTable];
RedshiftRange = If[DecaySwitch,{Min[DecRedshiftZ],Max[DecRedshiftZ]},{Min[AnnRedshiftZ],Max[AnnRedshiftZ]}];
LogVsLin[redshift_] := If[DecaySwitch, 1/redshift,1]

Example: Effective f-curves

(*Compute effective f curve*)
EffectiveFGrid= Flatten[Table[{{energy[[j]], Log[10,outputredshift[[i]]]},Log[10,(conversion[[i]]/DecaySwitchFunc[outputredshift[[i]]])Sum[fractions[[k,j,i]] RescaleFunction[inputredshift[[k]]] DecaySwitchFunc[inputredshift[[k]]]/conversion[[k]],{k,1,Length[inputredshift]}]]},{j,1,Length[energy]},{i,1,Length[outputredshift]}],1];
EffectiveFInterpolation=Interpolation[EffectiveFGrid];
EffectiveFCurve[redshift_,mwimp_] := NIntegrate[eparam Spectrum[eparam,mwimp]*10^EffectiveFInterpolation[Log[10,eparam], Log[10,redshift]],{eparam,1,mwimp}]/mwimp

(*Make effective f curve for a specified WIMP mass*)
mWIMP = 10^9 ; (*WIMP mass in eV*)
ApproxFCurve = Quiet[FunctionInterpolation[EffectiveFCurve[x,mWIMP],{x,Min[outputredshift],Max[outputredshift]}]]
LogLogPlot[ApproxFCurve[redshift],{redshift,Min[outputredshift],Max[outputredshift]}, PlotRange -> All, AxesOrigin -> {0,0}, LabelStyle -> 16, FrameLabel ->{"1+z","f(z)"},LabelStyle->16,  Frame -> {True,True,True, True}]

Example: Principal component decomposition and limits

(*Compute coefficients of principal components for chosen f(z) curve*)
CoeffList = Table[Quiet[NIntegrate[ApproxFCurve[redshift]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs}]

(*Constraints from WMAP7*)
WMAPSig[n_] := If[ExptLabel == 1,Sqrt[Sum[CoeffList[[i]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,n}]], "Please switch ExptLabel to 1 and rerun"]
If[ExptLabel == 1,ListPlot[Table[{n,WMAPSig[n]},{n,1,10}],LabelStyle->16, FrameLabel ->{"E (GeV)","dN/dlnE"},  Frame -> {True,True,False,False}, Axes -> None, Filling -> Axis, Joined -> True, Mesh -> Full]]
WMAPLinearLimit = 2*Epsilon[[1]]/WMAPSig[nPCs];If[ExptLabel == 1, "WMAP Limit from Linear Calculation =" <> ToString[WMAPLinearLimit, TraditionalForm]]
WMAPDecayDensityLimit = If[DecaySwitch,If[ExptLabel == 1, "Limit on Decaying Species Fraction =" <> ToString[TimeScale*WMAPLinearLimit, TraditionalForm], "Please switch ExptLabel to 1 and rerun"]]
WMAPLimit = If[!DecaySwitch,If[ExptLabel ==1,"WMAP Limit from First PC + Likelihood Analysis =" <>ToString[1.2 10^(-26) / CoeffList[[1]],TraditionalForm]]]

(*Constraints from Planck*)
PlanckSig[n_] := If[ExptLabel == 2,Sqrt[Sum[CoeffList[[i]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,n}]], "Please switch ExptLabel to 2 and rerun"]
If[ExptLabel == 2,ListPlot[Table[{n,PlanckSig[n]},{n,1,nPCs}],LabelStyle->16, FrameLabel ->{"Number of PCs","Significance of test curve (sigma)"},  Frame -> {True,True,False,False}, Axes -> None, Filling -> Axis, Joined -> True, Mesh -> Full]]
PlanckLinearLimit = 2*Epsilon[[1]]/PlanckSig[nPCs];If[ExptLabel == 2, "Planck Limit from Linear Calculation =" <> ToString[PlanckLinearLimit, TraditionalForm], "Please switch ExptLabel to 2 and rerun"]
PlanckDecayDensityLimit = If[DecaySwitch,If[ExptLabel == 2, "Limit on Decaying Species Fraction =" <> ToString[TimeScale*PlanckLinearLimit, TraditionalForm]]]
PlanckLimit = If[!DecaySwitch,If[ExptLabel ==2,"Planck Limit from First PC + Likelihood Analysis =" <>ToString[2.2 10^(-27) / CoeffList[[1]],TraditionalForm]]]
SigmaMC = {1.1, 2.4, 4.1} *10^(-27);
PlanckThreePCLimit = If[!DecaySwitch,If[ExptLabel ==2,"Planck Limit from Three PCs + Likelihood Analysis =" <>ToString[2/Sqrt[Sum[ CoeffList[[i]]^2 /SigmaMC[[i]]^2,{i,1,3}]],TraditionalForm]]]

(*Constraints from CVL*)
CVLSig[n_] := If[ExptLabel == 3,Sqrt[Sum[CoeffList[[i]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,n}]], "Please switch ExptLabel to 3 and rerun"]
If[ExptLabel == 3,ListPlot[Table[{n,CVLSig[n]},{n,1,nPCs}],LabelStyle->16, FrameLabel ->{"E (GeV)","dN/dlnE"},  Frame -> {True,True,False,False}, Axes -> None, Filling -> Axis, Joined -> True, Mesh -> Full]]
CVLLinearLimit = 2*Epsilon[[1]]/CVLSig[nPCs];If[DecaySwitch,If[ExptLabel == 3, "CVL Limit from Linear Calculation =" <> ToString[CVLLinearLimit, TraditionalForm]]]
CVLDecayDensityLimit = If[DecaySwitch,If[ExptLabel == 3, "Limit on Decaying Species Fraction =" <> ToString[TimeScale*CVLLinearLimit, TraditionalForm]]]
CVLLimit = If[!DecaySwitch,If[ExptLabel ==3,"CVL Limit from First PC + Likelihood Analysis =" <>ToString[10^(-27) / CoeffList[[1]],TraditionalForm]]]
SigmaMC = {0.5,1.1, 1.8,2.5, 3.4} *10^(-27);
CVLFivePCLimit = If[!DecaySwitch,If[ExptLabel ==3,"CVL Limit from Five PCs + Likelihood Analysis =" <>ToString[2/Sqrt[Sum[ CoeffList[[i]]^2 /SigmaMC[[i]]^2,{i,1,5}]],TraditionalForm]]]

Example: Scan of constraints vs decay lifetime

(*Set lifetime grid, compute effective f grid*)
LifetimeGrid = Flatten[{Table[10^logt,{logt,12,15,0.2}],Table[10^logt,{logt,16,28,2}]}];
EffectiveFGridDecay = Table[{{energy[[j]], Log[10,outputredshift[[i]]], Log[10,LifetimeGrid[[l]]]},Log[10,(conversion[[i]]/DecaySwitchFunc[outputredshift[[i]]])Sum[fractions[[k,j,i]] DecayFunction[inputredshift[[k]],LifetimeGrid[[l]]] DecaySwitchFunc[inputredshift[[k]]]/conversion[[k]],{k,1,Length[inputredshift]}]]},{j,1,Length[energy]},{i,1,Length[outputredshift]},{l,1,Length[LifetimeGrid]}];

EffectiveFGridDecayPhot = Table[{{energyphot[[j]], Log[10,outputredshiftphot[[i]]], Log[10,LifetimeGrid[[l]]]},Log[10,(conversionphot[[i]]/DecaySwitchFunc[outputredshiftphot[[i]]])Sum[fractionsphot[[k,j,i]] DecayFunction[inputredshiftphot[[k]],LifetimeGrid[[l]]] DecaySwitchFunc[inputredshiftphot[[k]]]/conversionphot[[k]],{k,1,Length[inputredshiftphot]}]]},{j,1,Length[energyphot]},{i,1,Length[outputredshiftphot]},{l,1,Length[LifetimeGrid]}];

LogEffectiveFCurveDecay = Interpolation[Flatten[EffectiveFGridDecay,2]]
LogEffectiveFCurveDecayPhot = Interpolation[Flatten[EffectiveFGridDecayPhot,2]]

(*Compute coefficients of principal components for the derived f(z) curves*)
PlanckCoeffListDecay = If[ExptLabel ==2,Table[Quiet[NIntegrate[10^LogEffectiveFCurveDecay[energy[[j]],Log[10,redshift],Log[10,LifetimeGrid[[k]]]]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs},{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"];

PlanckCoeffListDecayPhot = If[ExptLabel ==2,Table[Quiet[NIntegrate[10^LogEffectiveFCurveDecayPhot[energyphot[[j]],Log[10,redshift],Log[10,LifetimeGrid[[k]]]]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs},{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"];

(*Constraints from Planck*)
PlanckSigTable = If[ExptLabel == 2,Table[Sqrt[Sum[PlanckCoeffListDecay[[i,j,k]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,nPCs}]],{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}], "Please switch ExptLabel to 2 and rerun"];
PlanckLinearLimitTable = 2*Epsilon[[1]]/PlanckSigTable;
PlanckDecayDensityLimitTable = If[DecaySwitch,If[ExptLabel == 2,Table[LifetimeGrid[[k]]*PlanckLinearLimitTable[[j,k]],{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"],"Please switch DecaySwitch to True and rerun"];

PlanckSigTablePhot = If[ExptLabel == 2,Table[Sqrt[Sum[PlanckCoeffListDecayPhot[[i,j,k]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,nPCs}]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}], "Please switch ExptLabel to 2 and rerun"];
PlanckLinearLimitTablePhot = 2*Epsilon[[1]]/PlanckSigTablePhot;
PlanckDecayDensityLimitTablePhot = If[DecaySwitch,If[ExptLabel == 2,Table[LifetimeGrid[[k]]*PlanckLinearLimitTablePhot[[j,k]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"],"Please switch DecaySwitch to True and rerun"];

WMAPCoeffListDecay = If[ExptLabel == 1,Table[Quiet[NIntegrate[10^LogEffectiveFCurveDecay[energy[[j]],Log[10,redshift],Log[10,LifetimeGrid[[k]]]]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs},{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}] ,"Please switch ExptLabel to 2 and rerun"];

WMAPCoeffListDecayPhot = If[ExptLabel == 1,Table[Quiet[NIntegrate[10^LogEffectiveFCurveDecayPhot[energyphot[[j]],Log[10,redshift],Log[10,LifetimeGrid[[k]]]]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs},{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}] ,"Please switch ExptLabel to 2 and rerun"];

(*Constraints from WMAP*)
WMAPSigTable = If[ExptLabel == 1,Table[Sqrt[Sum[WMAPCoeffListDecay[[i,j,k]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,nPCs}]],{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}], "Please switch ExptLabel to 1 and rerun"];
WMAPLinearLimitTable = 2*Epsilon[[1]]/WMAPSigTable;
WMAPDecayDensityLimitTable = If[DecaySwitch,If[ExptLabel == 1,Table[LifetimeGrid[[k]]*WMAPLinearLimitTable[[j,k]],{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"],"Please switch DecaySwitch to 2 and rerun"];

WMAPSigTablePhot= If[ExptLabel == 1,Table[Sqrt[Sum[WMAPCoeffListDecayPhot[[i,j,k]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,nPCs}]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}], "Please switch ExptLabel to 1 and rerun"];
WMAPLinearLimitTablePhot = 2*Epsilon[[1]]/WMAPSigTablePhot;
WMAPDecayDensityLimitTablePhot = If[DecaySwitch,If[ExptLabel == 1,Table[LifetimeGrid[[k]]*WMAPLinearLimitTablePhot[[j,k]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"],"Please switch DecaySwitch to 2 and rerun"];

WMAPDensityLimitTable = Table[{energy[[j]],Log[10,LifetimeGrid[[i]]],Log[10,WMAPDecayDensityLimitTable[[j,i]]]},{i,1,Length[LifetimeGrid]-1}, {j,1,Length[energy]}];
thisPlot = Show[{ListContourPlot[Flatten[WMAPDensityLimitTable,1],  Contours ->  {-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1}]},FrameLabel ->{"Energy (eV)","Lifetime (s)"},LabelStyle->16, FrameTicks->tickst,Epilog->{clabel[3.5,14.1,PrettyExp[-10]],clabel[3.5,15.3,PrettyExp[-9]],clabel[3.5,16.3,PrettyExp[-8]],clabel[3.5,17.3,PrettyExp[-7]], clabel[3.5,18.3,PrettyExp[-6]],clabel[3.5,19.3,PrettyExp[-5]],clabel[3.5,20.3,PrettyExp[-4]],clabel[3.5,21.3,PrettyExp[-3]],clabel[3.5,22.3,PrettyExp[-2]],clabel[3.5,23.3,PrettyExp[-1]],clabel[3.5,24.3,PrettyExp[0]],clabel[3.5,25.3,10]}]
Export[PlotFolder<>"lifetimelimit_WMAP_contour.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]
WMAPDensityLimitTablePhot= Table[{energy[[j]],Log[10,LifetimeGrid[[i]]],Log[10,WMAPDecayDensityLimitTablePhot[[j,i]]]},{i,1,Length[LifetimeGrid]-1}, {j,1,Length[energyphot]}];
thisPlot = Show[{ListContourPlot[Flatten[WMAPDensityLimitTablePhot,1],  Contours ->  {-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1}]},FrameLabel ->{"Energy (eV)","Lifetime (s)"},LabelStyle->16, FrameTicks->tickst,Epilog->{clabel[3.5,14.1,PrettyExp[-10]],clabel[3.5,15.3,PrettyExp[-9]],clabel[3.5,16.3,PrettyExp[-8]],clabel[3.5,17.4,PrettyExp[-7]], clabel[3.5,18.4,PrettyExp[-6]],clabel[3.5,19.4,PrettyExp[-5]],clabel[3.5,20.4,PrettyExp[-4]],clabel[3.5,21.4,PrettyExp[-3]],clabel[3.5,22.4,PrettyExp[-2]],clabel[3.5,23.4,PrettyExp[-1]],clabel[3.5,24.4,PrettyExp[0]],clabel[3.5,25.4,10]}]
Export[PlotFolder<>"lifetimelimit_WMAP_phot_contour.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

PlanckDensityLimitTable = Table[{energy[[j]],Log[10,LifetimeGrid[[i]]],Log[10,PlanckDecayDensityLimitTable[[j,i]]]},{i,1,Length[LifetimeGrid]-1}, {j,1,Length[energy]}];
thisPlot = Show[{ListContourPlot[Flatten[PlanckDensityLimitTable,1],  Contours ->  {-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0}]},FrameLabel ->{"Energy (eV)","Lifetime (s)"},LabelStyle->16, FrameTicks->tickst,Epilog->{clabel[3.5,15.0,PrettyExp[-10]],clabel[3.5,16.1,PrettyExp[-9]],clabel[3.5,17.1,PrettyExp[-8]],clabel[3.5,18.1,PrettyExp[-7]], clabel[3.5,19.1,PrettyExp[-6]],clabel[3.5,20.1,PrettyExp[-5]],clabel[3.5,21.1,PrettyExp[-4]],clabel[3.5,22.1,PrettyExp[-3]],clabel[3.5,23.1,PrettyExp[-2]],clabel[3.5,24.1,PrettyExp[-1]],clabel[3.5,25.1,PrettyExp[0]]}]
Export[PlotFolder<>"lifetimelimit_Planck_contour.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]
PlanckDensityLimitTablePhot= Table[{energy[[j]],Log[10,LifetimeGrid[[i]]],Log[10,PlanckDecayDensityLimitTablePhot[[j,i]]]},{i,1,Length[LifetimeGrid]-1}, {j,1,Length[energyphot]}];
thisPlot = Show[{ListContourPlot[Flatten[PlanckDensityLimitTablePhot,1],  Contours ->  {-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0}]},FrameLabel ->{"Energy (eV)","Lifetime (s)"},LabelStyle->16, FrameTicks->tickst,Epilog->{clabel[3.5,15.0,PrettyExp[-10]],clabel[3.5,16.1,PrettyExp[-9]],clabel[3.5,17.1,PrettyExp[-8]],clabel[3.5,18.1,PrettyExp[-7]], clabel[3.5,19.1,PrettyExp[-6]],clabel[3.5,20.1,PrettyExp[-5]],clabel[3.5,21.1,PrettyExp[-4]],clabel[3.5,22.1,PrettyExp[-3]],clabel[3.5,23.1,PrettyExp[-2]],clabel[3.5,24.1,PrettyExp[-1]],clabel[3.5,25.1,PrettyExp[0]]}]
Export[PlotFolder<>"lifetimelimit_Planck_phot_contour.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

ListLogLogPlot[Flatten[{Table[Table[{LifetimeGrid[[i]],WMAPDecayDensityLimitTable[[j,i]]},{i,1,Length[LifetimeGrid]}], {j,1,32}],Table[Table[{LifetimeGrid[[i]],PlanckDecayDensityLimitTable[[j,i]]},{i,1,Length[LifetimeGrid]}], {j,1,32}]},1], Joined -> True, PlotRange -> {10^(-11),1},Frame->True,FrameLabel -> {"Decay lifetime (s)","Maximum allowed mass fraction of (Subscript[Ω, DM])(h^2)"} ,LabelStyle -> 16]
ListLogLogPlot[Flatten[{Table[Table[{LifetimeGrid[[i]],PlanckDecayDensityLimitTablePhot[[j,i]]},{i,1,Length[LifetimeGrid]}], {j,1,Length[energyphot]}],Table[Table[{LifetimeGrid[[i]],WMAPDecayDensityLimitTablePhot[[j,i]]},{i,1,Length[LifetimeGrid]}], {j,1,Length[energyphot]}]},1],PlotRange -> {10^(-11),1}, Joined -> True ,Frame->True,FrameLabel -> {"Decay lifetime (s)","Maximum allowed mass fraction of (Subscript[Ω, DM])(h^2)"}, LabelStyle -> 16]

WMAPMaxCurve = Table[Max[Flatten[{Table[WMAPDecayDensityLimitTablePhot[[j,i]],{j,1,Length[energyphot]}],Table[WMAPDecayDensityLimitTable[[j,i]],{j,1,Length[energy]}]}]],{i,1,Length[LifetimeGrid]}];
WMAPMinCurve = Table[Min[Flatten[{Table[WMAPDecayDensityLimitTablePhot[[j,i]],{j,1,Length[energyphot]}],Table[WMAPDecayDensityLimitTable[[j,i]],{j,1,Length[energy]}]}]],{i,1,Length[LifetimeGrid]}];
PlanckMaxCurve = Table[Max[Flatten[{Table[PlanckDecayDensityLimitTablePhot[[j,i]],{j,1,Length[energyphot]}],Table[PlanckDecayDensityLimitTable[[j,i]],{j,1,Length[energy]}]}]],{i,1,Length[LifetimeGrid]}];
PlanckMinCurve = Table[Min[Flatten[{Table[PlanckDecayDensityLimitTablePhot[[j,i]],{j,1,Length[energyphot]}],Table[PlanckDecayDensityLimitTable[[j,i]],{j,1,Length[energy]}]}]],{i,1,Length[LifetimeGrid]}];
thisPlot = ListLogLogPlot[{Table[{LifetimeGrid[[i]],WMAPMaxCurve[[i]]},{i,1,Length[LifetimeGrid]}],Table[{LifetimeGrid[[i]],WMAPMinCurve[[i]]},{i,1,Length[LifetimeGrid]}], Table[{LifetimeGrid[[i]],PlanckMaxCurve[[i]]},{i,1,Length[LifetimeGrid]}],Table[{LifetimeGrid[[i]],PlanckMinCurve[[i]]},{i,1,Length[LifetimeGrid]}]}, Joined -> True, PlotRange -> {10^(-11),10},Frame->True,FrameLabel -> {"Decay lifetime (s)","Maximum allowed mass fraction of (Subscript[Ω, DM])(h^2)"} ,LabelStyle -> 16, Filling -> {1 -> {2}, 3 -> {4}}, ImageSize -> 500]
Export[PlotFolder<>"lifetimelimit_bands.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

Example: Scan of constraints vs annihilation timescale (asymmetric DM)

(*Set lifetime grid, compute effective f grid*)
LifetimeGrid = Flatten[{Table[10^logt,{logt,12,15,0.2}],Table[10^logt,{logt,16,18,1}]}];
EffectiveFGridOsc = Table[{{energy[[j]], Log[10,outputredshift[[i]]], Log[10,LifetimeGrid[[l]]]},(conversion[[i]]/DecaySwitchFunc[outputredshift[[i]]])Sum[fractions[[k,j,i]] OscillationFunction[inputredshift[[k]],LifetimeGrid[[l]]] DecaySwitchFunc[inputredshift[[k]]]/conversion[[k]],{k,1,Length[inputredshift]}]},{j,1,Length[energy]},{i,1,Length[outputredshift]},{l,1,Length[LifetimeGrid]}];

EffectiveFGridOscPhot = Table[{{energyphot[[j]], Log[10,outputredshiftphot[[i]]], Log[10,LifetimeGrid[[l]]]},(conversionphot[[i]]/DecaySwitchFunc[outputredshiftphot[[i]]])Sum[fractionsphot[[k,j,i]] OscillationFunction[inputredshiftphot[[k]],LifetimeGrid[[l]]] DecaySwitchFunc[inputredshiftphot[[k]]]/conversionphot[[k]],{k,1,Length[inputredshiftphot]}]},{j,1,Length[energyphot]},{i,1,Length[outputredshiftphot]},{l,1,Length[LifetimeGrid]}];

EffectiveFCurveOsc = Interpolation[Flatten[EffectiveFGridOsc,2]]
EffectiveFCurveOscPhot = Interpolation[Flatten[EffectiveFGridOscPhot,2]]

(*Compute coefficients of principal components for the derived f(z) curves*)
PlanckCoeffListOsc = If[ExptLabel ==2,Table[Quiet[NIntegrate[EffectiveFCurveOsc[energy[[j]],Log[10,redshift],Log[10,LifetimeGrid[[k]]]]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs},{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"];

PlanckCoeffListOscPhot = If[ExptLabel ==2,Table[Quiet[NIntegrate[EffectiveFCurveOscPhot[energyphot[[j]],Log[10,redshift],Log[10,LifetimeGrid[[k]]]]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs},{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"];

(*Constraints from Planck*)
PlanckSigTable = If[ExptLabel == 2,Table[Sqrt[Sum[PlanckCoeffListOsc[[i,j,k]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,nPCs}]],{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}], "Please switch ExptLabel to 2 and rerun"];
PlanckLinearLimitTable = 2*Epsilon[[1]]/PlanckSigTable;
PlanckEpsilonLimitTable = If[!DecaySwitch,If[ExptLabel == 2,Table[2*PlanckLinearLimitTable[[j,k]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"],"Please switch DecaySwitch to False and rerun"];

PlanckSigTablePhot = If[ExptLabel == 2,Table[Sqrt[Sum[PlanckCoeffListOscPhot[[i,j,k]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,nPCs}]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}], "Please switch ExptLabel to 2 and rerun"];
PlanckLinearLimitTablePhot = 2*Epsilon[[1]]/PlanckSigTablePhot;
PlanckEpsilonLimitTablePhot = If[!DecaySwitch,If[ExptLabel == 2,Table[2*PlanckLinearLimitTablePhot[[j,k]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"],"Please switch DecaySwitch to False and rerun"];

WMAPCoeffListOsc = If[ExptLabel == 1,Table[Quiet[NIntegrate[EffectiveFCurveOsc[energy[[j]],Log[10,redshift],Log[10,LifetimeGrid[[k]]]]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs},{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}] ,"Please switch ExptLabel to 2 and rerun"];

WMAPCoeffListOscPhot = If[ExptLabel == 1,Table[Quiet[NIntegrate[EffectiveFCurveOscPhot[energyphot[[j]],Log[10,redshift],Log[10,LifetimeGrid[[k]]]]*(PCTable[[i]])[redshift]*LogVsLin[redshift],{redshift,Min[RedshiftRange],Max[RedshiftRange]}]]/RescaleTable[[i]],{i,1,nPCs},{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}] ,"Please switch ExptLabel to 2 and rerun"];

(*Constraints from WMAP*)
WMAPSigTable = If[ExptLabel == 1,Table[Sqrt[Sum[WMAPCoeffListOsc[[i,j,k]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,nPCs}]],{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}], "Please switch ExptLabel to 1 and rerun"];
WMAPLinearLimitTable = 2*Epsilon[[1]]/WMAPSigTable;
WMAPEpsilonLimitTable = If[!DecaySwitch,If[ExptLabel == 1,Table[2*WMAPLinearLimitTable[[j,k]],{j,1,Length[energy]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"],"Please switch DecaySwitch to False and rerun"];

WMAPSigTablePhot= If[ExptLabel == 1,Table[Sqrt[Sum[WMAPCoeffListOscPhot[[i,j,k]]^2*EigenvalueTable[[ExptLabel,i]],{i,1,nPCs}]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}], "Please switch ExptLabel to 1 and rerun"];
WMAPLinearLimitTablePhot = 2*Epsilon[[1]]/WMAPSigTablePhot;
WMAPEpsilonLimitTablePhot = If[!DecaySwitch,If[ExptLabel == 1,Table[2*WMAPLinearLimitTablePhot[[j,k]],{j,1,Length[energyphot]},{k,1,Length[LifetimeGrid]}],"Please switch ExptLabel to 2 and rerun"],"Please switch DecaySwitch to False and rerun"];

WMAPXSecLimitTable = Table[{energy[[j]],Log[10,LifetimeGrid[[i]]],Log[10,WMAPEpsilonLimitTable[[j,i]]]},{i,1,Length[LifetimeGrid]-1}, {j,1,Length[energy]}];
thisPlot = Show[{ListContourPlot[Flatten[WMAPXSecLimitTable,1],  Contours ->  {-26,-25.5,-25,-24.5,-24,-23.5,-23,-22.5}]},FrameLabel ->{"Energy (eV)","Timescale for annihilation to begin (s)"},LabelStyle->16, FrameTicks->tickst,Epilog->{clabel[3.5,13.6,PrettyExp[-26]],clabel[3.6,14.3,PrettyExp[-25.5]],clabel[3.5,14.8,PrettyExp[-25]],clabel[3.6,15.3,PrettyExp[-24.5]], clabel[3.5,15.8,PrettyExp[-24]],clabel[3.6,16.4,PrettyExp[-23.5]],clabel[3.5,16.8,PrettyExp[-23]],clabel[12.0,16.8,PrettyExp[-22.5]]}]
Export[PlotFolder<>"xseclimit_WMAP_contour.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]
WMAPXSecLimitTablePhot= Table[{energy[[j]],Log[10,LifetimeGrid[[i]]],Log[10,WMAPEpsilonLimitTablePhot[[j,i]]]},{i,1,Length[LifetimeGrid]-1}, {j,1,Length[energyphot]}];
thisPlot = Show[{ListContourPlot[Flatten[WMAPXSecLimitTablePhot,1],  Contours ->  {-26,-25.5,-25,-24.5,-24,-23.5,-23,-22.5}]},FrameLabel ->{"Energy (eV)","Timescale for annihilation to begin (s)"},LabelStyle->16, FrameTicks->tickst,Epilog->{clabel[5.5,13.6,PrettyExp[-26]],clabel[5.6,14.3,PrettyExp[-25.5]],clabel[5.5,14.8,PrettyExp[-25]],clabel[5.6,15.3,PrettyExp[-24.5]], clabel[5.5,15.8,PrettyExp[-24]],clabel[5.6,16.4,PrettyExp[-23.5]],clabel[5.5,16.8,PrettyExp[-23]],clabel[12.0,16.8,PrettyExp[-22.5]]}]
Export[PlotFolder<>"xseclimit_WMAP_phot_contour.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

PlanckXSecLimitTable = Table[{energy[[j]],Log[10,LifetimeGrid[[i]]],Log[10,PlanckEpsilonLimitTable[[j,i]]]},{i,1,Length[LifetimeGrid]-1}, {j,1,Length[energy]}];
thisPlot = Show[{ListContourPlot[Flatten[PlanckXSecLimitTable,1],  Contours ->   {-27,-26.5,-26,-25.5,-25,-24.5,-24,-23.5}]},FrameLabel ->{"Energy (eV)","Timescale for annihilation to begin (s)"},LabelStyle->16, FrameTicks->tickst,Epilog->{clabel[6.0,12.4,PrettyExp[-27]],clabel[3.6,13.9,PrettyExp[-26.5]],clabel[3.5,14.4,PrettyExp[-26]],clabel[3.6,14.9,PrettyExp[-25.5]],clabel[3.5,15.4,PrettyExp[-25]], clabel[3.6,16.0,PrettyExp[-24.5]],clabel[3.5,16.4,PrettyExp[-24]],clabel[12.0,16.5,PrettyExp[-23.5]]}]
Export[PlotFolder<>"xseclimit_Planck_contour.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]
PlanckXSecLimitTablePhot= Table[{energy[[j]],Log[10,LifetimeGrid[[i]]],Log[10,PlanckEpsilonLimitTablePhot[[j,i]]]},{i,1,Length[LifetimeGrid]-1}, {j,1,Length[energyphot]}];
thisPlot = Show[{ListContourPlot[Flatten[PlanckXSecLimitTablePhot,1],  Contours ->   {-27,-26.5,-26,-25.5,-25,-24.5,-24,-23.5,-23}]},FrameLabel ->{"Energy (eV)","Timescale for annihilation to begin (s)"},LabelStyle->16, FrameTicks->tickst,Epilog->{clabel[3.5,12.4,PrettyExp[-27]],clabel[5.6,13.8,PrettyExp[-26.5]],clabel[5.5,14.5,PrettyExp[-26]],clabel[5.6,15.0,PrettyExp[-25.5]],clabel[5.5,15.5,PrettyExp[-25]], clabel[5.6,16.0,PrettyExp[-24.5]],clabel[5.5,16.4,PrettyExp[-24]],clabel[12.0,16.5,PrettyExp[-23.5]]}]
Export[PlotFolder<>"xseclimit_Planck_phot_contour.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

ListLogLogPlot[Flatten[{Table[Table[{LifetimeGrid[[i]],WMAPEpsilonLimitTable[[j,i]]},{i,1,Length[LifetimeGrid]}], {j,1,32}],Table[Table[{LifetimeGrid[[i]],PlanckEpsilonLimitTable[[j,i]]},{i,1,Length[LifetimeGrid]}], {j,1,32}]},1], Joined -> True, PlotRange -> {{10^12,10^18},{10^(-28),10^(-22)}},Frame->True,FrameLabel -> {"Timescale for annihilation to begin (s)","<σv>/(Subscript[m, DM]) ((cm^3)/s/GeV)"} ,LabelStyle -> 16]
ListLogLogPlot[Flatten[{Table[Table[{LifetimeGrid[[i]],PlanckEpsilonLimitTablePhot[[j,i]]},{i,1,Length[LifetimeGrid]}], {j,1,Length[energyphot]}],Table[Table[{LifetimeGrid[[i]],WMAPEpsilonLimitTablePhot[[j,i]]},{i,1,Length[LifetimeGrid]}], {j,1,Length[energyphot]}]},1],PlotRange -> {{10^12,10^18},{10^(-28),10^(-22)}}, Joined -> True ,Frame->True,FrameLabel -> {"Timescale for annihilation to begin (s)","<σv>/(Subscript[m, DM]) ((cm^3)/s/GeV)"}, LabelStyle -> 16]

WMAPMaxCurve = Table[Max[Flatten[{Table[WMAPEpsilonLimitTablePhot[[j,i]],{j,1,Length[energyphot]}],Table[WMAPEpsilonLimitTable[[j,i]],{j,1,Length[energy]}]}]],{i,1,Length[LifetimeGrid]}];
WMAPMinCurve = Table[Min[Flatten[{Table[WMAPEpsilonLimitTablePhot[[j,i]],{j,1,Length[energyphot]}],Table[WMAPEpsilonLimitTable[[j,i]],{j,1,Length[energy]}]}]],{i,1,Length[LifetimeGrid]}];
PlanckMaxCurve = Table[Max[Flatten[{Table[PlanckEpsilonLimitTablePhot[[j,i]],{j,1,Length[energyphot]}],Table[PlanckEpsilonLimitTable[[j,i]],{j,1,Length[energy]}]}]],{i,1,Length[LifetimeGrid]}];
PlanckMinCurve = Table[Min[Flatten[{Table[PlanckEpsilonLimitTablePhot[[j,i]],{j,1,Length[energyphot]}],Table[PlanckEpsilonLimitTable[[j,i]],{j,1,Length[energy]}]}]],{i,1,Length[LifetimeGrid]}];
thisPlot = ListLogLogPlot[{Table[{LifetimeGrid[[i]],WMAPMaxCurve[[i]]},{i,1,Length[LifetimeGrid]}],Table[{LifetimeGrid[[i]],WMAPMinCurve[[i]]},{i,1,Length[LifetimeGrid]}], Table[{LifetimeGrid[[i]],PlanckMaxCurve[[i]]},{i,1,Length[LifetimeGrid]}],Table[{LifetimeGrid[[i]],PlanckMinCurve[[i]]},{i,1,Length[LifetimeGrid]}]}, Joined -> True, PlotRange -> {{10^12,10^18},{10^(-28),10^(-22)}},Frame->True,FrameLabel -> {"Timescale for annihilation to begin (s)","<σv>/(Subscript[m, DM]) ((cm^3)/s/GeV)"} ,LabelStyle -> 16, Filling -> {1 -> {2}, 3 -> {4}}, ImageSize -> 500]
Export[PlotFolder<>"xseclimit_bands.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

Example: Treatment of DM structure formation

DMHaloBoost[redshift_,zh_,fh_] := 1+fh (1 - Erf[redshift/(1+zh)])/redshift^3;
(*The boost to the signal at injection redshift z from the presence of DM halos. zh is the characteristic redshift of halo formation, fh is a normalization parameter. These should be extracted from simulations. The parametric form is taken from 1209.0247.*)
zhval = 20;
fhval = 10^9;

EffectiveFWithHaloBoost = Table[{energy[[j]], Log[10,outputredshift[[i]]],Log[10,conversion[[i]]Sum[DMHaloBoost[inputredshift[[k]], zhval, fhval] fractions[[k,j,i]]/conversion[[k]],{k,1,Length[inputredshift]}]]},{j,1,Length[energy]},{i,1,Length[outputredshift]}];
thisPlot = Show[{ListContourPlot[Flatten[EffectiveFWithHaloBoost,1],  Contours -> Log[10, {0.2,0.4,0.6,0.8, 1.0,10,100,1000,10^4,10^5}],PlotRange -> {All,Log[10,{10,2000}],All}]},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, Epilog->{clabel[7.5,2.1,1],clabel[11.2,1.4,10],clabel[11.2,2.2,0.4],clabel[10.5,2.6,0.6], clabel[10.0,1.4,100],clabel[10.5,3,0.8],clabel[7.0,1.48,10^3],clabel[7.0,1.3,10^4],clabel[7.0,1.15,10^5]}]
Export[PlotFolder<>"fDM_electrons_ann_halos.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]
EffectiveFPhotonsWithHaloBoost = Table[{energyphot[[j]], Log[10,outputredshiftphot[[i]]],Log[10,conversionphot[[i]]Sum[DMHaloBoost[inputredshift[[k]], zhval, fhval]fractionsphot[[k,j,i]]/conversionphot[[k]],{k,1,Length[inputredshiftphot]}]]},{j,1,Length[energyphot]},{i,1,Length[outputredshiftphot]}];
thisPlot = Show[{ListContourPlot[Flatten[EffectiveFPhotonsWithHaloBoost,1],  Contours ->  Log[10,{0.2,0.4,0.6,0.8, 1.0,10,100,1000,10^4,10^5}],PlotRange -> {All,Log[10,{10,2000}],All}]},FrameLabel ->{"Energy (eV)","Redshift (1+z)"},LabelStyle->16, FrameTicks->ticksb, Epilog->{clabel[5.5,1.5,10],clabel[6.5,1.3,100],clabel[12,1.8,0.2],clabel[10.5,2.4,0.4],clabel[9.5,2.7,0.6], clabel[3.5,1.8,1],clabel[8.2,3,0.8]}]
Export[PlotFolder<>"fDM_photons_ann_halos.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

Needs["PlotLegends`"];
thisPlot = ListLogLogPlot[{Table[{10^EffectiveFWithHaloBoost[[25,i,2]],DMHaloBoost[10^EffectiveFWithHaloBoost[[25,i,2]],zhval,fhval]},{i,1,60}],Table[{10^EffectiveFWithHaloBoost[[13,i,2]],10^EffectiveFWithHaloBoost[[13,i,3]]},{i,1,60}],Table[{10^EffectiveFWithHaloBoost[[25,i,2]],10^EffectiveFWithHaloBoost[[25,i,3]]},{i,1,60}],Table[{10^EffectiveFWithHaloBoost[[37,i,2]],10^EffectiveFWithHaloBoost[[37,i,3]]},{i,1,60}]}, PlotRange -> {0.1,10^6}, Joined -> True,FrameLabel -> {"Redshift (1+z)","Energy deposition"},Frame -> True, PlotStyle -> Thick,PlotLegend -> {Style["On-the-spot",16],Style["1 MeV electrons",16],Style["1 GeV electrons",16],Style["1 TeV electrons",16]}, LegendShadow -> None,ImageSize -> 500, LegendPosition -> {-0.15,0.1},LegendSize -> {0.9,0.4},LabelStyle->16]
Export[PlotFolder<>"fDM_halos_comparison.pdf",thisPlot/.{EdgeForm[],r_?(MemberQ[{RGBColor,Hue,CMYKColor,GrayLevel},Head[#]]&),i___}:>{EdgeForm[r],r,i}]

Instructions for use of this notebook

You should first open and execute "Cosmological Parameters" and "Import Data". You will need to modify the values of DataFolder and PlotFolder for your own system, according to where you placed the .hdf files, and where you want to output any plots.

The first cell in "Import Data" will read in the .hdf files: if that's all you needed, congratulations! The files contain six arrays: "output_redshift", "log10(energy/eV)", "input_redshift", "deposition_fractions", "f_check", and "conversion_factor". Note the "redshift" arrays are 1+z, not z. deposition_fractions[i,j,k] describes the fraction of energy deposited in the timestep at redshift output_redshift[i], by a particle injected at energy 10^(log10(energy/ev)[j]) eV and redshift input_redshift[k]. The last two arrays are for convenience and hold no additional information; f_check describes the effective f-curve for annihilation computed directly from the code (so can be used as a check that you're computing the effective f-curve correctly from the deposition_fractions table), while conversion_factor is just the weighting factor applied to deposition_fractions to compute the f(z) curves, H(z)/(1+z)^3 (see the letter for details). If this is unclear, please email me at tslatyer@gmail.com for clarification.

If you wish to reproduce the plots from the letter, go through "Make Basic Plots"; this may also help with building intuition and understanding the data files. "Redshift-dependence functions" has already-implemented basic functions giving the redshift dependence of energy injection for late-decaying species and a sample asymmetric DM scenario with late-time annihilations; if you are working with a model where the annihilations and decays have non-trivial redshift dependence, you can define it here. "Spectrum for electrons from DM annihilation/decay to muons" is just the spectrum of electrons produced by DM annihilation to a μ+ μ- pair, for an example slightly less trivial than a delta function in electron energy; you can similarly insert your own favorite electron/photon spectrum.

In the "Choose and implement" module, you will set the parameters you want for your study. DecaySwitch determines which set of principal components you are using, and the energy injection function that multiplies the effective f-curve; if you are studying a model that involves neither annihilation nor decay, pick whichever is closer (in terms of redshift dependence). Provided you keep enough principal components, it should make no difference to the results.
The Spectrum function gives the spectrum of electrons/positrons/photons you wish to study; define it as the appropriate spectrum for your model.
RescaleFunction describes the redshift dependence of the energy injection, multiplying the basic "decay" or "annihilation" redshift dependence; if you are interested in conventional WIMP annihilation or decay with a lifetime much longer than the age of the universe, you can set this to 1.
ExptSwitch and nPCs are relevant if you are computing estimated CMB constraints.
Executing the rest of this group of cells picks out the appropriate PCs and eigenvalues (according to DecaySwitch), and normalizes them.

The first example generates an effective f-curve by integrating over your energy injection history (defined by RescaleFunction) and spectrum of injected particles (defined by Spectrum).

The second example shows how to compute the estimated CMB constraints. If an effective f-curve has been generated, it computes the dot product of this curve with the principal components (up to nPCs). Executing the cell relevant to your preferred experiment returns a plot showing how the estimated significance of the signal varies with the number of included PCs, the estimated limit on the epsilon parameter (see the letter), and for decays with a specified lifetime, the resulting limit on the fraction of the DM energy density liberated by decays.

The third and fourth examples are scans over a wide range of models, and take considerable time to run. They are included primarily as self-contained examples; none of the previous examples are necessary to execute them. I recommend using a program other than Mathematica for large scans.

The fifth and final example is also self-contained, and includes a simple parameterization for the boost to annihilation from DM structure formation (see the letter for further details).