Basic Neutrino Oscillation Functionality

Project.toml

@@ -4,8 +4,14 @@ authors = ["Johannes Schumann"]
 version = "0.1.0"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
+IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 julia = "1"

docs/src/index.md

@@ -12,6 +12,94 @@
 The main focus of the package lies on atmospheric neutrino flux and the neutrino
 propagation through earth.
 
+### Usage
+First of all the basic vacuum properties have to be defined via creating a
+`OscillationParameters` struct with fixed number of neutrino flavours of the 
+considered model.
+```julia
+julia> using Neurthino
+
+julia> osc = OscillationParameters(3);
+```
+The values of the mixing angles, mass squared differences and CP phases are 
+initialised to 0 and have to be set individually.
+```
+julia> osc.mixing_angles[1,2] = 0.59;
+
+julia> osc.mixing_angles[1,3] = 0.15;
+
+julia> osc.mixing_angles[2,3] = 0.84;
+
+julia> osc.cp_phases[1,3] = 3.86;
+```
+The mass squared differences are defined as <img src="https://render.githubusercontent.com/render/math?math=\Delta_{ij}=m_i^2-m_j^2"> and
+within the package the convention <img src="https://render.githubusercontent.com/render/math?math=\forall%20i%3Cj:m_i%3Cm_j"> is kept.
+```
+julia> osc.mass_squared_diff[1,3] = -2.523e-3;
+
+julia> osc.mass_squared_diff[2,3] = -2.523e-3;
+
+julia> osc.mass_squared_diff[1,2] = -7.39e-5;
+```
+These oscillation parameters can now be used in order to calculate the transition
+probabilities between the flavour states. 
+```
+julia> transprob(osc, 1, 10000)
+3×3 Array{Float64,2}:
+ 0.684582  0.0964727  0.218946
+ 0.059416  0.790557   0.150027
+ 0.256002  0.11297    0.631027
+```
+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 `transprob` function can also be called in the following way:
+```
+julia> U = PMNSMatrix(osc)
+3×3 Array{Complex,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 SparseArrays.SparseVector{Float64,Int64} with 3 stored entries:
+  [1]  =  -0.000865633
+  [2]  =  -0.000816367
+  [3]  =  0.001682
+
+julia> transprob(U, H, 1, 10000)
+3×3 Array{Float64,2}:
+ 0.684582  0.0964727  0.218946
+ 0.059416  0.790557   0.150027
+ 0.256002  0.11297    0.631027
+```
+For the neutrino propagation through matter a modified PMNS-Matrix and Hamiltonian
+has to be determined. The matter is parametrised by its density. 
+```
+julia> H_mat, U_mat = MatterOscillationMatrices(U, H, 13);
+
+julia> H_mat
+3-element Array{Complex{Float64},1}:
+ -0.0008318913110425209 + 1.402405683460369e-20im 
+  0.0009755131119987117 + 4.535370547007611e-21im 
+  0.0018408194174787428 - 2.1947559170628515e-20im
+
+julia> U_mat
+3×3 Array{Complex{Float64},2}:
+ 0.0117205-2.18724e-5im    0.8666+0.0im       -0.374878+0.32914im   
+ -0.665633+0.00279569im  0.287503+0.2445im     0.643807+0.0im       
+  0.746182+0.0im         0.241939+0.219159im   0.580206-0.00274678im
+```
+The transition probabilities can then be calculated using the `transprob` function
+again.
+```
+julia> transprob(U_mat, H_mat, 1, 10000)
+3×3 Array{Float64,2}:
+ 0.252869  0.428127  0.319004
+ 0.414751  0.306335  0.278914
+ 0.33238   0.265538  0.402082
+```
+
 <!-- ```@index -->
 <!-- ``` -->
 <!--  -->

src/Neurthino.jl

@@ -1,7 +1,23 @@
 module Neurthino
 
-import Polynomials: Poly
+using LinearAlgebra
+using SparseArrays
+using PhysicalConstants
+using Unitful
+using Polynomials
+
+import Base
+export transprob, OscillationParameters, PMNSMatrix, Hamiltonian, MatterOscillationMatrices
+
+# TEMPORARY UNTIL PhysicalConstants.jl GETS UPDATED
+G_F = 1.1663787e-5u"GeV^-2"
+G_F = G_F * (PhysicalConstants.CODATA2018.SpeedOfLightInVacuum * PhysicalConstants.CODATA2018.ReducedPlanckConstant)^3
+G_F = uconvert(u"eV*cm^3", G_F)
+
 
 include("PREM.jl")
+include("Oscillation.jl")
+
+
 
 end # module

src/Oscillation.jl

@@ -0,0 +1,212 @@
+struct OscillationParameters
+    dim::Integer
+    mixing_angles::UnitUpperTriangular{T, <: AbstractMatrix{T}} where {T <: Real}
+    mass_squared_diff::UnitUpperTriangular{T, <: AbstractMatrix{T}} where {T <: Real}
+    cp_phases::UnitUpperTriangular{T, <: AbstractMatrix{T}} where {T <: Real}
+
+    OscillationParameters(dim::Integer) = begin
+            new(dim,
+                UnitUpperTriangular(zeros(Float64, (dim, dim))),
+                UnitUpperTriangular(zeros(Float64, (dim, dim))),
+                UnitUpperTriangular(zeros(Float64, (dim, dim))))
+    end
+end
+
+function _generate_ordered_index_pairs(n::Integer)
+    number_of_angles = number_mixing_angles(n)
+    indices = Vector{Pair{Int64,Int64}}(undef, number_of_angles)
+    a = 1
+    for i in 1:n 
+        for j in 1:i-1
+            indices[a] = Pair(j,i)
+            a += 1
+        end
+    end
+    indices
+end
+
+"""
+    MatterOscillationMatrices(osc_params::OscillationParameters, matter_density)
+
+Create modified oscillation parameters for neutrino propagation through matter
+
+# Arguments
+- `osc_vacuum::OscillationParameters`: Oscillation parameters in vacuum
+- `matter_density`: Matter density in g*cm^-3 
+
+"""
+function MatterOscillationMatrices(osc_vacuum::OscillationParameters, matter_density)
+    H_vacuum = Diagonal(Hamiltonian(osc_vacuum)) 
+    U_vacuum = PMNSMatrix(osc_vacuum)
+    return MatterOscillationMatrices(H_vacuum, U_vacuum, matter_density)
+end
+
+"""
+    MatterOscillationMatrices(U_vac, H_vac, matter_density)
+
+Create modified oscillation parameters for neutrino propagation through matter
+
+# Arguments
+- `U_vac`: Vacuum PMNSMatrix
+- `H_vac`: Vacuum Hamiltonian
+- `matter_density`: Matter density in g*cm^-3 
+
+"""
+function MatterOscillationMatrices(U_vac, H_vac, matter_density)
+    H_flavour = U_vac * Diagonal(H_vac) * adjoint(U_vac)
+    A = 2 * sqrt(2) * ustrip(G_F) * ustrip(PhysicalConstants.CODATA2018.AvogadroConstant) * 1e9
+    A *= matter_density
+    H_flavour[1,1] += A  
+    U_matter = eigvecs(H_flavour)
+    H_matter = eigvals(H_flavour)
+    return H_matter, U_matter
+end
+
+function PMNSMatrix(osc_params::OscillationParameters)
+"""
+    PMNSMatrix(osc_params::OscillationParameters)
+
+Create rotation matrix (PMNS) based on the given oscillation parameters
+
+# Arguments
+- `osc_params::OscillationParameters`: Oscillation parameters
+
+"""
+    pmns = Matrix{Complex}(1.0I, osc_params.dim, osc_params.dim) 
+    indices = _generate_ordered_index_pairs(osc_params.dim)
+    for (i, j) in indices
+        rot = sparse((1.0+0im)I, osc_params.dim, osc_params.dim) 
+        mixing_angle = osc_params.mixing_angles[i, j]
+        c, s = cos(mixing_angle), sin(mixing_angle)
+        rot[i, i] = c
+        rot[j, j] = c
+        rot[i, j] = s
+        rot[j, i] = -s
+        if CartesianIndex(i, j) in findall(!iszero, osc_params.cp_phases)
+            cp_phase = osc_params.cp_phases[i, j]
+            cp_term = exp(-1im * cp_phase)
+            rot[i, j] *= cp_term
+            rot[j, i] *= conj(cp_term)
+        end
+        pmns = rot * pmns 
+    end
+    pmns
+end
+
+function Hamiltonian(osc_params::OscillationParameters)
+"""
+    Hamiltonian(osc_params::OscillationParameters)
+
+Create modified hamiltonian matrix consisting of the squared mass differences
+based on the given oscillation parameters
+
+# Arguments
+- `osc_params::OscillationParameters`: Oscillation parameters
+
+"""
+    Hamiltonian(osc_params, zeros(Float64, osc_params.dim))
+end 
+
+function Hamiltonian(osc_params::OscillationParameters, lambda)
+"""
+    Hamiltonian(osc_params::OscillationParameters)
+
+Create modified hamiltonian matrix consisting of the squared mass differences
+based on the given oscillation parameters
+
+# Arguments
+- `osc_params::OscillationParameters`:  Oscillation parameters
+- `lambda`:                             Decay parameters for each mass eigenstate
+
+"""
+    H = zeros(Complex, osc_params.dim)
+    for i in 1:osc_params.dim
+        for j in 1:osc_params.dim
+            if i < j
+                H[i] += osc_params.mass_squared_diff[i,j]
+            elseif j < i
+                H[i] -= osc_params.mass_squared_diff[j,i]
+            end
+        end
+        H[i] += 1im * lambda[i]
+    end
+    H /= osc_params.dim
+    H
+end
+
+
+"""
+    transprob(U, H, energy, baseline)
+
+Calculate the transistion probability between the neutrino flavours
+
+# Arguments
+- `U`:          Unitary transition matrix
+- `H`:          Energy eigenvalues
+- `energy`:     Baseline [km]
+- `baseline`:   Energy [GeV]
+
+"""
+function transprob(U, H, energy, baseline)  
+    H_diag = 2.534 * Diagonal(H) * baseline / energy 
+    A = U * exp(-1im * H_diag) * adjoint(U)
+    P = abs.(A) .^ 2
+end
+
+"""
+    transprob(osc_params::OscillationParameters, energy, baseline)
+
+Calculate the transistion probability between the neutrino flavours
+
+# Arguments
+- `osc_params::OscillationParameters`:  Oscillation parameters
+- `energy`:                             Baseline [km]
+- `baseline`:                           Energy [GeV]
+
+"""
+function transprob(osc_params::OscillationParameters, energy, baseline)  
+    H = Hamiltonian(osc_params)
+    U = PMNSMatrix(osc_params)
+    H_diag = 2.534 * Diagonal(H) * baseline / energy 
+    A = U * exp(-1im * H_diag) * adjoint(U)
+    P = abs.(A) .^ 2
+end
+
+"""
+    number_cp_phases(n::Unsigned)
+
+Returns the number of CP violating phases at given number of neutrino types
+
+# Arguments
+- `n::Unsigned`: number of neutrino types in the supposed model
+
+# Examples
+```julia-repl
+julia> Neurthino.number_cp_phases(3)
+1
+```
+"""
+function number_cp_phases(n::T) where {T <: Integer}
+    if (n < 1) return 0 end
+    cp_phases = div( (n-1)*(n-2) , 2 )
+end
+
+"""
+    number_mixing_angles(n::Unsigned)
+
+Returns the number of mixing angles at given number of neutrino types
+
+# Arguments
+- `n::Unsigned`: number of neutrino types in the supposed model
+
+# Examples
+```julia-repl
+julia> Neurthino.number_mixing_phases(3)
+3
+```
+"""
+function number_mixing_angles(n::T) where {T <: Integer}
+    if (n < 1) return 0 end
+    mixing_angles = div( n*(n-1) , 2 )
+end
+

src/PREM.jl

@@ -2,7 +2,7 @@ struct Layer
     name::AbstractString
     r_min
     r_max
-    density::Poly
+    density::Polynomial
 end
 
 
@@ -13,17 +13,17 @@ end
 
 PREM = EarthModel(
     [
-        Layer("Inner Core", 0.0, 1221.5, Poly([13.0885, 0.0, -8.8381])),
-        Layer("Outer Core", 1221.5, 3480.0, Poly([12.5815, -1.2638, -3.6426, -5.5281])),
-        Layer("Lower Mantle", 3480.0, 5701.0, Poly([7.9565, -6.4761, 5.5283, -3.0807])),
-        Layer("Transition Zone 1", 5701.0, 5771.0, Poly([5.3197, -1.4836])),
-        Layer("Transition Zone 2", 5771.0, 5971.0, Poly([11.2494, -8.0298])),
-        Layer("Transition Zone 3", 5971.0, 6151.0, Poly([7.1089, -3.8045])),
-        Layer("LVZ", 6151.0, 6291.0, Poly([2.6910, 0.6924])),
-        Layer("LID", 6291.0, 6346.6, Poly([2.6910, 0.6924])),
-        Layer("Crust 1", 6346.6, 6356.0, Poly([2.9])),
-        Layer("Crust 2", 6356.0, 6368.0, Poly([2.6])),
-        Layer("Ocean", 6368.0, 6371.0, Poly([1.020]))
+        Layer("Inner Core", 0.0, 1221.5, Polynomial([13.0885, 0.0, -8.8381])),
+        Layer("Outer Core", 1221.5, 3480.0, Polynomial([12.5815, -1.2638, -3.6426, -5.5281])),
+        Layer("Lower Mantle", 3480.0, 5701.0, Polynomial([7.9565, -6.4761, 5.5283, -3.0807])),
+        Layer("Transition Zone 1", 5701.0, 5771.0, Polynomial([5.3197, -1.4836])),
+        Layer("Transition Zone 2", 5771.0, 5971.0, Polynomial([11.2494, -8.0298])),
+        Layer("Transition Zone 3", 5971.0, 6151.0, Polynomial([7.1089, -3.8045])),
+        Layer("LVZ", 6151.0, 6291.0, Polynomial([2.6910, 0.6924])),
+        Layer("LID", 6291.0, 6346.6, Polynomial([2.6910, 0.6924])),
+        Layer("Crust 1", 6346.6, 6356.0, Polynomial([2.9])),
+        Layer("Crust 2", 6356.0, 6368.0, Polynomial([2.6])),
+        Layer("Ocean", 6368.0, 6371.0, Polynomial([1.020]))
     ]
 )