# Incorporating matter effects on active-sterile neutrino oscillation with mathematica

I recently wanted to study the effect of a constant matter potential on a simple 2-flavor active-sterile neutrino oscillation system using Wolfram Mathematica, but I encountered a big problem in the result and I don't know what I did wrong. Maybe anyone of you also encountered a problem like this.

First, I defined the $$2 \times 2$$ mixing matrix with the mixing angle $$\Theta$$ and wrote down the vacuum Hamiltonian in flavor space:

U2[\[CapitalTheta]_] := {{Cos[\[CapitalTheta]],
Sin[\[CapitalTheta]]}, {-Sin[\[CapitalTheta]],
Cos[\[CapitalTheta]]}};
Mass2nu[dm_] := {{0, 0}, {0, dm}};
Hvac2nu[En_, \[CapitalTheta]_, dm_] :=
1/(2 En)*U2[\[CapitalTheta]].Mass2nu[dm].ConjugateTranspose[
U2[\[CapitalTheta]]];


Then I defined the different matter potentials (for derivation see https://arxiv.org/abs/hep-ph/0504264) as functions of the electron, neutron and proton densities $$N_e, N_n, N_p$$.

GF = 1.1664 * 10^-23;(*in eV^-2, because Subscript[G, F] = 1.1664 * 10^-5 GeV^-2*)
\[CapitalTheta]w = ArcSin[Sqrt[0.231]];
VMatter2nuCC[Nn_, Np_, Ne_] := {{Sqrt GF Ne, 0}, {0, 0}};
VMatter2nuNC[Nn_, Np_, Ne_] := {{-(GF (1 - 4 Sin[\[CapitalTheta]w]^2) (Ne - Np))/
Sqrt - (GF Nn)/Sqrt,
0}, {0, -(GF (1 - 4 Sin[\[CapitalTheta]w]^2) (Ne - Np))/
Sqrt - (GF Nn)/Sqrt}};
VMatter2nu[Nn_, Np_, Ne_] :=
VMatter2nuCC[Nn, Np, Ne] + VMatter2nuNC[Nn, Np, Ne];


,where $$\Theta_W$$ is the Weinberg angle and $$G_F$$ is the Fermi constant. Now I add the whole matter potential to the vacuum Hamiltonian

HMatter2nu[En_, \[CapitalTheta]_, dm_, Nn_, Np_, Ne_] :=
Hvac2nu[En, \[CapitalTheta], dm] + VMatter2nu[Nn, Np, Ne];


and define the oscillation probabilities as matrix exponentials of that Hamiltonian:

Pas2nuVac[L_, En_, \[CapitalTheta]_, dm_] :=
Abs[MatrixExp[-I*Hvac2nu[En, \[CapitalTheta], dm]*L][[2, 1]]]^2;
Paa2nuVac[L_, En_, \[CapitalTheta]_, dm_] :=
Abs[MatrixExp[-I*Hvac2nu[En, \[CapitalTheta], dm]*L][[1, 1]]]^2;
Pss2nuVac[L_, En_, \[CapitalTheta]_, dm_] :=
Abs[MatrixExp[-I*Hvac2nu[En, \[CapitalTheta], dm]*L][[2, 2]]]^2;
Pas2nuMatter[L_, En_, \[CapitalTheta]_, dm_, Nn_, Np_, Ne_] :=
Abs[MatrixExp[-I*HMatter2nu[En, \[CapitalTheta], dm, Nn, Np, Ne]*
L][[2, 1]]]^2;
Paa2nuMatter[L_, En_, \[CapitalTheta]_, dm_, Nn_, Np_, Ne_] :=
Abs[MatrixExp[-I*HMatter2nu[En, \[CapitalTheta], dm, Nn, Np, Ne]*
L][[1, 1]]]^2;
Pss2nuMatter[L_, En_, \[CapitalTheta]_, dm_, Nn_, Np_, Ne_] :=
Abs[MatrixExp[-I*HMatter2nu[En, \[CapitalTheta], dm, Nn, Np, Ne]*
L][[2, 2]]]^2;


Now we have everything to make plots. I chose the following test-variables:

(*For the Earth's core*)
ZAIron = 0.466;
DensityCore = 11* 10^3*  4.5*10^15(* (kg/(m^3))~4.5*10^15 eV^4 *);
NeCore = ZAIron *DensityCore;
NpCore = NeCore;
NnCore = NpCore;
EZero = 10^-1;
EMaximum = 10^10;
LTest = 10^2;
\[CapitalTheta]2nu = ArcSin[Sqrt[0.1]];
dm2nu = 1.2;


and used this code for plotting

LogLinearPlot[{Pas2nuVac[LTest, En, \[CapitalTheta]2nu, dm2nu],
Paa2nuVac[LTest, En, \[CapitalTheta]2nu, dm2nu],
Pas2nuVac[LTest, En, \[CapitalTheta]2nu, dm2nu] +
Paa2nuVac[LTest, En, \[CapitalTheta]2nu, dm2nu]}, {En, EZero,
EMaximum},
PlotLegends -> {"\!$$\*SubscriptBox[\(P$$, $$as$$]\)",
"\!$$\*SubscriptBox[\(P$$, $$aa$$]\)",
"\!$$\*SubscriptBox[\(\[CapitalSigma]$$, \
$$i$$]\)\!$$\*SubscriptBox[\(P$$, $$ai$$]\)"}, PlotPoints -> 500,
PlotLabel ->
StringForm["active/sterile oscillation probabilities in vacuum"],
PlotRange -> {{EZero, EMaximum}, {0, 1}}]
LogLinearPlot[{Pas2nuMatter[LTest, En, \[CapitalTheta]2nu, dm2nu,
NeCore, NeCore, NeCore],
Paa2nuMatter[LTest, En, dm2nu, \[CapitalTheta]2nu, NeCore, NeCore,
NeCore],
Pas2nuMatter[LTest, En, \[CapitalTheta]2nu, dm2nu, NeCore, NeCore,
NeCore] +
Paa2nuMatter[LTest, En, dm2nu, \[CapitalTheta]2nu, NeCore, NeCore,
NeCore]}, {En, EZero, EMaximum},
PlotLegends -> {"\!$$\*SubscriptBox[\(P$$, $$as$$]\)",
"\!$$\*SubscriptBox[\(P$$, $$aa$$]\)",
"\!$$\*SubscriptBox[\(\[CapitalSigma]$$, \
$$i$$]\)\!$$\*SubscriptBox[\(P$$, $$ai$$]\)"}, PlotPoints -> 500,
PlotLabel ->
StringForm["active/sterile oscillation probabilities in matter"],
PlotRange -> {{EZero, EMaximum}, {0, 1}}]


My program then spits out these plots

You can see, that in the vacuum case everything is as expected, whereas in the matter case something goes horribly wrong, since probability is not conserved here (as you can see in the yellow graph). Can anybody tell me what I did wrong here?

• From physics.stackexchange.com/help/on-topic: “Some kinds of questions should not be asked here: Implementation details of computational tasks While computational physics is on topic, we are not a programming site. If your question is about implementing computational code - in particular, if it's about writing, compiling, debugging or optimizing code, or about a specific language or library - then it is off topic. It may be suitable for Computational Science or Stack Overflow, however. – G. Smith Sep 3 '19 at 15:47
• @G.Smith my question is not really about optimization or the programming itself, but about the implementation of the physics problem. I don't think Comp. Science or Stack Overflow would be suitable sites for this question. – DomDoe Sep 4 '19 at 9:22