# -*- coding: utf-8 -*-
r"""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
:math:`2\times2` Hermitian Hamiltonian, using the SU(2) exponential
expansion described in [1]_.
The Hamiltonian is expanded in the basis of Pauli matrices,
.. math:: H = h_0 \mathbb{1} + h_k \sigma^k ,
and the time-evolution operator in the same basis,
.. math:: U_2(L) = u_0 \mathbb{1} + i u_k \sigma^k .
The term :math:`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
:math:`H L` is dimensionless. Elsewhere in **NuOscProbExact** the
Hamiltonian is in eV and the baseline in eV\ :sup:`-1`; the module
:mod:`globaldefs` provides ``CONV_KM_TO_INV_EV`` to convert a baseline
in km into eV\ :sup:`-1`.
Routine listings
----------------
* hamiltonian_2nu_coefficients - Returns the :math:`h_k`
* modulus - Returns the modulus :math:`|h|` of a vector
* evolution_operator_2nu_u_coefficients - Returns the :math:`u_k`
* evolution_operator_2nu - Returns the evolution operator :math:`U_2`
* probabilities_2nu - Returns the oscillation probabilities
References
----------
.. [1] Mauricio Bustamante, "Exact neutrino oscillation probabilities
with arbitrary time-independent Hamiltonians", arXiv:1904.12391.
"""
__author__ = "Mauricio Bustamante"
__email__ = "mbustamante@gmail.com"
__all__ = ['CHECK_HERMITICITY', 'SMALL_BATCH',
'hamiltonian_2nu_coefficients', 'modulus',
'evolution_operator_2nu_u_coefficients', 'evolution_operator_2nu',
'probabilities_2nu']
from typing import List, Tuple, Union
import math
import numpy as np
try:
import fastkernels
except ImportError: # pragma: no cover
# Copying this file on its own into another project is a supported way
# to use **NuOscProbExact**, and is documented as such --- but it stopped
# working in 1.6.0, when the optional compiled backend was added and
# imported unconditionally. A lone copy raised ImportError on the first
# line that mattered, which is a poor return for a promise the
# documentation makes twice.
#
# The backend is optional by design, so its absence is answered the same
# way its being switched off is: `worthwhile` says no, and the NumPy path
# runs. Nothing else in this module touches it.
class _NoFastKernels:
r"""Stands in for :mod:`fastkernels` when it is not importable."""
HAVE_NUMBA = False
USE_NUMBA = False
@staticmethod
def available():
r"""Returns False: there is no compiled backend here."""
return False
@staticmethod
def worthwhile(n_flavors, size):
r"""Returns False: no stack is worth a backend that is absent."""
return False
fastkernels = _NoFastKernels()
SMALL_BATCH = 11
r"""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 :data:`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.
"""
CHECK_HERMITICITY = True
r"""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
:mod:`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 :math:`2 \times 10^{-17}` relative, against a
tolerance of :math:`10^{-12}`.
.. versionadded:: 1.11.0
"""
_HERMITICITY_TOL = 1.e-12
r"""float: Module-level constant.
Relative tolerance for `CHECK_HERMITICITY`, measured against the largest
entry of the Hamiltonian.
"""
def _check_hermitian(h_matrix: np.ndarray, caller: str) -> None:
r"""Raises unless `h_matrix` is Hermitian to `_HERMITICITY_TOL`.
Compares the 1 independent pairs and the imaginary parts of the
diagonal, rather than forming ``H - H^dagger``, which would allocate
a temporary the size of the whole stack. Real and imaginary parts
are compared separately, on views rather than copies: the condition
is :math:`\mathrm{Re}\,H_{ij} = \mathrm{Re}\,H_{ji}` together
with :math:`\mathrm{Im}\,H_{ij} = -\mathrm{Im}\,H_{ji}`, which
needs no complex arithmetic and, unlike :func:`numpy.abs` on a
complex array, no square root per element. That is worth about a
factor of three on a large stack, where this check would otherwise
cost several times the evaluation it guards.
That describes the stack. A single matrix takes its own branch and
compares its entries as Python complex numbers, because the
reductions above are all fixed cost at one element and would
otherwise dominate a scalar probability.
Parameters
----------
h_matrix : numpy.ndarray
Hamiltonian, or stack of them, of shape ``(..., 2, 2)``.
caller : str
Name of the calling routine, used in the error message.
Returns
-------
None
Nothing; the routine either returns or raises.
Raises
------
ValueError
If any element of the stack is not Hermitian.
"""
if h_matrix.size == 0:
return
non_finite = (
'%s: the Hamiltonian has a non-finite entry, so it is neither '
'Hermitian nor usable. Set %s.CHECK_HERMITICITY = False to '
'skip this check.')
complaint = (
'%s: the Hamiltonian is not Hermitian%s. The expansion assumes '
'Hermiticity, and without it the probabilities returned are '
'meaningless even though they still sum to one. Set '
'oscprob2nu.CHECK_HERMITICITY = False to skip this check.')
# A single matrix takes its own path, because everything below is
# fixed cost at this size. The reductions run about sixty
# microseconds on one matrix against half a microsecond per element
# on a stack of two thousand, and a scalar probability costs eight
# microseconds in total --- so when this check arrived in 1.11.0 it
# became nine tenths of the work, and made short stacks slower than
# the batched path they exist to avoid. Comparing the entries as
# Python complex numbers is one conversion and a few scalar
# operations, and reaches the same verdict.
if h_matrix.ndim == 2:
entries = h_matrix.tolist()
# Tested per entry rather than by letting a non-finite value
# reach `scale`: `max` keeps its running value when the
# comparison is with a nan, since every comparison against one is
# false, so a nan would never arrive and `isfinite` would pass.
# The array path has no such hole, `np.max` propagating nan --- so
# the first draft of this branch returned probabilities for a
# Hamiltonian the batched path refuses, which is the divergence
# this check exists to prevent.
scale = 0.0
for row in entries:
for entry in row:
real, imaginary = abs(entry.real), abs(entry.imag)
if not (math.isfinite(real) and math.isfinite(imaginary)):
raise ValueError(non_finite % (caller, 'oscprob2nu'))
scale = max(scale, real, imaginary)
tolerance = (_HERMITICITY_TOL*scale if scale > 0.0
else _HERMITICITY_TOL)
for i in range(2):
if abs(entries[i][i].imag) > tolerance:
raise ValueError(complaint % (
caller, ' --- the diagonal entry (%d, %d) has a non-zero '
'imaginary part' % (i, i)))
for j in range(i+1, 2):
upper, lower = entries[i][j], entries[j][i]
if (abs(upper.real - lower.real) > tolerance
or abs(upper.imag + lower.imag) > tolerance):
raise ValueError(complaint % (
caller, ' --- entry (%d, %d) is not the complex '
'conjugate of entry (%d, %d)' % (i, j, j, i)))
return
real, imaginary = h_matrix.real, h_matrix.imag
# `np.abs(...).max()` allocates a float array the size of the stack,
# and is still the quickest way to the largest entry: replacing it with
# four reductions over `real` and `imaginary`, which allocate nothing,
# was measured 1.4x *slower* on two hundred thousand elements, because
# those views are strided over the complex array while `np.abs` reads
# it contiguously. Tried, measured, reverted.
scale = max(float(np.max(np.abs(real))), float(np.max(np.abs(imaginary))))
# A non-finite entry has to be caught here rather than left to
# propagate. It would otherwise make `scale` infinite or nan, hence
# `tolerance` infinite or nan, and every comparison below false ---
# so a Hamiltonian that is both non-finite *and* non-Hermitian would
# pass a check whose whole purpose is to refuse the second.
if not np.isfinite(scale):
raise ValueError(non_finite % (caller, 'oscprob2nu'))
tolerance = _HERMITICITY_TOL*scale if scale > 0.0 else _HERMITICITY_TOL
for i in range(2):
if np.any(np.abs(imaginary[..., i, i]) > tolerance):
raise ValueError(complaint % (
caller, ' --- the diagonal entry (%d, %d) has a non-zero '
'imaginary part' % (i, i)))
for j in range(i+1, 2):
if (np.any(np.abs(real[..., i, j] - real[..., j, i]) > tolerance)
or np.any(np.abs(imaginary[..., i, j]
+ imaginary[..., j, i]) > tolerance)):
raise ValueError(complaint % (
caller, ' --- entry (%d, %d) is not the complex '
'conjugate of entry (%d, %d)' % (i, j, j, i)))
[docs]
def hamiltonian_2nu_coefficients(
hamiltonian_matrix: Union[list, np.ndarray]
) -> List[float]:
r"""Returns the :math:`h_k` of the SU(2) expansion of the Hamiltonian.
Computes the coefficients :math:`h_1, h_2, h_3` of the SU(2)
expansion :math:`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 :math:`h_0` contributes
only an overall phase to the evolution operator and is not returned.
.. versionadded:: 1.0.0
.. versionchanged:: 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.
.. versionchanged:: 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:
:func:`numpy.real`, :func:`numpy.imag`, :obj:`numpy.arccos`,
:func:`numpy.clip` and :obj:`numpy.sqrt` on one number give way
to attribute access and the :mod:`math` module.
Parameters
----------
hamiltonian_matrix : array_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 :math:`|h|` of the returned vector.
Examples
--------
.. jupyter-execute::
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))
"""
H11 = hamiltonian_matrix[0][0]
H12 = hamiltonian_matrix[0][1]
H22 = hamiltonian_matrix[1][1]
# h0 = (H11+H22)/2.0 is not returned: it multiplies the identity and
# so contributes only an overall phase to U2, which cancels in the
# oscillation probabilities.
h1 = H12.real
h2 = -H12.imag
h3 = (H11-H22).real/2.0
return [float(h1), float(h2), float(h3)]
[docs]
def modulus(h_coeffs: Union[list, np.ndarray]) -> float:
r"""Returns the modulus :math:`|h|` of the vector of coefficients.
Returns the modulus of the vector of coefficients :math:`h_k` of the
SU(2) expansion of the two-neutrino Hamiltonian,
:math:`|h| = \sqrt{|h_1|^2 + |h_2|^2 + |h_3|^2}`.
.. versionadded:: 1.0.0
.. versionchanged:: 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 :mod:`math` rather than NumPy.
Parameters
----------
h_coeffs : array_like
Three-component vector of coefficients :math:`h_k`, as returned
by `hamiltonian_2nu_coefficients`.
Returns
-------
float
The modulus :math:`|h|`.
Examples
--------
.. jupyter-execute::
import oscprob2nu
print('%.6f' % oscprob2nu.modulus([0.0, -2.0, -1.0]))
"""
return math.sqrt(sum([abs(h)**2.0 for h in h_coeffs]))
def _hamiltonian_2nu_coefficients_batch(h_matrix: np.ndarray) -> np.ndarray:
r"""Returns the :math:`h_k` for a stack of Hamiltonians.
The vectorised counterpart of `hamiltonian_2nu_coefficients`.
Parameters
----------
h_matrix : numpy.ndarray
Hamiltonians, of shape ``(..., 2, 2)``.
Returns
-------
numpy.ndarray
The coefficients, of shape ``(3, ...)``.
Notes
-----
The component index is the *first* axis, so that each :math:`h_k` is
a contiguous array and the batch axes are the trailing ones, which
is what lets NumPy right-align them against the baselines. This
mirrors :func:`oscprob3nu._hamiltonian_3nu_coefficients_batch`.
"""
return np.stack([
h_matrix[..., 0, 1].real,
-h_matrix[..., 0, 1].imag,
(h_matrix[..., 0, 0]-h_matrix[..., 1, 1]).real/2.0,
], axis=0)
def _u_coefficients_2nu_batch(h: np.ndarray, L: np.ndarray) -> np.ndarray:
r"""Returns the four :math:`u_k` for a stack of Hamiltonians.
Parameters
----------
h : numpy.ndarray
The coefficients :math:`h_k`, of shape ``(3, ...)``, real.
L : numpy.ndarray
Baselines, of shape ``(...)``, broadcastable against `h`.
Returns
-------
numpy.ndarray
The coefficients, of shape ``(4, ...)``, real.
"""
# Pad the Hamiltonian's batch axes so that they right-align against
# the baselines; reshape returns a view, so this costs nothing
full = np.broadcast_shapes(h.shape[1:], np.shape(L))
extra = len(full) - (h.ndim - 1)
if extra > 0:
h = h.reshape(h.shape[:1] + (1,)*extra + h.shape[1:])
# |h| depends on the Hamiltonian alone; the baselines enter only in
# the trigonometric factors below
h_abs = np.sqrt((h*h).sum(0))
positive = h_abs > 0.0
safe_h_abs = np.where(positive, h_abs, 1.0)
phase = h_abs*L
u0 = np.broadcast_to(np.cos(phase), full)
# The limit of -sin(|h|L)/|h| as |h| -> 0 is -L
ss = np.where(positive, -np.sin(phase)/safe_h_abs,
-np.broadcast_to(L, full))
return np.concatenate([u0[None], np.broadcast_to(h*ss, (3,)+full)],
axis=0)
def _evolution_operator_2nu_batch(
h_matrix: Union[list, np.ndarray],
L: Union[int, float, list, np.ndarray]
) -> np.ndarray:
r"""Returns :math:`U_2(L)` for a stack of Hamiltonians and baselines.
Parameters
----------
h_matrix : array_like
Hamiltonians, of shape ``(..., 2, 2)``.
L : array_like
Baselines, broadcastable against the leading axes of `h_matrix`.
Returns
-------
numpy.ndarray
The evolution operators, of shape ``(..., 2, 2)``, complex.
"""
h_matrix = np.asarray(h_matrix, dtype=complex)
L = np.asarray(L, dtype=float)
# Check that the two broadcast against each other, and fail here with
# a clear message rather than deep inside the expansion
np.broadcast_shapes(h_matrix.shape[:-2], L.shape)
if CHECK_HERMITICITY:
_check_hermitian(h_matrix, 'evolution_operator_2nu')
u0, u1, u2, u3 = _u_coefficients_2nu_batch(
_hamiltonian_2nu_coefficients_batch(h_matrix), L)
return np.stack([
np.stack([u0+1.j*u3, 1.j*u1+u2], axis=-1),
np.stack([1.j*u1-u2, u0-1.j*u3], axis=-1),
], axis=-2)
def _probabilities_2nu_batch(
h_matrix: Union[list, np.ndarray],
L: Union[int, float, list, np.ndarray]
) -> np.ndarray:
r"""Returns the four probabilities for a stack, without forming U.
For a Hermitian :math:`2\times2` Hamiltonian the coefficients
:math:`u_k` are real, so
.. math::
|U_{ee}|^2 = u_0^2 + u_3^2 , \qquad
|U_{\mu e}|^2 = |U_{e\mu}|^2 = u_1^2 + u_2^2 ,
which leaves only two distinct numbers, and unitarity makes the
second the complement of the first. So neither the evolution
operator nor the coefficients need to be built: the transition
probability follows from the Hamiltonian directly, exactly as on
the scalar path.
Parameters
----------
h_matrix : array_like
Hamiltonians, of shape ``(..., 2, 2)``.
L : array_like
Baselines, broadcastable against the leading axes of `h_matrix`.
Returns
-------
numpy.ndarray
The probabilities, of shape ``(..., 4)``, ordered
``(Pee, Pem, Pme, Pmm)``.
"""
h_matrix = np.asarray(h_matrix, dtype=complex)
L = np.asarray(L, dtype=float)
batch = np.broadcast_shapes(h_matrix.shape[:-2], L.shape)
size = int(np.prod(batch, dtype=np.int64))
if CHECK_HERMITICITY:
_check_hermitian(h_matrix, 'probabilities_2nu')
if size > 0 and fastkernels.worthwhile(2, size):
return fastkernels.probabilities_2nu_kernel(
np.broadcast_to(h_matrix, batch+(2, 2)),
np.broadcast_to(L, batch))
if size <= SMALL_BATCH:
flat_h = np.broadcast_to(h_matrix, batch+(2, 2)).reshape(-1, 2, 2)
flat_l = np.broadcast_to(L, batch).reshape(-1)
out = np.empty((flat_l.shape[0], 4))
for n in range(flat_l.shape[0]):
# .tolist() gives Python complex numbers, on which the scalar
# path is quicker than on NumPy scalars, by more than the
# conversion costs
out[n] = probabilities_2nu(flat_h[n].tolist(), float(flat_l[n]))
return out.reshape(batch+(4,))
h = _hamiltonian_2nu_coefficients_batch(h_matrix)
h1, h2, h3 = h
h_sq = h1*h1 + h2*h2 + h3*h3
positive = h_sq > 0.0
safe_h_sq = np.where(positive, h_sq, 1.0)
sin_phase = np.sin(np.sqrt(safe_h_sq)*L)
# A Hamiltonian proportional to the identity drives no transitions
p_em = np.where(positive, (h1*h1 + h2*h2)/safe_h_sq*sin_phase*sin_phase,
0.0)
p_ee = 1.0 - p_em
return np.stack([p_ee, p_em, p_em, p_ee], axis=-1)
def _is_batched(
hamiltonian_matrix: Union[list, np.ndarray],
L: Union[int, float, list, np.ndarray]
) -> bool:
r"""Returns whether the arguments describe a stack of problems.
A single Hamiltonian is an ``n``-by-``n`` matrix and a single
baseline is a scalar, so anything with more axes than that is a
stack, and the vectorised path applies.
This runs on every scalar call, so it is written to be cheap: an
exact type check short-circuits the common case, and
``numpy.ndim`` --- which would convert a nested list to an array
every time --- is reached only for an argument that is neither a
plain Python number nor a NumPy array.
Parameters
----------
hamiltonian_matrix : array_like
Hamiltonian, or stack of them.
L : int or float or array_like
Baseline, or array of baselines.
Returns
-------
bool
Whether the vectorised path applies.
"""
if type(L) is not float and type(L) is not int:
if np.ndim(L) > 0:
return True
if type(hamiltonian_matrix) is np.ndarray:
return hamiltonian_matrix.ndim > 2
return isinstance(hamiltonian_matrix[0][0], (list, tuple, np.ndarray))
[docs]
def evolution_operator_2nu_u_coefficients(
hamiltonian_matrix: Union[list, np.ndarray],
L: Union[int, float]
) -> List[float]:
r"""Returns the coefficients :math:`u_0, \ldots, u_3`.
Returns the four coefficients :math:`u_0, \ldots, u_3` of the
two-neutrino time-evolution operator :math:`U_2(L)` in its SU(2)
exponential expansion,
:math:`U_2 = u_0 \mathbb{1} + i u_k \sigma^k`.
.. versionadded:: 1.0.0
.. versionchanged:: 1.1.0
Degenerate Hamiltonians are handled exactly instead of returning
NaN, by taking the limit :math:`\sin(|h|L)/|h| \to L`.
.. versionchanged:: 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:
:func:`numpy.real`, :func:`numpy.imag`, :obj:`numpy.arccos`,
:func:`numpy.clip` and :obj:`numpy.sqrt` on one number give way
to attribute access and the :mod:`math` module.
Parameters
----------
hamiltonian_matrix : array_like
Two-flavor Hermitian Hamiltonian, given as the nested list
``[[H11, H12], [H21, H22]]``.
L : float
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 :math:`i` that
multiplies :math:`u_k` is part of the expansion, not of the
coefficients.
Notes
-----
When :math:`|h| = 0` the Hamiltonian is proportional to the
identity, there is no flavor evolution, and the limit
:math:`\sin(|h| L)/|h| \to L` is used.
Examples
--------
.. jupyter-execute::
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))
"""
# [h1, h2, h3]
h_coeffs = hamiltonian_2nu_coefficients(hamiltonian_matrix)
# h_abs = |h|
h_abs = modulus(h_coeffs)
phase = h_abs*L
u0 = math.cos(phase)
# The limit of -sin(|h|L)/|h| as |h| -> 0 is -L
ss = -L if h_abs == 0.0 else -math.sin(phase)/h_abs
uk = [h_coeffs[k]*ss for k in range(0, 3)]
# [u0, u1, u2, u3]
return [u0]+uk
[docs]
def evolution_operator_2nu(
hamiltonian_matrix: Union[list, np.ndarray],
L: Union[int, float, list, np.ndarray]
) -> Union[List[List[complex]], np.ndarray]:
r"""Returns the two-neutrino time-evolution operator.
Returns the two-neutrino time-evolution operator :math:`U_2(L)` in
its SU(2) exponential expansion
:math:`U_2(L) = u_0 \mathbb{1} + i u_k \sigma^k`. This is a
:math:`2\times2` unitary matrix.
.. versionadded:: 1.0.0
.. versionchanged:: 1.1.0
Degenerate Hamiltonians are handled exactly instead of returning
NaN, by taking the limit :math:`\sin(|h|L)/|h| \to L`.
.. versionchanged:: 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.
.. versionchanged:: 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:
:func:`numpy.real`, :func:`numpy.imag`, :obj:`numpy.arccos`,
:func:`numpy.clip` and :obj:`numpy.sqrt` on one number give way
to attribute access and the :mod:`math` module.
.. versionchanged:: 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 :math:`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_matrix : array_like
Two-flavor Hermitian Hamiltonian, given as the nested list
``[[H11, H12], [H21, H22]]``, or a stack of them, of shape
``(..., 2, 2)``.
L : 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 time-evolution
operator :math:`U_2(L)` --- a :math:`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
--------
.. jupyter-execute::
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]))
"""
if _is_batched(hamiltonian_matrix, L):
return _evolution_operator_2nu_batch(hamiltonian_matrix, L)
if CHECK_HERMITICITY:
_check_hermitian(np.asarray(hamiltonian_matrix, dtype=complex),
'evolution_operator_2nu')
u0, u1, u2, u3 = \
evolution_operator_2nu_u_coefficients(hamiltonian_matrix, L)
return [
[u0+1.j*u3, 1.j*u1+u2],
[1.j*u1-u2, u0-1.j*u3]
]
[docs]
def probabilities_2nu(
hamiltonian_matrix: Union[list, np.ndarray],
L: Union[int, float, list, np.ndarray]
) -> Union[Tuple[float, float, float, float], np.ndarray]:
r"""Returns the two-neutrino oscillation probabilities.
Returns the two-neutrino flavor-transition probabilities
:math:`P_{ee}, P_{e\mu}, P_{\mu e}, P_{\mu\mu}`, where
:math:`P_{\alpha\beta} \equiv P(\nu_\alpha \to \nu_\beta)`.
.. versionadded:: 1.0.0
.. versionchanged:: 1.1.0
The :math:`h_2` contribution was restored. The transition
probability is :math:`|U_{\mu e}|^2 = u_1^2 + u_2^2`, but the
routine computed only :math:`|h_1|^2/|h|^2 \sin^2(|h|L)`. Since
:math:`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.
.. versionchanged:: 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.
.. versionchanged:: 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:
:func:`numpy.real`, :func:`numpy.imag`, :obj:`numpy.arccos`,
:func:`numpy.clip` and :obj:`numpy.sqrt` on one number give way
to attribute access and the :mod:`math` module. A scalar
two-flavor probability is 3.3x quicker, measured best-of-seven
interleaved against 1.3.0.
.. versionchanged:: 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 :math:`|U_{ee}|^2 = u_0^2 + u_3^2` and
:math:`|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 :math:`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.
.. versionchanged:: 1.6.0
Two dispatch decisions, neither changing the result. A stack of
at most :data:`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 :mod:`fastkernels` for the measured
thresholds and why the two-flavor one is high.
Parameters
----------
hamiltonian_matrix : array_like
Two-flavor Hermitian Hamiltonian, given as the nested list
``[[H11, H12], [H21, H22]]``, or a stack of them, of shape
``(..., 2, 2)``.
L : 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 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 :func:`oscprob3nu.probabilities_3nu` for the two scans
that broadcast naturally.
The transition probability is
.. math::
P_{e\mu} = \frac{|h_1|^2 + |h_2|^2}{|h|^2} \sin^2(|h| L) ,
i.e. :math:`|U_{\mu e}|^2 = u_1^2 + u_2^2`. Both :math:`h_1` and
:math:`h_2` contribute; :math:`h_2` vanishes only when the
off-diagonal entry of the Hamiltonian is real.
Examples
--------
.. jupyter-execute::
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))
"""
if _is_batched(hamiltonian_matrix, L):
return _probabilities_2nu_batch(hamiltonian_matrix, L)
if CHECK_HERMITICITY:
_check_hermitian(np.asarray(hamiltonian_matrix, dtype=complex),
'probabilities_2nu')
# [h1, h2, h3]
h_coeffs = hamiltonian_2nu_coefficients(hamiltonian_matrix)
# h_abs = |h|
h_abs = modulus(h_coeffs)
if h_abs == 0.0:
# The Hamiltonian is proportional to the identity: no flavor
# transitions occur, whatever the baseline.
Pem = 0.0
else:
h1, h2 = h_coeffs[0], h_coeffs[1]
sin_phase = math.sin(h_abs*L)
Pem = (h1*h1 + h2*h2)/(h_abs*h_abs) * sin_phase*sin_phase
Pme = Pem
Pee = 1.0-Pem
Pmm = 1.0-Pme
return Pee, Pem, Pme, Pmm