API reference

Generated from the docstrings. Every Examples block below is executed when this page is built, so the results shown are what the code returns rather than numbers written beside it. The regression suite runs the same blocks on every supported Python.

Core modules

The three core modules are self-contained: they need only numpy, and they accept any Hermitian Hamiltonian. Copying one into your own project is a supported way to use NuOscProbExact.

All three are unit-agnostic: they require only that the Hamiltonian and the baseline be given in reciprocal units, so that \(H L\) is dimensionless.

oscprob2nu

Compute the two-neutrino flavor-transition probabilities.

This module contains the routines needed to compute two-neutrino flavor-transition probabilities for an arbitrary time-independent \(2\times2\) Hermitian Hamiltonian, using the SU(2) exponential expansion described in [1].

The Hamiltonian is expanded in the basis of Pauli matrices,

\[H = h_0 \mathbb{1} + h_k \sigma^k ,\]

and the time-evolution operator in the same basis,

\[U_2(L) = u_0 \mathbb{1} + i u_k \sigma^k .\]

The term \(h_0\) contributes only an overall phase and is dropped; all routines therefore work with the traceless part of the Hamiltonian, which leaves the oscillation probabilities unchanged.

evolution_operator_2nu and probabilities_2nu accept either a single Hamiltonian and baseline or a stack of them, in which case the whole stack is evaluated at once.

Units

The routines are unit-agnostic: they require only that the Hamiltonian and the baseline be given in reciprocal units, so that the product \(H L\) is dimensionless. Elsewhere in NuOscProbExact the Hamiltonian is in eV and the baseline in eV^-1; the module globaldefs provides CONV_KM_TO_INV_EV to convert a baseline in km into eV^-1.

Routine listings

• hamiltonian_2nu_coefficients - Returns the \(h_k\)

• modulus - Returns the modulus \(|h|\) of a vector

• evolution_operator_2nu_u_coefficients - Returns the \(u_k\)

• evolution_operator_2nu - Returns the evolution operator \(U_2\)

• probabilities_2nu - Returns the oscillation probabilities

References

[1]

Mauricio Bustamante, “Exact neutrino oscillation probabilities with arbitrary time-independent Hamiltonians”, arXiv:1904.12391.

oscprob2nu.CHECK_HERMITICITY = True

bool: Module-level switch.

Whether to verify that the Hamiltonian handed in is Hermitian before evaluating anything. It is on by default, and the reason is that the failure it catches is silent: a non-Hermitian matrix does not raise, and does not produce obviously broken output either — the probabilities it returns still sum to one, so the usual sanity check a caller would apply cannot tell that anything went wrong. The expansion assumes Hermiticity; without it the result is meaningless rather than merely inaccurate.

Checking is not free, and the cost is stated here rather than buried, because it is larger than one might expect. Validating a stack is a pass over it, which is the same order of work as evaluating it — and the compiled kernel has made evaluating it fast, so on a large stack the check dominates. Measured by interleaving the two settings and taking the best of fifteen rounds each:

Stack	2 flavors	4 flavors
2 000 points	1.5x	1.3x
200 000 points	5.7x	3.2x

Three flavors sits between them, at 1.8x and 3.9x. So a scan that the compiled backend was installed to speed up can spend most of its time here instead. On a single matrix the check costs about 1.35x, and it takes its own branch to manage that: the reductions above are all fixed cost at one element, so without the branch one probability spent nine tenths of its time here. For production scans whose Hamiltonians come from a construction already known to be Hermitian — everything hamiltonians2nu builds is Hermitian to round-off, as the table below records — set this to False.

The default is nevertheless True, because the alternative is a library that silently returns meaningless numbers to anyone who makes this mistake once, and finding out costs far more than the check does.

The tolerance is relative to the largest entry of the Hamiltonian, so a matrix assembled in floating point passes: the ones this library builds are Hermitian to about \(2 \times 10^{-17}\) relative, against a tolerance of \(10^{-12}\).

Added in version 1.11.0.

oscprob2nu.SMALL_BATCH = 11

int: Module-level constant.

Stacks with at most this many elements are evaluated one at a time through the scalar path, whose fixed cost is lower than the array machinery’s. The measured crossover is twelve elements, and it governs every two-flavor stack below fastkernels.MIN_BATCH, which is fifty thousand here — so unlike oscprob3nu.SMALL_BATCH this one is consulted whether or not the compiled backend is installed.

The threshold was 6, measured before CHECK_HERMITICITY existed and made a scalar call several times dearer; with the check given a path for a single matrix it was re-measured. It sits close to the three-flavor crossover despite that expansion doing much more work per element, because what is being amortised is the array machinery’s fixed cost rather than the arithmetic.

oscprob2nu.hamiltonian_2nu_coefficients(hamiltonian_matrix: list | ndarray) → List[float]

Returns the \(h_k\) of the SU(2) expansion of the Hamiltonian.

Computes the coefficients \(h_1, h_2, h_3\) of the SU(2) expansion \(H = h_0 \mathbb{1} + h_k \sigma^k\) of the two-flavor Hamiltonian hamiltonian_matrix, which is assumed to be given in the flavor basis. The coefficient \(h_0\) contributes only an overall phase to the evolution operator and is not returned.

Added in version 1.0.0.

Changed in version 1.1.0: Returns real floats. The coefficients of a Hermitian Hamiltonian are real by construction, but the routine previously returned a mixture of floats and complex numbers.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module.

Parameters:

hamiltonian_matrixarray_like

Two-flavor Hamiltonian, given as the nested list [[H11, H12], [H21, H22]]. It must be Hermitian, i.e. H21 == conj(H12) and H11, H22 real.

Returns:

list of float

The three coefficients [h1, h2, h3]. They are real, because the Hamiltonian is Hermitian.

See also

modulus

Returns the modulus \(|h|\) of the returned vector.

Examples

import oscprob2nu

hamiltonian_matrix = [[1.0+0.0j, 0.0+2.0j],
                      [0.0-2.0j, 3.0+0.0j]]
h_coeffs = oscprob2nu.hamiltonian_2nu_coefficients(hamiltonian_matrix)
print('%.6f  %.6f  %.6f' % tuple(h_coeffs))

0.000000  -2.000000  -1.000000

oscprob2nu.modulus(h_coeffs: list | ndarray) → float

Returns the modulus \(|h|\) of the vector of coefficients.

Returns the modulus of the vector of coefficients \(h_k\) of the SU(2) expansion of the two-neutrino Hamiltonian, \(|h| = \sqrt{|h_1|^2 + |h_2|^2 + |h_3|^2}\).

Added in version 1.0.0.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The square root is taken with math rather than NumPy.

Parameters:

h_coeffsarray_like

Three-component vector of coefficients \(h_k\), as returned by hamiltonian_2nu_coefficients.

Returns:

float

The modulus \(|h|\).

Examples

import oscprob2nu

print('%.6f' % oscprob2nu.modulus([0.0, -2.0, -1.0]))

2.236068

oscprob2nu.evolution_operator_2nu_u_coefficients(hamiltonian_matrix: list | ndarray, L: int | float) → List[float]

Returns the coefficients \(u_0, \ldots, u_3\).

Returns the four coefficients \(u_0, \ldots, u_3\) of the two-neutrino time-evolution operator \(U_2(L)\) in its SU(2) exponential expansion, \(U_2 = u_0 \mathbb{1} + i u_k \sigma^k\).

Added in version 1.0.0.

Changed in version 1.1.0: Degenerate Hamiltonians are handled exactly instead of returning NaN, by taking the limit \(\sin(|h|L)/|h| \to L\).

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module.

Parameters:

hamiltonian_matrixarray_like

Two-flavor Hermitian Hamiltonian, given as the nested list [[H11, H12], [H21, H22]].

Lfloat

Baseline, in units reciprocal to those of the Hamiltonian.

Returns:

list of float

The four coefficients [u0, u1, u2, u3]. They are real, because the Hamiltonian is Hermitian; the factor \(i\) that multiplies \(u_k\) is part of the expansion, not of the coefficients.

Notes

When \(|h| = 0\) the Hamiltonian is proportional to the identity, there is no flavor evolution, and the limit \(\sin(|h| L)/|h| \to L\) is used.

Examples

import oscprob2nu

hamiltonian_matrix = [[1.0+0.0j, 0.0+2.0j],
                      [0.0-2.0j, 3.0+0.0j]]
u_coeffs = oscprob2nu.evolution_operator_2nu_u_coefficients(hamiltonian_matrix,
                                                 1.0)
print('%.6f  %.6f  %.6f  %.6f' % tuple(u+0.0 for u in u_coeffs))

-0.617273  0.000000  0.703690  0.351845

oscprob2nu.evolution_operator_2nu(hamiltonian_matrix: list | ndarray, L: int | float | list | ndarray) → List[List[complex]] | ndarray

Returns the two-neutrino time-evolution operator.

Returns the two-neutrino time-evolution operator \(U_2(L)\) in its SU(2) exponential expansion \(U_2(L) = u_0 \mathbb{1} + i u_k \sigma^k\). This is a \(2\times2\) unitary matrix.

Added in version 1.0.0.

Changed in version 1.1.0: Degenerate Hamiltonians are handled exactly instead of returning NaN, by taking the limit \(\sin(|h|L)/|h| \to L\).

Changed in version 1.2.0: Accepts a stack of Hamiltonians of shape (..., n, n), an array of baselines, or both broadcast against each other, returning an array with the broadcast leading axes. A single Hamiltonian with a scalar baseline returns exactly what it returned before.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module.

Changed in version 1.5.0: Faster, with identical results; the probabilities agree with 1.4.0 to 1.6e-13 across every code path. The batched coefficients are laid out with the component index first, so that each \(h_k\) is contiguous; every downstream step works one component at a time, and reading those from a strided view cost about a third more.

Parameters:

hamiltonian_matrixarray_like

Two-flavor Hermitian Hamiltonian, given as the nested list [[H11, H12], [H21, H22]], or a stack of them, of shape (..., 2, 2).

Lfloat or array_like

Baseline, in units reciprocal to those of the Hamiltonian, or an array of baselines broadcastable against the leading axes of hamiltonian_matrix.

Returns:

list of list of complex or numpy.ndarray

For a single Hamiltonian and baseline, the time-evolution operator \(U_2(L)\) — a \(2\times2\) unitary complex matrix — as a nested list. If either argument is a stack, an array of shape (..., 2, 2).

See also

probabilities_2nu

Returns the probabilities directly, more cheaply.

Examples

import oscprob2nu

hamiltonian_matrix = [[1.0+0.0j, 0.0+2.0j],
                      [0.0-2.0j, 3.0+0.0j]]
U2 = oscprob2nu.evolution_operator_2nu(hamiltonian_matrix, 1.0)
for row in U2:
    print('  '.join(['%+.6f%+.6fj' % (z.real+0.0, z.imag+0.0)
                     for z in row]))

-0.617273+0.351845j  +0.703690+0.000000j
-0.703690+0.000000j  -0.617273-0.351845j

oscprob2nu.probabilities_2nu(hamiltonian_matrix: list | ndarray, L: int | float | list | ndarray) → Tuple[float, float, float, float] | ndarray

Returns the two-neutrino oscillation probabilities.

Returns the two-neutrino flavor-transition probabilities \(P_{ee}, P_{e\mu}, P_{\mu e}, P_{\mu\mu}\), where \(P_{\alpha\beta} \equiv P(\nu_\alpha \to \nu_\beta)\).

Added in version 1.0.0.

Changed in version 1.1.0: The \(h_2\) contribution was restored. The transition probability is \(|U_{\mu e}|^2 = u_1^2 + u_2^2\), but the routine computed only \(|h_1|^2/|h|^2 \sin^2(|h|L)\). Since \(h_2 = -\mathrm{Im}(H_{12})\), this affected every Hamiltonian with a complex off-diagonal entry; oscillations in vacuum and in matter of constant density were unaffected.

Changed in version 1.2.0: Accepts a stack of Hamiltonians of shape (..., n, n), an array of baselines, or both broadcast against each other, returning an array with the broadcast leading axes. A single Hamiltonian with a scalar baseline returns exactly what it returned before.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module. A scalar two-flavor probability is 3.3x quicker, measured best-of-seven interleaved against 1.3.0.

Changed in version 1.5.0: Faster, with identical results; the probabilities agree with 1.4.0 to 1.6e-13 across every code path. The routine no longer builds the evolution operator, squares it, then transposes and reshapes the result. The coefficients of a Hermitian 2x2 Hamiltonian are real, so \(|U_{ee}|^2 = u_0^2 + u_3^2\) and \(|U_{\mu e}|^2 = u_1^2 + u_2^2\) are two numbers, the second the complement of the first by unitarity: neither the operator nor the coefficients are needed. The batched coefficients are laid out with the component index first, so that each \(h_k\) is contiguous; every downstream step works one component at a time, and reading those from a strided view cost about a third more. The scalar path is 3.3x quicker and a 2000-point scan 4.6x.

Changed in version 1.6.0: Two dispatch decisions, neither changing the result. A stack of at most SMALL_BATCH elements is evaluated one at a time through the scalar path, because a batched call carries a fixed cost that a handful of points does not amortise. Above fastkernels.MIN_BATCH[2] elements, and only if the optional numba extra is installed, the stack is evaluated by a compiled kernel instead; see fastkernels for the measured thresholds and why the two-flavor one is high.

Parameters:

hamiltonian_matrixarray_like

Two-flavor Hermitian Hamiltonian, given as the nested list [[H11, H12], [H21, H22]], or a stack of them, of shape (..., 2, 2).

Lfloat or array_like

Baseline, in units reciprocal to those of the Hamiltonian, or an array of baselines broadcastable against the leading axes of hamiltonian_matrix.

Returns:

tuple of float or numpy.ndarray

For a single Hamiltonian and baseline, the probabilities (Pee, Pem, Pme, Pmm) as a tuple. If either argument is a stack, an array of shape (..., 4) in the same order.

Notes

Passing arrays evaluates the whole stack at once, which on a 2000-point scan is about ninety times faster than calling this routine in a Python loop — the figure measured in README.md and guarded by tests/test_documented_figures.py; see the notes on oscprob3nu.probabilities_3nu() for the two scans that broadcast naturally.

The transition probability is

\[P_{e\mu} = \frac{|h_1|^2 + |h_2|^2}{|h|^2} \sin^2(|h| L) ,\]

i.e. \(|U_{\mu e}|^2 = u_1^2 + u_2^2\). Both \(h_1\) and \(h_2\) contribute; \(h_2\) vanishes only when the off-diagonal entry of the Hamiltonian is real.

Examples

import oscprob2nu

hamiltonian_matrix = [[1.0+0.0j, 0.0+2.0j],
                      [0.0-2.0j, 3.0+0.0j]]
Pee, Pem, Pme, Pmm = oscprob2nu.probabilities_2nu(hamiltonian_matrix, 1.0)
print('%.6f  %.6f  %.6f  %.6f' % (Pee, Pem, Pme, Pmm))

0.504821  0.495179  0.495179  0.504821

oscprob3nu

Compute the three-neutrino flavor-transition probabilities.

This module contains the routines needed to compute three-neutrino flavor-transition probabilities for an arbitrary time-independent \(3\times3\) Hermitian Hamiltonian, using the SU(3) exponential expansion described in [2].

The Hamiltonian is expanded in the basis of Gell-Mann matrices, whose structure constants and \(d\) tensor follow the conventions of [1],

\[H = h_0 \mathbb{1} + h_k \lambda^k ,\]

and the time-evolution operator in the same basis,

\[U_3(L) = u_0 \mathbb{1} + i u_k \lambda^k .\]

The term \(h_0\) contributes only an overall phase and is dropped; all routines therefore work with the traceless part of the Hamiltonian, which leaves the oscillation probabilities unchanged.

evolution_operator_3nu and probabilities_3nu accept either a single Hamiltonian and baseline or a stack of them, in which case the whole stack is evaluated at once. See the notes on probabilities_3nu.

Units

The routines are unit-agnostic: they require only that the Hamiltonian and the baseline be given in reciprocal units, so that the product \(H L\) is dimensionless. Elsewhere in NuOscProbExact the Hamiltonian is in eV and the baseline in eV^-1; the module globaldefs provides CONV_KM_TO_INV_EV to convert a baseline in km into eV^-1.

Routine listings

• hamiltonian_3nu_coefficients - Returns the \(h_k\)

• tensor_d - Returns the SU(3) tensor \(d_{ijk}\)

• star - Returns the SU(3) star product \((h \star h)_i\)

• su3_invariants - Returns the SU(3) invariants \(|h|^2, \langle h \rangle\)

• psi_roots - Returns the roots of the characteristic equation

• evolution_operator_3nu_u_coefficients - Returns the \(u_k\)

• evolution_operator_3nu - Returns the evolution operator \(U_3\)

• probabilities_3nu - Returns the oscillation probabilities

References

[1]

A.J. MacFarlane, A. Sudbery, and P.H. Weisz, “On Gell-Mann’s \(\lambda\)-matrices, \(d\)- and \(f\)-tensors, octets, and parametrizations of SU(3)”, Commun. Math. Phys. 11, 77 (1968).

[2]

Mauricio Bustamante, “Exact neutrino oscillation probabilities with arbitrary time-independent Hamiltonians”, arXiv:1904.12391.

oscprob3nu.CHECK_HERMITICITY = True

bool: Module-level switch.

Whether to verify that the Hamiltonian handed in is Hermitian before evaluating anything. It is on by default, and the reason is that the failure it catches is silent: a non-Hermitian matrix does not raise, and does not produce obviously broken output either — the probabilities it returns still sum to one, so the usual sanity check a caller would apply cannot tell that anything went wrong. The expansion assumes Hermiticity; without it the result is meaningless rather than merely inaccurate.

Checking is not free, and the cost is stated here rather than buried, because it is larger than one might expect. Validating a stack is a pass over it, which is the same order of work as evaluating it — and the compiled kernel has made evaluating it fast, so on a large stack the check dominates. Measured by interleaving the two settings and taking the best of fifteen rounds each:

Stack	2 flavors	4 flavors
2 000 points	1.5x	1.3x
200 000 points	5.7x	3.2x

Three flavors sits between them, at 1.8x and 3.9x. So a scan that the compiled backend was installed to speed up can spend most of its time here instead. On a single matrix the check costs about 1.35x, and it takes its own branch to manage that: the reductions above are all fixed cost at one element, so without the branch one probability spent nine tenths of its time here. For production scans whose Hamiltonians come from a construction already known to be Hermitian — everything hamiltonians3nu builds is Hermitian to round-off, as the table below records — set this to False.

The default is nevertheless True, because the alternative is a library that silently returns meaningless numbers to anyone who makes this mistake once, and finding out costs far more than the check does.

The tolerance is relative to the largest entry of the Hamiltonian, so a matrix assembled in floating point passes: the ones this library builds are Hermitian to about \(2 \times 10^{-17}\) relative, against a tolerance of \(10^{-12}\).

Added in version 1.11.0.

oscprob3nu.SMALL_BATCH = 12

int: Module-level constant.

Stacks with at most this many elements are evaluated one at a time through the scalar path. A batched call carries a fixed cost of a couple of hundred microseconds — allocating and reducing a dozen small arrays — which for a short stack exceeds what the scalar path spends on the whole job. Measured crossover: thirteen elements.

This governs the NumPy path only. With the compiled backend installed fastkernels.worthwhile() sends every three- and four-flavor stack to the kernel before this is consulted, because fastkernels.MIN_BATCH is one there.

The threshold was 10, measured before CHECK_HERMITICITY existed. The check then made a scalar call nine times dearer, which moved the crossover below one — batching won at every size, and this constant was sending stacks the slow way round. With the check given a path for a single matrix the scalar route is cheap again, and the crossover re-measured at thirteen. The two-flavor expansion does less work per element; see oscprob2nu.SMALL_BATCH.

oscprob3nu.hamiltonian_3nu_coefficients(hamiltonian_matrix: list | ndarray) → List[float]

Returns the \(h_k\) of the SU(3) expansion of the Hamiltonian.

Computes the coefficients \(h_1, \ldots, h_8\) of the SU(3) expansion \(H = h_0 \mathbb{1} + h_k \lambda^k\) of the three-flavor Hamiltonian hamiltonian_matrix, which is assumed to be given in the flavor basis. The coefficient \(h_0\) contributes only an overall phase to the evolution operator and is not returned.

Added in version 1.0.0.

Changed in version 1.1.0: Returns real floats. The coefficients of a Hermitian Hamiltonian are real by construction, but the routine previously returned a mixture of floats and complex numbers.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module.

Parameters:

hamiltonian_matrixarray_like

Three-flavor Hamiltonian, given as the nested list [[H11, H12, H13], [H21, H22, H23], [H31, H32, H33]]. It must be Hermitian.

Returns:

list of float

The eight coefficients [h1, h2, h3, h4, h5, h6, h7, h8]. They are real, because the Hamiltonian is Hermitian.

Examples

import oscprob3nu

hamiltonian_matrix = [[1.0+0.0j, 0.0+2.0j, 0.0-1.0j],
                      [0.0-2.0j, 3.0+0.0j, 3.0+0.0j],
                      [0.0+1.0j, 3.0-0.0j, -5.0+0.0j]]
h_coeffs = oscprob3nu.hamiltonian_3nu_coefficients(hamiltonian_matrix)
print('  '.join(['%.6f' % (h+0.0) for h in h_coeffs]))

0.000000  -2.000000  -1.000000  0.000000  1.000000  3.000000  0.000000  4.041452

oscprob3nu.tensor_d(i: int, j: int, k: int) → float

Returns the tensor \(d_{ijk}\) of the SU(3) algebra.

Returns the totally symmetric SU(3) tensor \(d_{ijk} = \frac{1}{4}\mathrm{Tr} (\{\lambda_i, \lambda_j\} \lambda_k)\), defined in [1].

Added in version 1.0.0.

Changed in version 1.1.0: Validates its indices and raises IndexError on one outside 0-7. The dispatch previously fell off the end of an elif chain and returned None, so the failure surfaced far from its cause.

Parameters:

iint

First index, in the range 0–7 (i.e., \(d_{ijk}\) is indexed from zero, so that i = 0 corresponds to \(d_{1jk}\)).

jint

Second index, in the range 0–7.

kint

Third index, in the range 0–7.

Returns:

float

The value of \(d_{ijk}\).

Raises:

IndexError

If any index lies outside the range 0–7.

References

[1]

A.J. MacFarlane, A. Sudbery, and P.H. Weisz, “On Gell-Mann’s \(\lambda\)-matrices, \(d\)- and \(f\)-tensors, octets, and parametrizations of SU(3)”, Commun. Math. Phys. 11, 77 (1968).

Examples

import oscprob3nu

print('%.6f' % oscprob3nu.tensor_d(0, 0, 7))
print('%.6f' % oscprob3nu.tensor_d(0, 1, 2))

0.577350
0.000000

oscprob3nu.star(i: int, h_coeffs: list | ndarray) → float

Returns the SU(3) star product \((h \star h)_i\).

Returns the SU(3) star product \((h \star h)_i = d_{ijk} h^j h^k\), summed over repeated indices.

Added in version 1.0.0.

Parameters:

iint

Index of the star product, in the range 0–7.

h_coeffsarray_like

Eight-component vector of coefficients \(h_k\), as returned by hamiltonian_3nu_coefficients.

Returns:

float

The star product \((h \star h)_i\).

Examples

import oscprob3nu

print('%.6f' % oscprob3nu.star(0, [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]))
print('%.6f' % oscprob3nu.star(7, [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]))

0.000000
0.577350

oscprob3nu.su3_invariants(h_coeffs: list | ndarray) → Tuple[float, float]

Returns the two SU(3) invariants, \(|h|^2\) and \(\langle h \rangle\).

Returns the two invariants of the SU(3) expansion, \(|h|^2 = h_i h_i\) and \(\langle h \rangle = d_{ijk} h_i h_j h_k\). They equal, respectively, \(\mathrm{Tr}(H_0^2)/2\) and \(\mathrm{Tr}(H_0^3)/2\), with \(H_0\) the traceless part of the Hamiltonian.

Added in version 1.0.0.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The star product is computed once and passed on, rather than once to form \(\langle h \rangle\) and again inside the expansion, and the sparse expansion of \(d_{ijk} h_j h_k\) replaces the dense table for the eight scalar components.

Parameters:

h_coeffsarray_like

Eight-component vector of coefficients \(h_k\), as returned by hamiltonian_3nu_coefficients.

Returns:

h2float

The SU(3) invariant \(|h|^2\).

h3float

The SU(3) invariant \(\langle h \rangle\).

Examples

import oscprob3nu

h2, h3 = oscprob3nu.su3_invariants([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
print('%.6f  %.6f' % (h2, h3))

1.000000  0.000000

oscprob3nu.psi_roots(h2: int | float, h3: int | float) → List[float]

Returns the three latent roots \(\psi\).

Returns the three latent roots \(\psi\) of the characteristic equation \(\psi^3 - |h|^2 \psi - \frac{2}{3}\langle h \rangle = 0\), which are the eigenvalues of minus the traceless part of the Hamiltonian. The roots are independent of the baseline.

Added in version 1.0.0.

Changed in version 1.1.0: No longer returns NaN when the traceless part of the Hamiltonian vanishes, and the arc-cosine argument is clipped to \([-1, 1]\) so that round-off cannot make the roots complex and the evolution operator non-unitary.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module.

Changed in version 1.5.0: Faster, with identical results; the probabilities agree with 1.4.0 to 1.6e-13 across every code path. The prefactor takes \(\sqrt{|h|^2}\) once rather than three times. The arc-cosine argument still forms \(|h|^2\) to the power of -1.5, which is what the code does and what an earlier version of this note wrongly described as a division.

Parameters:

h2float

The SU(3) invariant \(|h|^2\).

h3float

The SU(3) invariant \(\langle h \rangle\).

Returns:

list of float

The three roots [psi1, psi2, psi3]. They are real, because the Hamiltonian is Hermitian.

Notes

For a Hermitian Hamiltonian the argument of the arc cosine lies in \([-1, 1]\); it is clipped to that interval so that round-off cannot produce spurious complex roots, which would spoil the unitarity of the evolution operator. When \(|h|^2 = 0\) the Hamiltonian is proportional to the identity and all three roots vanish.

Examples

import oscprob3nu

psi = oscprob3nu.psi_roots(1.0, 0.0)
print('  '.join(['%.6f' % (round(p, 9)+0.0) for p in sorted(psi)]))

-1.000000  0.000000  1.000000

oscprob3nu.evolution_operator_3nu_u_coefficients(hamiltonian_matrix: list | ndarray, L: int | float) → List[complex]

Returns the coefficients \(u_0, \ldots, u_8\).

Returns the nine coefficients \(u_0, \ldots, u_8\) of the three-neutrino time-evolution operator \(U_3(L)\) in its SU(3) exponential expansion, \(U_3 = u_0 \mathbb{1} + i u_k \lambda^k\).

Added in version 1.0.0.

Changed in version 1.1.0: Degenerate Hamiltonians are handled exactly instead of returning NaN. Lagrange interpolation over the latent roots divides by \(3\psi_m^2 - |h|^2\), which vanishes at a repeated root; the two degenerate cases are now taken in their confluent limit.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The sum over the latent roots forms its \(k\)-independent factors once instead of inside the loop over the eight \(k\), and the star product is computed once and passed in rather than recomputed.

Changed in version 1.5.0: Faster, with identical results; the probabilities agree with 1.4.0 to 1.6e-13 across every code path. Around the latent roots, \(\sqrt{|h|^2}\) is taken once rather than three times, the arc-cosine argument is a division rather than a power of -1.5, and the degeneracy test uses two minima instead of stacking three gap arrays. The scalar exponentials use cmath.rect().

Parameters:

hamiltonian_matrixarray_like

Three-flavor Hermitian Hamiltonian, given as the nested list [[H11, H12, H13], [H21, H22, H23], [H31, H32, H33]].

Lfloat

Baseline, in units reciprocal to those of the Hamiltonian.

Returns:

list of complex

The nine coefficients [u0, u1, ..., u8].

Notes

The general expression divides by \(3\psi_m^2 - |h|^2\), which vanishes when two latent roots coincide. Two degenerate cases are therefore handled separately, and exactly:

• \(|h|^2 = 0\), when the Hamiltonian is proportional to the identity and \(U_3 = \mathbb{1}\);

• a doubly degenerate root \(\psi_a = \psi_b \neq \psi_c\), when the spectral decomposition reduces to a single projector and \(U_3 = e^{i\psi_a L}\mathbb{1} + (e^{i\psi_c L} - e^{i\psi_a L}) P_c\), with \(P_c = (h_k\lambda^k + \psi_a \mathbb{1})/(\psi_a - \psi_c)\).

Examples

import oscprob3nu

hamiltonian_matrix = [[1.0+0.0j, 0.0+2.0j, 0.0-1.0j],
                      [0.0-2.0j, 3.0+0.0j, 3.0+0.0j],
                      [0.0+1.0j, 3.0-0.0j, -5.0+0.0j]]
u_coeffs = oscprob3nu.evolution_operator_3nu_u_coefficients(hamiltonian_matrix,
                                                 1.0)
print('%+.6f%+.6fj' % (u_coeffs[0].real, u_coeffs[0].imag))

+0.621522-0.047327j

oscprob3nu.evolution_operator_3nu(hamiltonian_matrix: list | ndarray, L: int | float | list | ndarray) → List[List[complex]] | ndarray

Returns the three-neutrino time-evolution operator.

Returns the three-neutrino time-evolution operator \(U_3(L)\) in its SU(3) exponential expansion \(U_3(L) = u_0 \mathbb{1} + i u_k \lambda^k\). This is a \(3\times3\) unitary matrix.

Added in version 1.0.0.

Changed in version 1.1.0: Degenerate Hamiltonians are handled exactly instead of returning NaN. Lagrange interpolation over the latent roots divides by \(3\psi_m^2 - |h|^2\), which vanishes at a repeated root; the two degenerate cases are now taken in their confluent limit.

Changed in version 1.2.0: Accepts a stack of Hamiltonians of shape (..., n, n), an array of baselines, or both broadcast against each other, returning an array with the broadcast leading axes. A single Hamiltonian with a scalar baseline returns exactly what it returned before.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The sum over the three latent roots is rewritten to form its \(k\)-independent factors once, removing an (N, 3, 8) intermediate from the batched path.

Changed in version 1.5.0: Faster, with identical results; the probabilities agree with 1.4.0 to 1.6e-13 across every code path. The batched star product uses the sparse expansion rather than an numpy.einsum() contraction of the dense \(d\) tensor. The batched coefficients are laid out with the component index first, so that each \(h_k\) is contiguous; every downstream step works one component at a time, and reading those from a strided view cost about a third more.

Parameters:

hamiltonian_matrixarray_like

Three-flavor Hermitian Hamiltonian, given as the nested list [[H11, H12, H13], [H21, H22, H23], [H31, H32, H33]], or a stack of them, of shape (..., 3, 3).

Lfloat or array_like

Baseline, in units reciprocal to those of the Hamiltonian, or an array of baselines broadcastable against the leading axes of hamiltonian_matrix.

Returns:

list of list of complex or numpy.ndarray

For a single Hamiltonian and baseline, the time-evolution operator \(U_3(L)\) — a \(3\times3\) unitary complex matrix — as a nested list. If either argument is a stack, an array of shape (..., 3, 3).

See also

probabilities_3nu

Returns the probabilities built from this matrix.

Examples

import oscprob3nu

hamiltonian_matrix = [[1.0+0.0j, 0.0+2.0j, 0.0-1.0j],
                      [0.0-2.0j, 3.0+0.0j, 3.0+0.0j],
                      [0.0+1.0j, 3.0-0.0j, -5.0+0.0j]]
U3 = oscprob3nu.evolution_operator_3nu(hamiltonian_matrix, 1.0)
for row in U3:
    print('  '.join(['%+.6f%+.6fj' % (z.real+0.0, z.imag+0.0)
                     for z in row]))

+0.546090-0.496423j  -0.600964-0.114920j  -0.278885+0.056655j
+0.600964+0.114920j  +0.430462+0.614384j  -0.171381+0.183031j
+0.278885-0.056655j  -0.171381+0.183031j  +0.888015-0.259943j

oscprob3nu.probabilities_3nu(hamiltonian_matrix: list | ndarray, L: int | float | list | ndarray) → Tuple[float, ...] | ndarray

Returns the three-neutrino oscillation probabilities.

Returns the three-neutrino flavor-transition probabilities \(P_{ee}, P_{e\mu}, P_{e\tau}, P_{\mu e}, P_{\mu\mu}, P_{\mu\tau}, P_{\tau e}, P_{\tau\mu}, P_{\tau\tau}\), where \(P_{\alpha\beta} \equiv P(\nu_\alpha \to \nu_\beta) = |[U_3]_{\beta\alpha}|^2\).

Added in version 1.0.0.

Changed in version 1.1.0: Degenerate Hamiltonians are handled exactly instead of returning NaN. Lagrange interpolation over the latent roots divides by \(3\psi_m^2 - |h|^2\), which vanishes at a repeated root; the two degenerate cases are now taken in their confluent limit.

Changed in version 1.2.0: Accepts a stack of Hamiltonians of shape (..., n, n), an array of baselines, or both broadcast against each other, returning an array with the broadcast leading axes. A single Hamiltonian with a scalar baseline returns exactly what it returned before.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The sum over the three latent roots is rewritten so that the two factors independent of \(k\) are formed once rather than inside the loop over the eight \(k\), which also removes an (N, 3, 8) intermediate from the batched path; the star product is computed once per call rather than twice; and \(|z|^2\) skips the square root that squaring undoes. A scalar three-flavor probability is 3.2x quicker.

Changed in version 1.5.0: Faster, with identical results; the probabilities agree with 1.4.0 to 1.6e-13 across every code path. The batched star product no longer contracts the dense 8x8x8 \(d\) tensor through numpy.einsum(), which with no path plan walked the whole table for every element and was 70% of a 2000-point energy scan; the sparse expansion the scalar path already used vectorises unchanged. The routine also forms the nine entries and squares them directly instead of building and transposing the evolution operator. The batched coefficients are laid out with the component index first, so that each \(h_k\) is contiguous; every downstream step works one component at a time, and reading those from a strided view cost about a third more. The scalar path is 3.2x quicker, a 2000-point energy scan 3.5x, and an oscillogram 1.8x.

Changed in version 1.6.0: Two dispatch decisions, neither changing the result. A stack of at most SMALL_BATCH elements is evaluated one at a time through the scalar path, because a batched call carries a fixed cost that a handful of points does not amortise. Larger stacks are evaluated by a compiled kernel instead, if the optional numba extra is installed; for three flavors that kernel wins at every size. See fastkernels.

Parameters:

hamiltonian_matrixarray_like

Three-flavor Hermitian Hamiltonian, given as the nested list [[H11, H12, H13], [H21, H22, H23], [H31, H32, H33]], or a stack of them, of shape (..., 3, 3).

Lfloat or array_like

Baseline, in units reciprocal to those of the Hamiltonian, or an array of baselines broadcastable against the leading axes of hamiltonian_matrix.

Returns:

tuple of float or numpy.ndarray

For a single Hamiltonian and baseline, the nine probabilities (Pee, Pem, Pet, Pme, Pmm, Pmt, Pte, Ptm, Ptt) as a tuple, ordered with the initial flavor varying slowest. If either argument is a stack, an array of shape (..., 9) in the same order, with the leading axes given by broadcasting the two arguments together.

Notes

Passing arrays evaluates the whole stack at once, which is between one and two orders of magnitude faster than calling this routine in a Python loop. The two common scans both broadcast naturally:

• versus baseline, with one Hamiltonian and an array of baselines. The characteristic equation depends only on the Hamiltonian, so it is solved once for the whole scan;

• versus energy, with an array of Hamiltonians — one per energy, since \(H \propto 1/E\) — and a single baseline.

An oscillogram is the outer combination of the two, obtained by giving the Hamiltonians and the baselines separate axes, e.g. probabilities_3nu(H[:, None, :, :], L[None, :]), which returns an array of shape (len(H), len(L), 9).

Examples

import oscprob3nu

hamiltonian_matrix = [[1.0+0.0j, 0.0+2.0j, 0.0-1.0j],
                      [0.0-2.0j, 3.0+0.0j, 3.0+0.0j],
                      [0.0+1.0j, 3.0-0.0j, -5.0+0.0j]]
prob = oscprob3nu.probabilities_3nu(hamiltonian_matrix, 1.0)
print('  '.join(['%.6f' % p for p in prob[0:3]]))
print('  '.join(['%.6f' % p for p in prob[3:6]]))
print('  '.join(['%.6f' % p for p in prob[6:9]]))

0.544650  0.374364  0.080986
0.374364  0.562764  0.062872
0.080986  0.062872  0.856142

oscprob4nu

The \(n = 4\) member, and the last one that exists in closed form; see Why the method stops at four. Its accuracy on stiff 3+1 spectra deserves a look before use: Stiff spectra, and what they cost.

Compute the four-neutrino flavor-transition probabilities.

This module contains the routines needed to compute four-neutrino flavor-transition probabilities for an arbitrary time-independent \(4\times4\) Hermitian Hamiltonian, using the SU(4) exponential expansion. It is the \(n = 4\) member of the family whose \(n = 2\) and \(n = 3\) members are oscprob2nu and oscprob3nu, following the method of [1]. It is the last one: see Why the method stops at four in the methodology page. The Cayley-Hamilton route at \(n = 4\) in constant matter was set out in [2].

The Hamiltonian is expanded in the basis of the fifteen generalized Gell-Mann matrices,

\[H = h_0 \mathbb{1} + h_a \lambda^a , \qquad a = 1, \ldots, 15 ,\]

and the time-evolution operator in the same basis,

\[U_4(L) = u_0 \mathbb{1} + i u_a \lambda^a .\]

The term \(h_0\) contributes only an overall phase and is dropped; all routines therefore work with the traceless part of the Hamiltonian, which leaves the oscillation probabilities unchanged.

What is new at four flavors

Three ingredients have no counterpart in oscprob3nu:

• A third invariant. SU(4) has rank three, so the traceless part carries three independent invariants rather than two,

  \[I_2 = \tfrac12 \mathrm{Tr}\,\tilde{H}^2 , \qquad I_3 = \tfrac12 \mathrm{Tr}\,\tilde{H}^3 , \qquad I_4 = \tfrac12 \left(\mathrm{Tr}\,\tilde{H}^4 - I_2^2\right) ,\]

  and the cubic characteristic equation of the three-flavor case becomes the quartic

  \[\psi^4 - I_2 \psi^2 - \tfrac23 I_3 \psi + \tfrac14 \left(I_2^2 - 2 I_4\right) = 0 .\]

• A quartic that still solves in closed form. Euler’s method reduces it to a resolvent cubic whose three roots are real and non-negative because \(\tilde{H}\) is Hermitian, so the same trigonometric formula that oscprob3nu.psi_roots() uses solves it. The SU(3) machinery is literally nested inside the SU(4) solution.

• A longer star-product tower. The three-flavor identity \((h \star h) \star h = \tfrac13 |h|^2 h\) is a Cayley-Hamilton accident of \(n = 3\) and is false for SU(4) — it is off by some tens of per cent on a random Hamiltonian — so \(((h \star h) \star h)_a\) enters as independent data.

Accuracy

For a generic Hamiltonian the expansion is exact to round-off, like its two- and three-flavor counterparts. A stiff spectrum is the case that needs care, and the physically interesting 3+1 scenario is stiff: with \(\Delta m^2_{41} \sim 1\) eV² the eigenvalues span four orders of magnitude, three of them clustering, and the information needed to separate the cluster is destroyed when \(I_2, I_3, I_4\) are formed in double precision. No amount of care in solving the quartic recovers it.

The roots are therefore refined against the matrix, by one Newton step on \(\chi(\psi) = \det(\psi \mathbb{1} - \tilde{H})\), which uses the Hamiltonian entries rather than the three compressed invariants. That restores the roots to round-off and the probabilities from about \(5 \times 10^{-7}\) to \(10^{-9}\).

Neither figure is anywhere near a measurable effect — probabilities are confronted with data at the per-cent level — so this matters for the exactness claim, for error accumulating when slabs and earth compose operators across layers, and for a regression suite tight enough to catch a real mistake. POLISH_ROOTS has the measured comparison against the alternatives, including why LAPACK and extended precision both lose.

evolution_operator_4nu and probabilities_4nu accept either a single Hamiltonian and baseline or a stack of them, in which case the whole stack is evaluated at once, exactly as at two and three flavors.

Units

The routines are unit-agnostic: they require only that the Hamiltonian and the baseline be given in reciprocal units, so that the product \(H L\) is dimensionless. Elsewhere in NuOscProbExact the Hamiltonian is in eV and the baseline in eV^-1; the module globaldefs provides CONV_KM_TO_INV_EV to convert a baseline in km into eV^-1.

Routine listings

• generators_su4 - Returns the fifteen generalized Gell-Mann matrices

• hamiltonian_4nu_coefficients - Returns the \(h_a\)

• su4_invariants - Returns the invariants \(I_2, I_3, I_4\)

• psi_roots_4nu - Returns the roots of the characteristic equation

• evolution_operator_4nu_u_coefficients - Returns the \(u_a\)

• evolution_operator_4nu - Returns the evolution operator \(U_4\)

• probabilities_4nu - Returns the oscillation probabilities

References

[1]

Mauricio Bustamante, “Exact neutrino oscillation probabilities with arbitrary time-independent Hamiltonians”, arXiv:1904.12391.

[2]

S. Kamo et al., “Matter enhanced transitions of active and sterile neutrinos”, Eur. Phys. J. C 28, 211 (2003), which applies the Cayley-Hamilton approach at \(n = 4\) in constant matter.

oscprob4nu.CHECK_HERMITICITY = True

bool: Module-level switch.

Whether to verify that the Hamiltonian handed in is Hermitian before evaluating anything. It is on by default, and the reason is that the failure it catches is silent: a non-Hermitian matrix does not raise, and does not produce obviously broken output either — the probabilities it returns still sum to one, so the usual sanity check a caller would apply cannot tell that anything went wrong. The expansion assumes Hermiticity; without it the result is meaningless rather than merely inaccurate.

Checking is not free, and the cost is stated here rather than buried, because it is larger than one might expect. Validating a stack is a pass over it, which is the same order of work as evaluating it — and the compiled kernel has made evaluating it fast, so on a large stack the check dominates. Measured by interleaving the two settings and taking the best of fifteen rounds each:

Stack	2 flavors	4 flavors
2 000 points	1.5x	1.3x
200 000 points	5.7x	3.2x

Three flavors sits between them, at 1.8x and 3.9x. So a scan that the compiled backend was installed to speed up can spend most of its time here instead. On a single matrix the check costs about 1.35x, and it takes its own branch to manage that: the reductions above are all fixed cost at one element, so without the branch one probability spent nine tenths of its time here. For production scans whose Hamiltonians come from a construction already known to be Hermitian — everything hamiltonians4nu builds is Hermitian to round-off, as the table below records — set this to False.

The default is nevertheless True, because the alternative is a library that silently returns meaningless numbers to anyone who makes this mistake once, and finding out costs far more than the check does.

The tolerance is relative to the largest entry of the Hamiltonian, so a matrix assembled in floating point passes: the ones this library builds are Hermitian to about \(2 \times 10^{-17}\) relative, against a tolerance of \(10^{-12}\).

Added in version 1.11.0.

oscprob4nu.LAMBDA_SU4 = array([[[ 0.        +0.j,  1.        +0.j,  0.        +0.j,           0.        +0.j],         [ 1.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j, -0.        -1.j,  0.        +0.j,           0.        +0.j],         [ 0.        +1.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  1.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 1.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j, -0.        -1.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +1.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,           1.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 1.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,          -0.        -1.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +1.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  1.        +0.j,           0.        +0.j],         [ 0.        +0.j,  1.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j, -0.        -1.j,           0.        +0.j],         [ 0.        +0.j,  0.        +1.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           1.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  1.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,          -0.        -1.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +1.j,  0.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           1.        +0.j],         [ 0.        +0.j,  0.        +0.j,  1.        +0.j,           0.        +0.j]],         [[ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,          -0.        -1.j],         [ 0.        +0.j,  0.        +0.j,  0.        +1.j,           0.        +0.j]],         [[ 1.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j, -1.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.57735027+0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.57735027+0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j, -1.15470054+0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j]],         [[ 0.40824829+0.j,  0.        +0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.40824829+0.j,  0.        +0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.40824829+0.j,           0.        +0.j],         [ 0.        +0.j,  0.        +0.j,  0.        +0.j,          -1.22474487+0.j]]])

numpy.ndarray: Module-level constant.

The fifteen generalized Gell-Mann matrices, of shape (15, 4, 4), tabulated once at import time. At three flavors the analogous table is the \(d\) tensor of oscprob3nu.tensor_d(); here the generators themselves are tabulated, because every quantity this module needs is a trace against them and none needs \(d_{abc}\) explicitly.

oscprob4nu.POLISH_ROOTS = True

bool: Module-level switch.

Whether to refine the latent roots against the Hamiltonian matrix, by one Newton step on \(\chi(\psi) = \det(\psi \mathbb{1} - \tilde{H})\), after solving the quartic in closed form.

This is on by default and should stay on. The closed-form roots are limited by the conditioning of \(I_2, I_3, I_4\), not by the quartic solver: perturbing the three invariants at the \(10^{-16}\) level moves the roots of a stiff 3+1 spectrum by \(6 \times 10^{-11}\) relative, which is what the unrefined closed form achieves and what no better root-finder can beat. The Newton step reads the matrix entries instead, and is not subject to that floor.

Measured against mpmath at fifty decimal digits, on stiff 3+1 Hamiltonians, with cost quoted for a 200 000-point scan:

Strategy for the latent roots	Rel. error	Cost
Closed form alone	8.3e-11	0.17 s
Closed form + one Newton step	1.1e-16	0.41 s
numpy.linalg.eigvalsh	7.4e-16	0.17 s
Closed form in numpy.longdouble	4.5e-11	0.43 s

Note the second row beats the third: the Newton step is some seven times more accurate than LAPACK, because eigvalsh reduces by similarity transforms that each carry a backward error of order \(\epsilon\|H\|\), while this converges onto the root of \(\det(\psi\mathbb{1} - \tilde{H})\) for the matrix it was given. Extended precision was rejected for buying under a digit, being slower, and silently being float64 on Apple Silicon and Windows.

The refined figure has been confirmed from outside: against nuSQuIDS, which integrates the density matrix numerically, the four-flavor probabilities agree to \(4 \times 10^{-16}\) on a benign spectrum and \(3 \times 10^{-10}\) on the stiffest one tested. See notebook 17.

In probabilities the difference is \(5 \times 10^{-7}\) unrefined against \(10^{-9}\) refined. Both are orders of magnitude below what any experiment resolves; the reasons to want the smaller one are the library’s claim to exactness, error accumulating when slabs and earth compose operators across many layers, and a regression suite with no room for a bug to hide in.

A second Newton step changes nothing — one already reaches the floor — so exactly one is taken. It costs about 40% of the runtime, which brings the four-flavor closed form to parity with a batched eigh rather than ahead of it.

The step is applied unconditionally rather than only where a spectrum looks stiff, and that is a measured decision rather than a lazy one. Two skip criteria were tried on 6300 Hamiltonians: the gap-based amplification that perturbation theory suggests, which misjudges doubly paired spectra by four orders of magnitude, and a matrix residual comparing \(\prod_m \psi_m\) with \(\det \tilde{H}\), which is one constraint on four roots and misses errors that cancel in the product. Neither can safely skip any elements at all. The reason is structural: a criterion complete enough to certify four roots must evaluate \(\chi\) at four roots, which is this refinement — the check and the fix are the same computation. See the methodology page.

Set to False to skip it, which is useful only for reproducing the unrefined figures or for spectra known to be well separated.

oscprob4nu.SMALL_BATCH = 10

int: Module-level constant.

Nothing reads this. At two and three flavors the constant of the same name selects between a scalar path and a batched one; four flavors has no separate scalar closed form to select, so every stack goes through the array path whatever its length, and there is no short-stack case for a threshold to name. It is exported and documented as being what it is, rather than quietly kept.

It also no longer matches oscprob3nu.SMALL_BATCH, which it once did and was documented as doing: that one was re-measured when the scalar path was given a cheap Hermiticity check, and this one, governing nothing, was not.

oscprob4nu.generators_su4() → ndarray

Returns the fifteen generalized Gell-Mann matrices of SU(4).

The generators are traceless, Hermitian, and normalized as \(\mathrm{Tr}(\lambda^a \lambda^b) = 2 \delta^{ab}\), which is the convention oscprob3nu uses at \(n = 3\).

They are ordered as six symmetric and six antisymmetric off-diagonal matrices, one pair for each of the six index pairs \((j, k)\) with \(j < k\) in lexicographic order, followed by the three diagonal matrices that span the Cartan subalgebra.

Added in version 1.9.0.

Returns:

numpy.ndarray

Complex array of shape (15, 4, 4), the generators \(\lambda^1, \ldots, \lambda^{15}\).

Examples

import numpy as np

import oscprob4nu

lam = oscprob4nu.generators_su4()
gram = np.einsum('aij,bji->ab', lam, lam).real/2.0

print('shape:', lam.shape)
print('orthonormal to %.1e' % np.max(np.abs(gram - np.eye(15))))

shape: (15, 4, 4)
orthonormal to 1.1e-16

oscprob4nu.hamiltonian_4nu_coefficients(hamiltonian_matrix: list | ndarray) → List[float]

Returns the \(h_a\) of the SU(4) expansion of the Hamiltonian.

Computes the coefficients \(h_1, \ldots, h_{15}\) of the SU(4) expansion \(H = h_0 \mathbb{1} + h_a \lambda^a\) of the four-flavor Hamiltonian hamiltonian_matrix, which is assumed to be given in the flavor basis. The coefficient \(h_0\) contributes only an overall phase to the evolution operator and is not returned.

They follow from the normalization of the generators as \(h_a = \tfrac12 \mathrm{Tr}(H \lambda^a)\), and are real for a Hermitian Hamiltonian.

Added in version 1.9.0.

Parameters:

hamiltonian_matrixarray_like

Four-flavor Hermitian Hamiltonian, given as a nested list or an array of shape (4, 4).

Returns:

list of float

The fifteen coefficients \(h_1, \ldots, h_{15}\).

Examples

import numpy as np

import oscprob4nu

hamiltonian = np.diag([1.0, 2.0, 3.0, -6.0]).astype(complex)
h = oscprob4nu.hamiltonian_4nu_coefficients(hamiltonian)

print('%d coefficients' % len(h))
print('the three Cartan ones: %.4f, %.4f, %.4f' % tuple(h[12:]))

15 coefficients
the three Cartan ones: -0.5000, -0.8660, 4.8990

oscprob4nu.su4_invariants(hamiltonian_matrix: list | ndarray) → Tuple[float, float, float]

Returns the three SU(4) invariants of the Hamiltonian.

Returns \(I_2 = |h|^2\), \(I_3 = \langle h \rangle\) and \(I_4 = |h \star h|^2\) of the traceless part \(\tilde{H}\) of hamiltonian_matrix, computed from traces of its powers,

\[I_2 = \tfrac12 \mathrm{Tr}\,\tilde{H}^2 , \qquad I_3 = \tfrac12 \mathrm{Tr}\,\tilde{H}^3 , \qquad I_4 = \tfrac12 \left(\mathrm{Tr}\,\tilde{H}^4 - I_2^2\right) .\]

SU(4) has rank three, so there are three of them; SU(3) has rank two, which is why oscprob3nu.su3_invariants() returns two. Taking them from traces avoids ever building the \(d\) tensor of SU(4), a \(15 \times 15 \times 15\) table that nothing else here needs.

Added in version 1.9.0.

Parameters:

hamiltonian_matrixarray_like

Four-flavor Hermitian Hamiltonian, given as a nested list or an array of shape (4, 4).

Returns:

tuple of float

The invariants (I2, I3, I4).

Examples

import numpy as np

import oscprob4nu

hamiltonian = np.diag([1.0, 2.0, 3.0, -6.0]).astype(complex)
i2, i3, i4 = oscprob4nu.su4_invariants(hamiltonian)

print('I2 = %.4f' % i2)
print('I3 = %.4f' % i3)
print('I4 = %.4f' % i4)

I2 = 25.0000
I3 = -90.0000
I4 = 384.5000

oscprob4nu.psi_roots_4nu(invariant_2: int | float, invariant_3: int | float, invariant_4: int | float) → List[float]

Returns the roots of the four-flavor characteristic equation.

Returns the four real roots \(\psi_m\) of

\[\psi^4 - I_2 \psi^2 - \tfrac23 I_3 \psi + \tfrac14 \left(I_2^2 - 2 I_4\right) = 0 ,\]

the characteristic equation of the traceless part of the Hamiltonian, obtained by Euler’s solution of the quartic: the resolvent cubic

\[z^3 - 2 I_2 z^2 + 2 I_4 z - \tfrac49 I_3^2 = 0\]

has roots \(z_i = (\psi_i + \psi_j)^2\), real and non-negative because the Hamiltonian is Hermitian, and then

\[\psi_m = \tfrac12 \left(s_1 \sqrt{z_1} + s_2 \sqrt{z_2} + s_3 \sqrt{z_3}\right) ,\]

with signs \(s_i = \pm1\) fixed by \(s_1 s_2 s_3 \sqrt{z_1 z_2 z_3} = \tfrac23 I_3\).

These roots carry only the accuracy that \(I_2, I_3, I_4\) carry. For a stiff spectrum that is not enough, which is why the routines that use them refine them against the matrix; see POLISH_ROOTS.

Added in version 1.9.0.

Parameters:

invariant_2int or float

The invariant \(I_2 = |h|^2\).

invariant_3int or float

The invariant \(I_3 = \langle h \rangle\).

invariant_4int or float

The invariant \(I_4 = |h \star h|^2\).

Returns:

list of float

The four roots, in ascending order.

Examples

import numpy as np

import oscprob4nu

hamiltonian = np.diag([1.0, 2.0, 3.0, -6.0]).astype(complex)
i2, i3, i4 = oscprob4nu.su4_invariants(hamiltonian)
psi = oscprob4nu.psi_roots_4nu(i2, i3, i4)

print('roots    :', np.round(psi, 6))
print('reference:', np.round(np.linalg.eigvalsh(hamiltonian) - 0.0, 6))

roots    : [-6.  1.  2.  3.]
reference: [-6.  1.  2.  3.]

oscprob4nu.evolution_operator_4nu_u_coefficients(hamiltonian_matrix: list | ndarray, L: int | float) → List[complex]

Returns the \(u_a\) of the SU(4) expansion of \(U_4\).

Returns the sixteen coefficients \(u_0, u_1, \ldots, u_{15}\) of the expansion \(U_4(L) = u_0 \mathbb{1} + i u_a \lambda^a\) of the time-evolution operator, in the same convention as oscprob3nu.evolution_operator_3nu_u_coefficients().

The four-flavor analogue of Eqs. (10)-(11) of arXiv:1904.12391 is

\[u_0 = \frac14 \sum_m e^{-i \psi_m L} , \qquad i u_a = \sum_m e^{-i \psi_m L}\, \frac{\left(\psi_m^2 - \tfrac12 I_2\right) h_a + \psi_m (h \star h)_a + ((h \star h) \star h)_a} {\chi'(\psi_m)} ,\]

with \(\chi'(\psi_m) = 4\psi_m^3 - 2 I_2 \psi_m - \tfrac23 I_3\). Note the third term in the numerator: at three flavors the tower closes on itself and no such term appears.

Added in version 1.9.0.

Changed in version 1.10.1: Results changed for a spectrum with two nearly coincident latent roots. The Newton refinement of POLISH_ROOTS divides by a product of gaps, and was guarded only against a gap of exactly zero; a pair separated by one unit in the last place passed that guard and was thrown across the spectrum. The step is now refused whenever it would carry a root more than halfway to its nearest neighbour. Nothing else moved: ordinary and stiff spectra are bit-for-bit what 1.10.0 gave.

Parameters:

hamiltonian_matrixarray_like

Four-flavor Hermitian Hamiltonian, given as a nested list or an array of shape (4, 4).

Lint or float

Baseline, in units reciprocal to those of the Hamiltonian.

Returns:

list of complex

The sixteen coefficients [u0, u1, ..., u15].

Examples

import numpy as np

import oscprob4nu

hamiltonian = np.diag([1.0, 2.0, 3.0, -6.0]).astype(complex)
u = oscprob4nu.evolution_operator_4nu_u_coefficients(hamiltonian, 1.0)

print('%d coefficients' % len(u))
print('u0 = %.6f%+.6fj' % (u[0].real, u[0].imag))

16 coefficients
u0 = 0.023583-0.542826j

oscprob4nu.evolution_operator_4nu(hamiltonian_matrix: list | ndarray, L: int | float | list | ndarray) → List[List[complex]] | ndarray

Returns the four-neutrino time-evolution operator \(U_4\).

Returns \(U_4(L) = e^{-i \tilde{H} L}\), with \(\tilde{H}\) the traceless part of hamiltonian_matrix. The discarded trace contributes only the overall phase \(e^{-i h_0 L}\), which cancels in the probabilities.

Added in version 1.9.0.

Changed in version 1.10.1: Results changed for a spectrum with two nearly coincident latent roots. The Newton refinement of POLISH_ROOTS divides by a product of gaps, and was guarded only against a gap of exactly zero; a pair separated by one unit in the last place passed that guard and was thrown across the spectrum. The step is now refused whenever it would carry a root more than halfway to its nearest neighbour. Nothing else moved: ordinary and stiff spectra are bit-for-bit what 1.10.0 gave.

Parameters:

hamiltonian_matrixarray_like

Four-flavor Hermitian Hamiltonian, given as a nested list of shape (4, 4), or a stack of them, of shape (..., 4, 4).

Lint or float or array_like

Baseline, in units reciprocal to those of the Hamiltonian, or an array of baselines broadcastable against the leading axes of hamiltonian_matrix.

Returns:

list of list of complex or numpy.ndarray

For a single Hamiltonian and baseline, the \(4\times4\) evolution operator as a nested list. If either argument is a stack, an array of shape (..., 4, 4).

Examples

import numpy as np

import oscprob4nu

hamiltonian = np.diag([1.0, 2.0, 3.0, -6.0]).astype(complex)
operator = oscprob4nu.evolution_operator_4nu(hamiltonian, 1.0)

unitarity = np.asarray(operator).conj().T @ np.asarray(operator)
print('|U| unitary to %.1e' % np.max(np.abs(unitarity - np.eye(4))))

|U| unitary to 3.3e-16

oscprob4nu.probabilities_4nu(hamiltonian_matrix: list | ndarray, L: int | float | list | ndarray) → Tuple[float, ...] | ndarray

Returns the four-neutrino oscillation probabilities.

Returns the sixteen flavor-transition probabilities \(P_{\alpha\beta} \equiv P(\nu_\alpha \to \nu_\beta) = |[U_4]_{\beta\alpha}|^2\), ordered with the initial flavor varying slowest, exactly as at two and three flavors. With the fourth state read as sterile, the flavor order is \((\nu_e, \nu_\mu, \nu_\tau, \nu_s)\).

Added in version 1.9.0.

Changed in version 1.10.1: Results changed for a spectrum with two nearly coincident latent roots. The Newton refinement of POLISH_ROOTS divides by a product of gaps, and was guarded only against a gap of exactly zero; a pair separated by one unit in the last place passed that guard and was thrown across the spectrum. The step is now refused whenever it would carry a root more than halfway to its nearest neighbour. Nothing else moved: ordinary and stiff spectra are bit-for-bit what 1.10.0 gave.

Parameters:

hamiltonian_matrixarray_like

Four-flavor Hermitian Hamiltonian, given as a nested list of shape (4, 4), or a stack of them, of shape (..., 4, 4).

Lint or float or array_like

Baseline, in units reciprocal to those of the Hamiltonian, or an array of baselines broadcastable against the leading axes of hamiltonian_matrix.

Returns:

tuple of float or numpy.ndarray

For a single Hamiltonian and baseline, the sixteen probabilities as a tuple, ordered \(P_{ee}, P_{e\mu}, P_{e\tau}, P_{es}, P_{\mu e}, \ldots, P_{ss}\). If either argument is a stack, an array of shape (..., 16) in the same order.

Notes

Passing arrays evaluates the whole stack at once, which is far faster than calling this routine in a Python loop, exactly as at two and three flavors; see Scanning: pass arrays.

The accuracy deserves a word, because four flavors is where it stops being automatic. Against numpy.linalg.eigh(), a generic Hamiltonian agrees to about \(2 \times 10^{-13}\). A stiff 3+1 spectrum, with \(\Delta m^2_{41}\) four orders of magnitude above \(\Delta m^2_{21}\), agrees to about \(10^{-9}\) — limited not by this expansion but by what double precision retains when the invariants are formed. See POLISH_ROOTS, which is what keeps that figure from being \(5 \times 10^{-7}\).

Examples

import numpy as np

import oscprob4nu

hamiltonian = np.diag([1.0, 2.0, 3.0, -6.0]).astype(complex)
prob = oscprob4nu.probabilities_4nu(hamiltonian, 1.0)

print('%d probabilities' % len(prob))
print('P_ee = %.6f' % prob[0])
print('they sum to %.6f' % sum(prob[0:4]))

16 probabilities
P_ee = 1.000000
they sum to 1.000000

Sample Hamiltonians

These build the Hamiltonian for a number of standard scenarios. They are examples, not a limit on what can be computed: the core routines accept any Hermitian matrix.

Energies are in eV, mass-squared differences in eV², baselines in eV^-1, and potentials in eV. See Sign conventions for the sign the vacuum Hamiltonians adopt, and why it matters once a matter potential is added.

hamiltonians2nu

Compute two-neutrino Hamiltonians for selected scenarios.

This module contains routines that build the two-neutrino Hamiltonian for a number of standard scenarios — oscillations in vacuum, in matter of constant density, in matter with non-standard interactions (NSI), and in a CPT-odd Lorentz invariance-violating (LIV) background — together with the textbook oscillation formulas for vacuum and matter, which serve to validate the exact SU(2) computation performed by oscprob2nu and described in [1].

The Hamiltonians built here are meant to be passed to oscprob2nu.probabilities_2nu(). They are examples: the exact computation accepts any Hermitian \(2\times2\) Hamiltonian.

The routines that take a neutrino energy also accept an array of energies, and then return one Hamiltonian per energy, stacked along a leading axis. That stack is exactly what oscprob2nu.probabilities_2nu() expects, so a whole energy scan is two calls rather than a loop.

Units

Throughout this module,

Quantity	Units
Mass-squared difference	eV2
Neutrino energy	eV
Baseline	eV-1
Matter potential	eV
LIV eigenvalues and scale	eV

The routine hamiltonian_2nu_vacuum_energy_independent() returns the energy-independent part of the vacuum Hamiltonian, i.e. it has units of eV² and must still be divided by the neutrino energy. The module globaldefs provides CONV_KM_TO_INV_EV to convert a baseline in km into eV^-1.

Sign convention

The vacuum Hamiltonian is

\[\begin{split}H_{\rm vac} = \frac{\Delta m^2}{4E} \begin{pmatrix} -\cos 2\theta & \sin 2\theta \\ \sin 2\theta & \cos 2\theta \end{pmatrix} ,\end{split}\]

i.e. the mass eigenstate with the larger mass-squared value is the second one. This sign matters: it fixes the sign of the matter potential relative to the vacuum term, and hence whether the routines describe neutrinos (as they do) or antineutrinos. An overall sign flip of the vacuum Hamiltonian alone is invisible in vacuum but moves the Mikheyev-Smirnov-Wolfenstein resonance from neutrinos to antineutrinos.

Routine listings

• mixing_matrix_2nu - Returns the \(2\times2\) rotation matrix

• hamiltonian_2nu_vacuum_energy_independent - Returns \(H_{\rm vac}\) without the \(1/E\)

• probabilities_2nu_vacuum_std - Vacuum probabilities, standard formula

• hamiltonian_2nu_matter - Returns \(H_{\rm matter}\)

• probabilities_2nu_matter_std - Matter probabilities, standard formula

• hamiltonian_2nu_nsi - Returns \(H_{\rm NSI}\)

• hamiltonian_2nu_liv - Returns \(H_{\rm LIV}\)

References

[1]

Mauricio Bustamante, “Exact neutrino oscillation probabilities with arbitrary time-independent Hamiltonians”, arXiv:1904.12391.

hamiltonians2nu.mixing_matrix_2nu(sth: int | float) → List[List[float]]

Returns the \(2\times2\) rotation matrix.

Computes and returns the real \(2\times2\) rotation matrix parametrized by a single rotation angle \(\theta\).

Added in version 1.0.0.

Changed in version 1.1.0: \(\cos\theta\) is taken as \(\sqrt{1 - \sin^2\theta}\) rather than through the angle itself. The matrix is real, and is still returned as a nested list, as it always has been.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module.

Parameters:

sthfloat

\(\sin\theta\), with \(\theta\) in the first quadrant, so that \(\cos\theta = \sqrt{1-\sin^2\theta} \geq 0\).

Returns:

list of list of float

The rotation matrix [[cth, sth], [-sth, cth]], with cth = \(\cos\theta\).

Examples

import hamiltonians2nu

R = hamiltonians2nu.mixing_matrix_2nu(0.6)
print('%.6f  %.6f' % (R[0][0], R[0][1]))

0.800000  0.600000

hamiltonians2nu.hamiltonian_2nu_vacuum_energy_independent(sth: int | float, Dm2: int | float, compute_matrix_multiplication: bool = False) → ndarray

Returns the two-neutrino Hamiltonian for vacuum oscillations.

Computes and returns the energy-independent part of the real \(2\times2\) two-neutrino Hamiltonian for oscillations in vacuum, parametrized by a single mixing angle \(\theta\) and a single mass-squared difference \(\Delta m^2\). The Hamiltonian is \(H = \frac{1}{4} R M^2 R^T\), with \(R\) the rotation matrix and \(M^2 = \mathrm{diag}(-\Delta m^2, \Delta m^2)\) the traceless mass matrix. The multiplicative factor \(1/E\) is not applied.

Added in version 1.0.0.

Changed in version 1.1.0: The sign convention was corrected. The Hamiltonian was built from \(M^2 = \mathrm{diag}(\Delta m^2, -\Delta m^2)\), which yields the negative of the textbook Hamiltonian. In vacuum this is invisible, but it reverses the sign of the matter potential relative to the vacuum term, so results in matter, with NSI, or with LIV were the antineutrino ones. It also returns a complex numpy.ndarray rather than a nested list.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module.

Parameters:

sthfloat

\(\sin\theta\).

Dm2float

Mass-squared difference \(\Delta m^2\) [eV²].

compute_matrix_multiplicationbool, optional

If False (default), use the pre-computed closed-form expressions; if True, carry out the matrix multiplication \(R M^2 R^T\) explicitly. Both give the same result; the option exists as a cross-check.

Returns:

numpy.ndarray

The \(2\times2\) Hamiltonian [eV²], to be divided by the neutrino energy before use.

Notes

See the module-level Sign convention section: the mass eigenstate with the larger mass-squared value is the second one, so that adding a positive matter potential to the \(ee\) entry describes neutrinos.

Examples

import hamiltonians2nu

H = hamiltonians2nu.hamiltonian_2nu_vacuum_energy_independent(0.5, 1.0)
print('%.6f  %.6f' % (H[0][0].real, H[0][1].real))

-0.125000  0.216506

hamiltonians2nu.probabilities_2nu_vacuum_std(sth: int | float, Dm2: int | float, energy: int | float, L: int | float) → List[float]

Returns the 2nu vacuum probabilities, standard computation.

Returns the probabilities for two-neutrino oscillations in vacuum, computed with the standard analytical expression

\[P_{e\mu} = \sin^2 2\theta \sin^2\left(\frac{\Delta m^2 L}{4E}\right).\]

This routine exists to validate the exact SU(2) computation in oscprob2nu; the two agree to round-off.

Added in version 1.0.0.

Changed in version 1.1.0: The signature changed: the energy is now given in eV and the baseline in \(\mathrm{eV}^{-1}\), like the rest of the library, rather than in GeV and km. The rounded constants 1.27 and 2.54 that folded in the old conversion overstated every phase by 0.242%.

Parameters:

sthfloat

\(\sin\theta\).

Dm2float

Mass-squared difference \(\Delta m^2\) [eV²].

energyfloat

Neutrino energy [eV].

Lfloat

Baseline [eV^-1].

Returns:

list of float

The probabilities [Pee, Pem, Pme, Pmm].

See also

oscprob2nu.probabilities_2nu

The exact SU(2) computation.

Examples

import hamiltonians2nu

prob = hamiltonians2nu.probabilities_2nu_vacuum_std(0.5, 2.5e-3, 1.0e9, 5.0e12)
print('%.6f  %.6f' % (prob[0], prob[1]))

0.999794  0.000206

hamiltonians2nu.hamiltonian_2nu_matter(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, VCC: int | float | list | ndarray) → ndarray

Returns the two-neutrino Hamiltonian for matter oscillations.

Computes and returns the \(2\times2\) two-neutrino Hamiltonian for oscillations in matter of constant density, obtained by adding the charged-current matter potential to the \(ee\) entry of the vacuum Hamiltonian.

Added in version 1.0.0.

Changed in version 1.1.0: Results changed, following the sign-convention correction in hamiltonian_2nu_vacuum_energy_independent: this routine previously returned the antineutrino Hamiltonian when asked for the neutrino one, placing the MSW resonance on the wrong side. Returns a complex numpy.ndarray.

Changed in version 1.3.0: Accepts an array of energies, returning one Hamiltonian per energy stacked along a leading axis; the matter potential may be an array too. A scalar energy still returns a single matrix, and the results are bit-for-bit what the equivalent loop produced.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent part of the two-neutrino vacuum Hamiltonian [eV²], as returned by hamiltonian_2nu_vacuum_energy_independent. It is not modified.

energyfloat or array_like

Neutrino energy [eV], or an array of energies, in which case one Hamiltonian is returned per energy.

VCCfloat or array_like

Potential due to charged-current interactions of \(\nu_e\) with electrons [eV]. Positive for neutrinos, negative for antineutrinos. May be an array, to scan across a density profile alongside the energy.

Returns:

numpy.ndarray

The \(2\times2\) Hamiltonian [eV], of shape (2, 2) for a scalar energy and (..., 2, 2) for an array of energies.

Examples

import hamiltonians2nu

H_vac = hamiltonians2nu.hamiltonian_2nu_vacuum_energy_independent(0.5, 2.5e-3)
H = hamiltonians2nu.hamiltonian_2nu_matter(H_vac, 1.0e9, 1.0e-13)
print('%.6e' % H[0][0].real)

-2.125000e-13

hamiltonians2nu.probabilities_2nu_matter_std(sth: int | float, Dm2: int | float, VCC: int | float, energy: int | float, L: int | float) → List[float]

Returns the 2nu matter probabilities, standard computation.

Returns the probabilities for two-neutrino oscillations in matter of constant density, computed with the standard analytical expression in terms of the effective mixing angle \(\theta_m\) and effective mass-squared difference \(\Delta m^2_m\),

\[\sin^2 2\theta_m = \frac{\sin^2 2\theta} {\sin^2 2\theta + (\cos 2\theta - x)^2} , \quad x \equiv \frac{2 V_{\rm CC} E}{\Delta m^2} .\]

This routine exists to validate the exact SU(2) computation in oscprob2nu; the two agree to round-off.

Added in version 1.0.0.

Changed in version 1.1.0: The signature changed: the energy is now given in eV and the baseline in \(\mathrm{eV}^{-1}\), like the rest of the library, rather than in GeV and km. The rounded constants 1.27 and 2.54 that folded in the old conversion overstated every phase by 0.242%. The sign of \(\cos 2\theta\) is also kept, where it was previously discarded by computing it as \(\sqrt{1 - \sin^2 2\theta}\); for \(\theta > \pi/4\) that put the matter resonance on the wrong side.

Parameters:

sthfloat

\(\sin\theta\).

Dm2float

Mass-squared difference \(\Delta m^2\) [eV²].

VCCfloat

Potential due to charged-current interactions of \(\nu_e\) with electrons [eV].

energyfloat

Neutrino energy [eV].

Lfloat

Baseline [eV^-1].

Returns:

list of float

The probabilities [Pee, Pem, Pme, Pmm].

See also

oscprob2nu.probabilities_2nu

The exact SU(2) computation.

Notes

The resonance sits at \(x = \cos 2\theta\), which for \(\theta < \pi/4\) lies at positive energy, i.e. in the neutrino channel. Note that \(\cos 2\theta\) is signed: computing it as \(\sqrt{1 - \sin^2 2\theta}\) would lose the sign and misplace the resonance for \(\theta > \pi/4\).

Examples

import hamiltonians2nu

prob = hamiltonians2nu.probabilities_2nu_matter_std(0.5, 2.5e-3, 1.0e-13, 1.0e9,
                                    5.0e12)
print('%.6f  %.6f' % (prob[0], prob[1]))

0.985595  0.014405

hamiltonians2nu.hamiltonian_2nu_nsi(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, VCC: int | float | list | ndarray, eps: list | ndarray) → ndarray

Returns the two-neutrino Hamiltonian for oscillations with NSI.

Computes and returns the \(2\times2\) two-neutrino Hamiltonian for oscillations with non-standard interactions (NSI) in matter of constant density.

Added in version 1.0.0.

Changed in version 1.1.0: The imaginary part of a complex \(\epsilon_{e\mu}\) is no longer discarded: the vacuum Hamiltonian was real, so the array was float64 and the in-place addition truncated the value. Results also changed with the sign-convention correction described under hamiltonian_2nu_matter.

Changed in version 1.3.0: Accepts an array of energies, returning one Hamiltonian per energy stacked along a leading axis; the matter potential may be an array too. A scalar energy still returns a single matrix, and the results are bit-for-bit what the equivalent loop produced.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent part of the two-neutrino vacuum Hamiltonian [eV²], as returned by hamiltonian_2nu_vacuum_energy_independent. It is not modified.

energyfloat or array_like

Neutrino energy [eV], or an array of energies, in which case one Hamiltonian is returned per energy.

VCCfloat or array_like

Potential due to charged-current interactions of \(\nu_e\) with electrons [eV]. May be an array, to scan across a density profile alongside the energy.

epsarray_like

The three NSI strength parameters [eps_ee, eps_em, eps_mm], adimensional. The diagonal parameters eps_ee and eps_mm are real; the off-diagonal eps_em may be complex, and its complex conjugate is placed in the lower off-diagonal entry so that the Hamiltonian stays Hermitian.

Returns:

numpy.ndarray

The \(2\times2\) complex Hamiltonian [eV], of shape (2, 2) for a scalar energy and (..., 2, 2) for an array of energies.

Examples

import hamiltonians2nu

H_vac = hamiltonians2nu.hamiltonian_2nu_vacuum_energy_independent(0.5, 2.5e-3)
H = hamiltonians2nu.hamiltonian_2nu_nsi(H_vac, 1.0e9, 1.0e-13, [0.06, -0.06+0.03j,
                                                1.2])
print('%+.6e%+.6ej' % (H[0][1].real, H[0][1].imag))

+5.352659e-13+3.000000e-15j

hamiltonians2nu.hamiltonian_2nu_liv(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, sxi: int | float, b1: int | float, b2: int | float, Lambda: int | float) → ndarray

Returns the two-neutrino Hamiltonian for oscillations with LIV.

Computes and returns the \(2\times2\) two-neutrino Hamiltonian for oscillations in a CPT-odd Lorentz invariance-violating (LIV) background. The LIV term is \((E/\Lambda) R B_2 R^T\), with \(B_2 = \mathrm{diag}(b_1, b_2)\) and \(R\) the rotation by the angle \(\xi\) between the eigenvectors of \(B_2\) and the flavor states.

Added in version 1.0.0.

Changed in version 1.1.0: \(\cos\xi\) was computed as sqrt(1 - sxi - sxi) rather than sqrt(1 - sxi*sxi). For \(0 < \sin\xi < 1/2\) the LIV term was not a rotation at all, and for \(\sin\xi \geq 1/2\) the whole Hamiltonian became NaN. Results also changed with the sign-convention correction described under hamiltonian_2nu_matter.

Changed in version 1.3.0: Accepts an array of energies, returning one Hamiltonian per energy stacked along a leading axis. The LIV term scales with the energy rather than being added at constant strength, so it is formed per entry. A scalar energy still returns a single matrix, and the results are bit-for-bit what the equivalent loop produced.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent part of the two-neutrino vacuum Hamiltonian [eV²], as returned by hamiltonian_2nu_vacuum_energy_independent. It is not modified.

energyfloat or array_like

Neutrino energy [eV], or an array of energies, in which case one Hamiltonian is returned per energy.

sxifloat

\(\sin\xi\), with \(\xi\) the rotation angle between the space of the eigenvectors of \(B_2\) and the flavor states.

b1float

Eigenvalue \(b_1\) of the LIV operator \(B_2\) [eV].

b2float

Eigenvalue \(b_2\) of the LIV operator \(B_2\) [eV].

Lambdafloat

Energy scale \(\Lambda\) of the LIV operator \(B_2\) [eV].

Returns:

numpy.ndarray

The \(2\times2\) complex Hamiltonian [eV], of shape (2, 2) for a scalar energy and (..., 2, 2) for an array of energies.

Examples

import hamiltonians2nu

H_vac = hamiltonians2nu.hamiltonian_2nu_vacuum_energy_independent(0.5, 2.5e-3)
H = hamiltonians2nu.hamiltonian_2nu_liv(H_vac, 1.0e9, 0.6, 1.0e-9, 2.0e-9, 1.0e12)
print('%.6e' % H[0][0].real)

1.047500e-12

hamiltonians3nu

Compute three-neutrino Hamiltonians for selected scenarios.

This module contains routines that build the three-neutrino Hamiltonian for a number of standard scenarios — oscillations in vacuum, in matter of constant density, in matter with non-standard interactions (NSI), and in a CPT-odd Lorentz invariance-violating (LIV) background — together with the textbook oscillation formula in vacuum, which serves to validate the exact SU(3) computation performed by oscprob3nu and described in [1].

The Hamiltonians built here are meant to be passed to oscprob3nu.probabilities_3nu(). They are examples: the exact computation accepts any Hermitian \(3\times3\) Hamiltonian.

The routines that take a neutrino energy also accept an array of energies, and then return one Hamiltonian per energy, stacked along a leading axis. That stack is exactly what oscprob3nu.probabilities_3nu() expects, so a whole energy scan is two calls rather than a loop.

Units

Throughout this module,

Quantity	Units
Mass-squared differences	eV2
Neutrino energy	eV
Baseline	eV-1
Matter potential	eV
LIV eigenvalues and scale	eV
CP-violation phases	radian

The routine hamiltonian_3nu_vacuum_energy_independent() returns the energy-independent part of the vacuum Hamiltonian, i.e. it has units of eV² and must still be divided by the neutrino energy. The module globaldefs provides CONV_KM_TO_INV_EV to convert a baseline in km into eV^-1.

Sign convention

The vacuum Hamiltonian is \(H_{\rm vac} = U M^2 U^\dagger / (2E)\), with \(M^2 = \mathrm{diag}(0, \Delta m^2_{21}, \Delta m^2_{31})\) and \(U\) the PMNS matrix, so that adding a positive matter potential to the \(ee\) entry describes neutrinos, not antineutrinos.

Routine listings

• pmns_mixing_matrix - Returns the \(3\times3\) PMNS matrix

• hamiltonian_3nu_vacuum_energy_independent - Returns \(H_{\rm vac}\) without the \(1/E\)

• delta - Kronecker delta

• J - Product of four entries of the PMNS matrix

• probabilities_3nu_vacuum_std - Vacuum probabilities, standard formula

• hamiltonian_3nu_matter - Returns \(H_{\rm matter}\)

• hamiltonian_3nu_nsi - Returns \(H_{\rm NSI}\)

• hamiltonian_3nu_liv - Returns \(H_{\rm LIV}\)

References

[1]

Mauricio Bustamante, “Exact neutrino oscillation probabilities with arbitrary time-independent Hamiltonians”, arXiv:1904.12391.

hamiltonians3nu.pmns_mixing_matrix(s12: int | float, s23: int | float, s13: int | float, dCP: int | float) → List[List[complex]]

Returns the \(3\times3\) PMNS mixing matrix.

Computes and returns the complex \(3\times3\) PMNS mixing matrix, parametrized by the three rotation angles \(\theta_{12}\), \(\theta_{23}\), \(\theta_{13}\), and the CP-violation phase \(\delta_{\rm CP}\), in the standard PDG parametrization.

Added in version 1.0.0.

Changed in version 1.1.0: The entries are complex throughout. They were previously a mixture of complex and real numbers, depending on which of them the CP phase reached. The matrix is still returned as a nested list, as it always has been.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The scalar path stopped dispatching NumPy for single numbers: numpy.real(), numpy.imag(), numpy.arccos, numpy.clip() and numpy.sqrt on one number give way to attribute access and the math module.

Parameters:

s12float

\(\sin\theta_{12}\).

s23float

\(\sin\theta_{23}\).

s13float

\(\sin\theta_{13}\).

dCPfloat

CP-violation phase \(\delta_{\rm CP}\) [radian].

Returns:

list of list of complex

The \(3\times3\) PMNS mixing matrix, as a nested list.

Examples

import hamiltonians3nu

U = hamiltonians3nu.pmns_mixing_matrix(0.55, 0.76, 0.15, 0.0)
print('%.6f  %.6f' % (U[0][0].real, U[0][1].real))

0.825716  0.543777

hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(s12: int | float, s23: int | float, s13: int | float, dCP: int | float, D21: int | float, D31: int | float, compute_matrix_multiplication: bool = False) → ndarray

Returns the three-neutrino Hamiltonian for vacuum oscillations.

Computes and returns the energy-independent part of the complex \(3\times3\) three-neutrino Hamiltonian for oscillations in vacuum, parametrized by three mixing angles — \(\theta_{12}\), \(\theta_{23}\), \(\theta_{13}\) — one CP-violation phase — \(\delta_{\rm CP}\) — and two mass-squared differences — \(\Delta m^2_{21}\), \(\Delta m^2_{31}\). The Hamiltonian is \(H = \frac{1}{2} U M^2 U^\dagger\), with \(U\) the PMNS matrix and \(M^2 = \mathrm{diag}(0, \Delta m^2_{21}, \Delta m^2_{31})\) the mass matrix. The multiplicative factor \(1/E\) is not applied.

Added in version 1.0.0.

Changed in version 1.1.0: Returns a complex numpy.ndarray rather than a nested list.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The closed form rebuilt the CP phase fifteen times across its nine entries and recomputed two products in five places; hoisting them makes the routine 1.9x quicker.

Parameters:

s12float

\(\sin\theta_{12}\).

s23float

\(\sin\theta_{23}\).

s13float

\(\sin\theta_{13}\).

dCPfloat

CP-violation phase \(\delta_{\rm CP}\) [radian].

D21float

Mass-squared difference \(\Delta m^2_{21}\) [eV²].

D31float

Mass-squared difference \(\Delta m^2_{31}\) [eV²].

compute_matrix_multiplicationbool, optional

If False (default), use the pre-computed closed-form expressions; if True, carry out the matrix multiplication \(U M^2 U^\dagger\) explicitly. Both give the same result; the option exists as a cross-check.

Returns:

numpy.ndarray

The \(3\times3\) complex Hamiltonian [eV²], to be divided by the neutrino energy before use.

Examples

import hamiltonians3nu

H = hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(0.55, 0.76, 0.15,
                                              0.0, 7.4e-5, 2.5e-3)
print('%.6e' % H[0][0].real)

3.906567e-05

hamiltonians3nu.delta(a: int, b: int) → int

Returns the Kronecker delta \(\delta_{ab}\).

Added in version 1.0.0.

Parameters:

aint

First index.

bint

Second index.

Returns:

int

1 if a == b, 0 otherwise.

Examples

import hamiltonians3nu

print(hamiltonians3nu.delta(0, 0), hamiltonians3nu.delta(0, 1))

1 0

hamiltonians3nu.J(U: list | ndarray, alpha: int, beta: int, k: int, j: int) → complex

Returns the quartic product of PMNS matrix entries.

Returns \(J = U_{\alpha k}^* U_{\beta k} U_{\alpha j} U_{\beta j}^*\), with \(U\) the PMNS mixing matrix. This product appears in the standard expression for the three-neutrino oscillation probability in vacuum; its imaginary part is the Jarlskog invariant, up to a sign.

Added in version 1.0.0.

Parameters:

Uarray_like

The \(3\times3\) complex PMNS mixing matrix.

alphaint

Index of the initial flavor (0: \(e\), 1: \(\mu\), 2: \(\tau\)).

betaint

Index of the final flavor (0: \(e\), 1: \(\mu\), 2: \(\tau\)).

kint

First index of the sum over mass eigenstates (0, 1, or 2).

jint

Second index of the sum over mass eigenstates (0, 1, or 2).

Returns:

complex

The product \(U_{\alpha k}^* U_{\beta k} U_{\alpha j} U_{\beta j}^*\).

Examples

import hamiltonians3nu

U = hamiltonians3nu.pmns_mixing_matrix(0.55, 0.76, 0.15, 0.0)
print('%.6f' % hamiltonians3nu.J(U, 0, 1, 1, 0).real)

-0.097579

hamiltonians3nu.probabilities_3nu_vacuum_std(U: list | ndarray, D21: int | float, D31: int | float, energy: int | float, L: int | float) → List[float]

Returns the 3nu vacuum probabilities, standard computation.

Returns the probabilities for three-neutrino oscillations in vacuum, computed with the standard analytical expression

\[P_{\alpha\beta} = \delta_{\alpha\beta} - 4 \sum_{k>j} \mathrm{Re}(J_{kj}) \sin^2\left(\frac{\Delta m^2_{kj} L}{4E}\right) + 2 \sum_{k>j} \mathrm{Im}(J_{kj}) \sin\left(\frac{\Delta m^2_{kj} L}{2E}\right) .\]

This routine exists to validate the exact SU(3) computation in oscprob3nu; the two agree to round-off.

Added in version 1.0.0.

Changed in version 1.1.0: The signature changed: the energy is now given in eV and the baseline in \(\mathrm{eV}^{-1}\), like the rest of the library, rather than in GeV and km. The rounded constants 1.27 and 2.54 that folded in the old conversion overstated every phase by 0.242%.

Parameters:

Uarray_like

The \(3\times3\) complex PMNS mixing matrix, as returned by pmns_mixing_matrix.

D21float

Mass-squared difference \(\Delta m^2_{21}\) [eV²].

D31float

Mass-squared difference \(\Delta m^2_{31}\) [eV²].

energyfloat

Neutrino energy [eV].

Lfloat

Baseline [eV^-1].

Returns:

list of float

The nine probabilities [Pee, Pem, Pet, Pme, Pmm, Pmt, Pte, Ptm, Ptt], ordered with the initial flavor varying slowest.

See also

oscprob3nu.probabilities_3nu

The exact SU(3) computation.

Examples

import hamiltonians3nu

U = hamiltonians3nu.pmns_mixing_matrix(0.55, 0.76, 0.15, 0.0)
prob = hamiltonians3nu.probabilities_3nu_vacuum_std(U, 7.4e-5, 2.5e-3, 1.0e9,
                                    5.0e12)
print('%.6f  %.6f' % (prob[0], prob[1]))

0.992787  0.001981

hamiltonians3nu.hamiltonian_3nu_matter(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, VCC: int | float | list | ndarray) → ndarray

Returns the three-neutrino Hamiltonian for matter oscillations.

Computes and returns the \(3\times3\) three-neutrino Hamiltonian for oscillations in matter of constant density, obtained by adding the charged-current matter potential to the \(ee\) entry of the vacuum Hamiltonian.

Added in version 1.0.0.

Changed in version 1.1.0: Returns a complex numpy.ndarray rather than a nested list.

Changed in version 1.3.0: Accepts an array of energies, returning one Hamiltonian per energy stacked along a leading axis; the matter potential may be an array too. A scalar energy still returns a single matrix, and the results are bit-for-bit what the equivalent loop produced.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent part of the three-neutrino vacuum Hamiltonian [eV²], as returned by hamiltonian_3nu_vacuum_energy_independent. It is not modified.

energyfloat or array_like

Neutrino energy [eV], or an array of energies, in which case one Hamiltonian is returned per energy.

VCCfloat or array_like

Potential due to charged-current interactions of \(\nu_e\) with electrons [eV]. Positive for neutrinos, negative for antineutrinos. May be an array, to scan across a density profile alongside the energy.

Returns:

numpy.ndarray

The \(3\times3\) complex Hamiltonian [eV], of shape (3, 3) for a scalar energy and (..., 3, 3) for an array of energies.

Examples

import hamiltonians3nu

H_vac = hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(0.55, 0.76, 0.15,
                                                  0.0, 7.4e-5,
                                                  2.5e-3)
H = hamiltonians3nu.hamiltonian_3nu_matter(H_vac, 1.0e9, 1.0e-13)
print('%.6e' % H[0][0].real)

1.390657e-13

hamiltonians3nu.hamiltonian_3nu_nsi(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, VCC: int | float | list | ndarray, eps: list | ndarray) → ndarray

Returns the three-neutrino Hamiltonian for oscillations w/ NSI.

Computes and returns the \(3\times3\) three-neutrino Hamiltonian for oscillations with non-standard interactions (NSI) in matter of constant density.

Added in version 1.0.0.

Changed in version 1.1.0: Returns a complex numpy.ndarray rather than a nested list.

Changed in version 1.3.0: Accepts an array of energies, returning one Hamiltonian per energy stacked along a leading axis; the matter potential may be an array too. A scalar energy still returns a single matrix, and the results are bit-for-bit what the equivalent loop produced.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent part of the three-neutrino vacuum Hamiltonian [eV²], as returned by hamiltonian_3nu_vacuum_energy_independent. It is not modified.

energyfloat or array_like

Neutrino energy [eV], or an array of energies, in which case one Hamiltonian is returned per energy.

VCCfloat or array_like

Potential due to charged-current interactions of \(\nu_e\) with electrons [eV]. May be an array, to scan across a density profile alongside the energy.

epsarray_like

The six NSI strength parameters [eps_ee, eps_em, eps_et, eps_mm, eps_mt, eps_tt], adimensional. The diagonal parameters are real; the off-diagonal ones may be complex, and their complex conjugates are placed in the lower off-diagonal entries so that the Hamiltonian stays Hermitian.

Returns:

numpy.ndarray

The \(3\times3\) complex Hamiltonian [eV], of shape (3, 3) for a scalar energy and (..., 3, 3) for an array of energies.

Examples

import hamiltonians3nu

H_vac = hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(0.55, 0.76, 0.15,
                                                  0.0, 7.4e-5,
                                                  2.5e-3)
H = hamiltonians3nu.hamiltonian_3nu_nsi(H_vac, 1.0e9, 1.0e-13,
                        [0.06, -0.06+0.03j, 0.0, 1.2, 0.0, 0.0])
print('%+.6e%+.6ej' % (H[0][1].real, H[0][1].imag))

+1.445471e-13+3.000000e-15j

hamiltonians3nu.hamiltonian_3nu_liv(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, sxi12: int | float, sxi23: int | float, sxi13: int | float, dxiCP: int | float, b1: int | float, b2: int | float, b3: int | float, Lambda: int | float) → ndarray

Returns the three-neutrino Hamiltonian for oscillations w/ LIV.

Computes and returns the \(3\times3\) three-neutrino Hamiltonian for oscillations in a CPT-odd Lorentz invariance-violating (LIV) background. The LIV term is \((E/\Lambda) R B_3 R^\dagger\), with \(B_3 = \mathrm{diag}(b_1, b_2, b_3)\) and \(R\) a PMNS-like matrix built from the mixing angles \(\xi_{ij}\) and the phase \(\delta_{\xi,\rm CP}\) that relate the eigenvectors of \(B_3\) to the flavor states.

Added in version 1.0.0.

Changed in version 1.1.0: Returns a complex numpy.ndarray rather than a nested list.

Changed in version 1.3.0: Accepts an array of energies, returning one Hamiltonian per energy stacked along a leading axis. The LIV term scales with the energy rather than being added at constant strength, so it is formed per entry. A scalar energy still returns a single matrix, and the results are bit-for-bit what the equivalent loop produced.

Changed in version 1.4.0: Faster, with identical results — all 42 figures generated by run_testsuite.py are byte-for-byte those of 1.3.0. The gain is indirect, through pmns_mixing_matrix, which this routine builds the rotation from.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent part of the three-neutrino vacuum Hamiltonian [eV²], as returned by hamiltonian_3nu_vacuum_energy_independent. It is not modified.

energyfloat or array_like

Neutrino energy [eV], or an array of energies, in which case one Hamiltonian is returned per energy.

sxi12float

\(\sin\xi_{12}\), with \(\xi_{12}\) one of the mixing angles between the space of the eigenvectors of \(B_3\) and the flavor states.

sxi23float

\(\sin\xi_{23}\), likewise.

sxi13float

\(\sin\xi_{13}\), likewise.

dxiCPfloat

CP-violation phase of the LIV operator \(B_3\) [radian].

b1float

Eigenvalue \(b_1\) of the LIV operator \(B_3\) [eV].

b2float

Eigenvalue \(b_2\) of the LIV operator \(B_3\) [eV].

b3float

Eigenvalue \(b_3\) of the LIV operator \(B_3\) [eV].

Lambdafloat

Energy scale \(\Lambda\) of the LIV operator \(B_3\) [eV].

Returns:

numpy.ndarray

The \(3\times3\) complex Hamiltonian [eV], of shape (3, 3) for a scalar energy and (..., 3, 3) for an array of energies.

Examples

import hamiltonians3nu

H_vac = hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(0.55, 0.76, 0.15,
                                                  0.0, 7.4e-5,
                                                  2.5e-3)
H = hamiltonians3nu.hamiltonian_3nu_liv(H_vac, 1.0e9, 0.3, 0.4, 0.5, 0.7, 1.0e-9,
                        1.5e-9, 2.0e-9, 1.0e12)
print('%.6e' % H[0][0].real)

1.322816e-12

hamiltonians4nu

Sample four-neutrino Hamiltonians, for 3+1 scenarios.

This module builds the \(4\times4\) Hamiltonians that oscprob4nu evaluates, in the same spirit as hamiltonians3nu does at three flavors and as [1] describes: they are examples, not limitations. oscprob4nu.probabilities_4nu() takes any Hermitian \(4\times4\) matrix, so a scenario not built here is a matrix away.

The fourth state is written as sterile throughout, so the flavor basis is \((\nu_e, \nu_\mu, \nu_\tau, \nu_s)\), and the mixing matrix carries three extra angles \(\theta_{14}, \theta_{24}, \theta_{34}\) on top of the three standard ones. Nothing in oscprob4nu depends on that reading: a fourth active state, or any other four-level system, works identically.

Why 3+1 is in scope here

A 3+1 system is often described as “leaky” from the three-flavor point of view, because probability disappears from the active block into the sterile state. That is a statement about the \(3\times3\) subsystem, not about the physics: over all four states the evolution is closed and unitary, which is exactly the assumption the SU(4) expansion needs. Treating it at \(n = 4\) therefore brings it back inside the scope of an exact closed-form method.

The matter potential

Active neutrinos feel the charged-current potential \(V_{CC}\) (the electron flavor only) and the flavor-universal neutral-current potential \(V_{NC}\) (all three); a sterile state feels neither. A term proportional to the identity contributes only a global phase, so subtracting \(V_{NC}\mathbb{1}\) from all four states costs nothing and leaves

\[A_4 = \mathrm{diag}\left(V_{CC},\, 0,\, 0,\, -V_{NC}\right) ,\]

with \(-V_{NC} = +G_F n_n/\sqrt{2}\) positive. The sterile entry is therefore not zero, and getting it wrong is the four-flavor analogue of the antineutrino sign trap: the difference is invisible in vacuum and sets the position of the matter resonance.

Routine listings

• mixing_matrix_4nu - Returns the 3+1 mixing matrix

• hamiltonian_4nu_vacuum_energy_independent - Vacuum Hamiltonian

• hamiltonian_4nu_matter - Adds matter of constant density

• hamiltonian_4nu_nsi - Adds non-standard interactions

• hamiltonian_4nu_liv - Adds a Lorentz invariance-violating term

References

[1]

Mauricio Bustamante, “Exact neutrino oscillation probabilities with arbitrary time-independent Hamiltonians”, arXiv:1904.12391.

hamiltonians4nu.mixing_matrix_4nu(s12: int | float, s23: int | float, s13: int | float, s14: int | float, s24: int | float, s34: int | float, dCP: int | float, d14: int | float = 0.0, d24: int | float = 0.0) → ndarray

Returns the 3+1 lepton mixing matrix.

Built in the common 3+1 ordering

\[U = R_{34} R_{24}(\delta_{24}) R_{14}(\delta_{14}) R_{23} R_{13}(\delta_{CP}) R_{12} ,\]

which reduces to the standard PDG three-flavor matrix of hamiltonians3nu.pmns_mixing_matrix() in the upper-left block when the three new angles vanish.

All six mixing parameters are sines of the angles, not the angles, matching the convention used throughout NuOscProbExact.

Added in version 1.9.0.

Parameters:

s12int or float

Sine of \(\theta_{12}\).

s23int or float

Sine of \(\theta_{23}\).

s13int or float

Sine of \(\theta_{13}\).

s14int or float

Sine of \(\theta_{14}\).

s24int or float

Sine of \(\theta_{24}\).

s34int or float

Sine of \(\theta_{34}\).

dCPint or float

Standard Dirac CP-violating phase, in radian.

d14int or float, optional

Extra phase on the 1-4 rotation, in radian. Default: 0.

d24int or float, optional

Extra phase on the 2-4 rotation, in radian. Default: 0.

Returns:

numpy.ndarray

Complex array of shape (4, 4).

Examples

import numpy as np

import hamiltonians4nu

mixing = hamiltonians4nu.mixing_matrix_4nu(
    np.sqrt(0.310), np.sqrt(0.582), np.sqrt(2.240e-2),
    np.sqrt(0.10), np.sqrt(0.10), 0.0, 217.0/180.0*np.pi)

unitarity = mixing.conj().T @ mixing
print('unitary to %.1e' % np.max(np.abs(unitarity - np.eye(4))))

unitary to 1.7e-16

hamiltonians4nu.hamiltonian_4nu_vacuum_energy_independent(s12: int | float, s23: int | float, s13: int | float, s14: int | float, s24: int | float, s34: int | float, dCP: int | float, D21: int | float, D31: int | float, D41: int | float, d14: int | float = 0.0, d24: int | float = 0.0) → ndarray

Returns the energy-independent four-neutrino vacuum Hamiltonian.

Returns \(U M^2 U^\dagger / 2\), with \(M^2 = \mathrm{diag}(0, \Delta m^2_{21}, \Delta m^2_{31}, \Delta m^2_{41})\), so that the vacuum Hamiltonian at energy \(E\) is this matrix divided by \(E\).

The factor \(1/E\) is deliberately left out, so that the energy-independent part can be computed once and reused across an energy scan — the same arrangement as hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent().

Added in version 1.9.0.

Parameters:

s12int or float

Sine of \(\theta_{12}\).

s23int or float

Sine of \(\theta_{23}\).

s13int or float

Sine of \(\theta_{13}\).

s14int or float

Sine of \(\theta_{14}\).

s24int or float

Sine of \(\theta_{24}\).

s34int or float

Sine of \(\theta_{34}\).

dCPint or float

Standard Dirac CP-violating phase, in radian.

D21int or float

Mass-squared difference \(\Delta m^2_{21}\), in eV².

D31int or float

Mass-squared difference \(\Delta m^2_{31}\), in eV².

D41int or float

Mass-squared difference \(\Delta m^2_{41}\), in eV².

d14int or float, optional

Extra phase on the 1-4 rotation, in radian. Default: 0.

d24int or float, optional

Extra phase on the 2-4 rotation, in radian. Default: 0.

Returns:

numpy.ndarray

Complex array of shape (4, 4), in units of eV². Divide by the energy in eV to obtain a Hamiltonian in eV.

Examples

import numpy as np

import globaldefs as gd
import hamiltonians4nu
import oscprob4nu

h_vacuum = 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(
    h_vacuum/1.0e9, 1300.0*gd.CONV_KM_TO_INV_EV)

print('P_mumu = %.6f' % prob[5])
print('P_mus  = %.6f' % prob[7])

P_mumu = 0.407166
P_mus  = 0.011663

hamiltonians4nu.hamiltonian_4nu_matter(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, VCC: int | float | list | ndarray, VNC: int | float | list | ndarray) → ndarray

Returns the four-neutrino Hamiltonian in matter.

Adds to the vacuum term the matter potential

\[A_4 = \mathrm{diag}\left(V_{CC},\, 0,\, 0,\, -V_{NC}\right) ,\]

which is what remains after the flavor-universal neutral-current potential of the three active states is removed as a global phase. Because the sterile state does not feel \(V_{NC}\), removing it leaves \(-V_{NC}\) on the sterile entry rather than nothing.

VCC is positive for neutrinos. For antineutrinos, reverse the sign of both potentials and conjugate the vacuum term, exactly as at three flavors.

Added in version 1.9.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent four-flavor vacuum Hamiltonian, of shape (4, 4), in eV².

energyint or float or array_like

Neutrino energy, in eV, or an array of energies.

VCCint or float or array_like

Charged-current matter potential, in eV. Positive for neutrinos.

VNCint or float or array_like

Neutral-current matter potential, in eV. Negative for neutrinos, equal to \(-G_F n_n/\sqrt{2}\); see globaldefs.VNC_EARTH_CRUST.

Returns:

numpy.ndarray

Complex array of shape (4, 4) for a scalar energy, or (..., 4, 4) for an array of energies, in eV.

Examples

import numpy as np

import globaldefs as gd
import hamiltonians4nu
import oscprob4nu

h_vacuum = 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)

h_matter = hamiltonians4nu.hamiltonian_4nu_matter(
    h_vacuum, 1.0e9, gd.VCC_EARTH_CRUST, gd.VNC_EARTH_CRUST)

prob = oscprob4nu.probabilities_4nu(
    h_matter, 1300.0*gd.CONV_KM_TO_INV_EV)

print('P_ee   = %.6f' % prob[0])
print('P_mumu = %.6f' % prob[5])

P_ee   = 0.837799
P_mumu = 0.387694

hamiltonians4nu.hamiltonian_4nu_nsi(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, VCC: int | float | list | ndarray, VNC: int | float | list | ndarray, eps: list | ndarray) → ndarray

Returns the four-neutrino Hamiltonian with matter and NSI.

Adds non-standard interactions to the active block only. Sterile states have no standard-model interactions by construction, so they have no non-standard ones either: the sterile row and column of the NSI matrix are zero, and the sterile entry keeps the \(-V_{NC}\) of hamiltonian_4nu_matter().

The eps parameters follow hamiltonians3nu.hamiltonian_3nu_nsi(): six of them, with the three off-diagonal ones allowed to be complex.

Added in version 1.9.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent four-flavor vacuum Hamiltonian, of shape (4, 4), in eV².

energyint or float or array_like

Neutrino energy, in eV, or an array of energies.

VCCint or float or array_like

Charged-current matter potential, in eV.

VNCint or float or array_like

Neutral-current matter potential, in eV.

epsarray_like

The six NSI strength parameters [eps_ee, eps_em, eps_et, eps_mm, eps_mt, eps_tt], relative to VCC. The off-diagonal ones may be complex.

Returns:

numpy.ndarray

Complex array of shape (4, 4) for a scalar energy, or (..., 4, 4) for an array of energies, in eV.

Examples

import numpy as np

import globaldefs as gd
import hamiltonians4nu
import oscprob4nu

h_vacuum = 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)

h_nsi = hamiltonians4nu.hamiltonian_4nu_nsi(
    h_vacuum, 1.0e9, gd.VCC_EARTH_CRUST, gd.VNC_EARTH_CRUST,
    gd.EPS_3)

prob = oscprob4nu.probabilities_4nu(
    h_nsi, 1300.0*gd.CONV_KM_TO_INV_EV)

print('P_mue with NSI = %.6f' % prob[4])

P_mue with NSI = 0.027737

hamiltonians4nu.hamiltonian_4nu_liv(h_vacuum_energy_independent: list | ndarray, energy: int | float | list | ndarray, sxi12: int | float, sxi23: int | float, sxi13: int | float, sxi14: int | float, sxi24: int | float, sxi34: int | float, dxiCP: int | float, b1: int | float, b2: int | float, b3: int | float, b4: int | float, Lambda: int | float) → ndarray

Returns the four-neutrino Hamiltonian for oscillations w/ LIV.

The four-flavor counterpart of hamiltonians3nu.hamiltonian_3nu_liv(). The LIV term is \((E/\Lambda) R B_4 R^\dagger\), with \(B_4 = \mathrm{diag}(b_1, b_2, b_3, b_4)\) and \(R\) a mixing matrix of the same 3+1 form as mixing_matrix_4nu(), built from the angles \(\xi_{ij}\) and the phase \(\delta_{\xi,\rm CP}\) that relate the eigenvectors of \(B_4\) to the flavor states.

Nothing here privileges the fourth state: \(b_4\) is an eigenvalue like the others, so a sterile neutrino may couple to the LIV background whether or not it couples to matter. Setting the three new angles to zero and \(b_4\) equal to the trace-shifted remainder recovers the three-flavor term in the active block.

Added in version 1.11.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent four-flavor vacuum Hamiltonian, of shape (4, 4), in eV². It is not modified.

energyint or float or array_like

Neutrino energy, in eV, or an array of energies.

sxi12int or float

Sine of \(\xi_{12}\).

sxi23int or float

Sine of \(\xi_{23}\).

sxi13int or float

Sine of \(\xi_{13}\).

sxi14int or float

Sine of \(\xi_{14}\).

sxi24int or float

Sine of \(\xi_{24}\).

sxi34int or float

Sine of \(\xi_{34}\).

dxiCPint or float

CP-violation phase of the LIV operator, in radian.

b1int or float

Eigenvalue \(b_1\) of the LIV operator \(B_4\) [eV].

b2int or float

Eigenvalue \(b_2\) of the LIV operator \(B_4\) [eV].

b3int or float

Eigenvalue \(b_3\) of the LIV operator \(B_4\) [eV].

b4int or float

Eigenvalue \(b_4\) of the LIV operator \(B_4\) [eV].

Lambdaint or float

Energy scale \(\Lambda\) of the LIV operator [eV].

Returns:

numpy.ndarray

Complex array of shape (4, 4) for a scalar energy, or (..., 4, 4) for an array of energies, in eV.

Examples

import numpy as np

import globaldefs as gd
import hamiltonians4nu
import oscprob4nu

h_vacuum = 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)

h_liv = hamiltonians4nu.hamiltonian_4nu_liv(
    h_vacuum, 1.0e9, 0.3, 0.4, 0.5, 0.0, 0.0, 0.0, 0.7,
    1.0e-9, 1.5e-9, 2.0e-9, 2.5e-9, 1.0e12)

prob = oscprob4nu.probabilities_4nu(
    h_liv, 1300.0*gd.CONV_KM_TO_INV_EV)

print('P_ee with LIV = %.6f' % prob[0])

P_ee with LIV = 0.370795

Piecewise-constant matter

The exact expansions assume a Hamiltonian that does not change along the trajectory. A neutrino crossing the Earth does not have one, so its path is cut into slabs, each solved exactly and the results multiplied. Within a slab nothing is approximated; the only approximation is how finely a continuously varying profile is sliced, and that is an argument the caller controls.

slabs

Oscillation probabilities across a sequence of adjacent slabs.

NuOscProbExact computes the evolution operator exactly for a Hamiltonian that does not change along the trajectory. A great many interesting cases are piecewise constant instead: a neutrino crossing the Earth, a beam through a layered detector, a castle-wall profile built to enhance CP-violating effects. This module handles those by doing the only thing the exactness of the method allows — solving each slab exactly and multiplying the results.

For a trajectory divided into \(n\) slabs, the slab \(k\) having Hamiltonian \(H_k\) and width \(L_k\), the evolution operator is

\[U = U_n(L_n) \cdots U_2(L_2) U_1(L_1) ,\]

with the slab the neutrino meets first applied first — rightmost, since the operators act to the left on the initial state. Each \(U_k\) is the exact SU(2), SU(3) or SU(4) expansion of oscprob2nu, oscprob3nu or oscprob4nu, so the only approximation in the result is the one the caller makes in choosing how finely to slice a continuously varying profile. Within each slab there is none.

The per-slab operators are evaluated in a single batched call, so the cost of \(n\) slabs is one vectorised evaluation plus \(n-1\) small matrix products rather than \(n\) separate evaluations.

One thing carries over from the single-slab case and is easier to overlook here. The expansions return \(e^{-i H_0 L}\), with \(H_0\) the traceless part of the Hamiltonian, dropping the phase \(e^{-i h_0 L}\) that the trace contributes. Each slab therefore drops its own phase, and their product differs from \(\prod_k e^{-i H_k L_k}\) by the single scalar \(\exp(i \sum_k h_0^{(k)} L_k)\). That is still one overall phase, so every probability is unaffected — but a caller comparing the returned operator against an independent matrix exponential must compare against the traceless one, exactly as the single-slab tests do.

For the Earth specifically, earth builds the slabs for you from the Preliminary Reference Earth Model.

Routine listings

• evolution_operator_2nu_slabs - Two-flavor evolution operator

• evolution_operator_3nu_slabs - Three-flavor evolution operator

• evolution_operator_4nu_slabs - Four-flavor evolution operator

• probabilities_2nu_slabs - Two-flavor probabilities

• probabilities_3nu_slabs - Three-flavor probabilities

• probabilities_4nu_slabs - Four-flavor probabilities

slabs.evolution_operator_2nu_slabs(hamiltonian_matrices: list | ndarray, widths: list | ndarray) → ndarray

Returns the two-flavor evolution operator across adjacent slabs.

Returns the \(2\times2\) evolution operator \(U_2 = U_2^{(n)}(L_n) \cdots U_2^{(1)}(L_1)\) for a trajectory divided into slabs, each with its own constant Hamiltonian and width. Each slab is solved exactly by the SU(2) expansion of oscprob2nu.

Added in version 1.8.0.

Parameters:

hamiltonian_matricesarray_like

Hamiltonians, of shape (n, 2, 2), one per slab and ordered along the trajectory, in units of eV.

widthsarray_like

Slab widths, of shape (n,), in units of eV^-1. Use globaldefs.CONV_KM_TO_INV_EV to convert from km.

Returns:

numpy.ndarray

The evolution operator, of shape (2, 2), indexed (final, initial).

Raises:

ValueError

If the number of Hamiltonians and widths differ, if either is empty, if the Hamiltonians are not \(2\times2\), or if any width is negative.

Examples

import slabs

import numpy as np
H = np.array([[[0.0, 1.0], [1.0, 0.0]],
              [[0.0, 0.5], [0.5, 0.0]]])
U = slabs.evolution_operator_2nu_slabs(H, [0.3, 0.4])
print('%.6f' % abs(U[0][0]))

0.877583

slabs.evolution_operator_3nu_slabs(hamiltonian_matrices: list | ndarray, widths: list | ndarray) → ndarray

Returns the three-flavor evolution operator across adjacent slabs.

Returns the \(3\times3\) evolution operator \(U_3 = U_3^{(n)}(L_n) \cdots U_3^{(1)}(L_1)\) for a trajectory divided into slabs, each with its own constant Hamiltonian and width. Each slab is solved exactly by the SU(3) expansion of oscprob3nu.

Added in version 1.8.0.

Parameters:

hamiltonian_matricesarray_like

Hamiltonians, of shape (n, 3, 3), one per slab and ordered along the trajectory, in units of eV.

widthsarray_like

Slab widths, of shape (n,), in units of eV^-1. Use globaldefs.CONV_KM_TO_INV_EV to convert from km.

Returns:

numpy.ndarray

The evolution operator, of shape (3, 3), indexed (final, initial).

Raises:

ValueError

If the number of Hamiltonians and widths differ, if either is empty, if the Hamiltonians are not \(3\times3\), or if any width is negative.

Examples

import slabs

import numpy as np
H = np.array([np.diag([1.0, 0.0, -1.0]),
              np.diag([0.5, 0.0, -0.5])], dtype=complex)
U = slabs.evolution_operator_3nu_slabs(H, [0.2, 0.3])
print('%.6f' % abs(U[0][0]))

1.000000

slabs.evolution_operator_4nu_slabs(hamiltonian_matrices: list | ndarray, widths: list | ndarray) → ndarray

Returns the four-flavor evolution operator across adjacent slabs.

Returns the \(4\times4\) evolution operator \(U_4 = U_4^{(n)}(L_n) \cdots U_4^{(1)}(L_1)\) for a trajectory divided into slabs, each with its own constant Hamiltonian and width. Each slab is solved exactly by the SU(4) expansion of oscprob4nu.

This is what makes a 3+1 scenario propagable through layered matter: the sterile state’s matter entry is constant within a slab like every other, so nothing about the composition changes at four flavors. See hamiltonians4nu.hamiltonian_4nu_matter() for the entry itself, which is \(-V_{NC}\) rather than zero.

Added in version 1.11.0.

Parameters:

hamiltonian_matricesarray_like

Hamiltonians, of shape (n, 4, 4), one per slab and ordered along the trajectory, in units of eV.

widthsarray_like

Slab widths, of shape (n,), in units of eV^-1. Use globaldefs.CONV_KM_TO_INV_EV to convert from km.

Returns:

numpy.ndarray

The evolution operator, of shape (4, 4), indexed (final, initial).

Raises:

ValueError

If the number of Hamiltonians and widths differ, if either is empty, if the Hamiltonians are not \(4\times4\), or if any width is negative.

Examples

import slabs

import numpy as np
H = np.array([np.diag([1.0, 0.0, -0.5, -0.5]),
              np.diag([0.5, 0.0, -0.25, -0.25])], dtype=complex)
U = slabs.evolution_operator_4nu_slabs(H, [0.2, 0.3])
print('%.6f' % abs(U[0][0]))

1.000000

slabs.probabilities_2nu_slabs(hamiltonian_matrices: list | ndarray, widths: list | ndarray) → Tuple[float, float, float, float]

Returns the two-flavor probabilities across adjacent slabs.

Returns \(P_{ee}, P_{e\mu}, P_{\mu e}, P_{\mu\mu}\), where \(P_{\alpha\beta} \equiv P(\nu_\alpha \to \nu_\beta) = |[U_2]_{\beta\alpha}|^2\), for a trajectory divided into slabs.

Added in version 1.8.0.

Parameters:

hamiltonian_matricesarray_like

Hamiltonians, of shape (n, 2, 2), one per slab and ordered along the trajectory, in units of eV.

widthsarray_like

Slab widths, of shape (n,), in units of eV^-1.

Returns:

tuple of float

The probabilities \(P_{ee}, P_{e\mu}, P_{\mu e}, P_{\mu\mu}\).

Raises:

ValueError

If the slab sequence is malformed; see evolution_operator_2nu_slabs.

Examples

import slabs

import numpy as np
H = np.array([[[0.0, 1.0], [1.0, 0.0]],
              [[0.0, 0.5], [0.5, 0.0]]])
Pee, Pem, Pme, Pmm = slabs.probabilities_2nu_slabs(H, [0.3, 0.4])
print('%.6f  %.6f' % (Pee, Pem))

0.770151  0.229849

slabs.probabilities_3nu_slabs(hamiltonian_matrices: list | ndarray, widths: list | ndarray) → Tuple[float, float, float, float, float, float, float, float, float]

Returns the three-flavor probabilities across adjacent slabs.

Returns \(P_{ee}, P_{e\mu}, P_{e\tau}, P_{\mu e}, P_{\mu\mu}, P_{\mu\tau}, P_{\tau e}, P_{\tau\mu}, P_{\tau\tau}\), where \(P_{\alpha\beta} \equiv P(\nu_\alpha \to \nu_\beta) = |[U_3]_{\beta\alpha}|^2\), for a trajectory divided into slabs.

Added in version 1.8.0.

Parameters:

hamiltonian_matricesarray_like

Hamiltonians, of shape (n, 3, 3), one per slab and ordered along the trajectory, in units of eV.

widthsarray_like

Slab widths, of shape (n,), in units of eV^-1.

Returns:

tuple of float

The nine probabilities, with the initial flavor varying slowest.

Raises:

ValueError

If the slab sequence is malformed; see evolution_operator_3nu_slabs.

Examples

import slabs

import numpy as np
H = np.array([np.diag([1.0, 0.0, -1.0]),
              np.diag([0.5, 0.0, -0.5])], dtype=complex)
prob = slabs.probabilities_3nu_slabs(H, [0.2, 0.3])
print('%.6f  %.6f' % (prob[0], prob[1]))

1.000000  0.000000

slabs.probabilities_4nu_slabs(hamiltonian_matrices: list | ndarray, widths: list | ndarray) → Tuple[float, ...]

Returns the four-flavor probabilities across adjacent slabs.

Returns the sixteen probabilities \(P_{\alpha\beta} \equiv P(\nu_\alpha \to \nu_\beta) = |[U_4]_{\beta\alpha}|^2\), ordered with the initial flavor varying slowest, for a trajectory divided into slabs. With the fourth state read as sterile, the flavor order is \((\nu_e, \nu_\mu, \nu_\tau, \nu_s)\).

Added in version 1.11.0.

Parameters:

hamiltonian_matricesarray_like

Hamiltonians, of shape (n, 4, 4), one per slab and ordered along the trajectory, in units of eV.

widthsarray_like

Slab widths, of shape (n,), in units of eV^-1.

Returns:

tuple of float

The sixteen probabilities, with the initial flavor varying slowest.

Raises:

ValueError

If the slab sequence is malformed; see evolution_operator_4nu_slabs.

Examples

import slabs

import numpy as np
H = np.array([np.diag([1.0, 0.0, -0.5, -0.5]),
              np.diag([0.5, 0.0, -0.25, -0.25])], dtype=complex)
prob = slabs.probabilities_4nu_slabs(H, [0.2, 0.3])
print('%.6f  %.6f' % (prob[0], prob[1]))

1.000000  0.000000

earth

The Earth as a sequence of slabs: PREM, chord geometry, and probabilities.

A neutrino crossing the Earth sees a matter density that changes continuously along its path, so the Hamiltonian is not constant and the exact expansions of oscprob2nu and oscprob3nu do not apply to the trajectory as a whole. They apply to any piece of it over which the density is taken to be constant, which is what this module builds: the chord is cut into slabs, each slab is solved exactly by slabs, and the operators are multiplied.

The density comes from the Preliminary Reference Earth Model (PREM) [DA81], a piecewise-polynomial fit to seismological data, given as \(\rho(x)\) with \(x = r/R_\oplus\).

Where the slabs are cut matters, and two different things are going on.

Between shells the density jumps, so the chord is first split at every crossing of a PREM shell boundary: no amount of subdivision recovers a discontinuity that straddles a slab. This gives a set of chord segments. Note a segment is not a shell — a chord enters and leaves each shell it reaches, so a diametric chord has 19 segments across 10 shells.

Within a shell the density varies smoothly, since PREM gives it as a polynomial in \(x = r/R_\oplus\) rather than a constant, and a segment can be long: crossing the mantle, the density changes by 21% over a single 2200 km segment. Each segment is therefore divided further into n_slabs_per_segment equal sub-slabs, with the density taken at the midpoint of each.

Midpoint sampling is second-order accurate, so the result converges to the continuous answer as the sub-slabs are refined; the routines take that number as an argument so a caller can watch it converge rather than trust it. Sampling at the midpoint rather than an end also matters for the segment that straddles the closest approach: it enters and exits at the same radius, so its two ends have identical density while the interior differs from them by 2.5%.

Units follow the rest of the library: energies in eV, baselines in eV^-1, potentials in eV. The exceptions are the geometry routines, which work in km because that is how the Earth is described, and density_prem, which returns g cm^-3 because that is how PREM is stated. matter_potential is the bridge between the two.

The Earth is treated as a sphere of radius globaldefs.EARTH_RADIUS.

Routine listings

• dms_to_decimal - Degrees, minutes, seconds to decimal degrees

• coordinates_of_named_location - Coordinates of a named site

• density_prem - PREM density at a radius

• matter_potential - Charged-current potential from a density

• matter_potential_nc - Neutral-current potential, for a sterile state

• distance_traveled_inside_earth - Chord length for a given costhz

• earth_radial_distance_from_depth - Radius at a point on the chord

• prem_layer_edges_along_chord - Where a chord crosses PREM shells

• chord_length_inside_earth - Chord between two surface locations

• costhz_between_points_on_surface - Its zenith angle

• earth_slabs - Slab widths and densities along a chord

• probabilities_2nu_earth - Two-flavor probabilities across the Earth

• probabilities_3nu_earth - Three-flavor probabilities across the Earth

• probabilities_4nu_earth - Four-flavor probabilities across the Earth

• probabilities_2nu_between_locations - Between two named sites

• probabilities_3nu_between_locations - Between two named sites

• probabilities_4nu_between_locations - Between two named sites

earth.LOC_COORDS_DMS = {'baikal': {'lat': (51, 45, 54), 'lon': (104, 24, 54)}, 'cern': {'lat': (46, 14, 1.8), 'lon': (6, 3, 11.4)}, 'desy': {'lat': (53, 34, 19.79), 'lon': (9, 52, 27.59)}, 'ess': {'lat': (55, 44, 6), 'lon': (13, 15, 5.04)}, 'fermilab': {'lat': (41, 49, 55), 'lon': (-88, 15, 26)}, 'gran_sasso': {'lat': (42, 25, 15.8), 'lon': (13, 30, 58.43)}, 'homestake': {'lat': (44, 21, 5.76), 'lon': (-103, 45, 4.68)}, 'kamioka': {'lat': (36, 25, 50.05), 'lon': (137, 18, 41.15)}, 'km3net_arca': {'lat': (36, 16, 0), 'lon': (16, 6, 0)}, 'km3net_orca': {'lat': (42, 48, 0), 'lon': (6, 2, 0)}, 'north_pole': {'lat': (90, 0, 0), 'lon': (0, 0, 0)}, 'pyhaasalmi': {'lat': (63, 39, 31), 'lon': (26, 2, 28)}, 'snolab': {'lat': (46, 28, 18), 'lon': (-81, 11, 12)}, 'south_pole': {'lat': (-90, 0, 0), 'lon': (0, 0, 0)}, 'tokai': {'lat': (36, 27, 59), 'lon': (140, 36, 24)}}

dict: Predefined locations, in ISO 6709 convention.

North latitudes and East longitudes are positive; South and West are negative. Each entry gives lat and lon as (degree, minute, second) tuples. The same set of sites as the sibling Magnus package, so a trajectory named in one can be reproduced in the other.

earth.PREM_BOUNDARIES = array([1221.5, 3480. , 5701. , 5771. , 5971. , 6151. , 6346.6, 6356. ,        6368. ])

numpy.ndarray: Outer radius of each PREM shell but the last.

The last shell ends at globaldefs.EARTH_RADIUS. Units: [km].

earth.dms_to_decimal(degrees: int | float, minutes: int | float, seconds: int | float) → float

Returns a (degree, minute, second) coordinate in decimal degrees.

Added in version 1.8.0.

Parameters:

degreesint or float

Degree part of the coordinate. Carries the sign: a location at 5 degrees South is (-5, ...). For a coordinate between zero and one degree South or West, where the degree part cannot carry a sign, negate the minutes instead: 0 deg 30’ S is (0, -30, 0).

minutesint or float

Minute part of the coordinate. Normally positive; a negative value flips the sign of the whole coordinate, which is the only way to express a southern or western coordinate smaller than one degree.

secondsint or float

Second part of the coordinate, taken as positive.

Returns:

float

The coordinate in decimal degrees.

Examples

import earth

print('%.6f' % earth.dms_to_decimal(36, 25, 50.05))

36.430569

earth.coordinates_of_named_location(loc_name: str) → Tuple[Tuple[float, float, float], Tuple[float, float, float]]

Returns the coordinates of a predefined location.

Looks loc_name up in LOC_COORDS_DMS, case-insensitively and treating spaces as underscores.

Added in version 1.8.0.

Parameters:

loc_namestr

Name of the location, e.g. 'kamioka' or 'south pole'.

Returns:

tuple of tuple of float

The latitude and longitude, each as (degree, minute, second).

Raises:

ValueError

If the name is not one of the predefined locations.

Examples

import earth

lat, lon = earth.coordinates_of_named_location('South Pole')
print(lat, lon)

(-90, 0, 0) (0, 0, 0)

earth.density_prem(r: int | float | list | ndarray, tol: float = 1e-08) → float | ndarray

Returns the matter density inside the Earth, according to PREM.

Evaluates the Preliminary Reference Earth Model [DA81] at a radial distance measured from the centre of the Earth. Accepts a single radius or an array of radii, evaluated in one vectorised pass.

Added in version 1.8.0.

Parameters:

rint, float, list or numpy.ndarray

Radial distance(s) from the centre of the Earth, in units of km.

tolfloat, optional

Relative tolerance by which a radius may exceed globaldefs.EARTH_RADIUS before a ValueError is raised. Radii within the tolerance are clamped onto the surface, which is what makes a chord endpoint computed in floating point safe to pass in. Default: 1e-8.

Returns:

float or numpy.ndarray

The matter density, in units of g cm^-3.

Raises:

ValueError

If any radius exceeds globaldefs.EARTH_RADIUS by more than the relative tolerance, or if any radius is negative.

Examples

import earth

print('%.4f' % earth.density_prem(0.0))
print('%.4f' % earth.density_prem(6371.0))

13.0885
1.0200

earth.matter_potential(density: int | float | list | ndarray, electron_fraction: float = 0.5) → float | ndarray

Returns the charged-current matter potential for a density.

Returns \(V_{CC} = \sqrt{2} G_F n_e\), the potential that hamiltonians3nu.hamiltonian_3nu_matter and its two-flavor counterpart expect. It is positive for neutrinos; pass its negative for antineutrinos.

Added in version 1.8.0.

Parameters:

densityint, float, list or numpy.ndarray

Matter density, in units of g cm^-3.

electron_fractionfloat, optional

Electrons per nucleon. Default: globaldefs.ELECTRON_FRACTION_EARTH_CRUST, which is 0.5 and is a good approximation everywhere in the Earth.

Returns:

float or numpy.ndarray

The potential \(V_{CC}\), in units of eV.

Examples

import earth

print('%.4e' % earth.matter_potential(3.0))

1.1358e-13

earth.matter_potential_nc(density: int | float | list | ndarray, neutron_fraction: float | None = None, electron_fraction: float = 0.5) → float | ndarray

Returns the neutral-current matter potential for a density.

Returns \(V_{NC} = -G_F n_n/\sqrt{2}\), which is negative for neutrinos. It is felt equally by all three active flavors, so at two and three flavors it is proportional to the identity and drops out of the probabilities entirely — which is why matter_potential alone serves them.

It does not drop out once a sterile state is present, because the sterile state does not feel it. Removing it from all four states costs only a global phase and leaves \(-V_{NC}\) on the sterile entry; see hamiltonians4nu.hamiltonian_4nu_matter().

Added in version 1.11.0.

Parameters:

densityint, float, list or numpy.ndarray

Matter density, in units of g cm^-3.

neutron_fractionfloat, optional

Neutrons per nucleon. Default: 1 - electron_fraction, the isoscalar value, since a nucleon is either a proton — matched by an electron — or a neutron.

electron_fractionfloat, optional

Electrons per nucleon, used only to derive neutron_fraction when that is not given. Default: globaldefs.ELECTRON_FRACTION_EARTH_CRUST.

Returns:

float or numpy.ndarray

The potential \(V_{NC}\), in units of eV. Negative for neutrinos.

Examples

import earth

print('%.4e' % earth.matter_potential_nc(3.0))

-5.6791e-14

earth.distance_traveled_inside_earth(costhz: int | float) → float

Returns the chord length through the Earth for a given direction.

The neutrino is assumed to arrive at a detector on the surface, not underground, so the distance is zero for any down-going direction, costhz >= 0, which reaches the detector from the sky without entering the Earth at all.

Added in version 1.8.0.

Parameters:

costhzint or float

Cosine of the zenith angle of the neutrino direction. costhz = -1 is straight up through the centre of the Earth.

Returns:

float

The chord length, in units of km.

Raises:

ValueError

If costhz lies outside \([-1, 1]\), where it describes no direction.

Examples

import earth

print('%.1f' % earth.distance_traveled_inside_earth(-1.0))
print('%.1f' % earth.distance_traveled_inside_earth(0.5))

12742.0
0.0

earth.earth_radial_distance_from_depth(costhz: int | float, l: int | float | list | ndarray, tol: float = 1e-08) → float | ndarray

Returns the radius at a point along a chord through the Earth.

A neutrino with direction costhz travels from l = 0 at its point of entry to l = distance_traveled_inside_earth (costhz) at the detector. This returns its distance from the centre of the Earth at a given l.

Added in version 1.8.0.

Parameters:

costhzint or float

Cosine of the zenith angle of the neutrino direction.

lint, float, list or numpy.ndarray

Distance(s) along the chord from the point of entry, in units of km.

tolfloat, optional

Absolute tolerance, in km, by which l may exceed the chord length before a ValueError is raised; values within the tolerance are clamped onto the exit point. Default: 1e-8.

Returns:

float or numpy.ndarray

The radial distance from the centre of the Earth, in units of km.

Raises:

ValueError

If any l is negative or exceeds the chord length for this costhz by more than the tolerance.

Examples

import earth

print('%.1f' % earth.earth_radial_distance_from_depth(-1.0, 6371.0))

0.0

earth.prem_layer_edges_along_chord(costhz: int | float) → ndarray

Returns where a chord through the Earth crosses PREM boundaries.

The density is discontinuous across a PREM shell boundary, so a slab that straddles one cannot represent it however finely the rest of the chord is divided. These positions are therefore mandatory slab edges, and earth_slabs uses them as such.

The crossings solve \(r(l) = r_b\) for each boundary radius \(r_b\), a quadratic in \(l\): with \(u = d - l\) and \(d = -2 R \cos\theta_z\),

\[u^2 + 2 R \cos\theta_z\, u + \left(R^2 - r_b^2\right) = 0 .\]

Added in version 1.8.0.

Parameters:

costhzint or float

Cosine of the zenith angle of the neutrino direction. Crossings exist only for costhz < 0.

Returns:

numpy.ndarray

Sorted crossing positions along the chord, in units of km, each strictly inside (0, d). Empty if the chord crosses no boundary.

Examples

import earth

print(len(earth.prem_layer_edges_along_chord(-1.0)))
print(len(earth.prem_layer_edges_along_chord(0.5)))

18
0

earth.chord_length_inside_earth(lat1_dms: Tuple[float, float, float], lon1_dms: Tuple[float, float, float], lat2_dms: Tuple[float, float, float], lon2_dms: Tuple[float, float, float]) → float

Returns the straight-line distance between two surface locations.

Computes the chord — the straight line through the Earth’s interior, not the great-circle arc over its surface — between two points on a spherical Earth, via the haversine formula for the central angle.

Added in version 1.8.0.

Parameters:

lat1_dmstuple of float

Latitude of the first location, as (degree, minute, second).

lon1_dmstuple of float

Longitude of the first location, as (degree, minute, second).

lat2_dmstuple of float

Latitude of the second location, as (degree, minute, second).

lon2_dmstuple of float

Longitude of the second location, as (degree, minute, second).

Returns:

float

The chord length, in units of km.

Examples

import earth

lat1, lon1 = earth.coordinates_of_named_location('fermilab')
lat2, lon2 = earth.coordinates_of_named_location('homestake')
print('%.1f' % earth.chord_length_inside_earth(lat1, lon1, lat2, lon2))

1284.7

earth.costhz_between_points_on_surface(lat1_dms: Tuple[float, float, float], lon1_dms: Tuple[float, float, float], lat2_dms: Tuple[float, float, float], lon2_dms: Tuple[float, float, float]) → float

Returns the zenith angle of the chord between two locations.

The cosine of the zenith angle at which a neutrino must travel to reach the second location from the first through the Earth’s interior. Both are assumed to be on the surface, so the result is never positive.

Added in version 1.8.0.

Parameters:

lat1_dmstuple of float

Latitude of the first location, as (degree, minute, second).

lon1_dmstuple of float

Longitude of the first location, as (degree, minute, second).

lat2_dmstuple of float

Latitude of the second location, as (degree, minute, second).

lon2_dmstuple of float

Longitude of the second location, as (degree, minute, second).

Returns:

float

Cosine of the zenith angle of the connecting chord.

Examples

import earth

lat1, lon1 = earth.coordinates_of_named_location('cern')
lat2, lon2 = earth.coordinates_of_named_location('gran_sasso')
print('%.6f' % earth.costhz_between_points_on_surface(lat1, lon1,
                                                lat2, lon2))

-0.057179

earth.earth_slabs(costhz: int | float, n_slabs_per_segment: int = 8) → Tuple[ndarray, ndarray]

Returns the slab widths and densities along a chord.

Cuts the chord at every PREM shell boundary it crosses, divides each resulting segment into n_slabs_per_segment equal sub-slabs, and evaluates the density at the midpoint of each. The boundary cuts are what make the discretisation converge quickly: they keep every slab inside a single shell, where the density is smooth.

Added in version 1.8.0.

Parameters:

costhzint or float

Cosine of the zenith angle of the neutrino direction. Must be negative, so that the neutrino crosses the Earth at all.

n_slabs_per_segmentint, optional

Number of equal sub-slabs per chord segment. A segment runs between consecutive PREM boundary crossings; a chord crosses most shells twice, so there are more segments than shells. Default: 8.

Returns:

tuple of numpy.ndarray

The slab widths, in units of km, and the density in each slab, in units of g cm^-3, ordered along the trajectory.

Raises:

ValueError

If costhz >= 0, so that there is no path through the Earth, or if n_slabs_per_segment is not positive.

Examples

import earth

widths, densities = earth.earth_slabs(-1.0, n_slabs_per_segment=2)
print(len(widths), '%.1f' % sum(widths))

38 12742.0

earth.probabilities_2nu_earth(h_vacuum_energy_independent: list | ndarray, energy: int | float, costhz: int | float, n_slabs_per_segment: int = 8, electron_fraction: float = 0.5) → Tuple[float, float, float, float]

Returns the two-flavor probabilities across the Earth.

Builds the PREM slabs along the chord for the given direction and propagates through them exactly, slab by slab.

Added in version 1.8.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent two-flavor vacuum Hamiltonian, as returned by hamiltonians2nu.hamiltonian_2nu_vacuum_energy_independent.

energyint or float

Neutrino energy, in units of eV.

costhzint or float

Cosine of the zenith angle of the neutrino direction. Must be negative.

n_slabs_per_segmentint, optional

Number of equal sub-slabs per chord segment. A segment runs between consecutive PREM boundary crossings; a chord crosses most shells twice, so there are more segments than shells. Default: 8.

electron_fractionfloat, optional

Electrons per nucleon. Default: globaldefs.ELECTRON_FRACTION_EARTH_CRUST.

Returns:

tuple of float

The probabilities \(P_{ee}, P_{e\mu}, P_{\mu e}, P_{\mu\mu}\).

Raises:

ValueError

If costhz >= 0 or n_slabs_per_segment is not positive.

earth.probabilities_3nu_earth(h_vacuum_energy_independent: list | ndarray, energy: int | float, costhz: int | float, n_slabs_per_segment: int = 8, electron_fraction: float = 0.5) → Tuple[float, float, float, float, float, float, float, float, float]

Returns the three-flavor probabilities across the Earth.

Builds the PREM slabs along the chord for the given direction and propagates through them exactly, slab by slab. Raising n_slabs_per_segment and watching the result settle is the way to confirm the discretisation is fine enough for the energy in question; the number needed grows as the oscillation length falls.

Added in version 1.8.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent three-flavor vacuum Hamiltonian, as returned by hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent.

energyint or float

Neutrino energy, in units of eV.

costhzint or float

Cosine of the zenith angle of the neutrino direction. Must be negative.

n_slabs_per_segmentint, optional

Number of equal sub-slabs per chord segment. A segment runs between consecutive PREM boundary crossings; a chord crosses most shells twice, so there are more segments than shells. Default: 8.

electron_fractionfloat, optional

Electrons per nucleon. Default: globaldefs.ELECTRON_FRACTION_EARTH_CRUST.

Returns:

tuple of float

The nine probabilities, with the initial flavor varying slowest.

Raises:

ValueError

If costhz >= 0 or n_slabs_per_segment is not positive.

earth.probabilities_4nu_earth(h_vacuum_energy_independent: list | ndarray, energy: int | float, costhz: int | float, n_slabs_per_segment: int = 8, electron_fraction: float = 0.5) → Tuple[float, ...]

Returns the four-flavor probabilities across the Earth.

Builds the PREM slabs along the chord for the given direction and propagates through them exactly, slab by slab, exactly as at two and three flavors. The one thing that is new is the potential: a sterile state does not feel the neutral-current interaction, so \(V_{NC}\) no longer cancels between the flavors and is built per slab alongside \(V_{CC}\). This is what puts the sterile matter resonance where it belongs; see hamiltonians4nu.hamiltonian_4nu_matter().

Added in version 1.11.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent four-flavor vacuum Hamiltonian, as returned by hamiltonians4nu.hamiltonian_4nu_vacuum_energy_independent.

energyint or float

Neutrino energy, in units of eV.

costhzint or float

Cosine of the zenith angle of the neutrino direction. Must be negative.

n_slabs_per_segmentint, optional

Number of equal sub-slabs per chord segment. A segment runs between consecutive PREM boundary crossings; a chord crosses most shells twice, so there are more segments than shells. Default: 8.

electron_fractionfloat, optional

Electrons per nucleon. The neutron fraction is taken as its complement, which is what sets \(V_{NC}\). Default: globaldefs.ELECTRON_FRACTION_EARTH_CRUST.

Returns:

tuple of float

The sixteen probabilities, with the initial flavor varying slowest. With the fourth state read as sterile, the flavor order is \((\nu_e, \nu_\mu, \nu_\tau, \nu_s)\).

Raises:

ValueError

If costhz >= 0 or n_slabs_per_segment is not positive.

earth.probabilities_2nu_between_locations(h_vacuum_energy_independent: list | ndarray, energy: int | float, loc_name_1: str, loc_name_2: str, n_slabs_per_segment: int = 8, electron_fraction: float = 0.5) → Tuple[float, float, float, float]

Returns the two-flavor probabilities between two named locations.

Convenience wrapper: looks both locations up in LOC_COORDS_DMS, finds the zenith angle of the chord joining them, and evaluates probabilities_2nu_earth along it.

Added in version 1.8.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent two-flavor vacuum Hamiltonian.

energyint or float

Neutrino energy, in units of eV.

loc_name_1str

Name of the source location, e.g. 'fermilab'.

loc_name_2str

Name of the detector location, e.g. 'homestake'.

n_slabs_per_segmentint, optional

Number of equal sub-slabs per chord segment. A segment runs between consecutive PREM boundary crossings; a chord crosses most shells twice, so there are more segments than shells. Default: 8.

electron_fractionfloat, optional

Electrons per nucleon. Default: globaldefs.ELECTRON_FRACTION_EARTH_CRUST.

Returns:

tuple of float

The probabilities \(P_{ee}, P_{e\mu}, P_{\mu e}, P_{\mu\mu}\).

Raises:

ValueError

If either name is not predefined, or if the two locations coincide, so that there is no chord between them.

earth.probabilities_3nu_between_locations(h_vacuum_energy_independent: list | ndarray, energy: int | float, loc_name_1: str, loc_name_2: str, n_slabs_per_segment: int = 8, electron_fraction: float = 0.5) → Tuple[float, float, float, float, float, float, float, float, float]

Returns the three-flavor probabilities between two named locations.

Convenience wrapper: looks both locations up in LOC_COORDS_DMS, finds the zenith angle of the chord joining them, and evaluates probabilities_3nu_earth along it.

Added in version 1.8.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent three-flavor vacuum Hamiltonian.

energyint or float

Neutrino energy, in units of eV.

loc_name_1str

Name of the source location, e.g. 'cern'.

loc_name_2str

Name of the detector location, e.g. 'gran_sasso'.

n_slabs_per_segmentint, optional

Number of equal sub-slabs per chord segment. A segment runs between consecutive PREM boundary crossings; a chord crosses most shells twice, so there are more segments than shells. Default: 8.

electron_fractionfloat, optional

Electrons per nucleon. Default: globaldefs.ELECTRON_FRACTION_EARTH_CRUST.

Returns:

tuple of float

The nine probabilities, with the initial flavor varying slowest.

Raises:

ValueError

If either name is not predefined, or if the two locations coincide, so that there is no chord between them.

earth.probabilities_4nu_between_locations(h_vacuum_energy_independent: list | ndarray, energy: int | float, loc_name_1: str, loc_name_2: str, n_slabs_per_segment: int = 8, electron_fraction: float = 0.5) → Tuple[float, ...]

Returns the four-flavor probabilities between two named locations.

Convenience wrapper: looks both locations up in LOC_COORDS_DMS, finds the zenith angle of the chord joining them, and evaluates probabilities_4nu_earth along it.

Added in version 1.11.0.

Parameters:

h_vacuum_energy_independentarray_like

Energy-independent four-flavor vacuum Hamiltonian.

energyint or float

Neutrino energy, in units of eV.

loc_name_1str

Name of the source location, e.g. 'cern'.

loc_name_2str

Name of the detector location, e.g. 'gran_sasso'.

n_slabs_per_segmentint, optional

Number of equal sub-slabs per chord segment. A segment runs between consecutive PREM boundary crossings; a chord crosses most shells twice, so there are more segments than shells. Default: 8.

electron_fractionfloat, optional

Electrons per nucleon. Default: globaldefs.ELECTRON_FRACTION_EARTH_CRUST.

Returns:

tuple of float

The sixteen probabilities, with the initial flavor varying slowest.

Raises:

ValueError

If either name is not predefined, or if the two locations coincide, so that there is no chord between them.

Optional compiled backend

Optional Numba-compiled kernels for the batched evaluation paths.

NuOscProbExact needs only NumPy. If Numba happens to be installed, this module compiles the two-, three- and four-neutrino expansions into fused machine-code loops and oscprob2nu, oscprob3nu and oscprob4nu use them for large stacks; if it is not, HAVE_NUMBA is False, nothing here is defined, and the NumPy path is used instead. Nothing else in the library changes either way, and the results agree to round-off — see tests/test_fastkernels.py, which runs both paths against each other whichever is available.

Install the optional dependency with:

pip install "nuoscprobexact[fast]"

Why it is worth compiling

The NumPy path evaluates the expansion as a sequence of whole-array operations, so a stack of N Hamiltonians makes roughly fifteen passes over N-element arrays, each writing a temporary that the next pass reads back. The compiled kernel does the same arithmetic one element at a time, keeping every intermediate in registers, and spreads the elements over the available cores. Measured against the NumPy path on this library’s own benchmarks, best of seven runs with the two paths interleaved:

Stack	Speedup
200 000 energies, four flavors	~19x
20 000 energies, four flavors	~18x
200 000 energies, three flavors	~15x
20 000 energies, three flavors	~9x
100 x 100 oscillogram	~3.5x
200 000 baselines, two flavors	~1.5x

Four flavors gains the most, and not because the kernel is cleverer there: the NumPy path has the furthest to fall. Its expansion needs a quartic, a Newton refinement of the four roots against the matrix, and a Newton-form reconstruction, which as whole-array operations is some forty passes over the stack; done one element at a time none of it leaves the registers.

These are one machine on one day, and they move by tens of per cent between runs; read them as the shape of the gain rather than as constants. The figures quoted for 1.6.0 in CHANGELOG.md came from a different session and differ by up to a factor of two — which is why notebook 09 measures the comparison when it runs, on whatever machine is running it, rather than repeating a number from here.

Costs, so that the trade is visible

• importing Numba takes about 140 ms, against 65 ms for NumPy alone;

• the first call compiles, which takes a few seconds. The kernels are declared with cache=True, so that cost is paid once per machine and later runs load the compiled code from disk in milliseconds.

Both are why this is an optional extra rather than a dependency, and why the scalar path is deliberately left alone: a single probability takes about 8 microseconds, which is not worth a compilation pause.

Turning it off

Set fastkernels.USE_NUMBA = False to force the NumPy path even when Numba is installed — useful for checking that the two agree, which is what the test suite does.

Routine listings

• available - Whether the compiled kernels can be used at all

• worthwhile - Whether a stack is large enough to be worth compiling

• probabilities_2nu_kernel - Two-flavor probabilities for a stack

• probabilities_3nu_kernel - Three-flavor probabilities for a stack

• probabilities_4nu_kernel - Four-flavor probabilities for a stack

fastkernels.USE_NUMBA = True

bool: Module-level switch.

Set to False to force the NumPy path even when Numba is installed. available reports the two together.

fastkernels.MIN_BATCH = {2: 50000, 3: 1, 4: 1}

dict: Module-level constant.

The smallest stack for which the compiled kernel is worth using, by number of flavors. A backend that is sometimes slower than the path it replaces is worse than no backend, so these are measured rather than assumed.

For three flavors the kernel wins at every size, by between two and sixteen times, so the threshold is one. Four flavors is the same story only more so, and for a reason worth stating: the NumPy path there has no short-stack shortcut to fall back on — oscprob4nu has no separate scalar closed form, so even a stack of one pays for the whole array machinery, a batched determinant and all. Measured by alternating the two paths through oscprob4nu.probabilities_4nu() and taking the best of nine rounds each, the kernel leads by 15x at a single element, falls to 5x just below PARALLEL_THRESHOLD where it is still single-threaded, and settles at 18x once the threads are in use. It is never behind, so the threshold is one.

For two flavors it does not: that expansion reduces to a square root and a sine per element, which NumPy already does about as well as compiled code can, and the kernel additionally has to materialise the Hamiltonian stack — which for a scan over baselines is the same matrix repeated, costing 2.5 ms to copy at two hundred thousand points. Measured by alternating the two paths and taking the best of nine rounds each, the crossover sits between twenty and fifty thousand elements: at twenty thousand NumPy is still ahead by a few per cent, at fifty thousand the kernel leads by 1.3x and it grows slowly from there. The threshold is set at the first size where the kernel is unambiguously ahead, since the region around the crossover is broad and varies between machines.

fastkernels.PARALLEL_THRESHOLD = 256

int: Module-level constant.

Stacks with at least this many elements are spread over the available cores; smaller ones run in a single thread, because below roughly this size the cost of waking the thread pool exceeds what it saves.

fastkernels.available() → bool

Returns whether the compiled kernels can be used at all.

True when Numba was imported successfully and USE_NUMBA has not been turned off. Whether they are worth using for a given stack is a separate question; see worthwhile.

Added in version 1.6.0.

Returns:

bool

Whether probabilities_2nu_kernel, probabilities_3nu_kernel and probabilities_4nu_kernel may be called.

fastkernels.worthwhile(n_flavors: int, size: int) → bool

Returns whether the compiled kernel should be used for a stack.

The kernels are only used where they have been measured to win. Below the per-flavor threshold in MIN_BATCH the NumPy path is quicker, and using the kernel anyway would make installing the optional extra a pessimisation for those calls.

Added in version 1.6.0.

Parameters:

n_flavorsint

Number of neutrino flavors, 2, 3, or 4.

sizeint

Number of elements in the stack.

Returns:

bool

Whether to call the corresponding kernel.

Constants

Physical constants, unit-conversion factors, and the NuFit 4.0 best-fit oscillation parameters [EGGHC+19] for both mass orderings. The core modules do not need any of these. The sample non-standard-interaction strengths are deliberately large, so that the worked examples show a visible effect; the combination matter oscillations are sensitive to puts them in the LMA-D region [CDGG+17] rather than anywhere near a fit.

globaldefs

Physical constants, unit-conversion factors, and default parameters.

This module contains the values of the physical constants, the unit-conversion factors, and the default oscillation parameters used by the sample-Hamiltonian modules of NuOscProbExact, by earth, and by the notebooks and worked examples.

The core modules oscprob2nu, oscprob3nu and oscprob4nu do not need these constants: they accept an arbitrary Hermitian Hamiltonian in whatever units the user prefers.

Unless stated otherwise, quantities are expressed in natural units, in which energies are in eV, mass-squared differences in eV², and baselines in eV^-1.

Notes

The lepton mixing parameters are the best-fit values from NuFit 4.0 [1], including Super-Kamiokande atmospheric data, for both the normal ordering (suffix _NO_BF) and the inverted ordering (suffix _IO_BF). The non-standard-interaction strengths are taken from [2]. These constants support the method of [3].

References

[1]

I. Esteban et al., “Global analysis of three-flavour neutrino oscillations”, JHEP 01, 106 (2019), arXiv:1811.05487 (NuFit 4.0).

[2]

P. Coloma et al., “Curtailing the dark side in non-standard neutrino interactions”, arXiv:1805.04530.

[3]

Mauricio Bustamante, “Exact neutrino oscillation probabilities with arbitrary time-independent Hamiltonians”, arXiv:1904.12391.

globaldefs.CONV_KM_TO_INV_EV = 5067730000.0

float: Multiplicative conversion factor from km to eV^-1.

Units: [eV^-1 km^-1].

globaldefs.CONV_CM_TO_INV_EV = 50677.3

float: Multiplicative conversion factor from cm to eV^-1.

Units: [eV^-1 cm^-1].

globaldefs.CONV_INV_EV_TO_CM = 1.9732700834495916e-05

float: Multiplicative conversion factor from eV^-1 to cm.

Units: [eV cm].

globaldefs.CONV_EV_TO_G = 1.78266192e-33

float: Multiplicative conversion factor from eV to grams.

Converts a mass expressed in eV (i.e., eV/c²) into grams. Units: [g eV^-1].

Changed in version 1.11.0: Given to the precision of the CODATA value, 1.78266192e-33, rather than rounded to 1.783e-33. The rounded value was off by \(1.9 \times 10^{-4}\) relative, three orders of magnitude worse than every other constant in this module, which sit between \(10^{-7}\) and \(10^{-9}\). It propagates through NUM_DENSITY_E_EARTH_CRUST into VCC_EARTH_CRUST, and through earth.matter_potential() into every Earth crossing, so the matter potential moved by that much. Far below anything measurable, and still worth not carrying into a release.

globaldefs.CONV_G_TO_EV = 5.609588608927037e+32

float: Multiplicative conversion factor from grams to eV.

Units: [eV g^-1].

globaldefs.GF = 1.1663787e-23

float: Fermi constant.

Units: [eV^-2].

globaldefs.MASS_ELECTRON = 510998.9461

float: Electron mass.

Units: [eV].

globaldefs.MASS_PROTON = 938272046.0

float: Proton mass.

Units: [eV].

globaldefs.MASS_NEUTRON = 939565379.0

float: Neutron mass.

Units: [eV].

globaldefs.EARTH_RADIUS = 6371.0

float: Mean radius of the Earth.

The IUGG mean radius, which is what the Preliminary Reference Earth Model in earth is normalised against and what the chord geometry there assumes. The Earth is treated as a sphere throughout; the equatorial and polar radii differ from this by about 0.3%, which is far below the accuracy of any matter-density model.

Units: [km].

globaldefs.ELECTRON_FRACTION_EARTH_CRUST = 0.5

float: Electron fraction in the Earth’s crust.

Units: [adimensional].

globaldefs.DENSITY_MATTER_CRUST_G_PER_CM3 = 3.0

float: Average matter density in the Earth’s crust.

Units: [g cm^-3].

globaldefs.NUM_DENSITY_E_EARTH_CRUST = 6885791558.458149

float: Electron number density in the Earth’s crust.

Units: [eV³].

globaldefs.VCC_EARTH_CRUST = np.float64(1.1358172231000782e-13)

float: Charged-current matter potential in the Earth’s crust.

Equal to \(\sqrt{2} G_F n_e\). It is positive for neutrinos; use its negative for antineutrinos. Units: [eV].

globaldefs.NEUTRON_FRACTION_EARTH_CRUST = 0.5

float: Neutrons per nucleon in the Earth’s crust.

The crust is close to isoscalar, so with an electron fraction of one half there is about one neutron per electron. Units: [adimensional].

globaldefs.NUM_DENSITY_N_EARTH_CRUST = 6885791558.458149

float: Neutron number density in the Earth’s crust.

Units: [eV³].

globaldefs.VNC_EARTH_CRUST = np.float64(-5.6790861155003904e-14)

float: Neutral-current matter potential in the Earth’s crust.

Equal to \(-G_F n_n/\sqrt{2}\), and negative for neutrinos. It is felt equally by all three active flavors, so at three flavors it is proportional to the identity and drops out of the probabilities entirely — which is why hamiltonians3nu never needs it.

It matters as soon as a sterile state is present, because the sterile state does not feel it: subtracting it from all four states, which costs only a global phase, leaves \(-V_{NC}\) on the sterile entry. See hamiltonians4nu.hamiltonian_4nu_matter().

Added in version 1.9.0.

Units: [eV].

globaldefs.S12_NO_BF = np.float64(0.5567764362830022)

float: Lepton mixing angle \(\sin\theta_{12}\), normal ordering.

Units: [adimensional].

globaldefs.S23_NO_BF = np.float64(0.762889244910426)

float: Lepton mixing angle \(\sin\theta_{23}\), normal ordering.

Units: [adimensional].

globaldefs.S13_NO_BF = np.float64(0.14966629547095767)

float: Lepton mixing angle \(\sin\theta_{13}\), normal ordering.

Units: [adimensional].

globaldefs.DCP_NO_BF = 3.787364476827695

float: CP-violation phase \(\delta_{\rm CP}\), normal ordering.

Units: [radian].

globaldefs.D21_NO_BF = 7.39e-05

float: Mass-squared difference \(\Delta m^2_{21}\), normal ordering.

Units: [eV²].

globaldefs.D31_NO_BF = 0.002525

float: Mass-squared difference \(\Delta m^2_{31}\), normal ordering.

Units: [eV²].

globaldefs.S12_IO_BF = np.float64(0.5567764362830022)

float: Lepton mixing angle \(\sin\theta_{12}\), inverted ordering.

Units: [adimensional].

globaldefs.S23_IO_BF = np.float64(0.762889244910426)

float: Lepton mixing angle \(\sin\theta_{23}\), inverted ordering.

Units: [adimensional].

globaldefs.S13_IO_BF = np.float64(0.15043270920913443)

float: Lepton mixing angle \(\sin\theta_{13}\), inverted ordering.

Units: [adimensional].

globaldefs.DCP_IO_BF = 4.886921905584122

float: CP-violation phase \(\delta_{\rm CP}\), inverted ordering.

Units: [radian].

globaldefs.D21_IO_BF = 7.39e-05

float: Mass-squared difference \(\Delta m^2_{21}\), inverted ordering.

Units: [eV²].

globaldefs.D32_IO_BF = -0.002512

float: Mass-squared difference \(\Delta m^2_{32}\), inverted ordering.

Units: [eV²].

globaldefs.D31_IO_BF = -0.0024381

float: Mass-squared difference \(\Delta m^2_{31}\), inverted ordering.

Computed as \(\Delta m^2_{32} + \Delta m^2_{21}\). Units: [eV²].

globaldefs.EPS_EE = 0.06

float: NSI strength parameter \(\epsilon_{ee}\).

Units: [adimensional].

globaldefs.EPS_EM = -0.06

float: NSI strength parameter \(\epsilon_{e\mu}\).

May in general be complex. Units: [adimensional].

globaldefs.EPS_ET = 0.0

float: NSI strength parameter \(\epsilon_{e\tau}\).

May in general be complex. Units: [adimensional].

globaldefs.EPS_MM = 1.2

float: NSI strength parameter \(\epsilon_{\mu\mu}\).

Units: [adimensional].

globaldefs.EPS_MT = 0.0

float: NSI strength parameter \(\epsilon_{\mu\tau}\).

May in general be complex. Units: [adimensional].

globaldefs.EPS_TT = 0.0

float: NSI strength parameter \(\epsilon_{\tau\tau}\).

Units: [adimensional].

globaldefs.EPS_2 = [0.06, -0.06, 1.2]

list of float: NSI strengths for two-neutrino oscillations.

Ordered as [eps_ee, eps_em, eps_mm], ready to be passed to hamiltonians2nu.hamiltonian_2nu_nsi(). Units: [adimensional].

globaldefs.EPS_3 = [0.06, -0.06, 0.0, 1.2, 0.0, 0.0]

list of float: NSI strengths for three-neutrino oscillations.

Ordered as [eps_ee, eps_em, eps_et, eps_mm, eps_mt, eps_tt], ready to be passed to hamiltonians3nu.hamiltonian_3nu_nsi(). Units: [adimensional].

globaldefs.SXI12 = 0.0

float: LIV mixing angle \(\sin\xi_{12}\).

Units: [adimensional].

globaldefs.SXI23 = 0.0

float: LIV mixing angle \(\sin\xi_{23}\).

Units: [adimensional].

globaldefs.SXI13 = 0.0

float: LIV mixing angle \(\sin\xi_{13}\).

Units: [adimensional].

globaldefs.DXICP = 0.0

float: LIV CP-violation phase \(\delta_{\xi,\rm CP}\).

Units: [radian].

globaldefs.B1 = 1e-09

float: Eigenvalue \(b_1\) of the LIV operator.

Units: [eV].

globaldefs.B2 = 1e-09

float: Eigenvalue \(b_2\) of the LIV operator.

Units: [eV].

globaldefs.B3 = 2e-09

float: Eigenvalue \(b_3\) of the LIV operator.

Units: [eV].

globaldefs.LAMBDA = 1000000000000.0

float: Energy scale \(\Lambda\) of the LIV operator.

Units: [eV].