Simplified  Models for  Photohadronic Interactions in cosmic accelerators

Authors : S. Hümmer, M. Rüger, F. Spanier, W. Winter

This notebook illustrates how to reproduce some results for model Sim - B of arXiv : 1002.1310. It calculates pion, kaon, neutron muon and neutrino spectra.  Please cite arXiv: 1002.1310 if you use this notebook. Requires Mathematica 7.

General parameters

General note on units : All energies in this notebook are in GeV. All spectra are in units of (GeV^(-1)) (cm^(-3)) (steady) or (GeV^(-1)) (cm^(-3)) (s^(-1)) (injection) times the normalization factors, as given in App. C.

Input Proton spectrum as a function of energy (Eq. C2)

np[ep_]:=If[ep>1,(ep)^-2 Exp[-(ep/(6.9 10^8))^2],0]

Input Photon spectrum as a function of energy (Eq. C1)

ng[eg_]:=If[2 10^-10<eg<10^-6,(eg)^(-1),If[10^-6<eg<300 10^-6,10^-6 (eg)^-2,0]]

Energy range for which spectra are to be calculated

This notebook needs about 10 minutes to run through for the considered energy range chosen to be [1 ...10^11] GeV (logarithmically). The distance between the interpolation points is in this case deltae = 0.5 [GeV]. For smoother/better results, adjust this number to 0.1 (affects the computation time by a factor of five).

emin=0;emax=11;deltae=0.5;
energy=Table[{i},{i,emin,emax,deltae}];

Switches off the messages regarding slow converging and undeflows in computations.  For our accuracy, this treatment is sufficient.

Off[NIntegrate::slwcon, NIntegrate::ncvb, General::unfl]

Description of photohadronics

Here the response Functions (f) and the parameters are defined for the different interaction types.

Resonant production (see Sec. 4.2.2)

General version of f-functions in Eqs. (31) and (32):

funM[y_,emin_,emax_,σ_]:=If[2*y<emin,0,If[2*y<emax,σ/(2*y)^2*((2*y)^2-emin^2),σ/(2*y)^2*((emax)^2-emin^2)]]

Parameters from Table 4 in the form {emin, emax, σ, χ, (Subscript[M, (π^0)]), (Subscript[M, (SuperPlus[π])]), (Subscript[M, (SuperMinus[π])]), K, (Subscript[M, n])}  for different pion charges and interaction types {(LR; (π^0)),(LR; (SuperPlus[π])),(HR; (π^0)),(HR; (SuperPlus[π])),(HR; (SuperMinus[π]))}

RList={
{0.2,0.5,200,0.22,2/3,0,0,0.22, 1/3},
{0.2,0.5,200,0.22,0,1/3,0,0.22,0},
{0.5, 1.2, 90, 0.26, 0.47,0,0,0.39,0.43},
{0.5, 1.2, 90, 0.25, 0, 0.77, 0,0.39,0},
{0.5, 1.2, 90, 0.22, 0,0,0.34, 0.39,0}
};

Direct production (see Sec. 4.3)

Functions (I^T1) and (I^T2) from Eqs. (34) and (35):

IT2[y_]:=If[2 y<0.4,0,-3.4083+16.2864/(2 y)+40.7160 Log[2 y]];
IT1[y_]:=If[2 y<0.17,0,If[2 y<0.96,35.9533+84.0859 Log[10,y]+110.765 (Log[10,y])^2+102.728 (Log[10,y])^3+40.4699 (Log[10,y])^4,30.2004+40.5478 Log[10,y]+2.03074 (Log[10,y])^2-0.387884 (Log[10,y])^3+0.025044 (Log[10,y])^4 ]]

f-functions for single pion production (1) and 2-pion production (2) from Eq. (33):

funD1[y_,emin_,emax_]:=If[2*y<emin,0,If[2*y<emax,1/(2*y^2)*(IT1[y]-IT1[emin/2]),1/(2*y^2)*(IT1[emax/2]-IT1[emin/2])]];
funD2[y_,emin_,emax_]:=If[2*y<emin,0,If[2*y<emax,1/(2*y^2)*(IT2[y]-IT2[emin/2]),1/(2*y^2)*(IT2[emax/2]-IT2[emin/2])]];

Parameters from Table 5 {emin, emax,  χ, (Subscript[M, (π^0)]), (Subscript[M, (SuperPlus[π])]), (Subscript[M, (SuperMinus[π])]), K, (Subscript[M, n]), Z} for the different  interaction types {(Subscript[T, 1L]),(Subscript[T, 1M]),(Subscript[T, 1H]),(Subscript[T, 2aL]),(Subscript[T, 2aM]),(Subscript[T, 2aH]),(Subscript[T, 2b])};
here  Z=1 for single pion production and Z=2 for 2-pion production

DList={
{0.17, 0.56, 0.13, 0, 1,0, 0.13,1,1},
{0.56, 10, 0.05, 0,1, 0, 0.05,1,1},
{10, Infinity, 0.001, 0, 1, 0, 0.001,1,1},
{0.4,1.58,0.08, 0, 0.25,0.75, 0.28,1/6,2},
{1.58,100,0.02, 0, 0.25,0.75, 0.22, 1/6,2},
{100,Infinity,0.001, 0, 0.25,0.75, 0.201,1/6,2},
{0.4,Infinity,0.2, 1/6, 0.75,1/12, 0.2,0,2}
};

Multi-pion production (see Sec. 4.4.2)

f - function (see Eq. (40)) has the same general form as the one  for resonant production

Parameters from Table 6 {emin, emax, σ, χ, (Subscript[M, (π^0)]), (Subscript[M, (SuperPlus[π])]), (Subscript[M, (SuperMinus[π])]), K, (Subscript[M, n])} for the different interaction types (see first column of Table 6):

MList={
{0.5,0.9,60,0.1,0.32,0.34,0.04,0.27,0.31},
{0.5,0.9,60,0.4,0.17,0.29,0.05,0.27,0},
{0.9,1.5,85,0.15,0.42,0.31,0.07,0.34,0.31},
{0.9,1.5,85,0.35,0.19,0.35,0.08,0.34,0},
{1.5,5.0,120,0.15,0.59,0.57,0.3,0.39,0.31},
{1.5,5.0,120,0.35,0.16,0.21,0.13,0.39,0},
{5.0,50,120,0.07,1.38,1.37,1.11,0.49,0.31},
{5.0,50,120,0.35,0.16,0.25,0.23,0.49,0},
{50,500,120,0.02,3.01,2.86,2.64,0.45,0.31},
{50,500,120,0.5,0.2,0.21,0.14,0.45,0},
{500,5000,120,0.007,5.13,4.68,4.57,0.44,0.31},
{500,5000,120,0.5,0.27,0.29,0.12,0.44,0},
{5000,Infinity,120,0.002,7.59,6.8,6.65,0.44,0.31},
{5000,Infinity,120,0.6,0.26,0.27,0.13,0.44,0}
};

Kaon production (see Sec. 4.5)

f - function from Eq. (42):

funK[y_]:=If[2*y<1,0,If[2*y<1.2,2 (1-1/(2 y)^2),If[2*y<1.65,3.7 (1-0.3/(2 y)^2),2.7 (1-0.16/(2 y)^2)]]]

Computation of secondary spectra

Pion production; computed with Eq. (28) for the different interaction types; uses logarithmic integration and y as logarithmic variable

Pion spectra from resonant production

pi=Map[{#[[1]],NIntegrate[Map[3 10^-20 mp 10^(2 e)/(10^e) np[10^e/#[[4]]] Log[10] 10^y ng[10^y mp #[[4]]/10^e]   funM[10^y,#[[1]],#[[2]],#[[3]]]&,RList] /.e -> #[[1]] /. mp -> 0.938,{y,-1,10}]}&,energy];

Qpizerores=Map[{#[[1]],Total[#[[2]]*Map[#[[5]]&,RList]]}&,pi];
Qpiplusres=Map[{#[[1]],Total[#[[2]]*Map[#[[6]]&,RList]]}&,pi];
Qpiminusres=Map[{#[[1]],Total[#[[2]]*Map[#[[7]]&,RList]]}&,pi];

Pion spectra from direct production

pi=Map[{#[[1]],NIntegrate[Map[If[#[[9]]==1,3 10^-20 mp 10^(2 e)/(10^e) np[10^e/#[[3]]] Log[10] 10^y ng[10^y mp #[[3]]/10^e]   funD1[10^y,#[[1]],#[[2]]],3 10^-20 mp 10^(2 e)/(10^e) np[10^e/#[[3]]] Log[10] 10^y ng[10^y mp #[[3]]/10^e]   funD2[10^y,#[[1]],#[[2]]]]&,DList] /.e -> #[[1]] /. mp -> 0.938,{y,-2,10}]}&,energy];

Qpizerodir=Map[{#[[1]],Total[#[[2]]*Map[#[[4]]&,DList]]}&,pi];
Qpiplusdir=Map[{#[[1]],Total[#[[2]]*Map[#[[5]]&,DList]]}&,pi];
Qpiminusdir=Map[{#[[1]],Total[#[[2]]*Map[#[[6]]&,DList]]}&,pi];

Pion spectra from multi-pion production

pi=Map[{#[[1]],NIntegrate[Map[3 10^-20 mp 10^(2 e)/(10^e) np[10^e/#[[4]]] Log[10] 10^y ng[10^y mp #[[4]]/10^e]   funM[10^y,#[[1]],#[[2]],#[[3]]]&,MList] /.e -> #[[1]] /. mp -> 0.938,{y,-1,10}]}&,energy];

Qpizeromul=Map[{#[[1]],Total[#[[2]]*Map[#[[5]]&,MList]]}&,pi];
Qpiplusmul=Map[{#[[1]],Total[#[[2]]*Map[#[[6]]&,MList]]}&,pi];
Qpiminusmul=Map[{#[[1]],Total[#[[2]]*Map[#[[7]]&,MList]]}&,pi];

Total pion spectra

Qpizero=Table[{Qpizerores[[i,1]],Qpizerores[[i,2]]+Qpizerodir[[i,2]]+Qpizeromul[[i,2]]},{i,1,Length[Qpizerores]}];
Qpiplus=Table[{Qpiplusres[[i,1]],Qpiplusres[[i,2]]+Qpiplusdir[[i,2]]+Qpiplusmul[[i,2]]},{i,1,Length[Qpiplusres]}];
Qpiminus=Table[{Qpiminusres[[i,1]],Qpiminusres[[i,2]]+Qpiminusdir[[i,2]]+Qpiminusmul[[i,2]]},{i,1,Length[Qpiminusres]}];

Kaon production; see Eq. (28)

Qkaon=Map[{#[[1]],NIntegrate[3 10^-20 mp 10^(2 e)/(10^e) np[10^e/0.35] Log[10] 10^y ng[10^y mp 0.35/10^e] funK[10^y] /.e -> #[[1]] /. mp -> 0.938,{y,-1,10}]}&,energy];

Neutron production; see Eq. (B6)

Resonant production:

nr=Map[{#[[1]],Total[NIntegrate[Map[#[[9]] 3 10^-20 mp 10^(2 e)/(10^e) np[10^e/(1-#[[8]])] Log[10] 10^y ng[10^y mp (1-#[[8]])/10^e]   funM[10^y,#[[1]],#[[2]],#[[3]]]&,RList] /.e -> #[[1]] /. mp -> 0.938,{y,-1,10}]]}&,energy];

Direct production:

nd=Map[{#[[1]],Total[NIntegrate[Map[If[#[[9]]==1,#[[8]] 3 10^-20 mp 10^(2 e)/(10^e) np[10^e/(1-#[[7]])] Log[10] 10^y ng[10^y mp (1-#[[7]])/10^e] funD1[10^y,#[[1]],#[[2]]],#[[8]] 3 10^-20 mp 10^(2 e)/(10^e) np[10^e/(1-#[[7]])] Log[10] 10^y ng[10^y mp (1-#[[7]])/10^e] funD2[10^y,#[[1]],#[[2]]]]&,DList] /.e -> #[[1]] /. mp -> 0.938,{y,-1,10}]]}&,energy];

Multi-pion production:

nm=Map[{#[[1]],Total[NIntegrate[Map[#[[9]] 3 10^-20 mp 10^(2 e)/(10^e) np[10^e/(1-#[[8]])] Log[10] 10^y ng[10^y mp (1-#[[8]])/10^e]   funM[10^y,#[[1]],#[[2]],#[[3]]]&,MList] /.e -> #[[1]] /. mp -> 0.938,{y,-1,10}]]}&,energy];

Total neutron spectra :

Qneutron=Table[{nr[[i,1]],nr[[i,2]]+nd[[i,2]]+nm[[i,2]]},{i,1,Length[nr]}];

Figures for secondary spectra

Function constructing  interpolating functions (Parameters: data: table for which to define interpolating function,  xmin, xmax and ymax: to define ranges for the plots, func: interpolating function object to be returned)

part1[data_,xmin_,xmax_,ymax_,func_]:=
Module[{bdata,ndata,null,funcd},
ymax=Log[10,Max[Map[#[[2]]&,data]]];
bdata=Map[{#[[1]],If[#[[2]]<10^-160,0,10^(-2 #[[1]]) #[[2]]]}&,data]; (*Set all values below 10^-160 to zero*)
null=Position[Map[#[[2]]&,bdata],0]; (*Finds all entries with value zero*)
ndata=Delete[bdata,null]; (*Deletes entries with zero*)
If[Length[ndata]≠Length[bdata],AppendTo[ndata,{ndata[[Length[ndata],1]]+0.1,10^-160}]]; (*Adds value 10^-160 to the end*)
xmin=ndata[[1,1]]; (*Defines minimal value for Log[10,e]*)
xmax=ndata[[Length[ndata],1]];(*Defines maximal value for Log[10,e]*)
funcd=Interpolation[Map[{#[[1]],Log[10,#[[2]]]}&,ndata],InterpolationOrder -> 1]; (*Logarithmic interpolation*)
func[x_]:=If[Log[10,x]<xmin||Log[10,x]>xmax,0,10^funcd[Log[10,x]]]] (*Spectrum (not logarithmic in energy and function value)*)

Calculation of the interpoalting functions for pion, neutron and kaon spectra

part1[Qpiplusres,xminres,xmaxres,ymaxres,fQpiplusres];
part1[Qpiplusdir,xmindir,xmaxdir,ymaxdir,fQpiplusdir];
part1[Qpiplusmul,xminmul,xmaxmul,ymaxmul,fQpiplusmul];
part1[Qpiplus,xminpip,xmaxpip,ymaxpip,fQpiplus];
part1[Qpiminus,xminpim,xmaxpim,ymaxpim,fQpiminus];
part1[Qpizero,xminpin,xmaxpin,ymaxpin,fQpizero];
part1[Qneutron,xminn,xmaxn,ymaxn,fQneutron];
part1[Qkaon,xmink,xmaxk,ymaxk,fQkaon];

Needs["PlotLegends`"]

Contributions to the positive charged pion spectrum from resonant, direct and multi pion interactions (see Fig.5 top left corner)

ShowLegend[LogLogPlot[{x^2 fQpiplusres[x],x^2 fQpiplusdir[x],x^2 fQpiplusmul[x],x^2 fQpiplus[x]},{x,10^xminpip,10^xmaxpip}, PlotRange -> {10^(ymaxpip-5),10^(ymaxpip+1)},PlotStyle -> {Directive[RGBColor[1,0.2,0.2],AbsoluteThickness[1]],Directive[RGBColor[0.2,0.2,1], Dotted, AbsoluteThickness[1.5]],Directive[RGBColor[0.2,1,0.2],Dashed, AbsoluteThickness[1.5]],Directive[Black,AbsoluteThickness[1.5]]}, PlotLabel -> Style["(SuperPlus[π])",30],Frame -> True, AspectRatio -> 1,FrameTicksStyle-> 15,FrameLabel->{Style["(Subscript[E, π])/GeV",15],Style["(Subscript[E, π]^2) (Subscript[Q, π])/((Subscript[N, γ]) (Subscript[N, p]) GeV (cm^(-3)) (s^(-1)))",15]},Axes -> False, ImageSize -> 400,Frame -> True],{{{Graphics[{RGBColor[1,0.2,0.2],AbsoluteThickness[1],Line[{{0,0},{2,0}}]}],Style["Resonant ",15]},{Graphics[{RGBColor[0.2,0.2,1],Dotted, AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["Direct",15]},{Graphics[{RGBColor[0.2,1,0.2],Dashed, AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["Multi-π",15]},{Graphics[{Black,AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["Combined",15]}},LegendShadow -> None, LegendPosition -> {-0.2,-0.6},LegendSize -> {0.55,0.5}, LegendBorder -> None, LegendTextSpace -> 3.5}]

Pion and kaon spectra (see Fig. 6 left panel)

ShowLegend[LogLogPlot[{x^2 fQpiplus[x],x^2 fQpiminus[x],x^2 fQpizero[x],x^2 fQkaon[x]},{x,Min[10^xminpip,10^xminpin,10^xminpim,10^xmink],Max[10^xmaxpip,10^xmaxpim,10^xmaxpin,10^xmaxk]}, PlotRange -> {10^(Max[ymaxpip,ymaxpim,ymaxpin,ymaxk]-5),10^(Max[ymaxpip,ymaxpim,ymaxpin,ymaxk]+1)},PlotStyle -> {Directive[RGBColor[1,0.2,0.2],AbsoluteThickness[1]],Directive[RGBColor[0.2,0.2,1], Dotted,AbsoluteThickness[1.5]],Directive[RGBColor[0.2,1,0.2],Dashed,AbsoluteThickness[1.5]],Directive[RGBColor[0.6,0.6,0.6],AbsoluteThickness[1.5]]}, PlotLabel -> Style["π, (SuperPlus[K])",30],Frame -> True, AspectRatio -> 1,FrameTicksStyle-> 15,FrameLabel->{Style["(Subscript[E, π])/GeV",15],Style["(Subscript[E, π]^2) (Subscript[Q, π])/((Subscript[N, γ]) (Subscript[N, p]) GeV (cm^(-3)) (s^(-1)))",15]},Axes -> False, ImageSize -> 400],{{{Graphics[{RGBColor[1,0.2,0.2],AbsoluteThickness[1],Line[{{0,0},{2,0}}]}],Style["(SuperPlus[π])",15]},{Graphics[{RGBColor[0.2,0.2,1],Dotted,AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["(SuperMinus[π])",15]},{Graphics[{RGBColor[0.2,1,0.2],Dashed,AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["(π^0)",15]},{Graphics[{RGBColor[0.6,0.6,0.6],AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["(SuperPlus[K])",15]}},LegendShadow -> None, LegendPosition -> {0.05,-0.6},LegendSize -> {0.5,0.5}, LegendBorder -> None, LegendTextSpace -> 2}]

Description of Weak Decays

Scaling functions for the decays of pions and muons

famr[x_,r_]:=r (1-x)/((1-r)^2 x) (*see Eq.(46)*)
faml[x_,r_]:=(x-r)/((1-r)^2 x) (*see Eq.(47)*)
fpmn[r_]:=1/(1-r) (*see Eq.(14) Lipari et al.*)
fmn[x_,h_]:=5/3-3 x^2+4 x^3/3-h (1/3-3 x^2+8 x^3/3)(*see Eq.(48)*)
fen[x_,h_]:=2-6 x^2+4 x^3+h (2-12 x+18 x^2-8 x^3)(*see Eq.(49)*)
mpi=0.14;
mmu=0.105;
mk=0.494;

Pion decays - Muon spectra

Here we compute the decays of the pions into left- and right-handed muons and construct the interpolating functions for the muon spectra (for explanation see module part1)

scale[parent_,scalefunc_,rat_,xmin_,xmax_,ymax_,func_]:=
Module[{data,qamr,bdata,ndata,null,funcd},
qamr[emu_,x_]=parent[emu/x] 1/x scalefunc[x,rat]  ;
data=Table[{emu,NIntegrate[10^(2 emu) qamr[10^emu,x],{x,rat,1}]},{emu,emin,emax,deltae}]; (*see Eq.(2)*)
ymax=Log[10,Max[Map[#[[2]]&,data]]];
bdata=Map[{#[[1]],If[#[[2]]<10^-160,0,10^(-2 #[[1]]) #[[2]]]}&,data];
null=Position[Map[#[[2]]&,bdata],0];
ndata=Delete[bdata,null];
If[Length[ndata]≠Length[bdata],AppendTo[ndata,{ndata[[Length[ndata],1]]+0.1,10^-160}]];
xmin=ndata[[1,1]];
xmax=ndata[[Length[ndata],1]];
funcd=Interpolation[Map[{#[[1]],Log[10,#[[2]]]}&,ndata],InterpolationOrder -> 1];
func[x_]:=If[Log[10,x]<xmin||Log[10,x]>xmax,0,10^funcd[Log[10,x]]]]

scale[fQpiplus,famr,(mmu/mpi)^2,xminamr,xmaxamr,ymaxamr,fQmuplusr];
scale[fQpiplus,faml,(mmu/mpi)^2,xminaml,xmaxaml,ymaxaml,fQmuplusl];
scale[fQpiminus,faml,(mmu/mpi)^2,xminmr,xmaxmr,ymaxmr,fQmuminusr];
scale[fQpiminus,famr,(mmu/mpi)^2,xminml,xmaxml,ymaxml,fQmuminusl];

Pion decays - Neutrino spectra

Here we compute the decay of the pions into neutrinos and construct the interpolating functions for the neutrino spectra (for explanation see module part1)

scale2[parent_,scalefunc_,rat_,br_,xmin_,xmax_,ymax_,func_]:=
Module[{data,qamr,bdata,ndata,null,funcd},
qamr[emu_,x_]=parent[emu/x] 1/x scalefunc[rat]  ;
data=Table[{emu,NIntegrate[10^(2 emu)br  qamr[10^emu,x],{x,0,1-rat}]},{emu,emin,emax,deltae}];
ymax=Log[10,Max[Map[#[[2]]&,data]]];
bdata=Map[{#[[1]],If[#[[2]]<10^-160,0,10^(-2 #[[1]]) #[[2]]]}&,data];
null=Position[Map[#[[2]]&,bdata],0];
ndata=Delete[bdata,null];
If[Length[ndata]≠Length[bdata],AppendTo[ndata,{ndata[[Length[ndata],1]]+0.1,10^-160}]];
xmin=ndata[[1,1]];
xmax=ndata[[Length[ndata],1]];
funcd=Interpolation[Map[{#[[1]],Log[10,#[[2]]]}&,ndata],InterpolationOrder -> 1];
func[x_]:=If[Log[10,x]<xmin||Log[10,x]>xmax,0,10^funcd[Log[10,x]]]]

scale2[fQpiplus,fpmn,(mmu/mpi)^2,1,xminpmn,xmaxpmn,ymaxpmn,fQmuneufrompi];
scale2[fQpiminus,fpmn,(mmu/mpi)^2,1,xminapmn,xmaxapmn,ymaxapmn,fQantimuneufrompi];

Kaon decays - Neutrino spectrum

Here we compute the decay of the kaons into neutrinos and construct the interpolating functions for the neutrino spectrum

scale2[fQkaon,fpmn,(mmu/mk)^2,0.635,xminkmn,xmaxkmn,ymaxkmn,fQmuneufromk]

Neutron decays - Neutrino spectrum

Here we compute the decay of the neutrons into neutrinos and construct the interpolating functions for the neutrino spectrum (for explanation see module part1)

scale3[parent_,xmin_,xmax_,ymax_,func_]:=
Module[{data,qamr,bdata,ndata,null,funcd},
data=Table[{emu,10^(2 emu) parent[10^emu/(5.1 10^-4)] 1/(5.1 10^-4)},{emu,emin,emax,deltae}];
ymax=Log[10,Max[Map[#[[2]]&,data]]];
bdata=Map[{#[[1]],If[#[[2]]<10^-160,0,10^(-2 #[[1]]) #[[2]]]}&,data];
null=Position[Map[#[[2]]&,bdata],0];
ndata=Delete[bdata,null];
If[Length[ndata]≠Length[bdata],AppendTo[ndata,{ndata[[Length[ndata],1]]+0.1,10^-160}]];
xmin=ndata[[1,1]];
xmax=ndata[[Length[ndata],1]];
funcd=Interpolation[Map[{#[[1]],Log[10,#[[2]]]}&,ndata],InterpolationOrder -> 1];
func[x_]:=If[Log[10,x]<xmin||Log[10,x]>xmax,0,10^funcd[Log[10,x]]]]

scale3[fQneutron,xminnaen,xmaxnaen,ymaxnaen,fQantieneufromn]

Muon decays - Neutrino spectra

Here we compute the decays of the left - and right - handed muons into neutrinos and construct the interpolating functions for the neutrino spectra taking into account the helicity of the muons (for explanation see module part1).

scale4[parent_,parent2_,scalefunc_,rat_,xmin_,xmax_,ymax_,func_]:=
Module[{data,qamr,qaml,bdata,ndata,null,funcd},
qamr[emu_,x_]=parent[emu/x] 1/x scalefunc[x,rat]  ;
qaml[emu_,x_]=parent2[emu/x] 1/x scalefunc[x,-rat]  ;
data=Table[{emu,NIntegrate[(10^(2 emu) qamr[10^emu,x]+10^(2 emu) qaml[10^emu,x]),{x,0,1}]},{emu,emin,emax,deltae}];
ymax=Log[10,Max[Map[#[[2]]&,data]]];
bdata=Map[{#[[1]],If[#[[2]]<10^-160,0,10^(-2 #[[1]]) #[[2]]]}&,data];
null=Position[Map[#[[2]]&,bdata],0];
ndata=Delete[bdata,null];
If[Length[ndata]≠Length[bdata],AppendTo[ndata,{ndata[[Length[ndata],1]]+0.1,10^-160}]];
xmin=ndata[[1,1]];
xmax=ndata[[Length[ndata],1]];
funcd=Interpolation[Map[{#[[1]],Log[10,#[[2]]]}&,ndata],InterpolationOrder -> 1];
func[x_]:=If[Log[10,x]<xmin||Log[10,x]>xmax,0,10^funcd[Log[10,x]]]]

Helicity of the muons taken into account:

scale4[fQmuminusl,fQmuminusr,fmn,1,xminnm,xmaxnm,ymaxnm,fQmuneufrommu];
scale4[fQmuplusl,fQmuplusr,fmn,-1,xminamn,xmaxanm,ymaxamn,fQantimuneufrommu];
scale4[fQmuminusl,fQmuminusr,fen,1,xminaen,xmaxaen,ymaxaen,fQantieneufrommu];
scale4[fQmuplusl,fQmuplusr,fen,-1,xminen,xmaxen,ymaxen,fQeneufrommu];

Neglecting helicity of the muons:

scale4[fQmuminusl,fQmuminusr,fmn,0,xminnmh,xmaxnmh,ymaxnmh,fQmuneufrommuh];
scale4[fQmuplusl,fQmuplusr,fmn,0,xminamnh,xmaxanmh,ymaxamnh,fQantimuneufrommuh];
scale4[fQmuminusl,fQmuminusr,fen,0,xminaenh,xmaxaenh,ymaxaenh,fQantieneufrommuh];
scale4[fQmuplusl,fQmuplusr,fen,0,xminenh,xmaxenh,ymaxenh,fQeneufrommuh];

Figures for neutrino spectra

Muon neutrinos (see Fig. 10 left panel)

ShowLegend[LogLogPlot[{x^2 fQmuneufrommu[x],x^2 fQmuneufrompi[x],x^2 fQmuneufromk[x],x^2 (fQmuneufrommu[x]+fQmuneufrompi[x]+fQmuneufromk[x])},{x,Min[10^xminnm,10^xminkmn],Max[10^xmaxnm,10^xmaxkmn]},PlotRange -> {(10^ymaxnm+10^ymaxpmn+10^ymaxkmn) 10^-5,10^3 (10^ymaxnm+10^ymaxpmn+10^ymaxkmn)},PlotStyle -> {Directive[RGBColor[0.2,1,0.2],Dashed,AbsoluteThickness[1.5]],Directive[RGBColor[0.2,0.2,1], Dotted,AbsoluteThickness[1.5]],Directive[RGBColor[0.7,0.7,0.7],AbsoluteThickness[1]],Directive[RGBColor[0,0,0],AbsoluteThickness[1.5]],Directive[RGBColor[1,0.2,0.2],AbsoluteThickness[1]]}, PlotLabel -> Style["(Subscript[ν, μ])",30],Frame -> True, AspectRatio -> 1,FrameTicksStyle-> 15,FrameLabel->{Style["(Subscript[E, ν])/GeV",15],Style["(Subscript[E, ν]^2) (Subscript[Q, ,ν])/((Subscript[N, γ]) (Subscript[N, p]) GeV (cm^(-3))(s^(-1)))",15]},Axes -> False, ImageSize ->400],{{{Graphics[{RGBColor[0.2,1,0.2],Dashed,AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["From μ decay",14]},{Graphics[{RGBColor[0.2,0.2,1],Dotted,AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["From π decay",14]},{Graphics[{RGBColor[0.7,0.7,0.7],AbsoluteThickness[1],Line[{{0,0},{2,0}}]}],Style["From K decay",14]},{Graphics[{RGBColor[0,0,0],AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["Combined",14]}},LegendShadow -> None, LegendPosition -> {-0.5,0.15},LegendSize -> {0.6,0.5}, LegendBorder -> None, LegendTextSpace -> 5}]

Electron neutrinos  (see Fig. 11 left panel)

ShowLegend[LogLogPlot[{x^2 fQeneufrommu[x],x^2 fQantieneufrommu[x],x^2 fQantieneufromn[x]},{x,Min[10^xminen,10^xminaen,10^xminnaen],Max[10^xmaxen,10^xmaxaen,10^xmaxnaen]},PlotRange -> {Max[10^ymaxen,10^ymaxaen] 10^-5,10^3 Max[10^ymaxen,10^ymaxaen]},PlotStyle -> {Directive[RGBColor[0,0,0],Dashed,AbsoluteThickness[1]],Directive[RGBColor[0,0,0],AbsoluteThickness[1]],Directive[RGBColor[0.2,0.2,1],Dotted,AbsoluteThickness[1.5]]}, PlotLabel -> Style["(Subscript[ν, e]),(Subscript[(OverBar[ν]), e])",30],Frame -> True, AspectRatio -> 1,FrameTicksStyle-> 15,FrameLabel->{Style["(Subscript[E, ν])/GeV",15],Style["(Subscript[E, ν]^2) (Subscript[Q, ,ν])/((Subscript[N, γ]) (Subscript[N, p]) GeV (cm^(-3))(s^(-1)))",15]},Axes -> False, ImageSize ->400],{{{Graphics[{RGBColor[0,0,0],Dashed,AbsoluteThickness[1],Line[{{0,0},{2,0}}]}],Style["(Subscript[ν, e])",15]},{Graphics[{RGBColor[0,0,0],AbsoluteThickness[1],Line[{{0,0},{2,0}}]}],Style["(Subscript[(OverBar[ν]), e]) from μ decay",15]},{Graphics[{RGBColor[0.2,0.2,1],AbsoluteThickness[1.5],Dotted,Line[{{0,0},{2,0}}]}],Style["(Subscript[(OverBar[ν]), e]) from n decay",15]}},LegendShadow -> None, LegendPosition -> {-0.5,0.2},LegendSize -> {0.75,0.4}, LegendBorder -> None, LegendTextSpace -> 5}]

Flavor Ratio  (see Fig. 12 left panel)

ShowLegend[Plot[{(fQeneufrommu[10^x]+fQantieneufrommu[10^x])/(fQmuneufrommu[10^x]+fQmuneufrompi[10^x]+fQantimuneufrommu[10^x]+fQantimuneufrompi[10^x]),(fQeneufrommu[10^x]+fQantieneufrommu[10^x]+fQantieneufromn[10^x])/(fQmuneufrommu[10^x]+fQmuneufrompi[10^x]+fQantimuneufrompi[10^x]+fQantimuneufrommu[10^x]+fQmuneufromk[10^x]),(fQeneufrommuh[10^x]+fQantieneufrommuh[10^x])/(fQmuneufrommuh[10^x]+fQmuneufrompi[10^x]+fQantimuneufrompi[10^x]+fQantimuneufrommuh[10^x]+fQmuneufromk[10^x]),0.5},{x,Min[xminen,xminnm,xminkmn],Max[xmaxen,xmaxnm,xmaxkmn]},PlotRange -> {0.3,0.7},PlotStyle -> {Directive[RGBColor[0,0,0],AbsoluteThickness[1.5]],Directive[RGBColor[0.2,1,0.2], Dashed,AbsoluteThickness[1.5]],Directive[RGBColor[0.2,0.2,1],Dotted,AbsoluteThickness[1.5]]}, PlotLabel -> Style["((Subscript[(OverBar[ν]), e])+(Subscript[ν, e]))/((Subscript[(OverBar[ν]), μ])+(Subscript[ν, μ]))",30],Frame -> True, AspectRatio -> 1,FrameTicksStyle-> 15,FrameLabel->{Style["(Subscript[E, ν])/GeV",15],Style["Fluxratio: ((Subscript[(OverBar[ν]), e])+(Subscript[ν, e]))/((Subscript[(OverBar[ν]), μ])+(Subscript[ν, μ]))",15]},Axes -> False, ImageSize ->400],{{{Graphics[{RGBColor[0,0,0],AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["From π/μ decays",15]},{Graphics[{RGBColor[0.2,1,0.2],Dashed,AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["From π/μ/n/K decays",15]},{Graphics[{RGBColor[0.2,0.2,1],Dotted,AbsoluteThickness[1.5],Line[{{0,0},{2,0}}]}],Style["From π/μ dec., h=0",15]}},LegendShadow -> None, LegendPosition -> {-0.68,-0.7},LegendSize -> {1,0.5}, LegendBorder -> None, LegendTextSpace -> 6}]