QEdark-f2 analysis

set working directories

MyDirectory="My/Working/Directory/";(* working directory where files are stored *)
PlotDir=MyDirectory<>"plots/";

load factor files

Import data

(*in 900 bins of Δq=0.02 α (Subscript[m, e]) and 500 E bins of ΔE=0.1 eV*)
rawFFdataSi=Import[MyDirectory<>"C.Si.137.dat","table"];
rawFFdataGe=Import[MyDirectory<>"C.Ge.137.dat","table"];

(* put into square arrays: E across, q down*)
FFarraySi=Block[{Nqbins=900,NEbins=500},
ArrayReshape[rawFFdataSi,{NEbins,Nqbins}]]//Transpose;
FFarrayGe=Block[{Nqbins=900,NEbins=500},
ArrayReshape[rawFFdataGe,{NEbins,Nqbins}]]//Transpose;

(*Binning: 500 E bins of  ΔE = 0.1 eV = (0.1/13.6) (Subscript[R, y]) \[Cross] 900 bins of  Δq = 0.02 α me *)

unit conversions and constants

sec2year=60*60*24*365.25;
amu2kg=1.660538782 (10^(-27));
α=1/137.; (* fine structure constant *)
me=0.511 (10^6); (* electron mass [eV] *)
ρχ=0.4 (10^9); (* local DM density [eV/(cm^3)] *)

c=2.99792458 (10^10); (* speed of light [cm/s] *)
v0=230 (10^5); (* typical velocity [cm/s] *)
vE = 240 (10^5); (* velocity of Earth [cm/s] *)
vesc=600 (10^5); (* escape velocity [cm/s] *)
KK=(v0^3)(-2. Exp[-(vesc^2)/(v0^2)]π vesc/v0+(π^(3/2))Erf[vesc/v0]) (* (([cm/s])^3)*);

velocity distribution and (Subscript[F, DM]) definition

vmin[q_?NumberQ,Ee_?NumberQ,mχ_?NumberQ]:=c((q/(2 mχ))+(Ee/q));

(* Maxwell-Boltzmann distribution 
η((Subscript[v, min]))=∫(((d^3)v)/v)(Subscript[f, MB])(v)θ(v-(Subscript[v, min]))
*)
η[q_?NumberQ,Ee_?NumberQ,mχ_?NumberQ]:=(((v0^2)π)/(2 vE KK))*Piecewise[{{-4 Exp[-(vesc^2)/(v0^2)]vE+(Sqrt[π])v0(Erf[((vmin[q,Ee,mχ]+vE)/v0)]-Erf[((vmin[q,Ee,mχ]-vE)/v0)]),vmin[q,Ee,mχ]<vesc-vE},{-2 Exp[-(vesc^2)/(v0^2)](vE+vesc-vmin[q,Ee,mχ])+(Sqrt[π])v0(Erf[(vesc/v0)]-Erf[((vmin[q,Ee,mχ]-vE)/v0)]),vmin[q,Ee,mχ]<vesc+vE&&vesc-vE<vmin[q,Ee,mχ]},{0,vmin[q,Ee,mχ]>vesc+vE}}]

FDM[q_?NumberQ,n_?NumberQ]:=(((((α me)/q)))^n)

Rate

μχe[mχ_]:=((me mχ)/((me+ mχ)));
σetest[mχ_]:=1. (10^(-37));

Si

dRdESi is equivalent to (3.13)

(* (Subscript[f, crystal])=(...) Eq. (4.4) *)

dRdESi[mχ_,EE_?NumberQ,n_?NumberQ]:=With[{qunit=0.02 α me(* eV *),Eunit=0.1(* eV *),Eprefactor=2.0 (* eV *),wk=(2/137),MTarget=1 (*kg*),MCell=2*28.0855 amu2kg},  (c^2)sec2year ((ρχ)/(mχ))(MTarget/MCell) σetest[mχ] ((α (me^2))/((μχe[mχ])^2)) Sum[Eunit qunit((1 )/((qi*qunit)))η[qi*qunit,EE  ,mχ]((FDM[qi*qunit,n])^2)(Eprefactor (1/(qunit ))(1/Eunit)(wk/2)(wk/2)(1/wk)FFarraySi[[qi,Floor[EE/Eunit]]]),{qi,1,900}]]

RateSi[mχ_ (* MeV *),Emin_,n_]:=With[{Eunit=0.1},Sum[ dRdESi[mχ (10^6),EE,n],{EE,Emin,500 Eunit,Eunit}]]

plot (dR/dE)

With[{binsize=1.,Egap=1.1,Eunit=0.1,mχ=10. (10^6) (* eV *),iiFDM=0 (* FDM=1 *)},
tmptable=Table[{EE,dRdESi[mχ,EE,iiFDM]},{EE,0,50,binsize}];myf2SiPlot=ListLogPlot[tmptable,Joined->True,PlotRange->{{0,50},All},InterpolationOrder->0,PlotStyle->Black,Frame->True,FrameLabel->{"(Subscript[E, e]) [eV]","dR/dE [#/kg-year-eV]","(Subscript[(OverBar[σ]), e])= (10^(-37)) [(cm^2)], (Subscript[F, DM])(q)=1, mχ="<>ToString[mχ/(10^6)]<>" MeV"}]]

With[{binsize=1.,Egap=1.1,Eunit=0.1,mχ=1. (10^9) (* eV *),iiFDM=0 (* FDM=1 *)},
tmptable=Table[{EE,dRdESi[mχ,EE,iiFDM]},{EE,0,50,binsize}];myf2SiPlot=ListLogPlot[tmptable,Joined->True,PlotRange->{{0,50},All},InterpolationOrder->0,PlotStyle->Black,Frame->True,FrameLabel->{"(Subscript[E, e]) [eV]","dR/dE [#/kg-year-eV]","(Subscript[(OverBar[σ]), e])= (10^(-37)) [(cm^2)], (Subscript[F, DM])(q)=1, mχ="<>ToString[mχ/(10^6)]<>" MeV"}]]

calculate rate and (Subscript[(OverBar[σ]), e]) (* n.b. each mass takes ~100s to run for (Subscript[E, min])=(Subscript[E, gap])*)

nevents=3.6; (* events observed *)
years=1; (*exposure in years. 365.25 days = 1 year*)
kg=1.0; (* detector size in kg *)
exposure=years*kg;

With[{mX=10.,Emin=1.1,n=0},
tmp=RateSi[mX ,Emin,n];
Print["Rate = ",tmp," events/kg/year"];
Print["For (Subscript[m, χ]) = ",mX," MeV, (Subscript[(OverBar[σ]), e])=",(nevents/(tmp*exposure))σetest[mχ]," (cm^2)"] ]

With[{mX=1000.,Emin=1.1,n=0},
tmp=RateSi[mX ,Emin,n];
Print["Rate = ",tmp," events/kg/year"];
Print["For (Subscript[m, χ]) = ",mX," MeV, (Subscript[(OverBar[σ]), e])=",(nevents/(tmp*exposure))σetest[mχ]," (cm^2)"] ]

tmpSiTable=With[{Emin=1.1,n=0},Table[{mX,(nevents/(RateSi[mX ,Emin,n]*exposure))(10^(-37))},{mX,{0.1,0.5,1.0,2.0,3.0,5.0,7.0,10.,20.,50.,1000.}}]];

ListLogLogPlot[tmpSiTable,PlotStyle->Blend[{Blue,Gray}],FrameLabel->{"(Subscript[m, χ]) [MeV]","(Subscript[(OverBar[σ]), e])[\!(cm^2)]"}]

Ge

dRdEGe is equivalent to (3.13)

(* (Subscript[f, crystal])=(...) Eq. (4.4) *)

dRdEGe[mχ_,EE_?NumberQ,n_?NumberQ]:=With[{qunit=0.02 α me(* eV *),Eunit=0.1(* eV *),Eprefactor=1.8 (* eV *),wk=(2/137),MTarget=1 (*kg*),MCell=2*72.64 amu2kg},  (c^2)sec2year ((ρχ)/(mχ))(MTarget/MCell) σetest[mχ] ((α (me^2))/((μχe[mχ])^2)) Sum[Eunit qunit((1 )/((qi*qunit)))η[qi*qunit,EE  ,mχ]((FDM[qi*qunit,n])^2)(Eprefactor (1/(qunit ))(1/Eunit)(wk/2)(wk/2)(1/wk)FFarrayGe[[qi,Floor[EE/Eunit]]]),{qi,1,900}]]

RateGe[mχ_ (* MeV *),Emin_,n_]:=With[{Eunit=0.1},Sum[ dRdEGe[mχ (10^6),EE,n],{EE,Emin,500 Eunit,Eunit}]]

plot (dR/dE)

With[{binsize=1.,mχ=10. (10^6) (* eV *),iiFDM=0 (* FDM=1 *)},
tmptable=Table[{EE,dRdEGe[mχ,EE,iiFDM]},{EE,0,50,binsize}];myf2GePlot=ListLogPlot[tmptable,Joined->True,PlotRange->{{0,50},All},InterpolationOrder->0,PlotStyle->Black,Frame->True,FrameLabel->{"(Subscript[E, e]) [eV]","dR/dE [#/kg-year-eV]","(Subscript[(OverBar[σ]), e])= (10^(-37)) [(cm^2)], (Subscript[F, DM])(q)=1, mχ="<>ToString[mχ/(10^6)]<>" MeV"}]]

With[{binsize=1.,mχ=1. (10^9) (* eV *),iiFDM=0 (* FDM=1 *)},
tmptable=Table[{EE, dRdEGe[mχ,EE,iiFDM]},{EE,0,50,binsize}];myf2GePlot=ListLogPlot[tmptable,Joined->True,PlotRange->{{0,50},All},InterpolationOrder->0,PlotStyle->Black,Frame->True,FrameLabel->{"(Subscript[E, e]) [eV]","dR/dE [#/kg-year-eV]","(Subscript[(OverBar[σ]), e])= (10^(-37)) [(cm^2)], (Subscript[F, DM])(q)=1, mχ="<>ToString[mχ/(10^6)]<>" MeV"}]]

calculate rate and (Subscript[(OverBar[σ]), e]) (* n.b. each mass takes ~100s to run for (Subscript[E, min])=(Subscript[E, gap])*)

nevents=3.0; (* events observed *)
years=1; (*exposure in years. 365.25 days = 1 year*)
kg=1.0; (* detector size in kg *)
exposure=years*kg;

With[{mX=10.,Emin=0.67,n=0},
tmp=RateGe[mX ,Emin,n];
Print["Rate = ",tmp," events/kg/year"];
Print["For (Subscript[m, χ]) = ",mX," MeV, (Subscript[(OverBar[σ]), e])=",(nevents/(tmp *exposure))σetest[mχ]," (cm^2)"]]

With[{mX=1000.,Emin=0.67,n=0},
tmp=RateGe[mX ,Emin,n];
Print["Rate = ",tmp," events/kg/year"];
Print["For (Subscript[m, χ]) = ",mX," MeV, (Subscript[(OverBar[σ]), e])=",(nevents/(tmp*exposure))σetest[mχ]," (cm^2)"]
]

tmpGeTable=With[{Emin=0.67,n=0},Table[{mX,(nevents/(RateGe[mX ,Emin,n]*exposure))σetest[mX]},{mX,{0.1,0.5,1.0,2.0,3.0,5.0,7.0,10.,20.,50.,1000.}}]];

ListLogLogPlot[tmpGeTable,PlotStyle->Blend[{Blue,Gray}],FrameLabel->{"(Subscript[m, χ]) [MeV]","(Subscript[(OverBar[σ]), e])[\!(cm^2)]"}]