magnus.avgprob

avgprob.py

Contains the phase-averaged oscillation probabilities: the phase average over a relative energy spread, which average=True returns, and the exact \(L/E \to \infty\) limit it reduces to where every phase has decohered, as for high-energy astrophysical neutrinos.

Physical idea: a neutrino produced at a cosmological distance arrives with an oscillation phase \(\Delta m^2 L / 2E\) of order \(10^{15}\) or more, and neither the source distance, nor the production region, nor the detector’s energy resolution is known to anything close to that precision. Every oscillatory term is therefore averaged over many cycles and vanishes, leaving only the incoherent sum

\[P(\nu_\alpha \to \nu_\beta) = \sum_i |V_{\alpha i}|^2 |V_{\beta i}|^2 ,\]

where \(V\) diagonalizes the Hamiltonian. This is not an approximation to be refined: it is the exact limit, and it costs one matrix product rather than an integration. For standard vacuum oscillations the result does not depend on energy or baseline at all, so a single matrix serves an entire flux calculation.

Coherence is decided physically, not numerically

The formula above assumes every relative phase averages away. That is a statement about pairs of eigenvalues, not about the spectrum as a whole: the pair \((i,j)\) decoheres only if \((\lambda_i - \lambda_j) L\) sweeps through many cycles across the averaging window. Two eigenvalues that are close enough to keep their relative phase fixed stay coherent, and their cross term survives.

This module therefore groups the spectrum into blocks of mutually coherent eigenvalues and sums coherently inside each block,

\[P(\nu_\alpha \to \nu_\beta) = \sum_{b} \Big| \sum_{i \in b} V^*_{\alpha i} V_{\beta i} \Big|^2 ,\]

which reduces to the familiar expression when every block is a singleton. The distinction is not academic here: a sterile state with a small \(\Delta m^2_{41}\), or any degenerate spectrum, makes the naive sum quietly wrong.

The same per-pair phase decides whether the limit applies at all. A pair whose phase is neither much larger than \(2\pi\) (decohered) nor much smaller than one (coherent) sits in between, where the limit does not describe it; coherence_report() names those pairs.

The phase average

A measurement with a relative energy resolution \(\sigma\) averages a phase over the range it covers across that resolution, and does not average one that barely changes. The phase average keeps every interference term with its phase at the central energy and weights it by \(e^{-\sigma^2\phi'^2/2}\), \(\phi' = d\phi/d\ln E\): the limit above where every phase runs through many cycles, the oscillation probability where none does, and a smooth weighting between. Mixing and the eigenbases stay at the central energy, so a probability without interference is unchanged. phase_averaged_probabilities_constant_hamiltonian() and phase_averaged_probabilities_adiabatic() compute it, and they are what average=True in magnus.oscprob returns; the pairs coherence_report() names are where the result depends on \(\sigma\), and magnus.oscprob warns there. The functions that return the limit are unchanged.

This module stands apart from magnus.oscprob, so it can be applied to any Hermitian Hamiltonian of any dimension independently of the rest of the API. It depends on numpy and on magnus.adiabatic; everything except level_crossing_matrix() and averaged_probabilities_adiabatic() needs numpy alone.

Routine listings

coherence_blocks coherence_report averaged_probabilities_from_eigenbasis averaged_probabilities_constant_hamiltonian adiabatic_phase_differences level_crossing_matrix averaged_probabilities_adiabatic averaged_probabilities_numerically phase_averaged_probabilities_constant_hamiltonian phase_averaged_probabilities_adiabatic

Attributes

DECOHERENCE_PHASE_THRESHOLD	Module-level constant
COHERENCE_PHASE_THRESHOLD	Module-level constant
AVG_DEFAULT_ENERGY_SPREAD	Module-level constant
AVG_DEFAULT_N_SAMPLES	Module-level constant
SUDDEN_TRANSFER_THRESHOLD	Module-level constant
AVG_PHASE_SPREAD	Module-level constant
PHASE_AVERAGE_WINDOW_THRESHOLD	Module-level constant
PHASE_AVERAGE_PATCH_ATOL	Module-level constant
PHASE_AVERAGE_PHASE_TOL	Module-level constant
PHASE_SPREAD_SENSITIVITY_THRESHOLD	Module-level constant

Functions

coherence_blocks(→ List[List[int]])	Groups eigenvalues into blocks that stay mutually coherent.
coherence_report(→ Tuple[List[List[int]], ...)	Reports the coherence structure of a spectrum, and which pairs sit in
averaged_probabilities_from_eigenbasis(→ numpy.ndarray)	Phase-averaged oscillation probabilities from the eigenbasis of the
averaged_probabilities_constant_hamiltonian(...)	Phase-averaged oscillation probabilities for a constant Hamiltonian.
averaged_probabilities_numerically(...)	Averages a probability by sampling it across an energy window.
adiabatic_phase_differences(→ numpy.ndarray)	Relative phases accumulated between instantaneous eigenvalues.
level_crossing_matrix(→ Tuple[numpy.ndarray, ...)	Probability of ending on level \(j\) having started on level \(i\).
averaged_probabilities_adiabatic(...)	Phase-averaged probabilities for a position-dependent Hamiltonian.
phase_averaged_probabilities_constant_hamiltonian(...)	Phase-averaged probabilities for a constant Hamiltonian, from a flavor state at the start.
phase_averaged_probabilities_adiabatic(...)	Phase-averaged probabilities on a smooth position-dependent Hamiltonian.

Module Contents

magnus.avgprob.DECOHERENCE_PHASE_THRESHOLD

Module-level constant

Accumulated phase spread, in radians, above which a pair of eigenvalues is treated as fully decohered. One full cycle is the point at which the average of \(\cos\Delta\phi\) over the window has collapsed to a small fraction of its coherent value, and every further cycle only reduces it.

Added in version 1.0.0.

Type:

float

magnus.avgprob.COHERENCE_PHASE_THRESHOLD = 0.01

Module-level constant

Accumulated phase spread, in radians, below which a pair of eigenvalues is treated as fully coherent, so that its cross term is kept in full.

The gap between this and DECOHERENCE_PHASE_THRESHOLD is deliberate and is not a tolerance to be tightened away: a pair falling between the two is in neither limit, and no averaged expression describes it. Such a pair is still placed in a block – coherent below the decoherence threshold, decohered at or above it – so the accompanying number is a definite choice; what coherence_report() adds is that the choice is not made silently.

Added in version 1.0.0.

Type:

float

magnus.avgprob.coherence_blocks(eigenvalues: Sequence[float] | numpy.ndarray, phase_scale: float, decoherence_threshold: float | None = DECOHERENCE_PHASE_THRESHOLD) → List[List[int]]

Groups eigenvalues into blocks that stay mutually coherent.

Two eigenvalues belong to the same block when the phase they accumulate relative to each other, \(|\lambda_i - \lambda_j| \times\) phase_scale, stays below decoherence_threshold, so that their cross term in the probability is not averaged away.

Grouping is by transitive closure over that relation, which is the conservative choice: a chain of individually-close eigenvalues is kept in one block rather than split at an arbitrary point. A spectrum whose spacings are all comparable to the threshold therefore collapses into a single block, and is exactly the case coherence_report() flags as having no valid averaged limit.

Added in version 1.0.0.

Parameters:

• eigenvalues (list or np.ndarray) – Eigenvalues of the Hamiltonian [eV]. Need not be sorted.

• phase_scale (float) – Baseline over which the phase accumulates [\(\text{eV}^{-1}\)], so that (lambda_i - lambda_j)*phase_scale is a phase in radians.

• decoherence_threshold (float, optional) – Phase above which a pair is treated as decohered. Default: DECOHERENCE_PHASE_THRESHOLD.

Returns:

Indices of eigenvalues, grouped into blocks and sorted within each block. The blocks themselves are ordered by their smallest index, so the result is deterministic.

Return type:

list of list of int

Examples

A spectrum whose splittings are all large is fully decohered, one index per block; two eigenvalues sharing a value stay together.

import magnus.avgprob as ap

ap.coherence_blocks([0.0, 1.0, 2.0], phase_scale=1.0e3)

[[0], [1], [2]]

magnus.avgprob.coherence_report(eigenvalues: Sequence[float] | numpy.ndarray, phase_scale: float, decoherence_threshold: float | None = DECOHERENCE_PHASE_THRESHOLD, coherence_threshold: float | None = COHERENCE_PHASE_THRESHOLD) → Tuple[List[List[int]], List[Tuple[int, int, float]]]

Reports the coherence structure of a spectrum, and which pairs sit in neither limit.

Every pair of eigenvalues is in one of three regimes, set by the phase it accumulates relative to the others over phase_scale:

• far above decoherence_threshold, the cross term has averaged away and the pair contributes incoherently;

• far below coherence_threshold, the relative phase has barely advanced and the pair is still fully coherent;

• in between, neither statement holds, and no averaged expression is a valid description – the honest answer there is the full oscillation probability, not an average.

Added in version 1.0.0.

Parameters:

• eigenvalues (list or np.ndarray) – Eigenvalues of the Hamiltonian [eV].

• phase_scale (float) – Baseline over which the phase accumulates [\(\text{eV}^{-1}\)].

• decoherence_threshold (float, optional) – Phase above which a pair counts as decohered. Default: DECOHERENCE_PHASE_THRESHOLD.

• coherence_threshold (float, optional) – Phase below which a pair counts as fully coherent. Default: COHERENCE_PHASE_THRESHOLD.

Returns:

The coherence blocks, and the list of (i, j, phase) triples for pairs that are in neither limit. An empty second element means no pair sits between the two thresholds; the averaged result is then exact up to the residual the thresholds themselves allow, not exactly exact.

Return type:

(list of list of int, list of (int, int, float))

magnus.avgprob.averaged_probabilities_from_eigenbasis(eigenvectors: Sequence | numpy.ndarray, blocks: List[List[int]] | None = None) → numpy.ndarray

Phase-averaged oscillation probabilities from the eigenbasis of the Hamiltonian.

Computes

\[P_{\alpha\beta} = \sum_b \Big| \sum_{i \in b} V^*_{\alpha i} V_{\beta i} \Big|^2 ,\]

the sum over coherence blocks b of the squared modulus of the coherent amplitude within each block. With one index per block this is the familiar \(\sum_i |V_{\alpha i}|^2 |V_{\beta i}|^2\).

The result is symmetric, so the averaged probability is the same in both directions, and identical for neutrinos and antineutrinos: conjugating \(V\) leaves every term unchanged. CP violation does not survive the average, even though the mixing angles and phases do enter through \(|V_{\alpha i}|\).

Added in version 1.0.0.

Parameters:

• eigenvectors (list or np.ndarray) – Matrix whose columns are the eigenvectors of the Hamiltonian, shape (..., d, d). A leading batch axis is allowed and is broadcast over, so an array of energies costs one contraction.

• blocks (list of list of int, optional) – Coherence blocks, as returned by coherence_blocks(). If None (default), every eigenvalue is assumed to have decohered from every other, which is the astrophysical case.

Returns:

Averaged probability matrix, shape (..., d, d), with the initial flavor as the row index, so each row sums to one.

Return type:

np.ndarray

Examples

import numpy as np

import magnus.avgprob as ap
import magnus.hamiltonians as hams

U = hams.pmns_mixing_matrix(0.55, 0.68, 0.15, 3.7)
P = ap.averaged_probabilities_from_eigenbasis(U)
np.round(P, 4)

array([[0.5528, 0.2153, 0.2319],
       [0.2153, 0.4072, 0.3775],
       [0.2319, 0.3775, 0.3906]])

magnus.avgprob.averaged_probabilities_constant_hamiltonian(hamiltonian: Sequence | numpy.ndarray, baseline: float | None = None) → numpy.ndarray

Phase-averaged oscillation probabilities for a constant Hamiltonian.

Diagonalizes hamiltonian and applies averaged_probabilities_from_eigenbasis(). This covers every position-independent case – vacuum, matter of constant density, and their NSI and LIV variants – exactly, at the cost of one eigendecomposition.

Added in version 1.0.0.

Parameters:

• hamiltonian (list or np.ndarray) – Hermitian Hamiltonian [eV], shape (..., d, d). A leading batch axis (energies, say) is allowed.

• baseline (float, optional) – Baseline [\(\text{eV}^{-1}\)], used only to decide which eigenvalues have decohered from each other. If None (default), every pair is taken to be decohered, which is the astrophysical limit and makes the result independent of distance. Only for a single Hamiltonian: giving a baseline for a batch raises, since the coherence structure may differ from one entry to the next.

Returns:

Averaged probability matrix, shape (..., d, d), rows summing to one.

Return type:

np.ndarray

magnus.avgprob.AVG_DEFAULT_ENERGY_SPREAD = 0.1

Module-level constant

Half-width of the energy window, as a fraction of the energy, used when the averaged probability has to be obtained by sampling rather than in closed form.

Ten per cent is the order of a real detector’s energy resolution, and it is the smearing that does the averaging: the physical statement is that the oscillation phase varies by many cycles across whatever window the measurement integrates over. It is a default, not a property of the physics, so it is named here rather than buried, every use of it through the magnus.oscprob entry points is warned about, and callers with an actual resolution should pass theirs. Calling this module directly warns nobody: the width and the standard error come back in the result instead.

Added in version 1.0.0.

Type:

float

magnus.avgprob.AVG_DEFAULT_N_SAMPLES = 41

Module-level constant

Number of samples across the window used by averaged_probabilities_numerically().

The sampled phases are effectively independent when the accumulated phase is large, so the error of the mean falls only as \(1/\sqrt{N}\) – 41 samples give a few per cent. Raising it buys accuracy slowly and costs a full propagation each; the closed-form paths in this module exist precisely to avoid this trade.

Added in version 1.0.0.

Type:

int

magnus.avgprob.averaged_probabilities_numerically(prob_of_energy: Callable, energy: float, relative_spread: float | None = AVG_DEFAULT_ENERGY_SPREAD, n_samples: int | None = AVG_DEFAULT_N_SAMPLES) → Tuple[numpy.ndarray, float]

Averages a probability by sampling it across an energy window.

The fallback for cases with no closed form – a profile with discontinuities, say, where there is no instantaneous eigenbasis to decohere in. Unlike the closed forms in this module, this is not the \(L/E \to \infty\) limit: it is the average over a particular window, and the answer depends on that window. Its width is therefore an argument, and callers that leave it at the default should say so to their own callers.

Samples are uniform in \(1/E\), in which the oscillation phase is linear, so they are spread evenly in phase rather than bunched.

Added in version 1.0.0.

Parameters:

• prob_of_energy (Callable) – Returns the probability matrix at a given energy; called once per sample.

• energy (float) – Central energy [eV].

• relative_spread (float, optional) – Half-width of the window as a fraction of energy. Default: AVG_DEFAULT_ENERGY_SPREAD.

• n_samples (int, optional) – Number of samples. Default: AVG_DEFAULT_N_SAMPLES.

Returns:

The mean probability matrix, and the largest standard error of the mean across its entries – the honest uncertainty of the result, which a closed form would not have.

Return type:

(np.ndarray, float)

magnus.avgprob.adiabatic_phase_differences(H_func: Callable, l0: float, l1: float, n_points: int | None = 201) → numpy.ndarray

Relative phases accumulated between instantaneous eigenvalues.

In the adiabatic regime a neutrino stays on one level and accumulates the dynamical phase \(\int \lambda_i(l)\, dl\), so the phase that decides whether levels \(i\) and \(j\) still interfere is \(\Delta\phi_{ij} = \int_{l_0}^{l_1} [\lambda_i(l) - \lambda_j(l)]\, dl\). That integral, not the eigenvalue gap at any single point, is what the coherence tests in this module need for a position-dependent Hamiltonian.

Integrated with Simpson’s rule: the trapezoid leaves a residual here that is easily mistaken for a physical effect (the same error, in the same integral, once looked like a floor on the accuracy of adiabatic transport in magnus.adiabatic).

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian as a function of position, H_func(l) [eV].

• l0 (float) – Start and end of the trajectory [\(\text{eV}^{-1}\)].

• l1 (float) – Start and end of the trajectory [\(\text{eV}^{-1}\)].

• n_points (int, optional) – Number of sampling points. Truncated to an integer, raised to 3 if smaller, then raised to the next odd number for Simpson’s rule – all three silently. Default: 201.

Returns:

Matrix of accumulated phase differences, shape (d, d), antisymmetric.

Return type:

np.ndarray

magnus.avgprob.SUDDEN_TRANSFER_THRESHOLD = 0.001

Module-level constant

How much probability a feature must be able to move between levels before averaged_probabilities_adiabatic() stops trusting its own 200-point search to have seen it.

That search looks for non-adiabatic windows once, on a fixed probe grid. A density front narrower than the probe spacing falls between two probes and is never examined: no window opens, \(P^\text{cross}\) is the identity, and the answer is the fully adiabatic one, returned without a warning (issue #60). Such a front can be seen cheaply – one half of a probe interval carries nearly all of that interval’s change – but seeing it is not enough, because a solar-model table interpolated in log-density shows the same shape at every one of its grid points in the core, where nothing happens. What separates the two is whether the feature could move a neutrino at all. An instantaneous change from \(H(l_a)\) to \(H(l_b)\) moves at most \(\max_{i\ne j}|\langle v_i(l_a)|v_j(l_b)\rangle|^2\) between levels, and a monotone passage between the two positions moves less; when even that bound is below this threshold, the feature cannot change the averaged probability by more than the default tolerance, and today’s answer stands.

Measured (docs/dev/adversarial_batteries/sudden_transfer_sweep.py) as the largest bound over the intervals the probe grid finds concentrated:

population	largest bound	averaged answer today
BS05 solar model, cubic and linear, 1-30 MeV (core)	4.8e-07	escalating moves it <= 1.7e-16
issue #60’s shock ray, the 16 of 24 fronts it got wrong	0.20 - 0.74	off by 0.08 - 0.56
the same ray, broad fronts it got right	9.3e-03 - 0.70	right
supernova turbulence, 5-30 MeV	0.043 - 0.93	off by up to 0.25
Earth crust with undeclared layer edges, 5-30 MeV	7.3e-06 - 4.7e-03	not scored; 1 of 9 escalates

1e-3, the default tolerance, sits three orders above the solar model and two below the smallest front the engine gets wrong. Across the roughly 5,000 averaged calls the notebooks make, it escalates 16 pixels of paper Figure 5f and nothing else; three of those change, each to within 0.001 of a decohered reference.

Added in version 1.1.1.

Type:

float

magnus.avgprob.level_crossing_matrix(H_func: Callable, l0: float, l1: float, threshold: float | None = 0.1, n_probe: int | None = 200, fd_step_frac: float | None = 1e-06, magnus_exp_order: int | None = 6, integration_method: str | None = 'gl') → Tuple[numpy.ndarray, List[Tuple[float, float]], bool]

Probability of ending on level \(j\) having started on level \(i\).

Adiabatic evolution keeps a neutrino on the level it was produced on, so this matrix is the identity wherever the adiabatic approximation holds. It departs from the identity only across a non-adiabatic window – a resonance sharp enough for levels to exchange character faster than the state can follow – and it is exactly there that the averaged probability needs it.

The window is located with the Hellmann-Feynman diagnostic in magnus.adiabatic, and the transfer across it is computed with that module’s own convergence-checked Magnus patch rather than with a Landau-Zener formula, so it inherits an exact treatment of the crossing instead of an asymptotic approximation to it.

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian as a function of position, H_func(l) [eV].

• l0 (float) – Start and end of the trajectory [\(\text{eV}^{-1}\)].

• l1 (float) – Start and end of the trajectory [\(\text{eV}^{-1}\)].

• threshold (float, optional) – Adiabaticity threshold passed to magnus.adiabatic.find_nonadiabatic_windows(). Default: 0.1.

• n_probe (int, optional) – Density of the search grid for the same. Default: 200.

• fd_step_frac (float, optional) – Finite-difference step, as a fraction of the domain, for the same. Default: 1e-6.

• magnus_exp_order (int, optional) – Magnus order for the local patch. Default: 6.

• integration_method (str, optional) – Integration method for the local patch. Default: ‘gl’.

Returns:

The level-to-level probability matrix, with the starting level as the row index; the non-adiabatic windows found; and whether every local patch converged. A False in the last position means the crossing probabilities are not trustworthy, not that they are merely imprecise.

Return type:

(np.ndarray, list of (float, float), bool)

magnus.avgprob.averaged_probabilities_adiabatic(H_func: Callable, l0: float, l1: float, n_points: int | None = 201, threshold: float | None = 0.1, n_probe: int | None = 200, fd_step_frac: float | None = 1e-06, magnus_exp_order: int | None = 6, integration_method: str | None = 'gl') → Tuple[numpy.ndarray, dict]

Phase-averaged probabilities for a position-dependent Hamiltonian.

A neutrino produced at \(l_0\) decoheres in the eigenbasis there, is carried along the levels of the instantaneous Hamiltonian, and is detected in the eigenbasis at \(l_1\):

\[P_{\alpha\beta} = \sum_{ij} |V_{\alpha i}(l_0)|^2\, P^\text{cross}_{ij}\, |V_{\beta j}(l_1)|^2 ,\]

with \(P^\text{cross}\) from level_crossing_matrix() – the identity wherever the evolution is adiabatic. This is the standard MSW-plus-decoherence result, generalized to any number of levels and any number of crossings.

The windows come from a single search on n_probe points, which cannot see a front narrower than their spacing. So the profile is first checked for features that sharp and able to move probability (see SUDDEN_TRANSFER_THRESHOLD); where there is one, the windows are taken from magnus.adiabatic.hybrid_propagator() instead, which refines its search until it certifies, and which reports a profile it cannot resolve at all. Everywhere else the result is what it was, bit for bit.

Two things have to hold for the expression to mean anything, and both are checked rather than assumed. The levels must have decohered from each other by the time of detection, and if there is more than one crossing they must also have decohered between crossings, since otherwise composing the crossings as probabilities – rather than as amplitudes – discards interference that is still there. Both are reported.

Added in version 1.0.0.

Changed in version 1.1.1: Checks for features narrower than the probe spacing that could move probability, and takes the windows from magnus.adiabatic.hybrid_propagator() where it finds one (issue #60). The report gains 'escalated', 'resolved' and 'certified'.

Parameters:

• H_func (Callable) – Hamiltonian as a function of position, H_func(l) [eV].

• l0 (float) – Production and detection positions [\(\text{eV}^{-1}\)].

• l1 (float) – Production and detection positions [\(\text{eV}^{-1}\)].

• n_points (int, optional) – Sampling density for the accumulated-phase integrals. Default: 201.

• threshold (float, int, float, optional) – Passed to level_crossing_matrix(). Defaults: 0.1, 200 and 1e-6.

• n_probe (float, int, float, optional) – Passed to level_crossing_matrix(). Defaults: 0.1, 200 and 1e-6.

• fd_step_frac (float, int, float, optional) – Passed to level_crossing_matrix(). Defaults: 0.1, 200 and 1e-6.

• magnus_exp_order (int, optional) – Magnus order for the local patches. Default: 6.

• integration_method (str, optional) – Integration method for the local patches. Default: ‘gl’.

Returns:

The averaged probability matrix, rows summing to one, and a report with keys 'windows' (the non-adiabatic windows), 'patches_converged' (bool), 'undecided' (pairs that are in neither the coherent nor the decohered limit over the whole trajectory, as (i, j, phase) triples) and 'undecided_between_crossings' (every pair that has not decohered over an adiabatic stretch separating two crossings, coherent pairs included, since composing crossings as probabilities fails for those too; entries are (l_start, l_end, i, j, phase)). The returned matrix is always the fully decohered form, so these entries qualify a number that was computed regardless – unlike the constant-Hamiltonian route, which keeps coherent pairs coherent.

Three more keys say which search the windows came from. 'escalated' is True when the profile has a feature the n_probe grid cannot see and that could move probability. Then 'resolved' is whether magnus.adiabatic.hybrid_propagator() could resolve it – False means a discontinuity, the matrix is the unescalated one, and declaring the feature through t_breakpoints is the cure – and 'certified' is whether that refinement certified. Both are None when nothing escalated.

Return type:

(np.ndarray, dict)

magnus.avgprob.AVG_PHASE_SPREAD = 0.1

Module-level constant

Default relative energy spread \(\sigma\) of the phase average returned by phase_averaged_probabilities_constant_hamiltonian() and phase_averaged_probabilities_adiabatic(), and by average=True in magnus.oscprob.

The phase average keeps every interference term with its phase at the central energy and multiplies it by \(e^{-\sigma^2\phi'^2/2}\), where \(\phi' = d\phi/d\ln E\) is how fast that phase runs with energy. A term whose phase runs through many cycles across a spread \(\sigma\) is dropped, as the \(L/E \to \infty\) limit drops it; a term whose phase barely moves is kept with its real value; the ones in between are damped smoothly. Ten per cent is a typical resolution of neutrino detectors and telescopes. Mixing, crossing amplitudes and the eigenbases at the two ends of the path stay at the central energy: this averages phases, not probabilities, so a result without interference is returned unchanged.

Added in version 1.1.1.

Type:

float

magnus.avgprob.PHASE_AVERAGE_WINDOW_THRESHOLD = 0.01

Module-level constant

Adiabaticity threshold at which phase_averaged_probabilities_adiabatic() looks for non-adiabatic windows, lower than the 0.1 of averaged_probabilities_adiabatic().

Outside a window the evolution is carried as adiabatic, so whatever small transfer between levels happens there is lost; the decohered limit hides that loss, and the phase average does not, because it keeps the interference such a transfer carries. Measured against a brute-force average of the same definition on five solar chords from 10 GeV to 10 TeV: at 0.1 the error reaches 2.6e-03, at 0.03 and at 0.01 every chord is within 4.2e-05. On the two-level crossing of tests/test_phase_average.py 0.03 leaves 3.4e-04 and 0.01 leaves 5.1e-07, hence 0.01; the cost falls only on calls whose phases survive the spread, since the others never reach this search.

Added in version 1.1.1.

Type:

float

magnus.avgprob.PHASE_AVERAGE_PATCH_ATOL = 1e-05

Module-level constant

Tolerance on the elements of each window’s evolution operator in phase_averaged_probabilities_adiabatic(), looser than the 1e-7 of the Magnus patch the hybrid strategy uses. The slab count doubles until two successive operators agree to it, so the tolerance sets the cost of every window, at every node. Measured on three solar chords whose windows are long (10 GeV at \(b = 0.6\,R_\odot\), 100 GeV and 1 TeV at \(0.2\,R_\odot\), B16-GS98): at 1e-5 the probability moves by at most 1.4e-07 from its value at 1e-7, three orders below the 1e-4 at which magnus.oscprob returns the decohered limit instead, and the call is 5 to 16 times faster.

It is the default of phase_averaged_probabilities_adiabatic(), read at each call. The entry points of magnus.oscprob pass the tighter of their rtol and atol instead, 1e-3 by default (issue #65). With PHASE_AVERAGE_PHASE_TOL loosened with it, measured on twenty chords through the solar core (30 GeV to 3 TeV, \(b = 0.05\) to \(0.48\,R_\odot\)), the probability moves by at most 4.6e-06 from its value at 1e-5, and the call is 2.4 times faster at the median (0.9 to 6.3).

Added in version 1.1.1.

Type:

float

magnus.avgprob.PHASE_AVERAGE_PHASE_TOL = 1e-05

Module-level constant

Tolerance on the pair phases of each adiabatic stretch in phase_averaged_probabilities_adiabatic(), in rad. The stretch’s grid doubles until the pair phases that can still matter move by less than it, which moves a probability by at most as much. Like PHASE_AVERAGE_PATCH_ATOL, it is the default of that function, read at each call, and the entry points of magnus.oscprob pass the tighter of rtol and atol instead.

Added in version 1.1.1.

Type:

float

magnus.avgprob.PHASE_SPREAD_SENSITIVITY_THRESHOLD = 0.001

Module-level constant

Largest \(|\sigma\, \partial P / \partial\sigma|\) a phase-averaged probability may have before magnus.oscprob warns that it depends on the spread. That derivative is the change per e-fold of \(\sigma\); the threshold is the default tolerance of the package.

Added in version 1.1.1.

Type:

float

magnus.avgprob.phase_averaged_probabilities_constant_hamiltonian(hamiltonian: Sequence | numpy.ndarray, dH_dlnE: Sequence | numpy.ndarray, baseline: float | numpy.ndarray, spread: float | None = AVG_PHASE_SPREAD, dH_dlnE_step: float | None = None) → Tuple[numpy.ndarray, numpy.ndarray]

Phase-averaged probabilities for a constant Hamiltonian, from a flavor state at the start.

\[P_{\alpha\beta} = \sum_{ij} V^*_{\alpha i} V_{\beta i} V_{\alpha j} V^*_{\beta j}\, e^{-i\phi_{ij}}\, e^{-\sigma^2 \phi_{ij}'^2/2} , \qquad \phi_{ij} = (\lambda_i - \lambda_j) L ,\]

with \(\phi'_{ij} = d\phi_{ij}/d\ln E\) from the Hellmann-Feynman derivatives \(d\lambda_i/d\ln E = \langle v_i|\, dH/d\ln E\, |v_i\rangle\). At \(\sigma = 0\) this is the oscillation probability itself; for \(\sigma|\phi'| \gg 1\) on every pair it is the decohered sum \(\sum_i |V_{\alpha i}|^2 |V_{\beta i}|^2\). In vacuum \(\phi' = -\phi\), so a phase of \(2\pi\) keeps 82 per cent of its interference at \(\sigma = 10\%\), a phase of 30 rad about one per cent, and a phase of 40 rad \(3\times10^{-4}\). See AVG_PHASE_SPREAD.

Added in version 1.1.1.

Parameters:

• hamiltonian (list or np.ndarray) – Hermitian Hamiltonian [eV], shape (..., d, d); a leading batch axis is allowed.

• dH_dlnE (list or np.ndarray) – Its derivative with respect to \(\ln E\) [eV], same shape.

• baseline (float or np.ndarray) – Length of the path [\(\text{eV}^{-1}\)], broadcast against the batch axes.

• spread (float, optional) – Relative energy spread \(\sigma\). Default: AVG_PHASE_SPREAD.

• dH_dlnE_step (float, optional) – The step in \(\ln E\) of the finite difference dH_dlnE came from, if it came from one: it sets the round-off floor below which a pair’s slope is taken from its phase instead (see issue #61). None (default) means the derivative is exact.

Returns:

The probability matrix, shape (..., d, d), the initial flavor as the row index; and \(\max |\sigma\, \partial P/\partial\sigma|\) over its entries, shape (...).

Return type:

(np.ndarray, np.ndarray)

magnus.avgprob.phase_averaged_probabilities_adiabatic(H_func: Callable, dH_dlnE_func: Callable, l0: float, l1: float, spread: float | None = AVG_PHASE_SPREAD, threshold: float | None = PHASE_AVERAGE_WINDOW_THRESHOLD, windows: List[Tuple[float, float]] | None = None, n_probe: int | None = 200, fd_step_frac: float | None = 1e-06, magnus_exp_order: int | None = 6, integration_method: str | None = 'gl', dH_dlnE_step: float | None = None, patch_atol: float | None = None, phase_tol: float | None = None) → Tuple[numpy.ndarray, dict]

Phase-averaged probabilities on a smooth position-dependent Hamiltonian.

The neutrino starts decohered in the eigenbasis at \(l_0\), as in averaged_probabilities_adiabatic(), and is read out in the flavor basis at \(l_1\). In between, every interference term is kept with its phase and weighted by the spread of that phase across a relative energy spread \(\sigma\) (see AVG_PHASE_SPREAD): formally, the Gaussian average over \(u = \delta\ln E\) of the evolution under \(H + u\,D_\text{diag}\), with \(D_\text{diag}\) the part of \(dH/d\ln E\) diagonal in the instantaneous eigenbasis – an energy offset moves the eigenvalues and leaves the eigenvectors.

It is computed without sampling energies across the adiabatic stretches. Each non-adiabatic window is an amplitude matrix, evaluated at a few Gauss-Hermite nodes in \(u\) (a uniform grid when its internal phase runs too fast for that); each stretch between windows is a diagonal phase with an exact slope in \(u\). The density matrix is carried as terms labelled by accumulated slope, whose Gaussian average is analytic, and a term is dropped only once no later stretch can bring its slope back within reach. So the answer does not depend on where the windows are drawn: one window over a stretch or two windows with the stretch between them give the same number.

Where there is no window the evolution is adiabatic, a decohered start carries no interference, and the result is the decohered expression of averaged_probabilities_adiabatic().

Added in version 1.1.1.

Parameters:

• H_func (Callable) – Hamiltonian as a function of position, H_func(l) [eV]; arrays of positions are used where it accepts them.

• dH_dlnE_func (Callable) – Its derivative with respect to \(\ln E\), as a function of position [eV].

• l0 (float) – Production and detection positions [\(\text{eV}^{-1}\)].

• l1 (float) – Production and detection positions [\(\text{eV}^{-1}\)].

• spread (float, optional) – Relative energy spread \(\sigma\). Default: AVG_PHASE_SPREAD.

• threshold (float, optional) – Adiabaticity threshold of the window search. Default: PHASE_AVERAGE_WINDOW_THRESHOLD.

• windows (list of (float, float), optional) – Windows to use instead of searching. Default: None.

• n_probe (int, float, optional) – Passed to the window search. Defaults: 200 and 1e-6.

• fd_step_frac (int, float, optional) – Passed to the window search. Defaults: 200 and 1e-6.

• magnus_exp_order (int, optional) – Magnus order of the window patches. Default: 6.

• integration_method (str, optional) – Integration method of the window patches. Default: ‘gl’.

• dH_dlnE_step (float, optional) – See phase_averaged_probabilities_constant_hamiltonian(). Default: None.

• patch_atol (float, optional) – Tolerance on the elements of each window’s evolution operator. Default: None, which means PHASE_AVERAGE_PATCH_ATOL.

• phase_tol (float, optional) – Tolerance on the pair phases of each adiabatic stretch [rad]. Default: None, which means PHASE_AVERAGE_PHASE_TOL.

Returns:

The probability matrix, initial flavor as the row index; and a report with keys 'windows', 'escalated', 'resolved', 'certified' (as in averaged_probabilities_adiabatic()), 'patches_converged', 'phases_converged', 'n_nodes' and 'method' ('hermite', 'grid', or 'none' without windows), 'n_terms', and 'sigma_sensitivity', the largest \(|\sigma\, \partial P/\partial\sigma|\).

Return type:

(np.ndarray, dict)

Raises:

RuntimeError – If the number of terms would exceed an internal bound (_MAX_TERMS): many windows with many flavors whose phases never decohere.