Physics API

OscillationConfig{F, I, P, S, E}

Master configuration for the neutrino oscillation probability engine.

Each type parameter selects a different physical model via multiple dispatch. Construct with keyword arguments and pass to configure to obtain a ready-to-use Osc physics module.

Fields

• flavour::F = ThreeFlavour(): neutrino mass/flavour model (FlavourModel).
• interaction::I = Vacuum(): matter interaction model (InteractionModel).
• propagation::P = Basic(): propagation mechanism (PropagationModel).
• states::S = All(): eigenstate selection strategy (StateSelector).
• eigen_method::E = DefaultEigen(): eigendecomposition algorithm (EigenMethod).

Examples

# Standard 3-flavour vacuum oscillations (all defaults)
cfg = OscillationConfig()

# Inverted ordering with MSW matter effects
cfg = OscillationConfig(flavour=ThreeFlavour(ordering=:IO), interaction=SI())

# ADD model with decoherence and Barger eigensolver
cfg = OscillationConfig(
    flavour=ADD(N_KK=3),
    interaction=SI(),
    propagation=Decoherent(σₑ=0.05),
    eigen_method=BargerEigen()
)

Osc <: Newtrinos.Physics

Configured oscillation physics module, returned by configure.

This struct satisfies the Newtrinos.Physics interface and can be passed into any experiment's configure(physics=...) method.

Fields

• cfg::OscillationConfig: the full configuration used to build this module.
• params::NamedTuple: default oscillation parameter values (e.g. θ₁₂, Δm²₃₁, …).
• priors::NamedTuple: prior distributions for each parameter.
• matrices::Function: closure matrices(params) -> (U, h) returning mixing matrix and mass-squared eigenvalues.
• osc_prob::Function: closure computing oscillation probabilities. See get_osc_prob for the two call signatures.

configure(cfg::OscillationConfig=OscillationConfig()) -> Osc

Create a fully configured oscillation physics module from the given configuration.

Assembles default parameters, priors, mixing-matrix closure, and oscillation-probability closure into an Osc struct ready for use in experiment forward models.

Arguments

• cfg::OscillationConfig: oscillation configuration (defaults to standard 3-flavour vacuum).

Returns

An Osc instance.

Examples

using Newtrinos

# Default 3-flavour vacuum
physics = Newtrinos.osc.configure()

# With matter effects and inverted ordering
physics = Newtrinos.osc.configure(
    OscillationConfig(flavour=ThreeFlavour(ordering=:IO), interaction=SI())
)

ftype

Type alias for Float64, used as the default floating-point precision throughout the oscillation module.

Layer{T, U}

A spherical shell of matter with uniform proton and neutron number densities.

Used to represent concentric layers of the Earth in the PREM density model (see Newtrinos.earth_layers.EarthLayers in earth_layers.jl). A neutrino trajectory through the Earth is described by a sequence of Path segments, each referencing one layer.

Fields

• radius::T: outer radius of the shell [km].
• p_density::U: proton number density [mol/cm$^3$].
• n_density::U: neutron number density [mol/cm$^3$].

Path

A single segment of a neutrino propagation path through one matter Layer.

A full baseline is represented as a Vector{Path}, one entry per layer traversed. For vacuum oscillations the total baseline is simply sum(segment.length for segment in path).

Fields

• length::Float64: segment length [km].
• layer_idx::Int: index into the corresponding StructVector{Layer}.

get_osc_prob(cfg::OscillationConfig) -> Function

Construct the oscillation probability closure for the given configuration.

The returned function osc_prob has two methods:

Vacuum-style (baselines as distances):

osc_prob(E::AbstractVector, L::AbstractVector, params::NamedTuple; anti=false)

Matter-style (baselines as Earth layer paths):

osc_prob(E::AbstractVector, paths::VectorOfVectors{Path}, layers::StructVector{Layer},
         params::NamedTuple; anti=false)

Arguments

• E: neutrino energies [GeV].
• L: baselines [km] (vacuum method).
• paths: per-baseline layer traversals (matter method), from Newtrinos.earth_layers.EarthLayers.
• layers: Earth density layers (matter method).
• params::NamedTuple: oscillation parameters (mixing angles, mass splittings, etc.).
• anti::Bool = false: set true for antineutrinos (conjugates the PMNS matrix).

Returns

Array{T,4} of shape (n_E, n_L, n_flav, n_flav) where entry result[i, j, α, β] gives the transition probability $P(\nu_\alpha \to \nu_\beta)$ at energy E[i] and baseline/path index j (3rd index = input/source flavour, 4th index = output/detected flavour).

get_matrices(cfg::FlavourModel, eigen_method::EigenMethod=DefaultEigen()) -> Function

Construct a closure that computes the mixing matrix and mass-squared eigenvalues from oscillation parameters.

The returned function has signature matrices(params::NamedTuple) -> (U, h) where:

• U is the mixing matrix (unitary).
• h is a vector of mass-squared eigenvalues [eV²].

For ThreeFlavour, returns SMatrix{3,3} and SVector{3} (zero-allocation). For Sterile, returns 4×4 dense matrices. For ADD and dark-dimension models, returns $3(N_{KK}+1) \times 3(N_{KK}+1)$ dense matrices obtained by diagonalizing the full KK mass matrix via decompose. For NND and NNM returns the eigenvectors matrix 3Nx3N, a 3N vector of mass squared differences (from m1^2), a 3N vector of mass eigenvalues, 3 NxN eigenvector matrices (one per-flavour) Ve, Vm, V_t.

The closure is ForwardDiff-compatible: all intermediate computations preserve dual-number types through zero(T) / one(T) patterns and type promotion.

Arguments

• cfg::FlavourModel: flavour model instance (dispatches to the appropriate method).
• eigen_method::EigenMethod: eigendecomposition algorithm (only used by models that diagonalize a mass matrix, e.g. ADD).

See also get_PMNS, configure.

get_PMNS(params) -> SMatrix{3,3}

Construct the $3 \times 3$ PMNS leptonic mixing matrix from oscillation parameters.

The matrix is built as

\[U_{PMNS} = R_{23}(\theta_{23}) \;\cdot\; R_{13}(\theta_{13},\, \delta_{CP}) \;\cdot\; R_{12}(\theta_{12})\]

where each rotation matrix uses zero(T) / one(T) to preserve ForwardDiff.Dual number types for automatic differentiation.

Arguments

• params::NamedTuple: must contain θ₂₃, θ₁₃, θ₁₂, and δCP.

Returns

SMatrix{3,3,Complex{T}} — the PMNS mixing matrix, where T is inferred from the parameter types.

get_abs_masses(params) -> (m1, m2, m3)

Compute the absolute neutrino masses from the lightest mass and mass-squared splittings.

For normal ordering ($\Delta m^2_{31} > 0$): $m_1 = m_0$, $m_2 = \sqrt{\Delta m^2_{21} + m_0^2}$, $m_3 = \sqrt{\Delta m^2_{31} + m_0^2}$.

For inverted ordering ($\Delta m^2_{31} < 0$): $m_3 = m_0$, $m_1 = \sqrt{-\Delta m^2_{31} + m_0^2}$, $m_2 = \sqrt{\Delta m^2_{21} - \Delta m^2_{31} + m_0^2}$.

Arguments

• params::NamedTuple: must contain m₀, Δm²₂₁, and Δm²₃₁.

Returns

A tuple (m1, m2, m3) of absolute neutrino masses [eV].

FlavourModel

Abstract type selecting the neutrino mass/flavour model.

Each subtype defines its own set of oscillation parameters (mixing angles, mass splittings, and any BSM quantities) via get_params and get_priors dispatches, as well as a get_matrices method that returns the mixing matrix $U$ and mass-squared eigenvalues.

Standard and BSM subtypes:

• ThreeFlavour: standard three-neutrino oscillations.
• ThreeFlavourXYCP: three-flavour with shell parametrization of $\delta_{CP}$.
• Sterile: 3+1 sterile neutrino model.
• ADD: Arkani-Hamed–Dimopoulos–Dvali large extra dimensions.
• Darkdim_Lambda, Darkdim_Masses, Darkdim_cas: dark-dimension model variants.
• 'NND': N-Naturalness with Dirac neutrinos.
• 'NNM': N-Naturalness with Majorana neutrinos.

ThreeFlavour <: FlavourModel

Standard three-flavour neutrino oscillation model.

The PMNS mixing matrix is constructed as $U_{PMNS} = R_{23}(\theta_{23}) \cdot R_{13}(\theta_{13},\, \delta_{CP}) \cdot R_{12}(\theta_{12})$ and the oscillation parameters are:

Fields

• ordering::Symbol = :NO: neutrino mass ordering. :NO for normal ordering ($\Delta m^2_{31} > 0$), :IO for inverted ordering ($\Delta m^2_{31} < 0$).

ThreeFlavourXYCP <: FlavourModel

Three-flavour model with a shell (Cartesian) parametrization of $\delta_{CP}$.

Replaces the scalar δCP parameter with a 2-vector δCPshell to avoid discontinuities at the circular boundary $[0, 2\pi)$ during gradient-based fitting.

Fields

• three_flavour::ThreeFlavour = ThreeFlavour(): underlying three-flavour configuration.

See also ThreeFlavour.

Sterile <: FlavourModel

3+1 sterile neutrino model extending the standard three flavours.

Adds a fourth mass eigenstate with mass-squared splitting Δm²₄₁ and three new mixing angles θ₁₄, θ₂₄, θ₃₄. The PMNS matrix is extended to 4×4 via

\[U_{4\times4} = R_{34}\, R_{24}\, R_{14} \cdot \begin{pmatrix} U_{PMNS} & 0 \\ 0 & 1 \end{pmatrix}\]

Fields

• three_flavour::ThreeFlavour = ThreeFlavour(): underlying three-flavour configuration.

See also ThreeFlavour.

ADD <: FlavourModel

Arkani-Hamed–Dimopoulos–Dvali (ADD) large extra dimensions model.

Introduces a tower of $N_{KK}$ Kaluza-Klein excitations for each active neutrino, yielding a $3(N_{KK}+1) \times 3(N_{KK}+1)$ mass matrix. The Dirac mass matrix is constructed from the PMNS matrix and absolute neutrino masses, with off-diagonal couplings controlled by the extra-dimension compactification radius.

Additional parameters beyond ThreeFlavour:

Fields

• three_flavour::ThreeFlavour = ThreeFlavour(): underlying three-flavour configuration.
• N_KK::Int = 5: number of Kaluza-Klein modes to include.

See also Cut for excluding heavy KK states from coherent oscillation.

get_params(cfg::OscillationConfig) -> NamedTuple
get_params(cfg::FlavourModel) -> NamedTuple

Return the default oscillation parameter values for the given configuration or flavour model.

Delegates to the flavour-model-specific method. Each FlavourModel subtype defines its own set of parameters (mixing angles, mass splittings, and any BSM quantities).

Arguments

• cfg: an OscillationConfig or a FlavourModel subtype instance.

Returns

A NamedTuple of Symbol => Float64 (or Vector for shell parameters) with the nominal parameter values.

get_priors(cfg::OscillationConfig) -> NamedTuple
get_priors(cfg::FlavourModel) -> NamedTuple

Return prior distributions for each oscillation parameter.

Delegates to the flavour-model-specific method. Priors are Distributions.jl objects (typically Uniform or LogUniform) used for sampling and inference.

Arguments

• cfg: an OscillationConfig or a FlavourModel subtype instance.

Returns

A NamedTuple of Symbol => Distribution mapping each parameter to its prior.

InteractionModel

Abstract type selecting the neutrino-matter interaction model used when propagating through Earth layers. Subtypes:

• Vacuum: no matter effects.
• SI: Standard Interactions (MSW effect).
• NSI: Non-Standard Interactions.

Vacuum <: InteractionModel

Vacuum oscillations — no matter effects. Layer densities in Layer are ignored; the baseline is the sum of Path segment lengths.

SI <: InteractionModel

Standard Interactions — coherent forward scattering in matter (the MSW effect).

Adds the Wolfenstein charged-current potential $V_{CC} = \sqrt{2}\, G_F\, N_e$ to the $(\nu_e, \nu_e)$ element and a neutral-current potential to all diagonal elements of the effective Hamiltonian. The sign is flipped for antineutrinos.

NSI <: InteractionModel

Non-Standard Interactions in matter. Extends the matter Hamiltonian beyond the Standard Model Wolfenstein potential.

PropagationModel

Abstract type governing how neutrino quantum states evolve along the baseline.

Subtypes select different physical approximations for the oscillation amplitude:

• Basic: standard coherent quantum-mechanical propagation.
• Decoherent: density-matrix evolution with off-diagonal damping.
• Damping: amplitude-level low-pass filter with incoherent recovery.
• Spray: analytic ray-to-spray averaging over nearby energy/angle trajectories.

Basic <: PropagationModel

Standard coherent quantum-mechanical propagation (no decoherence effects).

The oscillation amplitude is computed as

\[A_{\alpha\beta} = \bigl(U \, \mathrm{diag}\bigl(e^{-i\,\Delta m^2_j\, L \,/\, 4E}\bigr) \, U^\dagger\bigr)_{\alpha\beta}\]

and the transition probability is $P_{\alpha\beta} = |A_{\alpha\beta}|^2$.

Decoherent <: PropagationModel

Propagation with energy-dependent decoherence via density-matrix evolution.

Off-diagonal elements of the density matrix in the mass eigenbasis are damped by

\[D_{ij} = \exp\!\bigl(-2\,|\Delta\phi_{ij}|\,\sigma_e^2\bigr), \qquad \Delta\phi_{ij} = \frac{\Delta m^2_{ij}\, L}{4E}\]

after coherent phase evolution at each layer crossing.

Fields

• σₑ::Float64 = 0.1: decoherence strength parameter (dimensionless).

Damping <: PropagationModel

Propagation with a low-pass filter that damps fast-oscillating amplitudes.

Each mass-eigenstate phase factor is multiplied by a decay $d_j = \exp(-2\,|\phi_j|\,\sigma_e^2)$ and the probability is corrected by an incoherent recovery term $|U|^2 \, \mathrm{diag}(1-d_j^2) \, |U|^{2\,\prime}$.

Fields

• σₑ::Float64 = 0.1: damping strength parameter (dimensionless).

StateSelector

Abstract type controlling which mass eigenstates contribute to the coherent oscillation calculation. Subtypes:

• All: include every eigenstate (default).
• Cut: exclude eigenstates above a mass-squared cutoff.

All <: StateSelector

Include all mass eigenstates in the oscillation calculation. This is the default.

Cut <: StateSelector

Exclude mass eigenstates whose $\sqrt{|m^2_j|}$ exceeds a cutoff threshold.

States above the cutoff are averaged incoherently, which is useful for ADD or dark-dimension models where heavy Kaluza-Klein states oscillate too rapidly to be resolved experimentally.

Fields

• cutoff::Float64 = Inf: mass eigenvalue cutoff [eV]. Default Inf keeps all states.

EigenMethod

Abstract type selecting the eigendecomposition algorithm for the effective matter Hamiltonian. Subtypes:

• DefaultEigen: Julia's built-in LinearAlgebra.eigen.
• BargerEigen (defined in barger_eigen.jl): fast analytic 3×3 decomposition.

DefaultEigen <: EigenMethod

Use Julia's built-in LinearAlgebra.eigen for Hermitian matrix diagonalization.

decompose(H::Hermitian, method::EigenMethod) -> Eigen

Eigendecompose the Hermitian Hamiltonian H using the selected EigenMethod.

This is the extension point for custom eigensolvers. To add a new method, define a subtype of EigenMethod and a corresponding decompose dispatch.

Arguments

• H::Hermitian: effective Hamiltonian matrix (mass-squared basis or matter-modified).
• method::EigenMethod: algorithm selector.

Returns

An Eigen factorization with fields vectors and values.

BargerEigen <: Newtrinos.osc.EigenMethod

Analytic 3×3 Hermitian eigendecomposition using closed-form expressions. Uses the method of Deledalle et al. 2017, which computes eigenvalues via the cubic formula and eigenvectors via closed-form expressions.

Only valid for 3×3 Hermitian matrices (i.e. three-flavour oscillations). Falls back to Julia's eigen for non-3×3 matrices.

fast_eigen(C::Hermitian{T, <:SMatrix{3,3,T}}; sortby=identity) -> Eigen

Analytic eigendecomposition of a 3×3 Hermitian matrix using closed-form expressions.

Implements the algorithm of Deledalle et al. (2017), hal-01501221. This avoids the iterative LAPACK eigen call and is significantly faster for the small matrices that arise in three-flavour neutrino oscillations with matter effects.

Called by Newtrinos.BargerEigen via Newtrinos.osc.decompose.

Eigenvalue algorithm: Computes the three roots of the characteristic polynomial via the Cardano trigonometric method. With $x_1$ and $x_2$ defined from the matrix elements, the discriminant $\Delta = 4x_1^3 - x_2^2 \geq 0$ and the roots are

\[\lambda_k = \frac{a+b+c}{3} + \frac{2\sqrt{x_1}}{3} \cos\!\left(\frac{\phi + 2\pi(k-1)}{3}\right), \quad k = 1,2,3\]

where $\phi = \operatorname{atan}(\sqrt{\Delta},\, x_2)$.

Eigenvector algorithm: Each eigenvector is obtained as the cross product of two rows of $(C - \lambda_k I)$. Three row-pair combinations are tried in order to handle degenerate or near-degenerate cases where a cross product vanishes.

Dispatches to specialized helpers _fast_eigen_efzero, _fast_eigen_edzero, _fast_eigen_dfzero when pairs of off-diagonal elements are zero (reducing to 2×2 subproblems).

The matrix elements are labelled as follows (a, b, c are real diagonal entries; d, e, f are complex lower-triangle entries):

\[C = \begin{pmatrix} a & d^* & f^* \\ d & b & e^* \\ f & e & c \end{pmatrix}\]

Arguments

• C::Hermitian{T, <:SMatrix{3,3,T}}: the 3×3 Hermitian matrix to decompose.
• sortby::F = identity: ordering function applied to eigenvalues (e.g. abs, real).

Returns

An Eigen object with values::SVector{3} and vectors::SMatrix{3,3}.

_fast_eigen_efzero(C; sortby=identity) -> Eigen

Edge-case eigendecomposition when C[3,1] == C[3,2] == 0 (off-diagonal elements e=f=0).

Row 3 decouples: $\lambda_3 = C_{33}$ with eigenvector $(0, 0, 1)$. The remaining two eigenvalues and vectors are obtained by calling Julia's eigen on the 2×2 upper-left submatrix $C_{1:2,1:2}$, then embedding the result back into 3×3.

_fast_eigen_edzero(C; sortby=identity) -> Eigen

Edge-case eigendecomposition when C[1,2] == C[2,3] == 0 (off-diagonal elements e=d=0).

Row 2 decouples: $\lambda_2 = C_{22}$ with eigenvector $(0, 1, 0)$. The remaining two eigenvalues are found from the 2×2 submatrix of rows/columns 1 and 3.

_fast_eigen_dfzero(C; sortby=identity) -> Eigen

Edge-case eigendecomposition when C[1,2] == C[1,3] == 0 (off-diagonal elements d=f=0).

Row 1 decouples: $\lambda_1 = C_{11}$ with eigenvector $(1, 0, 0)$. The remaining two eigenvalues are found from the 2×2 lower-right submatrix $C_{2:3,2:3}$.

_sorted_eigen(::Nothing, vals, vecs) -> Eigen

Passthrough overload that wraps eigenvalues and eigenvectors into an Eigen object without reordering. Called when sortby = nothing.

_sorted_eigen(by::F, vals, vecs) -> Eigen

Sort eigenvalues by the key function by and permute eigenvector columns accordingly.

Arguments

• by::F: ordering key function (e.g. identity, abs, real).
• vals: eigenvalue vector.
• vecs: eigenvector matrix (columns correspond to eigenvalues).

Returns

An Eigen object with eigenvalues sorted by by and eigenvectors permuted to match.

_sorted_eigen(by::Fun, F::Eigen) -> Eigen

Convenience overload that unpacks an Eigen object and delegates to _sorted_eigen(by, F.values, F.vectors).

osc_kernel(U, H, e, l) -> AbstractMatrix

Compute the coherent oscillation amplitude matrix for a single energy and baseline.

Evaluates $A = U \, \mathrm{diag}\!\bigl(e^{-i\, \phi_j}\bigr) \, U^\dagger$ where $\phi_j = \frac{\Delta m^2_j \, L}{4E}$ (in natural units via F_units).

Arguments

• U::AbstractMatrix{<:Number}: mixing matrix (unitary, $n \times n$).
• H::AbstractVector{<:Number}: mass-squared eigenvalues [eV²].
• e::Real: neutrino energy [GeV].
• l::Real: baseline length [km].

Returns

An $n \times n$ complex amplitude matrix $A_{\beta\alpha}$.

osc_kernel(U, H, e, l, σₑ) -> (A, decay)

Compute the oscillation amplitude matrix with low-pass damping.

Same as the coherent osc_kernel but each phase factor is additionally multiplied by a decay $d_j = \exp(-2\,|\phi_j|\,\sigma_e^2)$.

Arguments

• U::AbstractMatrix{<:Number}: mixing matrix.
• H::AbstractVector{<:Number}: mass-squared eigenvalues [eV²].
• e::Real: neutrino energy [GeV].
• l::Real: baseline length [km].
• σₑ::Real: damping strength parameter.

Returns

A tuple (A, decay) where A is the damped amplitude matrix and decay is the per-eigenstate decay vector.

compute_matter_matrices(H_eff, e, layer, anti, interaction::SI, eigen_method=DefaultEigen()) -> (U, h)

Add the MSW matter potential to the effective Hamiltonian and diagonalize.

Modifies H_eff by adding the Wolfenstein charged-current and neutral-current potentials for the given Layer density, then eigendecomposes the result via decompose. Two methods are provided: one for generic AbstractMatrix and an optimized one for SMatrix{3,3}.

Arguments

• H_eff: effective Hamiltonian in the flavour basis ($U\,\mathrm{diag}(h)\,U^\dagger$).
• e::Real: neutrino energy [GeV].
• layer::Layer: matter layer with proton and neutron densities.
• anti::Bool: true for antineutrinos (flips potential sign).
• interaction::SI: matter interaction selector.
• eigen_method::EigenMethod: eigendecomposition algorithm.

Returns

A tuple (U_m, h_m) of matter-modified eigenvectors and eigenvalues.

osc_reduce(matter_matrices, path, e, propagation) -> Matrix

Combine per-layer oscillation amplitudes along a multi-layer path into a single transition probability matrix for one energy.

Dispatches on PropagationModel:

• Basic: multiplies amplitude matrices across layers, then takes $|\cdot|^2$.
• Damping: multiplies damped amplitudes and adds an incoherent recovery term using a path-length-weighted average of $|U|^2$.

Arguments

• matter_matrices: vector of (U_m, h_m) tuples, one per layer (from compute_matter_matrices).
• path::Vector{Path}: sequence of path segments for this baseline.
• e::Real: neutrino energy [GeV].
• propagation::PropagationModel: propagation model instance.

Returns

An $n \times n$ real matrix of transition probabilities $P_{\beta\alpha}$.

matter_osc_per_e(H_eff, e, layers, paths, anti, propagation, interaction, eigen_method=DefaultEigen()) -> Array

Compute matter-affected oscillation probabilities at a single energy across all paths.

For each Layer, computes the matter-modified eigensystem via compute_matter_matrices, then reduces along each path via osc_reduce (for Basic/Damping) or via explicit density-matrix evolution (for Decoherent).

Arguments

• H_eff: effective Hamiltonian in the flavour basis.
• e::Real: neutrino energy [GeV].
• layers::StructVector{Layer}: Earth density layers.
• paths: per-baseline layer traversals (VectorOfVectors{Path}).
• anti::Bool: true for antineutrinos.
• propagation::PropagationModel: propagation model.
• interaction::InteractionModel: matter interaction model.
• eigen_method::EigenMethod: eigendecomposition algorithm.

Returns

An Array of shape (n_flav, n_flav, n_paths) with $P_{\beta\alpha}$ for each path.

select(U, h, cfg::StateSelector) -> (U_sel, h_sel, rest)

Partition mass eigenstates into coherent and incoherent subsets based on the StateSelector strategy.

• All: returns all states unchanged with rest = 0.0.
• Cut: states with $\sqrt{|m^2_j|}$ below the cutoff are kept for coherent oscillation; states above it are averaged incoherently and their contribution is returned as the rest matrix.

Arguments

• U: mixing matrix ($n_{flav} \times n_{states}$).
• h: mass-squared eigenvalues vector.
• cfg::StateSelector: selection strategy.

Returns

A tuple (U_sel, h_sel, rest) where:

• U_sel: mixing matrix columns for coherent states.
• h_sel: eigenvalues for coherent states.
• rest: incoherent probability correction matrix (scalar 0.0 or $n \times n$ matrix).

propagate(U, h, E, L, propagation) -> Array{T,4}
propagate(U, h, E, paths, layers, propagation, interaction, anti, ...) -> Array{T,4}

Compute oscillation probabilities over grids of energy and baseline/path.

Vacuum methods (dispatched by baseline vector L): Loops over all (E, L) pairs, evaluating osc_kernel and squaring amplitudes. Dispatches on PropagationModel for Basic, Damping, and Decoherent propagation.

Matter methods (dispatched by paths::VectorOfVectors{Path}):

• Vacuum interaction: collapses paths to total lengths and delegates to the vacuum method.
• SI/NSI: reconstructs the flavour-basis Hamiltonian and calls matter_osc_per_e for each energy.

Arguments

• U: mixing matrix.
• h: mass-squared eigenvalues [eV²].
• E: neutrino energies [GeV].
• L: baselines [km] (vacuum) or paths/layers (matter).
• propagation::PropagationModel: propagation model.
• interaction::InteractionModel: matter interaction model (matter methods only).
• anti::Bool: true for antineutrinos (matter methods only).

Returns

Array{T,4} of shape (n_flav, n_flav, n_E, n_L) with $P_{\beta\alpha}$ (note: the caller get_osc_prob transposes this to (n_E, n_L, n_flav, n_flav)).

DensityModel

Abstract type for Earth density profile models.

Each subtype provides a recipe for dividing the Earth into concentric shells of constant density. Currently the only implementations are PREM and VariableDensity.

PREM <: DensityModel

Preliminary Reference Earth Model (Dziewonski & Anderson, 1981).

Reads the tabulated PREM density profile from PREM_1s.csv and groups radial shells into zones defined by density boundaries.

Fields

• zones::Array{Float64} = [0., 4., 7.5, 12.5, 13.1]: density boundaries [g/cm³] defining the constant-density zones. Adjacent PREM rows whose density falls within the same bin are averaged into a single Newtrinos.osc.Layer.
• p_fractions::Vector{Float64} = [0.496, 0.494, 0.468, 0.466]: proton number fraction $Y_p = N_p / (N_p + N_n)$ for each density zone (one value per interval defined by zones), rather than a single value assumed constant across all layers.
• atm_heihgt::Float64 = 20.: atmospheric shell thickness [km] added above the Earth's surface (density = 0).

EarthLayers <: Newtrinos.Physics

Configured Earth density model, returned by configure.

Provides the layer structure and path-computation functions needed by the oscillation module's matter-effect calculations (see Newtrinos.osc.SI).

Fields

• cfg::DensityModel: the density model used to build this module.
• params::NamedTuple: oscillation parameters.
• priors::NamedTuple: prior distributions.
• compute_layers::Function: closure compute_layers() -> StructVector{Layer} returning the concentric density shells.
• compute_paths::Function: function compute_paths(cz, layers; r_detector) -> VectorOfVectors{Path} computing the layer traversal for each cosine-zenith value.

configure(cfg::PREM=PREM()) -> EarthLayers

Create a fully configured Earth density physics module.

Arguments

• cfg::PREM: density model (defaults to PREM).

Returns

An EarthLayers instance with compute_layers and compute_paths closures and empty parameters and priors.

Examples

using Newtrinos

earth = Newtrinos.earth_layers.configure()
layers = earth.compute_layers()
paths = earth.compute_paths([-1.0, -0.5, 0.0], layers)

configure(cfg::VariableDensity) -> EarthLayers

Create an Earth density physics module with a floatable overall density normalization.

Uses the same compute_layers/compute_paths closures as configure(::PREM) (built from cfg.prem), but additionally registers an electron_density_scale parameter (nominal 1.0, prior Normal(1.0, 0.068)). The scale is not applied by compute_layers itself — callers that want it applied must pass the parameter to scale_densities explicitly.

Arguments

• cfg::VariableDensity: density model wrapping the PREM profile to be scaled.

Returns

An EarthLayers instance with the electron_density_scale parameter and prior.

get_compute_layers(cfg::PREM) -> Function

Construct a closure that builds the Earth layer structure from the PREM density profile.

The returned function compute_layers() reads PREM_1s.csv, bins rows by the density boundaries in cfg.zones, and returns a StructVector{Layer} with one entry per zone plus an atmospheric shell.

Arguments

• cfg::PREM: PREM configuration with zone boundaries, per-zone proton fractions, and atmosphere height.

Returns

A zero-argument closure compute_layers() -> StructVector{Layer}.

compute_paths(cz::Number, layers, r_detector) -> StructArray{Path}
compute_paths(cz::AbstractArray, layers; r_detector=6369) -> VectorOfVectors{Path}

Compute the sequence of Newtrinos.osc.Path segments a neutrino traverses through the Earth.

For a single cosine-zenith value, determines which Newtrinos.osc.Layer shells are crossed using ray_circle_path_length, then builds an ordered list of segments. Layers below the detector are traversed twice (entry and exit), while layers above the detector are traversed once.

The array method broadcasts over multiple cosine-zenith values and returns a VectorOfVectors{Path}.

Arguments

• cz: cosine of the zenith angle (scalar or array). −1 is vertically upgoing.
• layers::StructVector{Layer}: Earth density layers from compute_layers().
• r_detector::Real: radial position of the detector [km] (default 6369, approximate IceCube depth).

Returns

• Scalar method: StructArray{Path} with length and layer_idx columns.
• Array method: VectorOfVectors{Path}, one Vector{Path} per zenith angle.

ray_circle_path_length(r, y, cz) -> Real

Compute the chord length of a ray through a circle of radius r.

The ray originates at radial distance y from the Earth's centre with direction cosine cz (cosine of the zenith angle). Returns zero if the ray does not intersect the circle or if the chord is shorter than 1 km (numerical noise filter).

Arguments

• r: radius of the spherical shell [km].
• y: radial position of the detector [km].
• cz: cosine of the zenith angle (−1 = vertically upgoing through the core).

Returns

Path length through the shell [km], or zero if no intersection.

NominalFluxModel

Abstract type for atmospheric neutrino flux predictions.

Each subtype provides a recipe for computing the unoscillated neutrino flux as a function of energy and zenith angle. Currently the only implementation is HKKM.

HKKM <: NominalFluxModel

Honda–Kajita–Kasahara–Midorikawa (HKKM) atmospheric neutrino flux model.

Reads a site-specific flux table (.d file) and constructs cubic-spline interpolations over $\log_{10}(E)$ and $\cos\theta_z$ for each neutrino flavour ($\nu_\mu$, $\bar\nu_\mu$, $\nu_e$, $\bar\nu_e$).

Fields

• fname::String = "spl-nu-20-01-000.d": filename of the HKKM flux table (relative to the physics/ data directory).

FluxSystematicsModel

Abstract type for atmospheric flux systematic uncertainty models.

Each subtype defines a set of nuisance parameters and a function that modifies the nominal flux to account for systematic uncertainties. Implementations are Barr and BarrEnergyBands.

Barr <: FluxSystematicsModel

Barr systematic uncertainty model for atmospheric neutrino fluxes.

Parametrizes flux uncertainties as energy- and zenith-dependent modifications following Barr et al., "Uncertainties in Atmospheric Neutrino Fluxes." The systematics cover:

All sigma parameters are centred at 0 with unit Gaussian priors (truncated at ±3).

See BarrEnergyBands for a variant that splits each flux-ratio uncertainty into independent sub-GeV/1-10 GeV/\>10 GeV energy-range parameters with additional coszen-dependent shaping.

AtmFluxConfig{F, S}

Configuration for the atmospheric neutrino flux module.

Fields

• nominal_model::F = HKKM(): nominal flux prediction model (NominalFluxModel).
• systematics_model::S = Barr(): systematic uncertainty model (FluxSystematicsModel).

AtmFlux <: Newtrinos.Physics

Configured atmospheric flux physics module, returned by configure.

Fields

• cfg::AtmFluxConfig: the configuration used to build this module.
• params::NamedTuple: default systematic parameter values.
• priors::NamedTuple: prior distributions for each systematic parameter.
• nominal_flux::Function: closure nominal_flux(energy, coszen) -> Table returning the unoscillated flux on a fine $(E, \cos\theta_z)$ grid.
• sys_flux::Function: closure sys_flux(flux, params) -> NamedTuple applying systematic modifications to the nominal flux and returning per-flavour flux arrays.

configure(cfg::AtmFluxConfig=AtmFluxConfig()) -> AtmFlux

Create a fully configured atmospheric flux physics module.

Arguments

• cfg::AtmFluxConfig: flux configuration (defaults to HKKM nominal model with Barr systematics).

Returns

An AtmFlux instance.

Examples

using Newtrinos

# Default South Pole flux
flux_physics = Newtrinos.atm_flux.configure()

# Custom site flux file
flux_physics = Newtrinos.atm_flux.configure(
    AtmFluxConfig(nominal_model=HKKM(fname="spl-nu-20-12-000.d"))
)

get_params(cfg::FluxSystematicsModel) -> NamedTuple

Return the default systematic parameter values for the given flux systematics model.

Arguments

• cfg::FluxSystematicsModel: a Barr or BarrEnergyBands instance.

Returns

A NamedTuple mapping parameter names to their nominal values (all default to 0.0).

get_params(cfg::BarrEnergyBands) -> NamedTuple

Return the default systematic parameter values for the BarrEnergyBands model.

Returns

A NamedTuple mapping parameter names to their nominal values (all default to 0.0), covering the _lo/_mid/_hi energy-range variants of the $\nu_e/\bar\nu_e$, $\nu_\mu/\bar\nu_\mu$, and $\nu_e/\nu_\mu$ ratio uncertainties, plus the spectral index, up/horizontal, and up/down parameters.

get_priors(cfg::FluxSystematicsModel) -> NamedTuple

Return prior distributions for each flux systematic parameter.

All sigma parameters use Truncated(Normal(0, 1), -3, 3) and the spectral index uses Truncated(Normal(0, 0.1), -0.3, 0.3).

Arguments

• cfg::FluxSystematicsModel: a Barr or BarrEnergyBands instance.

Returns

A NamedTuple of Symbol => Distribution priors.

get_priors(cfg::BarrEnergyBands) -> NamedTuple

Return prior distributions for each flux systematic parameter of the BarrEnergyBands model.

All sigma parameters (including the _lo/_mid/_hi energy-range variants) use Truncated(Normal(0, 1), -3, 3) and the spectral index uses Truncated(Normal(0, 0.1), -0.3, 0.3).

Returns

A NamedTuple of Symbol => Distribution priors.

get_nominal_flux(cfg::NominalFluxModel) -> Function

Construct the nominal flux closure for the given flux model.

For HKKM, returns a function nominal_flux(energy, coszen) -> Table that evaluates the HKKM spline interpolations on a meshgrid of energy [GeV] and coszen values. The returned Table contains columns true_energy, log10_true_energy, true_coszen, and per-flavour flux values (:numu, :numubar, :nue, :nuebar).

Arguments

• cfg::NominalFluxModel: a HKKM instance.

Returns

A closure nominal_flux(energy, coszen) -> Table.

get_sys_flux(cfg::FluxSystematicsModel) -> Function

Construct the systematic flux modification closure for the given systematics model.

For Barr, returns a function sys_flux(flux, params) -> NamedTuple that applies the following energy- and zenith-dependent corrections to the nominal flux:

1. Spectral tilt: $(E / E_\text{pivot})^{\Delta\gamma}$ with $E_\text{pivot} \approx 24.1$ GeV.
2. $\nu_e / \bar\nu_e$ ratio scaling via scale_flux.
3. $\nu_\mu / \bar\nu_\mu$ ratio scaling via scale_flux.
4. $\nu_e / \nu_\mu$ ratio scaling via scale_flux.
5. Up/down asymmetry via updown.
6. Up/horizontal anisotropy via uphorizontal (separate polynomial uncertainty fits for $\nu_e$ and $\nu_\mu$).

Arguments

• cfg::FluxSystematicsModel: a Barr instance.

Returns

A closure sys_flux(flux, params) -> NamedTuple{(:nue, :numu, :nuebar, :numubar)} of modified flux arrays.

get_sys_flux(cfg::BarrEnergyBands) -> Function

Construct the systematic flux modification closure for the BarrEnergyBands model.

Returns a function sys_flux(flux, params) -> NamedTuple that applies the following energy- and zenith-dependent corrections to the nominal flux:

1. Spectral tilt: $(E / E_\text{pivot})^{\Delta\gamma}$ with $E_\text{pivot} \approx 24.1$ GeV.
2. $\nu_e / \bar\nu_e$ ratio scaling via scale_flux, with the sigma selected per energy range (sub-GeV/1–10 GeV/>10 GeV via mask_lo/mask_mid/mask_hi) and its coszen dependence shaped by fun_nuenuebar.
3. $\nu_\mu / \bar\nu_\mu$ ratio scaling via scale_flux, with the sigma selected per energy range and its coszen dependence shaped by fun_numunumubar.
4. $\nu_e / \nu_\mu$ ratio scaling via scale_flux, with the sigma selected per energy range and its coszen dependence shaped by fun_numunue.
5. Up/down asymmetry via updown.
6. Up/horizontal anisotropy via uphorizontal (separate polynomial uncertainty fits for $\nu_e$ and $\nu_\mu$).

Returns

A closure sys_flux(flux, params) -> NamedTuple{(:nue, :numu, :nuebar, :numubar)} of modified flux arrays.

get_hkkm_flux(filename) -> OrderedDict{Symbol, Extrapolation}

Read an HKKM flux table and return cubic-spline interpolations for each flavour.

Parses the fixed-format .d file into 20 cosine-zenith chunks of 101 energy bins each, then builds 2D cubic-spline interpolations over $(\log_{10}(E/\text{GeV}),\, \cos\theta_z)$.

Arguments

• filename::String: absolute path to the HKKM flux data file.

Returns

An OrderedDict with keys :numu, :numubar, :nue, :nuebar, each mapping to a Interpolations.Extrapolation with linear extrapolation boundary conditions.

scale_flux(A, B, scale) -> (mod_A, mod_B)

Scale the ratio between two flux components while conserving their sum.

Given fluxes A and B, modifies their ratio by the factor scale while keeping $A + B$ constant: $A' / B' = (A / B) \cdot \text{scale}$.

Arguments

• A: flux array for the first component.
• B: flux array for the second component.
• scale: multiplicative factor applied to the $A/B$ ratio.

Returns

A tuple (mod_A, mod_B) of rescaled flux arrays.

uphorizontal(coszen, rel_error) -> Real

Compute the up/horizontal flux anisotropy correction factor.

Models the zenith-dependent flux distortion as the ratio of an ellipse to a circle, where the ellipse semi-axes are a = 1/rel_error and b = rel_error.

Arguments

• coszen: cosine of the zenith angle.
• rel_error: relative error parameter controlling the ellipticity.

Returns

A multiplicative correction factor for the flux.

updown(coszen, up_down_ratio) -> AbstractArray

Compute the up/down flux asymmetry correction factor.

Uses a smooth $\tanh(3\,\cos\theta_z)$ transition to interpolate between 1/up_down_ratio (downgoing) and up_down_ratio (upgoing).

Arguments

• coszen: cosine of the zenith angle (scalar or array).
• up_down_ratio: ratio of upgoing to downgoing flux modification.

Returns

A multiplicative correction factor array.

SNSFlux <: Newtrinos.Physics

Configured Spallation Neutron Source (SNS) neutrino flux module.

Supports two modes: data-driven (use_data=true), which loads 2D energy-time flux tables from CSV files, and analytic (use_data=false), which computes fluxes from pion-at-rest decay kinematics.

Fields

• params::NamedTuple: flux normalization parameters (flux_norm or sns_nu_per_POT (sns neutrinos per Proton On Target)).
• priors::NamedTuple: prior distributions for each parameter.
• assets::NamedTuple: precomputed energy grids, time bins, and raw flux arrays.
• flux::Function: closure flux(params) -> NamedTuple returning the scaled flux components and total flux.

configure(; exposure, distance, ecut, E_bin_width, tcut, tbins, use_data, flux_folder, beam_power, proton_energy) -> SNSFlux

Create a fully configured SNS flux physics module.

Arguments

• exposure: integrated exposure [GWh].
• distance: source-to-detector distance [m].
• ecut::Real = mmu/2: maximum neutrino energy [MeV] (default: $m_\mu / 2$).
• E_bin_width::Real = 0.5: energy bin width [MeV] (analytic mode only).
• tcut::Real = 6000.0: maximum time [ns] (data mode only).
• tbins = Nothing: explicit time-bin edges (data mode only). When given, overrides the default log-spaced time binning in get_assets.
• use_data::Bool = true: if true, load flux from CSV tables (default); if false, use analytic decay-at-rest spectra.
• flux_folder::String: path to the directory containing the SNS flux CSV files.
• beam_power::Real = 1.4: beam power [MW] (for time normalization).
• proton_energy::Real = 1.0: proton beam energy [GeV] (analytic mode only).

Returns

An SNSFlux instance.

get_params(use_data::Bool) -> NamedTuple

Return default flux normalization parameters.

• Data mode (use_data=true): returns (flux_norm=1.0, flux_onset=0.0).
• Analytic mode (use_data=false): returns (sns_nu_per_POT=0.09,).

Arguments

• use_data::Bool: flux mode (data/analytic) selector.

Returns

A NamedTuple with a single normalization parameter.

get_priors(use_data::Bool) -> NamedTuple

Return prior distributions for the flux normalization parameter.

• Data mode: Truncated(Normal(1.0, 0.1), 0.7, 1.3) for flux_norm and Uniform(-500.0, 500.0) for flux_onset.
• Analytic mode: Truncated(Normal(0.09, 0.009), 0.05, 0.15) for sns_nu_per_POT.

Arguments

• use_data::Bool: flux mode (data/analytic) selector.

Returns

A NamedTuple of Symbol => Distribution priors.

get_assets(; use_data, exposure, distance, beam_power, proton_energy, flux_folder, ecut, E_bin_width, tcut, tbins) -> NamedTuple

Precompute energy grids, time bins, and raw (unnormalized) flux arrays.

In data mode, reads 2D $(E, t)$ flux CSV files, applies energy and time cuts, and scales by the geometric factor $\eta = \text{exposure} \cdot s_{\text{per GWh}} \,/\, (4\pi d^2)$. Time is rebinned either into the explicit tbins edges, if given, or otherwise into a default set of logarithmically spaced bins.

In analytic mode, evaluates decay-at-rest spectra via flux_nu_mu, flux_nu_e, and flux_nu_mu_bar.

Arguments

• use_data::Bool: flux mode selector.
• exposure: integrated exposure [GWh].
• distance: source-to-detector distance [m].
• beam_power: beam power [MW].
• proton_energy: proton beam energy [GeV].
• flux_folder: path to CSV flux directory.
• ecut: energy cutoff [MeV].
• E_bin_width: energy bin width [MeV] (analytic mode).
• tcut: time cutoff [ns] (data mode, used only when tbins is not given).
• tbins: explicit time-bin edges (data mode only); Nothing (sentinel type, not a value) selects the default log-spaced binning.

Returns

A NamedTuple with fields E, T (data mode only), and per-flavour flux arrays.

get_flux(use_data, assets) -> Function

Construct the flux evaluation closure from precomputed assets.

The returned closure flux(params) -> NamedTuple applies the normalization parameter to the raw flux arrays and returns per-flavour fluxes and the total flux. In data mode, it additionally shifts each flavour's time histogram by params.flux_onset via shift_time_histogram before summing.

In data mode, the returned NamedTuple has fields: E, T, total_flux, flux_e, flux_mu, flux_mu_bar, flux_e_bar.

In analytic mode: E, total_flux, flux_e, flux_mu, flux_mu_bar.

Arguments

• use_data::Bool: flux mode (data/analytic) selector.
• assets::NamedTuple: precomputed energy grids and flux arrays from get_assets.

Returns

A closure flux(params) -> NamedTuple.

flux_nu_mu(E, E0, eta, bin_width) -> (E, flux)

Compute the $\nu_\mu$ flux from $\pi^+ \to \mu^+ \nu_\mu$ decay at rest.

The $\nu_\mu$ is monoenergetic at $E_0 = (m_\pi^2 - m_\mu^2)/(2 m_\pi)$, represented as a Gaussian-smeared delta function binned onto the energy grid.

Arguments

• E: energy grid [MeV].
• E0: pion decay-at-rest $\nu_\mu$ energy [MeV].
• eta: geometric normalization factor.
• bin_width: energy bin width [MeV].

Returns

A tuple (E, flux) where flux is the bin-integrated $\nu_\mu$ flux array.

flux_nu_e(E, eta, Emax, bin_width) -> (E, flux)

Compute the $\nu_e$ flux from $\mu^+ \to e^+ \nu_e \bar\nu_\mu$ decay at rest.

The spectrum follows the Michel distribution: $\frac{d\Phi}{dE} \propto 192\, \frac{E^2}{m_\mu^3} \left(\frac{1}{2} - \frac{E}{m_\mu}\right)$ for $0 \le E \le E_\text{max}$.

Arguments

• E: energy grid [MeV].
• eta: geometric normalization factor.
• Emax: maximum neutrino energy [MeV] ($\approx m_\mu / 2$).
• bin_width: energy bin width [MeV].

Returns

A tuple (E, flux) where flux is the bin-integrated $\nu_e$ flux array.

flux_nu_mu_bar(E, eta, Emax, bin_width) -> (E, flux)

Compute the $\bar\nu_\mu$ flux from $\mu^+ \to e^+ \nu_e \bar\nu_\mu$ decay at rest.

The spectrum follows: $\frac{d\Phi}{dE} \propto 64\, \frac{E^2}{m_\mu^3} \left(\frac{3}{4} - \frac{E}{m_\mu}\right)$ for $0 \le E \le E_\text{max}$.

Arguments

• E: energy grid [MeV].
• eta: geometric normalization factor.
• Emax: maximum neutrino energy [MeV] ($\approx m_\mu / 2$).
• bin_width: energy bin width [MeV].

Returns

A tuple (E, flux) where flux is the bin-integrated $\bar\nu_\mu$ flux array.

CevnsXsec <: Newtrinos.Physics

Configured Coherent Elastic neutrino-Nucleus Scattering (CEvNS) cross-section module.

Computes the differential cross section $d\sigma/dE_r$ for neutrino scattering off atomic nuclei, including optional beyond-Standard-Model deformations parameterized by cevns_xsec_a through cevns_xsec_d.

Fields

• params::NamedTuple: nominal values for $\sin^2\theta_W$, the nuclear radii Rn_* (one per isotope), and the four BSM deformation parameters.
• priors::NamedTuple: prior distributions for each parameter.
• diff_xsec::Function: closure returning the differential cross section. See get_diff_xsec for the call signature.

configure(isotopes, er_centers, enu_centers; ff_model=:helm, ff_kwargs=(;)) -> CevnsXsec

Create a CevnsXsec module from a list of isotope descriptors and energy grids.

Builds cross-section assets and assembles per-isotope parameters and priors via build_params_and_priors and get_assets.

Arguments

• isotopes: collection of isotope descriptors. Each entry must have the fields:
  • Rn_key::Symbol: parameter name for this isotope's nuclear radius (e.g. :Rn_Cs133).
  • Rn_nom::Float64: nominal nuclear radius $R_n$ [fm].
  • mass::Float64: nuclear mass [MeV/c²].
  • Z::Int: proton number.
  • N::Int: neutron number.
  • fraction::Float64: isotopic abundance fraction.
• er_centers: nuclear recoil energy grid [MeV].
• enu_centers: neutrino energy grid [MeV].
• ff_model::Symbol: nuclear form-factor model, one of :helm, :klein_nystrand, :sym_fermi (default :helm). See ffsq.
• ff_kwargs::NamedTuple: extra keyword arguments forwarded to the chosen form-factor model (e.g. (s_fm=0.9,)).

Returns

A CevnsXsec instance.

configure(params::NamedTuple, priors::NamedTuple) -> CevnsXsec

Create a CevnsXsec from externally supplied parameter and prior NamedTuples.

Use this constructor when the parameter set has already been assembled (e.g. by a previous call to build_params_and_priors) and no isotope list is available.

Arguments

• params::NamedTuple: parameter values (must include sin2thetaW, cevns_xsec_a–d, and one Rn_* entry per isotope).
• priors::NamedTuple: corresponding prior distributions.

Returns

A CevnsXsec instance.

build_params_and_priors(isotopes) -> (NamedTuple, NamedTuple)

Assemble default parameter values and prior distributions for a given set of isotopes.

Includes the weak mixing angle $\sin^2\theta_W$ and four BSM deformation parameters cevns_xsec_a through cevns_xsec_d that modify the SM differential cross section.

Per-isotope nuclear radius parameters Rn_* are appended automatically using the Rn_key field of each isotope descriptor.

Arguments

• isotopes: collection of isotope descriptors; each must have Rn_key::Symbol and Rn_nom::Float64.

Returns

A 2-tuple (params::NamedTuple, priors::NamedTuple).

get_assets(isotopes, er_centers, enu_centers; ff_model=:helm, ff_kwargs=(;)) -> NamedTuple

Package isotope data, energy grids, and the chosen form-factor model into an assets NamedTuple for use by get_diff_xsec.

Arguments

• isotopes: collection of isotope descriptors (same format as configure).
• er_centers: nuclear recoil energy grid [MeV].
• enu_centers: neutrino energy grid [MeV].
• ff_model::Symbol: nuclear form-factor model (default :helm). See ffsq.
• ff_kwargs::NamedTuple: extra keyword arguments forwarded to the form-factor model.

Returns

A NamedTuple with fields:

• isotopes: Dict{Symbol, NamedTuple} mapping each Rn_key to (mass, Z, N, fraction, Rn_nom).
• er_centers: the recoil energy grid.
• enu_centers: the neutrino energy grid.
• ff_model, ff_kwargs: the form-factor model choice and its extra keyword arguments.

ffsq(assets) -> (er,mn,rn) -> F(q)^2

Returns a closure that uses assets.ff_model and assets.ff_kwargs.

ffsq(er, mn, rn; model=:helm, ...)

Supported models (ASCII symbols):

• :helm
• :klein_nystrand
• :sym_fermi

ds(er, enu, params, nupar, Rn_key; ffsq_fn=ffsq) -> Matrix

Vectorized differential CEvNS cross section $d\sigma/dE_r$ on an $(N_{E_r} \times N_{E_\nu})$ grid.

Implements the BSM-extended SM formula (with $x = E_r/E_\nu$):

\[\frac{d\sigma}{dE_r} = \frac{G_F^2}{4\pi}\,M_N\,Q_W^2\,F^2(q^2)\; \max\!\Bigl(0,\; a + b\,\tfrac{M_N E_r}{E_\nu^2} + c\,x + d\,x^2\Bigr)\]

where $Q_W = N - (1 - 4\sin^2\theta_W)\,Z$ and a–d are the cevns_xsec_a–d deformation parameters (the SM cross section is recovered at their nominal values from build_params_and_priors). The bracket is additionally zeroed outside the kinematically allowed region $1 - M_N E_r / (2 E_\nu^2) \geq 0$.

Arguments

• er: recoil energy grid [MeV], shape (n_er,).
• enu: neutrino energy grid [MeV], shape (n_enu,).
• params::NamedTuple: must include sin2thetaW, cevns_xsec_a–d, and params[Rn_key].
• nupar: 3-tuple (mass, Z, N) — nuclear mass [MeV/c²], proton number, neutron number.
• Rn_key::Symbol: key used to look up the nuclear radius $R_n$ in params.
• ffsq_fn::Function: form-factor function, (er, mn, rn) -> F^2 (default ffsq).

Returns

Matrix of shape (n_er, n_enu) with units [MeV⁻¹].

get_diff_xsec(assets) -> Function

Build the differential cross-section closure over all isotopes from pre-computed assets.

The energy grids and isotope data are captured from assets at construction time. The returned function has signature:

diff_xsec(params::NamedTuple) -> Dict{Symbol, Matrix}

It iterates over all isotopes, calling ds for each, and returns a Dict mapping each Rn_key to its $(n_{E_r} \times n_{E_\nu})$ cross-section matrix.

Arguments

• assets: NamedTuple produced by get_assets.

Returns

A callable closure params -> Dict{Symbol, Matrix}.

get_diff_xsec_lar() -> Function

Return a closure that calls ds for the liquid-argon (LAr) target.

The returned function has signature:

diff_xsec_lar(er_centers, enu_centers, params, nupar, Rn_key) -> Matrix

get_diff_xsec_csi() -> Function

Return a closure that calls ds for the CsI target.

The returned function has signature:

diff_xsec_csi(er_centers, enu_centers, params, nupar, Rn_key) -> Matrix

XsecModel

Abstract type for neutrino cross-section models.

Each subtype defines how cross-section normalization factors are applied to event rates. Subtypes:

• SimpleScaling: global normalization for NC and $\nu_\tau$ CC channels.
• Differential_H2O: per-interaction-mode differential scaling for water Cherenkov detectors.
• H2O_PCA: GENIE-tune-based per-channel PCA shape and normalization systematics for water Cherenkov detectors.

SimpleScaling <: XsecModel

Simple cross-section model applying global normalization factors.

Provides two nuisance parameters:

• nc_norm: neutral-current cross-section scale factor.
• nutau_cc_norm: $\nu_\tau$ charged-current cross-section scale factor.

All other flavour/interaction combinations return unit scaling.

Differential_H2O <: XsecModel

Differential cross-section model for water (H₂O) targets.

Provides per-interaction-mode normalization parameters (cc1p1h_norm, cc2p2h_norm, cc1pi_norm, ccother_norm, ccdis_norm) in addition to nc_norm and nutau_cc_norm. Energy-dependent CC interaction fractions are interpolated from digitized data (T. Wester, Super-K PhD thesis, Figure 4.7) for $\nu_e$ and $\bar{\nu}_e$ separately.

Additional systematic parameters adjust the underlying shape of the interpolated curves: xsec_MA_QE/xsec_MA_Res (axial mass form-factor scaling for CCQE and resonance interactions), xsec_I12 (isospin-1/2 background normalization for CCother), xsec_fsi (final-state interaction strength, shifting weight between CC1p1h and CC1pi), an overall cc_norm, and nubar_ratio (the $\bar{\nu}/\nu$ cross-section ratio, applied in a normalization-conserving way).

Xsec <: Newtrinos.Physics

Configured cross-section physics module, returned by configure.

Fields

• cfg::XsecModel: the cross-section model used to build this module.
• params::NamedTuple: default normalization parameter values.
• priors::NamedTuple: prior distributions for each parameter.
• scale::Function: closure that computes the cross-section scale factor. See get_scale for the call signatures.
• dσdE::Union{Function, Nothing}: optional closure returning the absolute differential cross-section, when supported by cfg; nothing otherwise.
• scale_event::Union{Function, Nothing}: optional per-event reweighting closure using generator interaction codes, when supported by cfg; nothing otherwise.
• event_weights::Union{Function, Nothing}: optional precomputed fast-path per-event reweighting closure, when supported by cfg; nothing otherwise.
• grid_weights::Union{Function, Nothing}: optional precomputed fast-path per-energy-grid reweighting closure, when supported by cfg; nothing otherwise.

configure(cfg::XsecModel=SimpleScaling()) -> Xsec

Create a fully configured cross-section physics module.

Assembles default parameters, priors, and the scaling closure (plus any optional dσdE, per-event, or grid-reweighting closures supported by cfg, falling back to nothing when a model does not implement them) into an Xsec struct ready for use in experiment forward models.

Arguments

• cfg::XsecModel: cross-section model (defaults to SimpleScaling).

Returns

An Xsec instance.

Examples

using Newtrinos

# Default simple scaling
xsec_physics = Newtrinos.xsec.configure()

# Differential water cross-sections for Super-K
xsec_physics = Newtrinos.xsec.configure(Differential_H2O())

get_params(cfg::XsecModel) -> NamedTuple

Return the default cross-section normalization parameter values for the given model.

Arguments

• cfg::XsecModel: a SimpleScaling or Differential_H2O instance.

Returns

A NamedTuple mapping parameter names to their nominal values. Most normalization parameters default to 1.0; a few Differential_H2O systematic parameters (xsec_MA_QE, xsec_MA_Res, xsec_I12) default to their central physical values instead.

get_priors(cfg::XsecModel) -> NamedTuple

Return prior distributions for each cross-section normalization parameter.

Most priors are truncated normal distributions centred near 1.0; a few Differential_H2O systematic parameters (xsec_MA_QE, xsec_MA_Res, xsec_I12, xsec_fsi) use untruncated Normal priors centred on their nominal physical values instead.

Arguments

• cfg::XsecModel: a SimpleScaling or Differential_H2O instance.

Returns

A NamedTuple mapping parameter names to Distributions.UnivariateDistribution priors.

get_scale(cfg::XsecModel) -> Function

Construct the cross-section scaling closure for the given model.

SimpleScaling returns a function with signature:

scale(flav::Symbol, interaction::Symbol, params::NamedTuple) -> Real

Returns params.nc_norm for NC interactions, params.nutau_cc_norm for $\nu_\tau$ CC, and 1.0 for everything else.

Differential_H2O returns a function with signature:

scale(E::AbstractArray, flav::Symbol, interaction::Symbol, anti::Bool, params::NamedTuple) -> Real or AbstractArray

For NC interactions returns params.nc_norm. For CC interactions computes an energy-dependent weighted sum of per-mode normalizations using interpolated interaction fractions, further reweighted by axial-mass form-factor ratios (xsec_MA_QE, xsec_MA_Res), an FSI shift (xsec_fsi) between CC1p1h/CC1pi, an isospin-1/2 background scale (xsec_I12) for CCother, an overall cc_norm, and a normalization-conserving $\bar{\nu}/\nu$ ratio (nubar_ratio), with an additional $\nu_\tau$ CC factor when applicable.

Arguments

• cfg::XsecModel: cross-section model instance.

Returns

A closure computing the cross-section scale factor.

bin(df, out_centers; col=1)

Bins a list of events in DataFrame or Matrix df into out_centers using the specified column (col) for the observable. Returns a vector of counts per bin.

rebin(df, out_centers; var_col=1, count_col=2)

Rebins a table of (observable, count) pairs in DataFrame or Matrix df into out_centers. var_col specifies which column contains the observable to bin in. count_col specifies which column contains the counts. Returns a vector of summed counts per bin.