magnus.adiabatic

adiabatic.py

Contains the adiabatic-transport-plus-Magnus-patch (“hybrid”) propagator used as an alternative strategy to compute the neutrino evolution operator when the Hamiltonian is position-dependent and develops an extreme accumulated phase (e.g., low-energy solar neutrinos crossing an MSW resonance), the regime in which the plain slab-refinement Magnus engine in magnus.oscprob needs a very large number of slabs (and may raise ToleranceNotAchievedWarning).

Physical idea: away from an eigenvalue crossing (or narrowly-avoided crossing) of the instantaneous Hamiltonian \(H(l)\), the adiabatic theorem says the evolution operator is well approximated by transport in the instantaneous eigenbasis of \(H(l)\) – a dynamical phase (the integral of the eigenvalues) plus a geometric (Berry) phase, both cheap to compute on a coarse grid regardless of how large the accumulated phase is. Near a genuine MSW resonance, the adiabatic approximation breaks down over a narrow window, which is patched with an exact, short-baseline Magnus computation (magnus.magnus.magnus_expansion_multislab(), the package’s own, already-unitary integrator). The two pieces are stitched together with the exact composition law of quantum evolution, \(U(l_2, l_0) = U(l_2, l_1) U(l_1, l_0)\), so the result is exactly unitary regardless of the approximation’s accuracy.

Where a patch is needed is decided by an exact Hellmann-Feynman diagnostic (no finite-differenced eigenvectors, which are gauge- ambiguous), so this applies to any Hermitian Hamiltonian of any dimension, with any number of simultaneous or sequential resonances – see Adiabatic + Magnus Hybrid Strategy for the full derivation, validation, and worked examples.

This module is self-contained: it depends only on magnus.magnus (the Magnus-expansion core), not on magnus.oscprob, so it can be used directly on any Hamiltonian function, independently of the rest of the oscillation-probability API. magnus.oscprob calls hybrid_propagator() internally when strategy='hybrid' or strategy='auto' (the default) is passed to magnus.oscprob.osc_prob_matter_std_potential(), magnus.oscprob.osc_prob_matter_nsi(), or magnus.oscprob.osc_prob_liv() (and, transitively, to every osc_prob_*_sun/osc_prob_*_sun_nsi/osc_prob_*_sun_liv wrapper), and also when it is passed to the fully generic user-Hamiltonian entry points magnus.oscprob.osc_prob_sun() and magnus.oscprob.osc_prob_earth() (via magnus.oscprob._osc_prob_with_potential). For osc_prob_earth the hybrid path is normally declined, since a real Earth trajectory supplies PREM layer breakpoints; see Adiabatic + Magnus Hybrid Strategy.

Routine listings

• adiabatic_propagator - Evolution operator via pure adiabatic

  (instantaneous-eigenbasis) transport, no resonance patching

• find_hidden_features - Detects structure too narrow for any grid

  this package lays down to sample

• find_resonance_candidates - Locates every exact eigenvalue-gap

  critical point of H(l) via the Hellmann-Feynman theorem

• find_nonadiabatic_windows - Filters/grows/merges candidates into

  position windows that need a Magnus patch

• hybrid_propagator - Adiabatic transport with Magnus patches at

  non-adiabatic windows, self-certified against successive refinement of every internal tolerance knob

Attributes

GAMMA_TO_ERROR	Module-level constant
RESOLUTION_RATIO	Module-level constant
LOCAL_JUMP_RATIO	Module-level constant
N_LOCAL_CONFIRM	Module-level constant
MAX_LOCAL_CONFIRMATIONS	Module-level constant
HIDDEN_FEATURE_CONCENTRATION	Module-level constant
N_HIDDEN_FEATURE_SUBDIVISION	Module-level constant
THRESHOLD0_PROVENANCE	Module-level constant

Functions

find_hidden_features(→ Dict)	Looks for structure too narrow for any grid this package lays down to sample.
oscillation_sampling(→ Dict)	How finely does a scan sample the fastest oscillation on its trajectory?
adiabatic_propagator(→ numpy.ndarray)	Computes the evolution operator via pure adiabatic (instantaneous-eigenbasis) transport.
find_resonance_candidates(→ List[Dict])	Locates every exact eigenvalue-gap critical point of H_func between l0 and l1.
find_nonadiabatic_windows(→ Tuple[List[Tuple[float, ...)	Finds every position window along [l0, l1] where H_func needs a Magnus patch.
hybrid_propagator(→ Tuple[numpy.ndarray, ...)	Computes the evolution operator via adiabatic transport, with a Magnus patch at every

Module Contents

magnus.adiabatic.GAMMA_TO_ERROR = 0.85

Module-level constant

Probability error to budget per unit of the adiabaticity parameter \(\gamma\), when hybrid_propagator() decides whether a result with no non-adiabatic window may be certified. The rule is GAMMA_TO_ERROR * gamma_max <= atol + rtol.

Write the pure adiabatic answer’s error as \(|\Delta P| = k\,\gamma_\max\). Certifying an empty window list is safe exactly when \(k\,\gamma_\max \le\) tolerance, so this constant is an upper bound on \(k\) – and the only honest way to set it is to measure \(k\).

Measured over 149 configurations (docs/dev/adversarial_batteries/gamma_slack_sweep.py): resonance width swept over a decade, d = 2-5, 5-80 MeV, 0.5-2 density scale heights, scoring the pure adiabatic operator against solve_ivp/DOP853. Only the small-\(\gamma\) rows matter, because this rule governs the no-window case alone: once \(\gamma_\max\) exceeds the active threshold a window opens and the answer is patched exactly, whatever \(k\) is there.

population	n	max k	implied bound on gamma_max
gamma_max < 1e-2	76	0.812	<= 1.23 x tolerance
gamma_max < 3e-3	35	0.679	<= 1.47 x tolerance
gamma_max < 1e-3	12	0.502	<= 1.99 x tolerance
all rows, incl. patched	149	1.136	<= 0.88 x tolerance

\(k\) falls towards ~0.5 as \(\gamma\) shrinks, which is what the linear model predicts asymptotically. 0.85 covers the worst case in the governed regime (0.812) with a little margin.

This constant was wrong twice, in opposite directions. It began at 1.0 alongside a slack factor of 2.0, derived from five points that all happened to sit at \(\gamma_\max < 10^{-3}\) where \(k \approx 0.5\); that pair encoded \(k \le 0.5\), right for those five and optimistic by up to 1.6x elsewhere in the governed regime. Reading the unrestricted maximum (1.136) then argued for a far stricter bound – an over-correction, since those rows sit at \(\gamma_\max \sim 0.2\), open a window immediately, and are never decided by this rule. The slack factor is gone: it existed only to compensate for the mis-measured value, and with \(k\) measured in the regime that matters no fudge is needed.

Added in version 1.0.0.

Type:

float

magnus.adiabatic.RESOLUTION_RATIO = 0.7

Module-level constant

Threshold of the probe-scale resolution test in _profile_is_resolved, which decides whether H_func is sampled finely enough for this module’s finite-difference diagnostics to mean anything.

Within each probe interval, the test asks what fraction of the variation falls in one half. For a \(C^1\) Hamiltonian the two halves each carry about half, so the ratio tends to 0.5; a jump lands entirely inside one half, so the ratio tends to 1.0.

The two limits are 0.5 and 1.0, but the honest threshold is set by measurement rather than by their midpoint, because a jump comparable in size to the local smooth variation lands between them: a jump \(J\) on top of a smooth change \(S\) within one interval gives \((J + S/2)/(J + S)\). Measured:

population	statistic
192 smooth configurations (6 profile families plus
10 random Fourier sums, d = 2-5, 5-200 MeV)	<= 0.602
a jump 4.7x smaller than the steepest smooth step	0.773
15 random piecewise-constant profiles, d = 2 and 3	1.000

0.70 sits in the gap with margin on both sides – 16% above the smooth ceiling, 10% below the weakest genuine discontinuity. It also states what the test can and cannot catch: solving \((J + S/2)/(J + S) > t\) gives \(J/S > (2t-1)/(1-t)\), so at 0.70 a jump must be at least 1.33x the local smooth variation to be seen. A smaller one is genuinely indistinguishable from steep smooth behaviour at that sampling density.

Two earlier formulations were wrong, both caught by measurement. Comparing the global largest adjacent change at two grid densities masks any jump smaller than the largest smooth variation elsewhere on the path: a discontinuity 4.7x smaller than the steepest smooth step went undetected, and the answer came back wrong by 2.0e-02, silently. Comparing one half against the whole interval fixes that but false-positives at a smooth turning point, where the interval’s net change is near zero while each half is not. Comparing each half against the sum of the two is immune to both: it measures concentration of variation, which is what a jump is.

Added in version 1.0.0.

Type:

float

magnus.adiabatic.LOCAL_JUMP_RATIO = 0.5

Module-level constant

Threshold of the local confirmation in _profile_is_resolved: having flagged a probe interval whose variation is concentrated in one half, re-sample that interval alone on N_LOCAL_CONFIRM points and ask what fraction of its variation still falls in a single adjacent step.

The two limits are far apart and are set by arithmetic, not by taste. A jump is not diluted by refinement – one step still carries all of it – so the fraction tends to 1.0. A smooth feature is spread over the whole interval, so with 32 sub-steps the largest one carries roughly \(1/32\) of the variation. Measured (docs/dev/adversarial_batteries/resolution_fp.py), over the flagged intervals only, which is the population this constant decides:

population	n	statistic
10 smooth families x d = 2-5 x 3 energies x 12
sub-intervals (1440 configurations)	79	<= 0.087
3 piecewise-constant families, same sweep	348	>= 1.000

0.5 sits between them with a factor of 5.7 of margin on the smooth side and 2.0 on the jump side. Over the same sweep the completed test reports 0 of 1440 smooth configurations as unresolved, and every sub-interval that genuinely contains a jump as unresolved.

Why this exists. Without it the concentration test is a maximum over intervals, and one interval decides – so the interval containing a smooth turning point decides. There the two halves are genuinely asymmetric: one nearly cancels while the other does not, and the ratio is a draw in \([0.5, 1]\) that depends on where the extremum happens to fall inside its interval. Refining the grid does not remove it; it re-draws it, which is why the two-stage coarse-then-fine protocol could not separate the two either. Measured on a Gaussian bump of width \(10^{-2}\) of the trajectory – a profile the module answers to 3.5e-09 – the statistic hit 0.75 at the interval containing the peak, and 6 of 30 baselines of one ordinary scan were declared discontinuous. The local confirmation is the discriminator the concentration ratio alone does not have: it asks whether the concentration survives refinement, which is the one thing a jump does and a turning point does not.

Costs nothing on an ordinary call: it runs only on intervals the cheap test already flagged, and on a profile that flags none it is never reached.

Added in version 1.0.0.

Type:

float

magnus.adiabatic.N_LOCAL_CONFIRM = 33

Module-level constant

Points used to re-sample one flagged probe interval in _profile_is_resolved (see LOCAL_JUMP_RATIO). 32 sub-steps puts the smooth limit at \(1/32 \approx 0.03\), comfortably below LOCAL_JUMP_RATIO, and resolves a feature down to \(1/(199 \times 32) \approx 1.6\times10^{-4}\) of the trajectory – the same scale as the probe refinement ceiling max_n_probe = 6400, so the local test does not claim to see anything the caller’s own refinement could not.

Added in version 1.0.0.

Type:

int

magnus.adiabatic.MAX_LOCAL_CONFIRMATIONS = 8

Module-level constant

At most this many flagged intervals are confirmed locally, taken in decreasing order of how much variation they carry. A profile with hundreds of genuine jumps is found by the first one, and a bound is needed because the flagged set is unbounded in principle.

Added in version 1.0.0.

Type:

int

magnus.adiabatic.HIDDEN_FEATURE_CONCENTRATION = 0.3

Module-level constant

Threshold of find_hidden_features(), which looks for structure below the scale any grid in this package samples on. This is the one exposure the adversarial validation could not close: a Gaussian narrower than the probe spacing is invisible to the hybrid detector, to the general ladder’s slab grid and to the cumulative scan alike, so all three agree and all three are wrong – silently, by up to 2.9e-02 against a requested 1e-3.

The statistic is concentration, not size, and that choice is the whole design. Within each interval of a reference grid, compare the total variation a much denser grid sees inside it with the change its two endpoints show; the excess is variation hidden between reference samples. Report the largest such excess as a fraction of the profile’s total variation.

A first attempt used the cruder TV_dense/TV_reference and was wrong: a sinusoid at exactly the probe spacing sends the denominator to zero and that ratio to \(10^{13}\), while the package answers such a profile to ~1e-11. What separates the two is where the hidden variation sits. An aliased sinusoid hides some in every interval, so its share of the total is ~1/n_ref; one narrow bump hides all of it in one, so its share is ~1. The excess is summed over adjacent pairs of reference intervals, so that a feature landing on an interval boundary is not halved by the split.

Measured over 67 profiles the package serves – solar, multi-resonance, noisy, sinusoids at 1x, 2x and 1/2x the probe spacing, 400 crossings, a declared step, 30 random Fourier sums, and 30 Gaussian bumps of random width down to 1e-4 of the trajectory at random positions:

population	concentration
67 smooth/resolvable profiles	max 0.060 p99 0.047

and, over 60 random positions each, the detection rate for features in the band no grid here resolves:

feature width	detection at 0.2	detection at 0.3	detection at 0.5
3e-5	0.70	0.68	0.55
1e-5	0.90	0.90	0.90
3e-6	0.82	0.82	0.82
1e-6	0.73	0.73	0.73

0.3 gives zero false positives over all 67 smooth profiles – five times the measured ceiling – at the best detection the margin allows; 0.2 buys two points of detection for half the margin, and 0.5 costs thirteen.

The caller now varies the sampling density with the size of the request (see N_HIDDEN_FEATURE_SUBDIVISION), so the ceiling was re-measured at every density the dispatcher can choose rather than only at the default: 0.0597 at 8 sub-steps, 0.0601 at 16, 0.0602 at 32 – 0 of 67 in all three. The statistic is a fraction of the total variation, which is why it barely moves: refining the dense grid adds the same variation to numerator and denominator.

This detects most of the class, not all of it, and the shortfall is structural. A feature of width 3e-5 is right at the edge of what max_n_probe = 6400 can partially resolve, so its variation is partly visible to the reference grid and the statistic is diluted; a feature of width 1e-6 is far below the dense spacing (\((l_1-l_0)/51192 \approx 2\times10^{-5}\)), so whether a sample lands inside it is luck. Between those two limits detection is ~0.8-0.9. Against a prior state of zero detection and a silent 2.9e-02 error, that is the improvement on offer; it is not a guarantee, and the docstring says so rather than implying one.

Why the reference grid is the ceiling and not the starting density. Against n_probe0 = 200 the same statistic flags widths of 1e-3 and 1e-4 too – correctly, in that they are hidden at that density, and uselessly, in that the refinement resolves them. Keying it to max_n_probe asks the question that matters: is there structure left that no amount of refinement will reach?

Added in version 1.0.0.

Type:

float

magnus.adiabatic.N_HIDDEN_FEATURE_SUBDIVISION = 8

Module-level constant

Sub-steps per reference interval in find_hidden_features(), so 51 192 samples of the profile in total. Chosen on cost, because this runs on ordinary calls: the statistic is nearly independent of it (the separation is 4.7x at 2 sub-steps and 4.6x at 32), while the cost is not.

sub-steps	dense samples	scan cost
4	25 597	0.22 ms
8	51 193	0.37 ms
16	102 385	1.30 ms
32	204 769	2.85 ms

The jump past 8 is superlinear – the arrays stop fitting in cache – and 0.37 ms is about 3% of an ordinary 13 ms single-point call, where 2.85 ms would be 20% and fail the package’s own 10% performance criterion. What it costs in reach is the very narrowest features: the dense spacing is \((l_1-l_0)/51192\), and detection of anything below that is a matter of whether a sample lands inside it (measured 0.73 at a width of 1e-6, against 0.90 at 1e-5).

magnus.oscprob scales this with the number of requested points – 8 below four points, 16 below sixteen, 32 above – because the scan runs once per call whatever the point count, so its share of the work falls as the request grows. That holds it under about 7 % of the call at every size instead of spending 20 % of the cheapest one, and the false-positive rate was re-measured at each of those three densities (0 of 67 every time; see HIDDEN_FEATURE_CONCENTRATION).

Raise it if you have reason to think the profile hides something finer; the scan is find_hidden_features() and takes n_sub directly.

Added in version 1.0.0.

Type:

int

magnus.adiabatic.find_hidden_features(profile: Callable, l0: float, l1: float, n_ref: int | None = 6400, n_sub: int | None = None) → Dict

Looks for structure too narrow for any grid this package lays down to sample.

See HIDDEN_FEATURE_CONCENTRATION for the statistic, why it is a concentration rather than a size, and the measured separation.

Pass the scalar potential when there is one. For the separable Hamiltonians this package builds, \(H(l) = h_\mathrm{vac}/E + V_{CC}(l)\,P_{ee}\) is affine in \(V_{CC}\), so every difference is \(|\Delta V_{CC}|\) times a constant and the statistic is identical either way – verified bit-for-bit at d = 2, 3 and 5. It is also 18x cheaper (2.6 ms against 48 ms), because sampling a scalar avoids allocating a stack of \(2\times10^5\) matrices. This is a free choice, not an approximation.

Cost is a single vectorized evaluation of profile and two array passes; no eigendecomposition. It depends on the profile and the interval but not on energy, so a caller scanning many energies or baselines should run it once, not once per point.

Added in version 1.0.0.

Parameters:

• profile (Callable) – The position-dependent part of the problem: either the scalar potential (preferred) or H_func itself. Must accept an array of positions.

• l0 (float) – Interval to scan.

• l1 (float) – Interval to scan.

• n_ref (int, optional) – Reference grid, which should be the finest grid the caller’s machinery can reach – the hybrid strategy’s max_n_probe. Default: 6400.

• n_sub (int, optional) – Sub-steps per reference interval. Default: N_HIDDEN_FEATURE_SUBDIVISION.

Returns:

'concentration' (the statistic), 'hidden' (whether it exceeds HIDDEN_FEATURE_CONCENTRATION), 'l_lo'/'l_hi' (the reference interval carrying the excess) and 'l_centre'. On a constant profile, concentration 0.0.

Return type:

dict

magnus.adiabatic.oscillation_sampling(H_func: Callable, l0: float, l1: float, baselines: numpy.ndarray | None = None, n_probe: int | None = 8) → Dict

How finely does a scan sample the fastest oscillation on its trajectory?

A long trajectory carries an enormous oscillation phase, and a scan over it usually samples that phase far too coarsely to represent it – so the returned array is a set of correct values that must not be read as a curve. For a solar or supernova problem the observable is normally the phase-averaged probability anyway; see magnus.avgprob, whose magnus.avgprob.coherence_report() decides which pairs of eigenvalues have decohered, and which are in neither limit.

This reports; it never warns. A Nyquist criterion is objectively correct and fires on essentially everything: measured over the physical profile families, a scan would need 4400 points on a solar trajectory and 73 000 on a supernova ray to sample the fastest oscillation twice per cycle, so 44 of 45 realistic scan sizes are formally aliased (docs/dev/adversarial_batteries/alias_fp.py). A warning that fires on 98 % of calls is noise however correct each firing is, so the information is offered here and in strategy_info instead, and the caller decides.

Cost: n_probe evaluations of H_func and one eigvalsh each. Measured against the scan it describes: 5.5 % of the cheapest scan in the physical population (a 13 ms Earth chord), 1.8 % of a 21 ms one, and under 0.1 % of anything substantial (docs/dev/adversarial_batteries/alias_cost.py). Callers who do not ask for strategy_info pay nothing – oscprob only runs this when a report was requested. Eight samples are enough: the fastest oscillation comes from the largest eigenvalue spread, which tracks the matter potential and varies smoothly, and 8 samples agree with 4096 to 2 % over every family measured. (The smallest gap would need far denser sampling, because it has a sharp minimum at an MSW resonance – but that sets the slowest oscillation, not the fastest, and is not what aliases.)

Added in version 1.0.0.

Parameters:

• H_func (Callable) – The Hamiltonian as a function of position, H(l), at the energy of interest [eV].

• l0 (float) – Interval the scan traverses.

• l1 (float) – Interval the scan traverses.

• baselines (np.ndarray, optional) – The requested baselines. When given, the spacing and the aliasing verdict are computed from them; when omitted, only the trajectory-level quantities are returned.

• n_probe (int, optional) – Samples along the trajectory. Default: 8.

Returns:

'oscillation_length' (shortest on the trajectory [eV^-1]), 'cycles_over_trajectory', 'nyquist_points' (points a scan would need to sample the fastest oscillation twice per cycle), and, when baselines is given, 'spacing', 'cycles_per_step' and 'aliased'. Empty dict if the spectrum is degenerate or the interval has zero length.

Return type:

dict

Examples

import numpy as np

import magnus.globaldefs as gd
from magnus import adiabatic
from magnus.hamiltonians import hamiltonians3nu

p = gd.OSC_PARAMS_PREDEFINED['OSC_PARAMS_DEFAULT']
h_vac = np.asarray(hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
    p['s12'], p['s23'], p['s13'], p['dCP'], p['D21'], p['D31']))
e00 = np.diag([1.0, 0.0, 0.0])
energy = 10.0e6

def H(l):
    v = 1.0e-11*np.exp(-np.asarray(l, dtype=float)/gd.L_SCALE_SUN)
    return (1.0/energy)*h_vac + np.asarray(v)[..., None, None]*e00

print(adiabatic.oscillation_sampling(H, 0.0, 3.0*gd.L_SCALE_SUN, 100))

{'oscillation_length': 50129793295.43478, 'cycles_over_trajectory': 19989.20386274845, 'nyquist_points': 39980}

Reported, never warned about. A trajectory of twenty thousand cycles sampled at a hundred points is not necessarily wrong – an averaged observable may not care – but it is worth knowing before trusting an instantaneous one.

magnus.adiabatic.THRESHOLD0_PROVENANCE = 0.1

Module-level constant

The value hybrid_propagator() starts its adiabaticity-threshold ladder at, recorded here with its measurement because the measurement is what stopped it being changed.

What was measured. Over 3 profiles x d = 2, 3 x three requested tolerances (docs/dev/adversarial_batteries/constants_audit.py), sweeping threshold0 from 1 down to 1e-3 at a fixed baseline: accuracy is identical at every value in 16 of 18 rows, and at rtol <= 1e-3 a lower start is up to 6.5x cheaper (1.57 s to 0.24 s on a solar profile at rtol = 1e-5). That is because certifying an empty window list additionally requires \(\gamma\) to fit the tolerance (see GAMMA_TO_ERROR), so the ladder reaches whatever threshold the request needs regardless of where it starts; the start only decides how many iterations that takes, and each one re-runs the detector at doubled n_probe and the transport at doubled n_points. At rtol = 1e-2 the sign flips – there no window is needed at all, so starting low opens one that is not, and 0.21 s becomes 0.95 s at d = 3.

So the right value looked like a rule rather than a constant, and one was built: start at \((\mathrm{atol} + \mathrm{rtol})/\texttt{GAMMA\_TO\_ERROR}\), the exact \(\gamma\) at which an empty window list stops being certifiable, clipped to [min_threshold, 0.1]. At the default tolerance that is 2.4e-03 rather than 0.1.

It was then rejected, on evidence the sweep above could not produce. Scored against solve_ivp on the package’s bit-identity workloads, which include an energy scan the fixed-baseline sweep did not:

workload	t0=0.1	the rule	verdict
single point, solar	1.624e-06	1.184e-10	13711x better
sub-threshold scan, N = 8	3.220e-05	3.814e-05	1.2x worse
energy scan, fixed baseline	2.509e-05	4.954e-04	20x worse

All three stay inside the requested 1e-3, but 4.95e-04 spends half the budget where 2.5e-05 spent a fortieth. The mechanism is visible once looked for: starting low opens a window on the first iteration, and windows_next or windows_prev then short-circuits the \(\gamma\) check, so agreement can be accepted at a coarser transport grid than the old start would have forced. The saving and the loss are the same effect seen from two sides.

The lesson is the one this package keeps re-learning: the sweep that justified the rule ran at a fixed baseline, and the row that refuted it was an energy scan. A population that does not contain the workload you are about to change is not evidence about it – the same shape of mistake that made GAMMA_TO_ERROR wrong twice. 0.1 stays until a population that spans scans as well as points says otherwise.

Added in version 1.0.0.

Type:

float

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

Computes the evolution operator via pure adiabatic (instantaneous-eigenbasis) transport.

Diagonalizes H_func on a grid of n_points positions between l0 and l1, integrates each eigenvalue’s dynamical phase with Simpson’s rule (trapezoidal quadrature leaves a spurious residual that can look like a physics limit but is pure quadrature error), and reassembles the evolution operator in the original (flavor) basis:

\[U(l_1, l_0) \approx W(l_1)\, \mathrm{diag}\!\left(e^{-i\Phi_k}\right)\, W(l_0)^\dagger , \qquad \Phi_k = \int_{l_0}^{l_1} \lambda_k(l)\, dl ,\]

with \(W(l)\) the (parallel-transported; see _eigs_along) matrix of instantaneous eigenvectors of \(H(l)\) and \(\lambda_k(l)\) its eigenvalues. This is exact in the strict adiabatic limit (no eigenvalue crossing or narrowly-avoided crossing along the trajectory) and remains unitary by construction (a diagonal phase conjugated by unitary matrices) regardless of grid density – the only thing n_points controls is how well the quadrature/parallel-transport approximate the continuum limit, not whether the result is unitary. See hybrid_propagator() for what to do when the trajectory does cross a resonance.

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian, a function of position returning a square (Hermitian) matrix. May be real- or complex-valued.

• l0 (float) – Initial position.

• l1 (float) – Final position.

• n_points (int, optional) – Number of positions at which to diagonalize H_func between l0 and l1. Default: 201.

Returns:

The evolution operator, exactly unitary.

Return type:

np.ndarray

Examples

import numpy as np

import magnus.globaldefs as gd
from magnus import adiabatic
from magnus.hamiltonians import hamiltonians3nu

p = gd.OSC_PARAMS_PREDEFINED['OSC_PARAMS_DEFAULT']
h_vac = np.asarray(hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
    p['s12'], p['s23'], p['s13'], p['dCP'], p['D21'], p['D31']))
e00 = np.diag([1.0, 0.0, 0.0])
energy = 10.0e6

def H(l):
    v = 1.0e-11*np.exp(-np.asarray(l, dtype=float)/gd.L_SCALE_SUN)
    return (1.0/energy)*h_vac + np.asarray(v)[..., None, None]*e00

U = np.asarray(adiabatic.adiabatic_propagator(H, 0.0, 3.0*gd.L_SCALE_SUN))

print('shape', U.shape)
print('unitary to %.1e' % np.max(np.abs(U.conj().T @ U - np.eye(3))))

shape (3, 3)
unitary to 6.7e-16

Cheap however large the accumulated phase, because it transports in the instantaneous eigenbasis rather than resolving the oscillation. That is exactly why it needs a patch wherever the transport stops being adiabatic.

magnus.adiabatic.find_resonance_candidates(H_func: Callable, l0: float, l1: float, n_probe: int | None = 200, fd_step_frac: float | None = 1e-06, info: Dict | None = None) → List[Dict]

Locates every exact eigenvalue-gap critical point of H_func between l0 and l1.

For every pair of levels \((j, k)\), scans for sign changes of

\[f_{jk}(l) = \langle v_j(l)|\, dH/dl\, |v_j(l)\rangle - \langle v_k(l)|\, dH/dl\, |v_k(l)\rangle ,\]

refined by bisection to machine precision in position. By the Hellmann-Feynman theorem, \(d\lambda_k/dl = \langle v_k|\, dH/dl\, |v_k\rangle\) exactly (no eigenvector finite difference, which would be gauge-ambiguous and fragile), so a sign change of \(f_{jk}\) is an exact critical point of the gap \(\lambda_j - \lambda_k\) – a genuine crossing or near-crossing candidate, for any Hermitian H_func of any dimension, with no assumption of a separable or otherwise special structure. Every pair is scanned, so any number of simultaneous or sequential resonances (between any pair of levels) are all found.

A returned candidate is a structural fact about H_func (an extremum of that pair’s gap); whether it is actually non-adiabatic (whether it needs a Magnus patch) is a separate question, answered by find_nonadiabatic_windows().

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian, a function of position returning a square (Hermitian) matrix.

• l0 (float) – Initial position.

• l1 (float) – Final position.

• n_probe (int, optional) – Number of positions on the initial scan grid used to bracket sign changes. Default: 200.

• fd_step_frac (float, optional) –

  Finite-difference step for _dH_dl, as a fraction of l1 - l0. Default: 1e-6.

  Provenance. A central difference trades truncation error (\(\sim h^2\)) against subtractive cancellation (\(\sim \epsilon/h\)), so there is an optimum, and this default had never been checked against it. Measured (docs/dev/adversarial_batteries/constants_audit.py) against the analytic \(dH/dl\) – available because \(H = h_\mathrm{vac}/E + C\,n_e(l)\,P_{ee}\) and the profiles used have closed-form \(n_e'\), so the reference carries no error of its own – over d = 2, 3 and 10-50 MeV, as a fraction of the largest \(|dH/dl|\) on the path:

  profile	optimum	error at 1e-6	error at the optimum
  solar exponential	1e-5	6.8e-11	1.9e-11
  sinusoid, period span/7	1e-6	3.4e-10	3.4e-10
  Gaussian bump, w = 1e-2	1e-7 - 1e-8	2.5e-09	2.6e-11

  The optimum is not a single number – it moves with the profile’s shortest length scale, which is what the theory predicts. But the curve is flat enough that this does not matter: anywhere in \(10^{-8}\) to \(10^{-5}\) the error stays below \(3\times10^{-9}\) relative on every profile measured, which is six orders below anything that could move a probability at the tolerances this package works to. Outside that band it degrades fast in both directions – 1.6e-04 at \(10^{-12}\) (cancellation), 0.23 at \(10^{-2}\) (truncation) – so the band, not the value, is the thing to preserve.

• info (dict, optional) – If given, filled in place with the probe-grid quantities this function had to compute anyway – 'ls' (the grid), 'lam', 'W' (eigenvalues and eigenvectors, shapes (n, d) and (n, d, d)) and 'dH'. find_nonadiabatic_windows() sweeps \(\gamma\) on exactly this grid with exactly this finite-difference step, so without this it would recompute all of it: ~600 extra Hamiltonian evaluations and a second eigendecomposition, which measured as 1.4x on an ordinary single-point solar call. Default: None.

Returns:

One entry per candidate, with keys 'l' (position), 'j', 'k' (the level indices, j < k), and 'gap' (\(\lambda_k - \lambda_j\) at that position), sorted by position.

Return type:

list of dict

Examples

import numpy as np

import magnus.globaldefs as gd
from magnus import adiabatic
from magnus.hamiltonians import hamiltonians3nu

p = gd.OSC_PARAMS_PREDEFINED['OSC_PARAMS_DEFAULT']
h_vac = np.asarray(hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
    p['s12'], p['s23'], p['s13'], p['dCP'], p['D21'], p['D31']))
e00 = np.diag([1.0, 0.0, 0.0])
energy = 10.0e6

def H(l):
    v = 1.0e-11*np.exp(-np.asarray(l, dtype=float)/gd.L_SCALE_SUN)
    return (1.0/energy)*h_vac + np.asarray(v)[..., None, None]*e00

for c in adiabatic.find_resonance_candidates(H, 0.0, 3.0*gd.L_SCALE_SUN):
    print('levels %d-%d cross at l/l_scale = %.3f, gap %.2e eV'
          % (c['j'], c['k'], c['l']/gd.L_SCALE_SUN, c['gap']))

levels 0-1 cross at l/l_scale = 1.917, gap 3.46e-12 eV

A candidate is a critical point of a pairwise gap, found exactly through Hellmann-Feynman rather than by scanning for a minimum. Whether it is actually non-adiabatic is a separate question – find_nonadiabatic_windows() answers it.

magnus.adiabatic.find_nonadiabatic_windows(H_func: Callable, l0: float, l1: float, threshold: float | None = 0.1, n_probe: int | None = 200, fd_step_frac: float | None = 1e-06, info: Dict | None = None) → Tuple[List[Tuple[float, float]], List[Dict]]

Finds every position window along [l0, l1] where H_func needs a Magnus patch.

Calls find_resonance_candidates(), evaluates the adiabaticity parameter \(\gamma_{jk}\) (see _point_adiabaticity) at each candidate, grows a window around every candidate with \(\gamma_{jk} > \text{threshold}\) (see _estimate_window_bounds), and merges any windows that overlap or touch – so two (or more) resonances close enough together are correctly folded into a single patch, rather than either double-counted or (worse) silently dropped. This works for any number of simultaneous or sequential resonances, between any pair of levels.

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian, a function of position returning a square (Hermitian) matrix.

• l0 (float) – Initial position.

• l1 (float) – Final position.

• threshold (float, optional) – Adiabaticity parameter above which a candidate is treated as non-adiabatic. Default: 0.1.

• n_probe (int, optional) – Forwarded to find_resonance_candidates(). Default: 200.

• fd_step_frac (float, optional) – Forwarded to find_resonance_candidates(). Default: 1e-6.

• info (dict, optional) – If given, filled in place with diagnostics about this call, following the same out-parameter convention as convergence_info in magnus.oscprob.osc_prob(). Currently one key, 'gamma_max': the largest adiabaticity parameter seen anywhere on the probe grid or at any candidate, over every level pair. It is inf if some pair’s gap vanishes exactly. hybrid_propagator() uses it to decide whether an empty window list may be certified – without it, “no window opened” is indistinguishable from “no window was looked for hard enough”. Default: None.

Returns:

The merged, non-overlapping windows (each a (l_b, l_c) pair, sorted by position), and the candidate list from find_resonance_candidates(), each entry additionally carrying its evaluated 'gamma'.

Return type:

(list of (float, float), list of dict)

magnus.adiabatic.hybrid_propagator(H_func: Callable, l0: float, l1: float, rtol: float | None = 0.001, atol: float | None = 0.001, magnus_exp_order: int | None = 6, integration_method: str | None = 'gl', threshold0: float | None = 0.1, min_threshold: float | None = 1e-06, n_probe0: int | None = 200, max_n_probe: int | None = 6400, n_points0: int | None = 201, max_n_points: int | None = 12864, fd_step_frac: float | None = 1e-06, max_iters: int | None = 12, info: Dict | None = None) → Tuple[numpy.ndarray, List[Tuple[float, float]], bool]

Computes the evolution operator via adiabatic transport, with a Magnus patch at every non-adiabatic window, self-certified against successive refinement of every internal tolerance knob.

This is the main entry point of this module (see Adiabatic + Magnus Hybrid Strategy for the full derivation and validation). Given any Hermitian H_func of any dimension:

1. Locates every non-adiabatic window along [l0, l1] (see find_nonadiabatic_windows()).

2. If there are none, returns the pure adiabatic-transport operator (see adiabatic_propagator()).

3. Otherwise, composes adiabatic transport between windows with an exact local Magnus patch inside each window (see _local_evolution_operator), using the exact composition law of quantum evolution, \(U(l_2, l_0) = U(l_2, l_1)\, U(l_1, l_0)\), so the result is exactly unitary regardless of any approximation’s accuracy.

4. Self-certifies the result: a single fixed adiabaticity threshold (deciding which candidates count as non-adiabatic) is not safe in general – too loose, and a genuine resonance is patched too narrowly or missed; too tight, and windows are patched needlessly, at some (still usually modest) extra cost. Rather than trust one fixed value, the whole computation (window threshold, adiabatic-transport grid density, and the probe grid used to locate candidates) is repeated with the knobs tightened together (threshold divided by 3, n_points/n_probe doubled) until two successive results agree within rtol/atol, mirroring the successive-refinement discipline magnus.oscprob.osc_prob() already uses for the number of slabs.

   Each knob stops at its own ceiling (min_threshold, max_n_probe, max_n_points), which they reach at different iterations, so the later iterations tighten fewer knobs than the earlier ones. Once all of them have saturated, a further iteration would recompute bit-identical inputs and the agreement test would pass trivially, comparing a result with itself; the loop therefore stops at that point and reports certified=False rather than certifying on the strength of a comparison that carries no information.

Added in version 1.0.0.

Parameters:

• H_func (Callable) – Hamiltonian, a function of position returning a square (Hermitian) matrix. May be real- or complex-valued, of any dimension.

• l0 (float) – Initial position.

• l1 (float) – Final position.

• rtol (float, optional) – Relative tolerance on the agreement between successive refinement levels, and on the adiabaticity bound in _certified. Default: 1e-3. Like every tolerance in this package it is a stopping rule rather than a guaranteed accuracy: the loop halts when two successive levels agree, and no error of the returned operator is ever estimated. See the rtol entry of magnus.oscprob.osc_prob() for what that does and does not promise. certified is the flag that says whether the loop stopped because it agreed or because it ran out of room.

• atol (float, optional) – Absolute tolerance on the same agreement; see rtol. Default: 1e-3.

• magnus_exp_order (int, optional) – Magnus expansion order used for the local patch inside each non-adiabatic window. Default: 6.

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

• threshold0 (float, optional) –

  Starting adiabaticity threshold. Default: 0.1. See THRESHOLD0_PROVENANCE for what it was measured to do, and for why a tolerance-derived rule was built, tested and then rejected.

  Provenance, and a reframing. This constant used to decide correctness: if \(\gamma\) never crossed it, no window opened, successive refinements agreed with each other, and the answer was certified while wrong by up to 1.8e-02. Since certification of an empty window list additionally requires GAMMA_TO_ERROR * \(\gamma_\max\) to fit the tolerance, that failure mode is closed and this is now a cost knob: it sets where the threshold /= 3 ladder starts, not whether the result may be believed.

  Measured (docs/dev/adversarial_batteries/constants_audit.py) over 3 profiles x d = 2, 3 x three requested tolerances, sweeping threshold0 from 1 down to 1e-3:

  • Accuracy is identical at every value in 16 of 18 rows. The ladder reaches whatever threshold the tolerance requires regardless of where it starts.

  • At ``rtol <= 1e-3``, lower is monotonically cheaper – up to 6.5x (1.57 s to 0.24 s on a solar profile at rtol = 1e-5) – because every step of the ladder re-runs the detector at doubled n_probe and the transport at doubled n_points. Starting low skips those iterations.

  • At ``rtol = 1e-2`` the sign flips. There the tolerance does not require a window at all, so a low threshold opens one that is not needed: 0.21 s becomes 0.95 s at d = 3, buying an accuracy improvement (2.67e-03 to 1.44e-09) nobody asked for.

  So the brief’s hypothesis is confirmed: the right value is a rule, not a constant – low when the requested tolerance is tight, high when it is loose. 0.1 is not the optimum at the default rtol = 1e-3, where 0.01 is 2-3x cheaper at identical accuracy. It is deliberately left unchanged: three profiles at one energy is precisely the size of population that made GAMMA_TO_ERROR wrong twice, and retuning a default on it would repeat that mistake rather than learn from it. The measurement is recorded here so the next person starts from evidence.

• min_threshold (float, optional) –

  Floor below which the threshold is not tightened further. Default: 1e-6.

  Provenance. Swept over 1e-4 to 1e-8 across 18 ordinary workloads: the worst error is 4.49e-04 at every value, because the ladder stops long before reaching the floor. Even at rtol = atol = 1e-12 on a multi-resonance profile it converges in 7 iterations with a window open, and every value from 1e-4 to 1e-10 gives an identical answer in identical time.

  The regime this constant governs, found by construction rather than assumed: the floor is reached only when \(\gamma_\max\) is below it – so no window can ever open, however far the threshold falls – and the requested tolerance is tighter than GAMMA_TO_ERROR \(\times \gamma_\max\), so the \(\gamma\) rule cannot certify either. An almost-flat profile (\(\gamma_\max = 3\times10^{-7}\)) at rtol = atol = 1e-9 satisfies both:

  min_threshold	error	windows	iterations	time
  1e-4	8.49e-13	0	9	3.3 s
  1e-6	8.49e-13	0	13	7.4 s
  1e-8	3.24e-12	1	13	7.7 s
  1e-10	3.24e-12	1	13	7.9 s

  So it does change behaviour there – below \(\gamma_\max\) a window opens – but not usefully: the result is certified=False at every value, and the error is three orders inside the requested tolerance either way, with the window costing a factor of 2.4 in time and making the answer very slightly worse. The floor decides how much work is done in a regime where the answer is already good and known to be uncertified; it does not decide correctness anywhere measured.

• n_probe0 (int, optional) –

  Starting number of positions used to locate resonance candidates. Default: 200.

  Provenance. Swept over 50, 100, 200, 400, 800 across 18 workloads spanning single points, baseline scans and energy scans, 3 profile families, d = 2 and 3 (docs/dev/adversarial_batteries/constants_audit2.py). Worst error over all workloads: 4.98e-04, 5.22e-04, 4.49e-04, 3.38e-04, 3.36e-04 – flat within a factor of 1.6 and inside the requested 1e-3 everywhere. Not load-bearing: the refinement doubles it, so the starting value only shifts which iteration finds a given feature.

• max_n_probe (int, optional) – Ceiling on the probe grid density. Default: 6400. A cost ceiling rather than a calibration – reaching it is reported rather than absorbed – and it also sets what find_hidden_features() treats as resolvable.

• n_points0 (int, optional) –

  Starting number of positions used for adiabatic-transport quadrature. Default: 201.

  Provenance. Swept over 51, 101, 201, 401, 801 on the same 18 workloads: the worst error is 4.49e-04 at every value, identical to three digits. The refinement doubles this too, so the starting value is invisible in the answer; it buys only iterations.

• max_n_points (int, optional) – Ceiling on the adiabatic-transport grid density. Default: 12864. A cost ceiling, as max_n_probe.

• fd_step_frac (float, optional) – Finite-difference step for the Hellmann-Feynman diagnostics, as a fraction of l1 - l0. Default: 1e-6.

• max_iters (int, optional) – Maximum number of refinement iterations. Default: 12.

• info (dict, optional) – If given, filled in place with why this call ended as it did, following the same out-parameter convention as convergence_info in magnus.oscprob.osc_prob(). Keys: 'resolved' (whether H_func passed the probe-scale resolution test – see _profile_is_resolved), 'gamma_max', 'n_windows', 'iterations', and 'patches_converged'. certified=False on its own does not say which of these failed, and the cures are different: an unresolved profile wants t_breakpoints, an exhausted refinement wants a looser tolerance. magnus.oscprob uses 'resolved' to raise magnus.oscprob.UnmarkedDiscontinuityWarning on the hybrid path instead of declining in silence. Default: None.

Returns:

The evolution operator (exactly unitary regardless of certified), the non-adiabatic windows used in the last iteration, and whether the result is certified (True). certified is False if the refinement exhausted max_iters, if every knob reached its ceiling before two successive results agreed, if a local patch failed to converge within its own slab cap, or if H_func is not resolved at the probe scale (see _profile_is_resolved) – in all of these the returned operator is the best available estimate, still exactly unitary, but its accuracy is not certified to the requested tolerance.

Return type:

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

Notes

Two successive results agreeing is necessary but not sufficient, and step 4 above states the reason narrowly. When no window opens at all, successive iterations differ only in the adiabatic-transport grid: they converge to the same adiabatic limit and agree with each other whether or not that limit is the right answer. Certifying an empty window list therefore additionally requires the adiabaticity parameter itself to be small enough for the requested tolerance (see GAMMA_TO_ERROR); otherwise the loop keeps lowering the threshold until a window does open. Without that, a profile whose \(\gamma\) stays just below threshold0 everywhere is certified while wrong – measured at 1.8e-02 against a requested 1e-3.

Examples

import numpy as np

import magnus.globaldefs as gd
from magnus import adiabatic
from magnus.hamiltonians import hamiltonians3nu

p = gd.OSC_PARAMS_PREDEFINED['OSC_PARAMS_DEFAULT']
h_vac = np.asarray(hamiltonians3nu.hamiltonian_3nu_vacuum_energy_independent(
    p['s12'], p['s23'], p['s13'], p['dCP'], p['D21'], p['D31']))
e00 = np.diag([1.0, 0.0, 0.0])
energy = 10.0e6

def H(l):
    v = 1.0e-11*np.exp(-np.asarray(l, dtype=float)/gd.L_SCALE_SUN)
    return (1.0/energy)*h_vac + np.asarray(v)[..., None, None]*e00

U, windows, certified = adiabatic.hybrid_propagator(
    H, 0.0, 3.0*gd.L_SCALE_SUN)

print('non-adiabatic windows :', len(windows))
print('certified             :', certified)
print('unitary to %.1e'
      % np.max(np.abs(np.asarray(U).conj().T @ U - np.eye(3))))

non-adiabatic windows : 0
certified             : True
unitary to 5.6e-16

certified is the value to check: it says the result agreed with itself under tightening, not that it is correct. Zero windows means the whole trajectory was adiabatic, which for the real solar mixing angle is the usual answer.