## Neurthino.jl

Neutrino oscillation probability calculator
Author KM3NeT
Popularity
4 Stars
Updated Last
1 Year Ago
Started In
July 2020

# Neurthino.jl

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.

## Basic Usage

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
``````

## Oscillations in Vacuum

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.

## Neutrino Propagation 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]`:

### Required Packages

View all packages

### Used By Packages

No packages found.