# -*- coding: utf-8 -*-
# SPDX-License-Identifier: GPL-3.0-only
# Copyright (C) 2026 Mauricio Bustamante
r"""hamiltonians2nu.py
Compute two-neutrino Hamiltonians for selected scenarios.
This module contains the routines to compute the two-neutrino
Hamiltonians for the following scenarios: oscillations in vacuum, in
matter of constant density, in matter with non-standard interactions
(NSI), and in a Lorentz invariance-violating background (LIV).
Routine listings
----------------
* mixing_matrix_2nu - Returns 2x2 rotation matrix
* hamiltonian_2nu_vacuum_energy_independent - Returns H_vac (no 1/E)
* hamiltonian_2nu_vacuum_energy_independent_td - Returns H_vac (no
1/E), as a function of position
* hamiltonian_2nu_vacuum - Returns H_vac
* hamiltonian_2nu_vacuum_td - Returns H_vac, as a function of position
* hamiltonian_2nu_matter - Returns H_matter
* hamiltonian_2nu_matter_td - Returns H_matter, as a function of position
* hamiltonian_2nu_nsi - Returns H_NSI
* hamiltonian_2nu_nsi_td - Returns H_NSI, as a function of position
* hamiltonian_2nu_liv - Returns H_LIV
* hamiltonian_2nu_liv_energy_independent - Returns H_LIV (no energy
dependence)
"""
__author__ = "Mauricio Bustamante"
__email__ = "mbustamante@gmail.com"
import numpy as np
from magnus.hamiltonians import _angles
from magnus.hamiltonians import _broadcast
from typing import Optional, Callable
from magnus import _validate as _v
# from globaldefs import *
[docs]
def mixing_matrix_2nu(sth: float, angles: Optional[str]='sin') -> np.ndarray:
r"""Returns the 2x2 rotation matrix.
Computes and returns a 2x2 real rotation matrix parametrized by a single rotation angle theta.
.. versionadded:: 1.0.0
.. versionchanged:: 1.2.0
A mixing angle beyond 90 degrees, whose cosine is negative, is refused (issue #160 §11).
Parameters
----------
sth : float
Mixing angle :math:`\theta`, in whichever convention ``angles`` names -- by
default its sine.
angles : str, optional
How ``sth`` is stated: ``'sin'`` (default) its sine, ``'sin2'`` its sine
*squared* -- which is what global fits report -- ``'rad'`` the angle itself in
radians, or ``'deg'`` in degrees. Any other value raises.
Returns
-------
np.ndarray
Rotation matrix [[cth, sth], [-sth, cth]], with cth = cos(theta) and sth = sin(theta).
Examples
--------
.. jupyter-execute::
import numpy as np
from magnus.hamiltonians import hamiltonians2nu
U = np.asarray(hamiltonians2nu.mixing_matrix_2nu(np.sqrt(0.308)))
print(np.round(U, 6))
print('unitary to %.1e' % np.max(np.abs(U.conj().T @ U - np.eye(2))))
Note the argument is :math:`\sin\theta`, not :math:`\sin^2\theta`; fits
are usually quoted as the latter -- pass ``angles='sin2'`` and hand the fitted
value straight in, rather than taking the square root yourself.
"""
# Checked for calls from outside the package only (issue #160 §11); the wrappers
# reach this with arguments already validated.
if not _v._called_from_inside(2):
_v.check_physics_params('hamiltonians.mixing_matrix_2nu', {'sth': sth})
sines, _ = _angles.resolve('hamiltonians.mixing_matrix_2nu', angles, {'sth': sth})
sth = sines['sth']
cth = np.sqrt(1.0-sth*sth)
return np.array([[cth,sth],[-sth,cth]])
[docs]
def hamiltonian_2nu_vacuum_energy_independent(sth: float, Dm2: float,
compute_matrix_multiplication: Optional[bool]=False,
angles: Optional[str]='sin') -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for vacuum oscillations.
Computes and returns the 2x2 real two-neutrino Hamiltonian for oscillations in vacuum,
parametrized by a single mixing angle theta and a single mass-squared difference Dm2. The
Hamiltonian is H = (1/2)*R.M2.R^dagger, with R the 2x2 rotation matrix and M2 the mass matrix.
The multiplicative factor 1/E is not applied.
.. versionadded:: 1.0.0
.. versionchanged:: 1.2.0
A mixing angle beyond 90 degrees, whose cosine is negative, is refused (issue #160 §11).
Parameters
----------
sth : float
Mixing angle :math:`\theta`, in whichever convention ``angles`` names -- by
default its sine.
Dm2 : float
Mass-squared difference :math:`\Delta m^2`.
compute_matrix_multiplication : bool, optional
If False (default), use the pre-computed expressions; otherwise,
multiply R.M2.R^dagger live.
angles : str, optional
How ``sth`` is stated: ``'sin'`` (default) its sine, ``'sin2'`` its sine
*squared* -- which is what global fits report -- ``'rad'`` the angle itself in
radians, or ``'deg'`` in degrees. Any other value raises.
:func:`magnus.globaldefs.load_nufit_params` returns sines, so its output goes
with the default.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix.
Examples
--------
.. jupyter-execute::
import numpy as np
import magnus.globaldefs as gd
from magnus.hamiltonians import hamiltonians2nu
h = np.asarray(hamiltonians2nu.hamiltonian_2nu_vacuum_energy_independent(
np.sqrt(0.308), 7.49e-5))
print(np.round(h/1e-5, 6), ' [1e-5 eV^2]')
print('at 1 GeV, in eV: %.3e' % (h/(1.0*gd.UNIT_GEV))[0][0])
The energy is divided out separately, which is what makes this the piece to
build once and reuse across a scan over energies.
"""
# Checked for calls from outside the package only (issue #160 §11); the wrappers
# reach this with arguments already validated.
if not _v._called_from_inside(2):
_v.check_physics_params('hamiltonians.hamiltonian_2nu_vacuum_energy_independent', {'sth': sth, 'Dm2': Dm2, 'compute_matrix_multiplication': compute_matrix_multiplication})
sines, _ = _angles.resolve(
'hamiltonians.hamiltonian_2nu_vacuum_energy_independent', angles, {'sth': sth})
sth = sines['sth']
cth = np.sqrt(1.0-sth*sth)
c2th = cth*cth-sth*sth
s2th = 2.0*cth*sth
if not compute_matrix_multiplication:
# Standard convention, Dm2 = m2^2 - m1^2 > 0 for m2 > m1 (matches the
# three-neutrino convention M2 = diag(0, D21, D31) and the closed-form
# matter probability in oscprobstd): H = (Dm2/4E) [[-c2th, s2th], [s2th, c2th]]
return (Dm2/4.0)*np.array([[-c2th,s2th], [s2th,c2th]])
else:
# 2D mixing matrix
R = mixing_matrix_2nu(sth)
# Mass matrix
M2 = np.diag([-1.0, 1.0]) # (m1^2 - m2^2, m2^2 - m1^2)/Dm2, traceless form
# Hamiltonian
return (Dm2/4.0) * R @ M2 @ R.T # Use matrix multiplication operator
# H = (Dm2/4.0)*np.matmul(R, np.matmul(M2, np.transpose(R)))
# return H
[docs]
def hamiltonian_2nu_vacuum_energy_independent_td(l: float, sth: float, Dm2: float,
compute_matrix_multiplication: Optional[bool]=False,
angles: Optional[str]='sin') -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for vacuum oscillations, as a function of distance,
even if it does not depend on it.
Computes and returns the 2x2 real two-neutrino Hamiltonian for oscillations in vacuum, as a
function of distance, parametrized by a single mixing angle theta and a single mass-squared
difference Dm2. The Hamiltonian is H = (1/2)*R.M2.R^dagger, with R the 2x2 rotation matrix and
M2 the mass matrix. The multiplicative factor 1/E is not applied. The vacuum Hamiltonian does
not depend on distance in reality, but we include the dependence here as a way to validate the
routine to compute probabilities for time-dependent Hamiltonians.
.. versionadded:: 1.0.0
.. versionchanged:: 1.2.0
``l`` may be an array: the result is a stack of matrices, one per position (issue #155
§2).
Parameters
----------
l : float or array_like
Position at which the Hamiltonian is evaluated.
An array returns a stack of matrices, one per position.
sth : float
Mixing angle :math:`\theta`, in whichever convention ``angles`` names -- by
default its sine.
Dm2 : float
Mass-squared difference :math:`\Delta m^2`.
compute_matrix_multiplication : bool, optional
If False (default), use the pre-computed expressions; otherwise, multiply R.M2.R^dagger
live.
angles : str, optional
How ``sth`` is stated: ``'sin'`` (default) its sine, ``'sin2'`` its sine
*squared* -- which is what global fits report -- ``'rad'`` the angle itself in
radians, or ``'deg'`` in degrees. Any other value raises.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix, or a stack of them, shape ``(..., 2, 2)``, for array input.
"""
return _broadcast.over_positions(l, hamiltonian_2nu_vacuum_energy_independent(sth, Dm2,
compute_matrix_multiplication=compute_matrix_multiplication, angles=angles))
[docs]
def hamiltonian_2nu_vacuum(energy: float, sth: float, Dm2: float,
compute_matrix_multiplication: Optional[bool]=False,
angles: Optional[str]='sin') -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for vacuum oscillations.
Same as :func:`hamiltonian_2nu_vacuum_energy_independent`, but with the 1/E factor applied.
.. versionadded:: 1.0.0
.. versionchanged:: 1.2.0
``energy`` may be an array: the result is a stack of matrices, one per energy (issue #155
§2).
Parameters
----------
energy : float or array_like
Neutrino energy.
An array returns a stack of matrices, one per energy.
sth : float
Mixing angle :math:`\theta`, in whichever convention ``angles`` names -- by
default its sine.
Dm2 : float
Mass-squared difference :math:`\Delta m^2`.
compute_matrix_multiplication : bool, optional
If False (default), use the pre-computed expressions; otherwise,
multiply R.M2.R^dagger live.
angles : str, optional
How ``sth`` is stated: ``'sin'`` (default) its sine, ``'sin2'`` its sine
*squared* -- which is what global fits report -- ``'rad'`` the angle itself in
radians, or ``'deg'`` in degrees. Any other value raises.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix, or a stack of them, shape ``(..., 2, 2)``, for array input.
"""
# Runs at every quadrature node inside a user Hamiltonian, so only the energy is
# checked here, in one comparison (issue #160 §11); the other arguments are
# checked by the energy-independent builder.
if not (type(energy) is float and 0.0 < energy < _v._INF):
_v.check_physics_params('hamiltonians.hamiltonian_2nu_vacuum', {'energy': energy})
if type(energy) not in _broadcast.SCALARS and np.ndim(energy):
energy = _broadcast.stacked(energy)
return (1/energy)*hamiltonian_2nu_vacuum_energy_independent(sth, Dm2,
compute_matrix_multiplication=compute_matrix_multiplication, angles=angles)
[docs]
def hamiltonian_2nu_vacuum_td(l: float, energy: float, sth: float, Dm2: float,
compute_matrix_multiplication: Optional[bool]=False,
angles: Optional[str]='sin') -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for vacuum oscillations, as a function of distance,
even if it does not depend on it.
Computes and returns the 2x2 real two-neutrino Hamiltonian for oscillations in vacuum, as a
function of distance, parametrized by a single mixing angle theta and a single mass-squared
difference Dm2. The Hamiltonian is H = (1/2)*R.M2.R^dagger, with R the 2x2 rotation matrix and
M2 the mass matrix. The multiplicative factor 1/E is not applied. The vacuum Hamiltonian does
not depend on distance in reality, but we include the dependence here as a way to validate the
routine to compute probabilities for time-dependent Hamiltonians.
.. versionadded:: 1.0.0
.. versionchanged:: 1.2.0
``l`` may be an array: the result is a stack of matrices, one per position (issue #155
§2).
Parameters
----------
l : float or array_like
Position at which the Hamiltonian is evaluated.
An array returns a stack of matrices, one per position.
energy : float or array_like
Neutrino energy.
An array returns a stack of matrices, one per energy.
sth : float
Mixing angle :math:`\theta`, in whichever convention ``angles`` names -- by
default its sine.
Dm2 : float
Mass-squared difference :math:`\Delta m^2`.
compute_matrix_multiplication : bool, optional
If False (default), use the pre-computed expressions; otherwise,
multiply R.M2.R^dagger live.
angles : str, optional
How ``sth`` is stated: ``'sin'`` (default) its sine, ``'sin2'`` its sine
*squared* -- which is what global fits report -- ``'rad'`` the angle itself in
radians, or ``'deg'`` in degrees. Any other value raises.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix, or a stack of them, shape ``(..., 2, 2)``, for array input.
"""
return _broadcast.over_positions(l, hamiltonian_2nu_vacuum(energy, sth, Dm2,
compute_matrix_multiplication=compute_matrix_multiplication, angles=angles))
[docs]
def hamiltonian_2nu_matter(VCC: float) -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for matter oscillations.
Not validated: called on the hot path, at every quadrature node. Its inputs are
checked where they are set (issue #160 §11).
Computes and returns the 2x2 real two-neutrino Hamiltonian for oscillations in matter with
constant density.
.. versionadded:: 1.0.0
Parameters
----------
VCC : float or array_like
Potential due to charged-current interactions of nu_e with electrons.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix, or a stack of them, shape ``(..., 2, 2)``, for array input.
Examples
--------
.. jupyter-execute::
import numpy as np
from magnus.hamiltonians import hamiltonians2nu
print(np.asarray(hamiltonians2nu.hamiltonian_2nu_matter(1.0e-13)))
Only the electron-flavor entry is filled: matter is what the other
flavors do not feel.
"""
# The matter Hamiltonian is [[VCC,0],[0,0]]
# Built by broadcasting rather than np.diag so that VCC may be an array of
# positions: VCC[..., None, None] turns one potential per position into a
# stack of matrices, which is what lets a caller's H_func take the engine's
# vectorized path (see magnus.magnus.ScalarHamiltonianWarning). A scalar VCC
# still returns a plain (2, 2) matrix.
VCC = np.asarray(VCC, dtype=float)
e00 = np.zeros((2, 2))
e00[0, 0] = 1.0
return VCC[..., None, None] * e00
# h_matter = np.zeros((2,2))
# # Add the matter potential to the ee term to find the matter Hamiltonian
# h_matter[0][0] = VCC
# return h_matter
[docs]
def hamiltonian_2nu_matter_td(l: float, VCC_func: Callable) -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for matter oscillations, as a function of distance.
Not validated: called on the hot path, at every quadrature node. Its inputs are
checked where they are set (issue #160 §11).
Computes and returns the 2x2 real two-neutrino Hamiltonian for
oscillations in matter with a given density as a function of
position.
.. versionadded:: 1.0.0
Parameters
----------
l : float or array_like
Position at which the Hamiltonian is evaluated.
An array returns a stack of matrices, one per position.
VCC_func : Callable
Potential due to charged-current interactions of nu_e with electrons, as a function of
position, l.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix, or a stack of them, shape ``(..., 2, 2)``, for array input.
"""
return hamiltonian_2nu_matter(VCC_func(l))
[docs]
def hamiltonian_2nu_nsi(VCC: float, eps_aa: float, eps_ab: complex) -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for oscillations with NSI.
Computes and returns the 2x2 complex two-neutrino Hamiltonian for oscillations with
non-standard interactions (NSI) in matter with constant density.
Two flavors admit only one physically meaningful diagonal NSI degree of freedom: an overall
(flavor-universal) diagonal shift is proportional to the identity matrix, so it commutes with
everything, contributes only an unobservable global phase to the evolution operator, and cannot
affect any oscillation probability. ``eps_aa`` is therefore defined here as the non-universal
(flavor-off-diagonal-*difference*) coupling, following the convention eps_mumu = 0, i.e., it
parametrizes the coupling of :math:`\nu_e` alone, relative to :math:`\nu_\mu`. [Earlier
versions of this function placed eps_aa on *both* diagonal entries, making it a pure multiple
of the identity and therefore a no-op on every oscillation probability -- this was a bug, not a
convention choice, confirmed by direct calculation.]
.. versionadded:: 1.0.0
.. versionchanged:: 1.2.0
``VCC`` may be an array: the result is a stack of matrices, one per potential (issue #155
§2).
Parameters
----------
VCC : float or array_like
Potential due to charged-current interactions of nu_e with electrons.
eps_aa : float
Non-universal diagonal NSI coupling of nu_e (relative to nu_mu, whose diagonal coupling is
fixed to 0 by this convention).
eps_ab : complex
Flavor-off-diagonal (nu_e-nu_mu) NSI coupling.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix, or a stack of them, shape ``(..., 2, 2)``, for array input.
"""
# Runs at every quadrature node inside a user Hamiltonian, so only what a plain
# float cannot get wrong is left to check, in one comparison per coupling: a
# complex diagonal coupling makes H non-Hermitian (issue #160 §2).
if type(eps_aa) is not float:
_v.check_physics_params('hamiltonians.hamiltonian_2nu_nsi', {'eps_aa': eps_aa})
if type(VCC) not in _broadcast.SCALARS and np.ndim(VCC):
VCC = _broadcast.stacked(VCC)
return VCC * np.array([[eps_aa, eps_ab], [np.conj(eps_ab), 0.0]], dtype=np.complex128)
[docs]
def hamiltonian_2nu_nsi_td(l: float, VCC_func: Callable, eps_aa: float,
eps_ab: complex) -> np.ndarray:
r"""Returns the two-neutrino NSI Hamiltonian as a function of position.
Same as :func:`hamiltonian_2nu_nsi`, but evaluates the position-dependent matter potential
``VCC_func(l)`` first.
.. versionadded:: 1.0.0
Parameters
----------
l : float or array_like
Position at which the Hamiltonian is evaluated.
An array returns a stack of matrices, one per position.
VCC_func : Callable
Potential due to charged-current interactions of nu_e with electrons, as a function of
position, l.
eps_aa : float
Non-universal diagonal NSI coupling of nu_e; see :func:`hamiltonian_2nu_nsi`.
eps_ab : complex
Flavor-off-diagonal (nu_e-nu_mu) NSI coupling.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix, or a stack of them, shape ``(..., 2, 2)``, for array input.
"""
return hamiltonian_2nu_nsi(VCC_func(l), eps_aa, eps_ab)
[docs]
def hamiltonian_2nu_liv(energy: float, sxi: float, b1: float, b2: float, Lambda: float, n_liv: int,
nubar: Optional[bool]=False, angles: Optional[str]='sin') -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for oscillations with LIV.
Computes and returns the 2x2 real two-neutrino Hamiltonian for oscillations in a Lorentz
invariance-violating background. Same as
:func:`hamiltonian_2nu_liv_energy_independent`, but with the
:math:`E^{n_{\rm liv}}` energy dependence of the LIV operator applied.
.. versionadded:: 1.0.0
.. versionchanged:: 1.2.0
``energy`` may be an array: the result is a stack of matrices, one per energy (issue #155
§2).
``nubar`` takes effect: for antineutrinos, the term changes sign when ``n_liv`` is even
(a CPT-odd operator). It used to be ignored (issue #162).
Parameters
----------
energy : float or array_like
Neutrino energy.
An array returns a stack of matrices, one per energy.
sxi : float
Rotation angle :math:`\xi` between the space of the eigenvectors of B2 and the
flavor states, in whichever convention ``angles`` names -- by default its sine.
b1 : float
Eigenvalue b1 of the LIV operator B2.
b2 : float
Eigenvalue b2 of the LIV operator B2.
Lambda : float
Energy scale of the LIV operator B2.
n_liv : int
Power of the energy dependence of the LIV operator (dimension of the operator minus 3).
nubar : bool, optional
If True, compute the Hamiltonian for antineutrinos: when ``n_liv`` is even, the term
changes sign; the real 2-flavor rotation has nothing to conjugate. An operator of
dimension ``n_liv + 3`` is CPT-odd for odd dimension (even ``n_liv``) and
CPT-even for even dimension, as in the Standard-Model Extension (Kostelecky & Mewes, Phys.
Rev. D 85, 096005 (2012)). Default: False.
angles : str, optional
How ``sxi`` is stated: ``'sin'`` (default) its sine, ``'sin2'`` its sine
*squared*, ``'rad'`` the angle itself in radians, or ``'deg'`` in degrees. Any
other value raises.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix, or a stack of them, shape ``(..., 2, 2)``, for array input.
"""
# Runs at every quadrature node inside a user Hamiltonian, so only the energy is
# checked here, in one comparison (issue #160 §11); the other arguments are
# checked by the energy-independent builder.
if not (type(energy) is float and 0.0 < energy < _v._INF):
_v.check_physics_params('hamiltonians.hamiltonian_2nu_liv', {'energy': energy})
if type(energy) not in _broadcast.SCALARS and np.ndim(energy):
energy = _broadcast.stacked(energy)
return pow(energy, n_liv) * hamiltonian_2nu_liv_energy_independent(
sxi, b1, b2, Lambda, n_liv, nubar=nubar, angles=angles)
[docs]
def hamiltonian_2nu_liv_energy_independent(sxi: float, b1: float, b2: float,
Lambda: float, n_liv: int, nubar: Optional[bool]=False,
angles: Optional[str]='sin') -> np.ndarray:
r"""Returns the two-neutrino Hamiltonian for oscillations with LIV.
Computes and returns the 2x2 real two-neutrino Hamiltonian for oscillations in a Lorentz
invariance-violating background, without the energy-dependent prefactor.
.. versionadded:: 1.0.0
.. versionchanged:: 1.2.0
A LIV mixing angle beyond 90 degrees, whose cosine is negative, is refused (issue #160
§11).
Added ``nubar``: for antineutrinos, the term changes sign when ``n_liv`` is even (a
CPT-odd operator) (issue #162).
Parameters
----------
sxi : float
Rotation angle :math:`\xi` between the space of the eigenvectors of B2 and the
flavor states, in whichever convention ``angles`` names -- by default its sine.
b1 : float
Eigenvalue b1 of the LIV operator B2.
b2 : float
Eigenvalue b2 of the LIV operator B2.
Lambda : float
Energy scale of the LIV operator B2.
n_liv : int
Power of the energy dependence of the LIV operator (dimension of the operator minus 3).
nubar : bool, optional
If True, compute the Hamiltonian for antineutrinos: when ``n_liv`` is even, the term
changes sign; the real 2-flavor rotation has nothing to conjugate. An operator of
dimension ``n_liv + 3`` is CPT-odd for odd dimension (even ``n_liv``) and
CPT-even for even dimension, as in the Standard-Model Extension (Kostelecky & Mewes, Phys.
Rev. D 85, 096005 (2012)). Default: False.
angles : str, optional
How ``sxi`` is stated: ``'sin'`` (default) its sine, ``'sin2'`` its sine
*squared*, ``'rad'`` the angle itself in radians, or ``'deg'`` in degrees. Any
other value raises.
Returns
-------
np.ndarray
Hamiltonian 2x2 matrix.
"""
# Checked for calls from outside the package only (issue #160 §11); the wrappers
# reach this with arguments already validated.
if not _v._called_from_inside(2):
_v.check_physics_params('hamiltonians.hamiltonian_2nu_liv_energy_independent', {'sxi': sxi, 'b1': b1, 'b2': b2, 'Lambda': Lambda, 'n_liv': n_liv, 'nubar': nubar})
# Until now the LIV angle went through no guard at all: only the vacuum builders
# validated their sines, so an out-of-range sxi reached np.sqrt(1 - sxi*sxi) and
# became a NaN Hamiltonian in silence. resolve() closes that as a side effect of
# accepting the convention.
_sines, _ = _angles.resolve(
'hamiltonians.hamiltonian_2nu_liv_energy_independent', angles, {'sxi': sxi})
sxi = _sines['sxi']
# H = R . diag(b1, b2) . R^T, with R = mixing_matrix_2nu(sxi) -- the same convention used by
# every sibling Hamiltonian (2nu vacuum's slow path, and the 3/4/5nu LIV Hamiltonians). The
# off-diagonal sign below was previously flipped relative to this convention (a confirmed bug).
cxi = np.sqrt(1.0 - sxi * sxi)
delta_b = b2 - b1
# Antineutrinos (issue #162): the operator is transposed in flavor space, which leaves this
# real symmetric matrix as it is, and its CPT-odd part changes sign (Kostelecky & Mewes,
# Phys. Rev. D 85, 096005 (2012), Eq. 78). An operator of dimension n_liv + 3 is CPT-odd
# when that is odd, i.e. when n_liv is even.
sign = -1.0 if (nubar and n_liv % 2 == 0) else 1.0
return sign * pow(1.0 / Lambda, n_liv) * np.array([
[b1 * cxi * cxi + b2 * sxi * sxi, delta_b * cxi * sxi],
[delta_b * cxi * sxi, b2 * cxi * cxi + b1 * sxi * sxi]
], dtype=np.float64)