! ! Survival probability calculation for varying electron number density along neutrino path ! ! Based 'flavour basis' equations, section II, in ! ! Three flavour neutrino oscillations in matter: Flavour diagonal potentials, the adiabatic basis, and the CP phase ! James P. Kneller, Gail C. McLaughlin ! PHYSICAL REVIEW D 80, 053002 (2009) ! ! conventions ! ! mass eigenstates (vacuum) index imass = 1, 2 , 3 ! ! flavour eigenstates (e,mu,tau) index ialpha = 1, 2, 3 ! ! pseudo code outline ! ! acquire mixing parameters : mixing angles (3), phase angle (1) and known differences of squared mass (2) ! ! set flag to indicate normal / inverted mass hierachy and/or set angles omega(alpha) (to be used in equation 12) ! ! set values for arbritrary phases alpha, beta and epsilon (CP) ! ! acquire path information: distance travelled (n segments) and electron number density (n segments) ! ! calculate all variables which depend on mixing parameters: sin12, sin23, sin13, cos12, cos13, cos23, msq12, msq13, msq23 ! ! for all path length segments ! ! calculate elements of Hamiltonian matrix HF in flavour base (equations 6a to 6 f) ! ! calculate values of q(alpha) and r(alpha) (equations 10a to 10f) ! ! calculate values of TF, QF, RF (equations 8a to 8c or alternatively equations 9a to 10f) ! ! calculate eigenvalues eigen_value(alpha) (equation 12) ! ! find three eigenvectors eigen_vector(alpha) of Hamiltonian HF for given eigenvalues (use LINPACK, MKL, GSL, other ? ) ! ! Note that by definition : HF * eigen_vector(alpha) = eigen_value(alpha) * eigen_vector(alpha) ! Solution of Schroedinger equation : sum_over_alpha: ( a(alpha) * exp( -j * eigen_value(alpha) * time / hbar ) * eigen_vector(alpha) ) = sum_over_alpha: state(alpha,time) ! Boundary condition at start time: sum_over_alpha: ( a(alpha) * eigen_vector(alpha) ) = sum_over_alpha: previous_state(alpha) ( in practice a 3 x 3 system of equations ) ! ! find coefficients a(alpha) by numerically solving ( 3 x 3 ) linear system (question: are the coefficients real numbers?) ! ! calculate state at end of step: sum_over_alpha: state(alpha) = sum_over_alpha: a(alpha) * exp( -j * eigen_value(alpha) * time / hbar ) * eigen_vector(alpha) ! ! calculate survival probabilities P(beta) = ||**2 where alpha is state at end of last step ! ! enter next points P(beta,t) or P(beta,L) with L=c*t into a plot (histogram? points?) ! ! check: 0 < P(beta,L) < 1 and sum_over_beta: P(beta,L) = 1 ! ! if this is the end of the last step then exit loop ! ! generate plot, display, save plot ! ! Note(s): ! ! alternative algorythm: start with matrix HF and calculate both eigenvalues and eigenvectors numerically? Would this be any slower? ! ! TF, QF, RF or else eigen_value(alpha) plus eigen_vector(alpha) (total of 12 values) could be tabulated, and then recovered from a spline function ! ! test case: results should reproduce the vacuum case (zero matter density) ! ! test case: results should reproduce the constant matter density case (see journal paper) ! ! test case: results should reproduce the exponential matter density case (see journal paper) ! ! ! Questions: ! ! Implement as mixed language code? Fortran for complex numbers, vectors, matrices and if needed 128 bit precision? C++ to the use CERN root libraries for plots? ! ! Implement twice as two independent pieces of code as cross check? Is this needed as there are several test cases? ! ! Can precise results be approximated by average matter density over all path length? Or a low order polynominal mass density? ! ! Under which conditions can several neutrino path segments be combined into a single average density segment? ! ! What is the 'oscillation length scale' over which variations in abundance become significant? program neutrino_flux_calculation print*," Neutrino flux calculation" end program