Plot of neutrino oscillation probabilities using newest NuFit data.

oscillation.py

import numpy as np
# normal ordering values from NuFit 2.0, JHEP 11 (2014) 052 [arXiv:1409.5439]
θ_12 = np.deg2rad(33.48)
θ_23 = np.deg2rad(42.3)
θ_13 = np.deg2rad(8.50)
δ_CP = np.deg2rad(306)

Δm21 = np.sqrt(7.50e-5) # eV^2
Δm31 = np.sqrt(2.449e-3) # eV^2

masses = [0, Δm21, Δm31]

s12, c12 = np.sin(θ_12), np.cos(θ_12)
s23, c23 = np.sin(θ_23), np.cos(θ_23)
s13, c13 = np.sin(θ_13), np.cos(θ_13)

U = np.matrix([
    [c12*c13, s12*c13, s13*np.exp(-1j*δ_CP)],
    [-s12*c23 - c12*s23*s13*np.exp(1j*δ_CP), c12*c23 - s12*s23*s13*np.exp(1j*δ_CP), s23*c13],
    [s12*s23 - c12*c23*s13*np.exp(1j*δ_CP), -c12*s23 - s12*c23*s13*np.exp(1j*δ_CP), c23*c13]
])

def P(α, β, L, E):
    result = 0
    for i in range(3):
        # 1.267 accounts for factors of c and ħ
        result += U[α-1, i].conjugate() * U[β-1, i] * np.exp(-1j * 1.267 * masses[i]**2 * L/E)

    return abs(result)**2

import matplotlib.pyplot as plt

x = np.linspace(0, 70000, 2000)
E = 1
plt.plot(x, P(1, 1, x, E), label='$\\ell = e$')
plt.plot(x, P(1, 2, x, E), label='$\\ell = \\mu$')
plt.plot(x, P(1, 3, x, E), label='$\\ell = \\tau$')
plt.plot(x, P(1, 1, x, E) + P(1, 2, x, E) + P(1, 3, x, E), 'k--')
plt.ylabel('$P_{e\\to\\ell}$')
plt.xlabel('$\\frac{L}{E}\ /\ \\frac{\\mathrm{km}}{\\mathrm{GeV}}$')
plt.legend(loc='upper center')
plt.tight_layout()

plt.savefig('out.pdf')

out.pdf