Neurthino.jl is a package for calculating neutrino oscillation probabilities. The main focus of the package lies on atmospheric neutrino flux and the neutrino propagation through earth.
First of all the basic vacuum properties have to be defined by creating a
OscillationParameters struct with fixed number of neutrino flavours of the
considered model:
julia> using Neurthino
julia> osc = OscillationParameters(3);The values of the mixing angles (setθ!), mass squared differences (setΔm²)
and CP phases (setδ!) are initialised to 0 and have to be set individually:
julia> setθ!(osc, 1=>2, 0.59);
julia> setθ!(osc, 1=>3, 0.15);
julia> setθ!(osc, 2=>3, 0.84);
julia> setδ!(osc, 1=>3, 3.86);
The mass squared differences are defined as and
within the package the convention
is kept.
julia> setΔm²!(osc, 2=>3, -2.523e-3);
julia> setΔm²!(osc, 1=>2, -7.39e-5);
These oscillation parameters can now be used to calculate the oscillation probabilities between the flavour states:
julia> p = Pνν(osc, 1, 10000)
4-dimensional AxisArray{Float64,4,...} with axes:
:Energy, [1.0]
:Baseline, [10000.0]
:InitFlav, NeutrinoFlavour[Electron, Muon, Tau]
:FinalFlav, NeutrinoFlavour[Electron, Muon, Tau]
And data, a 1×1×3×3 Array{Float64,4}:
[:, :, 1, 1] =
0.40280077905806266
[:, :, 2, 1] =
0.24823028034134093
[:, :, 3, 1] =
0.348968940600596
[:, :, 1, 2] =
0.10025499082597984
[:, :, 2, 2] =
0.49250415138072934
[:, :, 3, 2] =
0.4072408577932906
[:, :, 1, 3] =
0.49694423011595723
[:, :, 2, 3] =
0.2592655682779296
[:, :, 3, 3] =
0.24379020160611306
The output is an AxisArray which provides intuitive indexing, e.g.
for P(νμ→ντ) at the given energy and baseline:
julia> p[Energy=1, Baseline=1, InitFlav=Muon, FinalFlav=Tau]
0.2592655682779296
The probabilities are calculated based on the transition matrix
(the so-called PMNS-Matrix) between flavour and mass eigenstates,
as well as the Hamiltonian in the mass eigenbasis. In order to calculating these
just once, the Pνν function can be utilised in the following way:
julia> U = PMNSMatrix(osc)
3×3 Array{Complex{Float64},2}:
0.82161+0.0im 0.550114+0.0im -0.112505+0.0983582im
-0.301737+0.0608595im 0.601232+0.0407488im 0.736282+0.0im
0.476688+0.0545516im -0.576975+0.0365253im 0.659968+0.0im
julia> H = Hamiltonian(osc)
3-element Array{Complex{Float64},1}:
-0.0008902666666666667 + 0.0im
-0.0008163666666666667 + 0.0im
0.0017066333333333333 + 0.0im
julia> Pνν(U, H, 1, 10000)
4-dimensional AxisArray{Float64,4,...} with axes:
:Energy, [1.0]
:Baseline, [10000.0]
:InitFlav, NeutrinoFlavour[Electron, Muon, Tau]
:FinalFlav, NeutrinoFlavour[Electron, Muon, Tau]
And data, a 1×1×3×3 Array{Float64,4}:
[:, :, 1, 1] =
0.40280077905806266
[:, :, 2, 1] =
0.24823028034134093
[:, :, 3, 1] =
0.348968940600596
[:, :, 1, 2] =
0.10025499082597984
[:, :, 2, 2] =
0.49250415138072934
[:, :, 3, 2] =
0.4072408577932906
[:, :, 1, 3] =
0.49694423011595723
[:, :, 2, 3] =
0.2592655682779296
[:, :, 3, 3] =
0.24379020160611306
For homogeneous matter with a fixed density, a modified PMNS-Matrix
and Hamiltonian can be determined and passed into Pνν, just like for
oscillations in vacuum. In order to determine the modified PMNS-Matrix and
Hamiltonian the neutrino energy and the matter density are required:
julia> U_mat, H_mat = MatterOscillationMatrices(U, H, 1, 13);
julia> H_mat
3-element Array{Complex{Float64},1}:
-0.0008404901318507502 - 2.5459232191294903e-20im
9.078126149399635e-5 - 1.75151351027943e-20im
0.0017419062876598283 - 1.8741859435908039e-19im
julia> U_mat
3×3 Array{Complex{Float64},2}:
0.0358018-0.000158113im 0.970863+0.0im -0.178275+0.156083im
-0.662778+0.00661213im 0.157174+0.116074im 0.722845+0.0im
0.74793+0.0im 0.0917808+0.104043im 0.649115-0.00104331im
The oscillation probabilities using the Pνν function, as described above:
julia> Pνν(U_mat, H_mat, 1, 10000)
4-dimensional AxisArray{Float64,4,...} with axes:
:Energy, [1]
:Baseline, [10000]
:InitFlav, NeutrinoFlavour[Electron, Muon, Tau]
:FinalFlav, NeutrinoFlavour[Electron, Muon, Tau]
And data, a 1×1×3×3 Array{Float64,4}:
[:, :, 1, 1] =
0.8340722296308641
[:, :, 2, 1] =
0.08290502782120308
[:, :, 3, 1] =
0.08302274254793415
[:, :, 1, 2] =
0.10825570726818898
[:, :, 2, 2] =
0.052976635020068
[:, :, 3, 2] =
0.8387676577117485
[:, :, 1, 3] =
0.05767206310094823
[:, :, 2, 3] =
0.8641183371587345
[:, :, 3, 3] =
0.07820959974032213
The second option is suitable for scenarios with more complex paths with sections of different densities. An example is shown in the next chapter, where we propagate neutrinos through the earth.
The Neurthino.jl package also includes features for the neutrino oscillation probabilities
through the Earth, i.e. it contains functions for generating a neutrino path based on the
PREM model. In the following example a neutrino oscillogram with a resolution of 200x200 bins
is determined. The zenith angles for up going neutrinos (cos(θ)ϵ[-1,0]) and
subsequently the neutrino paths are generated first:
julia> zenith = acos.(range(-1,stop=0,length=200));
julia> paths = Neurthino.prempath(zenith, 2.5, samples=100, discrete_densities=0:0.1:14);
The detector is assumed to be 2.5km under the earth's surface (a typical KM3NeT
detector block in the Mediterranean), which is a realistic scenario for
Water-Cherenkov-Detectors in sea or ice. Each path consists of 100 sections of
equal lengths while the matter density is taken from the PREM model.
If a vector of densities is passed as discrete_densities, the values are
clipped to the closest value.
julia> energies = 10 .^ range(0, stop=2, length=200);
julia> prob = Pνν(U, H, energies, paths);
The returned array prob is again of type AxisArray with an axis Path for the path index (instead of the Baseline axis).
P(νe→νe) is determined by prob[InitFlav=Electron, FinalFlav=Electron], which can be visualised by a heatmap:
and for P(νμ→νμ) or prob[InitFlav=Muon, FinalFlav=Muon]:
