Numerical recipes
What NuOscProbExact can compute, with the code that computes it.
Each recipe below is a few lines and a figure. The figures are the ones the notebooks produce, so the code shown here and the notebook linked beside it are the same calculation — there is no third version to drift out of step. Where a recipe is short enough to be worth running on the spot, it is executed when this page is built and its output is what you see.
One probability
The shortest useful thing the library does. Give it a Hermitian matrix and a baseline in reciprocal units, and it returns the exact probabilities.
import numpy as np
import globaldefs as gd
import hamiltonians3nu
import oscprob3nu
KM = gd.CONV_KM_TO_INV_EV
GEV = 1.0e9
h_vacuum = hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
gd.S12_NO_BF, gd.S23_NO_BF, gd.S13_NO_BF, gd.DCP_NO_BF,
gd.D21_NO_BF, gd.D31_NO_BF)
prob = oscprob3nu.probabilities_3nu(
np.asarray(h_vacuum)/(1.0*GEV), 1300.0*KM)
print('P_ee = %.6f' % prob[0])
print('P_emu = %.6f' % prob[1])
print('P_etau = %.6f' % prob[2])
P_ee = 0.927678
P_emu = 0.014323
P_etau = 0.057999
The nine probabilities come back with the initial flavor varying slowest. Full walk-through: notebook 01.
A scan, without a loop
Pass an array and the whole scan is one call. This is the single most useful thing to know about using the library well.
energies = np.logspace(-1.0, 1.5, 500)*GEV
stack = np.asarray(h_vacuum)/energies[:, None, None]
probabilities = oscprob3nu.probabilities_3nu(stack, 1300.0*KM)
print('shape returned:', probabilities.shape)
print('P_mue at the first three energies:',
np.round(probabilities[:3, 3], 6))
shape returned: (500, 9)
P_mue at the first three energies: [0.208266 0.30167 0.402009]
Note the shape: a batched call returns (..., 9), with the flavor index
last, so probabilities[:, 3] is \(P_{\mu e}\) along the scan. A
scalar call returns a tuple of nine instead.
Vacuum oscillations at a 1300 km baseline. Code: notebook 02.
Matter, and new physics
Matter, non-standard interactions and Lorentz-invariance violation are not special cases in the code. Each is a different Hermitian matrix handed to the same routine.
h_matter = hamiltonians3nu.hamiltonian_3nu_matter(
h_vacuum, energies, gd.VCC_EARTH_CRUST)
h_nsi = hamiltonians3nu.hamiltonian_3nu_nsi(
h_vacuum, energies, gd.VCC_EARTH_CRUST, gd.EPS_3)
p_matter = oscprob3nu.probabilities_3nu(h_matter, 1300.0*KM)
p_nsi = oscprob3nu.probabilities_3nu(h_nsi, 1300.0*KM)
print('largest difference NSI vs standard matter: %.4f'
% np.max(np.abs(p_nsi[:, 3] - p_matter[:, 3])))
largest difference NSI vs standard matter: 0.0600
The MSW resonance in constant-density matter. Code: notebook 03.
An oscillogram
A two-dimensional map of energy against baseline, in one call. Index the two arguments so they broadcast against each other and the grid falls out.
n_e, n_l = 240, 240
energies = np.logspace(-1.0, 1.5, n_e)*GEV
baselines = np.linspace(50.0, 12000.0, n_l)*KM
h_stack = hamiltonians3nu.hamiltonian_3nu_matter(
h_vacuum, energies, gd.VCC_EARTH_CRUST)
# (n_e, 1, 3, 3) against (1, n_l) -> an (n_e, n_l) grid
grid = oscprob3nu.probabilities_3nu(h_stack[:, None, :, :],
baselines[None, :])[:, :, 3]
print('grid shape:', grid.shape, '--', grid.size, 'probabilities')
print('P_mue runs from %.4f to %.4f' % (grid.min(), grid.max()))
grid shape: (240, 240) -- 57600 probabilities
P_mue runs from 0.0000 to 0.5822
57 600 probabilities, one call, no Python loop. Code: notebook 04.
CP violation
Plotting the neutrino appearance probability against the antineutrino one, as \(\delta_{CP}\) runs through \(2\pi\), traces an ellipse. Matter pushes it off the diagonal, which is what makes the measurement hard.
Antineutrinos need both changes: conjugate the vacuum Hamiltonian and reverse the sign of the potential.
h_nu = hamiltonians3nu.hamiltonian_3nu_matter(
h_vacuum, 1.0*GEV, gd.VCC_EARTH_CRUST)
h_nubar = hamiltonians3nu.hamiltonian_3nu_matter(
np.conj(h_vacuum), 1.0*GEV, -gd.VCC_EARTH_CRUST)
print('P(numu -> nue) = %.6f'
% oscprob3nu.probabilities_3nu(h_nu, 1300.0*KM)[3])
print('P(numubar -> nuebar) = %.6f'
% oscprob3nu.probabilities_3nu(h_nubar, 1300.0*KM)[3])
P(numu -> nue) = 0.025898
P(numubar -> nuebar) = 0.018414
Bi-probability ellipses in matter. Code: notebook 05, and notebook 13 for antineutrinos in full.
Through the Earth
The Earth’s density is not constant, so the expansions do not apply to a whole
trajectory. They apply to any piece of it over which the density is taken
constant, which is what earth builds from the Preliminary Reference
Earth Model [DA81].
import earth
print('chord at costhz = -1 : %.0f km'
% earth.distance_traveled_inside_earth(-1.0))
print('density at the centre: %.4f g/cm^3' % earth.density_prem(0.0))
probabilities = earth.probabilities_3nu_earth(
h_vacuum, 8.0*GEV, -0.8, n_slabs_per_segment=6)
print('P_mumu at 8 GeV, costhz = -0.8: %.6f' % probabilities[4])
chord at costhz = -1 : 12742 km
density at the centre: 13.0885 g/cm^3
P_mumu at 8 GeV, costhz = -0.8: 0.527807
The Preliminary Reference Earth Model. Code: notebook 06.
An Earth oscillogram, in energy and zenith angle. Code: notebook 07.
Between two places on the Earth
The chord between two named sites, and the probability along it.
earth.probabilities_3nu_between_locations() does the lookup, the
geometry and the PREM slabbing in one call.
for source, detector in (('cern', 'gran_sasso'),
('fermilab', 'homestake'),
('tokai', 'kamioka')):
lat1, lon1 = earth.coordinates_of_named_location(source)
lat2, lon2 = earth.coordinates_of_named_location(detector)
chord = earth.chord_length_inside_earth(lat1, lon1, lat2, lon2)
p_mue = earth.probabilities_3nu_between_locations(
h_vacuum, 1.0*GEV, source, detector, n_slabs_per_segment=6)[3]
print('%-22s %8.1f km P_mue = %.6f'
% (source + ' to ' + detector, chord, p_mue))
cern to gran_sasso 728.6 km P_mue = 0.059820
fermilab to homestake 1284.7 km P_mue = 0.020344
tokai to kamioka 294.7 km P_mue = 0.033812
Those are the baselines the experiments quote: CNGS is 730 km, T2K 295 km. Fermilab to Homestake comes out at 1285 km against DUNE’s quoted 1300, the difference being that DUNE quotes the distance to the detector hall rather than the surface chord. Code: notebook 07.
An arbitrary matter profile
slabs takes any sequence of widths and Hamiltonians, so a profile can be
built by hand. Castle-wall profiles are the interesting case: the arrangement
of the matter can change the answer even when the mean density does not.
The effect is resonant, not generic — at most energies the two agree closely, and near a particular one they do not.
import slabs
widths_km = np.full(24, 250.0)
castle = np.where(np.arange(24) % 2 == 0, 2.0, 8.0)
uniform = np.full(24, castle.mean())
def appearance(densities, energy):
h = hamiltonians3nu.hamiltonian_3nu_matter(
h_vacuum, energy, earth.matter_potential(densities))
return slabs.probabilities_3nu_slabs(h, widths_km*KM)[3]
print('mean density, both cases: %.1f g/cm^3' % castle.mean())
for energy_gev in (0.44, 3.0):
print('E = %4.2f GeV : castle %.4f uniform %.4f' %
(energy_gev,
appearance(castle, energy_gev*GEV),
appearance(uniform, energy_gev*GEV)))
mean density, both cases: 5.0 g/cm^3
E = 0.44 GeV : castle 0.0104 uniform 0.0028
E = 3.00 GeV : castle 0.0457 uniform 0.0459
At 3 GeV the two are indistinguishable; at 0.44 GeV the castle wall gives nearly four times the appearance probability of a uniform slab of the same mean density.
Four profiles, one mean density. Code: notebook 08.
Mass ordering and the octant
globaldefs carries the NuFit best fit for both orderings, so comparing them
needs no numbers typed in. Matter is what separates them: the potential enters
with a definite sign, so it enhances the resonance for one ordering and
suppresses it for the other.
def h_vacuum_3nu(ordering='NO', s23=None):
"""Energy-independent vacuum Hamiltonian, for either ordering."""
if ordering == 'NO':
pars = (gd.S12_NO_BF, gd.S23_NO_BF, gd.S13_NO_BF,
gd.DCP_NO_BF, gd.D21_NO_BF, gd.D31_NO_BF)
else:
pars = (gd.S12_IO_BF, gd.S23_IO_BF, gd.S13_IO_BF,
gd.DCP_IO_BF, gd.D21_IO_BF, gd.D31_IO_BF)
s12, s23_bf, s13, dcp, d21, d31 = pars
return hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
s12, s23_bf if s23 is None else s23, s13, dcp, d21, d31)
def p_matter(h, energy):
"""The nine probabilities in crust matter, at 1300 km."""
return oscprob3nu.probabilities_3nu(
hamiltonians3nu.hamiltonian_3nu_matter(
h, energy, gd.VCC_EARTH_CRUST), 1300.0*KM)
print('normal : Dm31 = %+.4e eV^2' % gd.D31_NO_BF)
print('inverted : Dm31 = %+.4e eV^2' % gd.D31_IO_BF)
for ordering in ('NO', 'IO'):
print(' %s : P_mue at 2.5 GeV = %.4f'
% (ordering, p_matter(h_vacuum_3nu(ordering), 2.5*GEV)[3]))
normal : Dm31 = +2.5250e-03 eV^2
inverted : Dm31 = -2.4381e-03 eV^2
NO : P_mue at 2.5 GeV = 0.0872
IO : P_mue at 2.5 GeV = 0.0393
The octant of \(\theta_{23}\) is the other open question, and it needs the appearance channel rather than the disappearance one:
for s23_squared in (0.45, 0.55):
p = p_matter(h_vacuum_3nu('NO', s23=np.sqrt(s23_squared)), 5.0*GEV)
print('sin^2(theta23) = %.2f : P_mumu = %.4f P_mue = %.4f'
% (s23_squared, p[4], p[3]))
sin^2(theta23) = 0.45 : P_mumu = 0.4816 P_mue = 0.0245
sin^2(theta23) = 0.55 : P_mumu = 0.4768 P_mue = 0.0299
Disappearance depends on \(\theta_{23}\) mainly through \(\sin^2 2\theta_{23}\), which is symmetric about maximal mixing, so the two values either side of it are nearly indistinguishable there — the octant degeneracy. Appearance carries \(\sin^2\theta_{23}\) instead and tells them apart.
Matter through the Earth separates the two orderings. Code: notebook 12.
A sterile neutrino
Four flavors is the same call with a bigger matrix. A 3+1 scenario is a closed four-state system, not a leak out of the three-flavor block, which is what brings it inside an exact method at all.
import numpy as np
import globaldefs as gd
import hamiltonians4nu
import oscprob4nu
h4 = hamiltonians4nu.hamiltonian_4nu_vacuum_energy_independent(
gd.S12_NO_BF, gd.S23_NO_BF, gd.S13_NO_BF,
np.sqrt(0.10), np.sqrt(0.10), 0.0,
gd.DCP_NO_BF, gd.D21_NO_BF, gd.D31_NO_BF, 1.0)
prob = oscprob4nu.probabilities_4nu(np.asarray(h4)/1.0e9,
1300.0*gd.CONV_KM_TO_INV_EV)
print('%d probabilities, initial flavor slowest' % len(prob))
print('P(nu_mu -> nu_mu) = %.5f' % prob[5])
print('P(nu_mu -> nu_s) = %.5f' % prob[7])
16 probabilities, initial flavor slowest
P(nu_mu -> nu_mu) = 0.40717
P(nu_mu -> nu_s) = 0.01166
In matter the sterile state changes the problem qualitatively. It feels neither potential, so the neutral-current term — which is proportional to the identity across the three active flavors, and therefore invisible at two and three flavors — no longer cancels. Removing it from all four states costs only a global phase and leaves \(-V_{NC}\) on the sterile entry, and that entry is what places the sterile matter resonance.
h4_matter = hamiltonians4nu.hamiltonian_4nu_matter(
h4, 1.0e9, gd.VCC_EARTH_CRUST, gd.VNC_EARTH_CRUST)
print('nu_e entry : %+.4e eV' % h4_matter[0][0].real)
print('sterile : %+.4e eV' % h4_matter[3][3].real)
nu_e entry : +5.0149e-11 eV
sterile : +4.0511e-10 eV
Through the Earth, the same PREM machinery applies:
earth.probabilities_4nu_earth() cuts the chord at every shell
boundary and builds both potentials per slab.
import earth
prob = earth.probabilities_4nu_earth(h4, 1.0e10, -0.8)
print('P(nu_mu -> nu_mu) through the Earth = %.5f' % prob[5])
P(nu_mu -> nu_mu) through the Earth = 0.73238
Full walk-through: notebook 16.
Where to go next
Quickstart — the shortest path to a first probability.
Methodology — what the SU(2), SU(3) and SU(4) expansions actually do, and the sign conventions that matter once a matter potential is added.
API reference — the full API reference, generated from the docstrings.
Notebook 17 — the same probabilities cross-checked against an independent external code, and against a published closed form, with the conventions that have to be matched first.
The notebooks — eighteen of them, carrying their figures inline.